DanteWatchers.C

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2010 by Olivier Stezowski                               *
 *   stezow(AT)ipnl.in2p3.fr                                               *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/


#include "DanteWatchers.h"
#include "TMath.h"

#include <fstream>

Short_t BaseDanteWatcher::gLUT[BaseDanteWatcher::gNumberOfItems][BaseDanteWatcher::gMaxNumberOfModules];

ClassImp(BaseDanteWatcher);

void MyDante::Calibrate () 
{
        Double_t xtmp = X, ytmp = Y;
        
        X = AX + BX*xtmp + CX*TMath::Power(xtmp,2) ;
        Y = AY + BY*ytmp + CY*TMath::Power(ytmp,2) ; 
        
        if ( X < threshold_x_left || X > threshold_x_right ) {
                has_fired = false;
        }
        if ( Y < threshold_y_down || Y > threshold_y_up ) {
                has_fired = false;
        }               
}

BaseDanteWatcher::BaseDanteWatcher(const char *name, const char *title,TDirectory *sp_dir, TDirectory *tag_dir): 
RancLegnaroWatcher(name,title,sp_dir,tag_dir)
{
        // to be sure the LUT is initialized
        InitLUT();
        fNumberOfModules = 0;
        
        if ( !LoadCalibration("Dante.conf") )
                fLog << warning << "Dante.conf not found in the current directory ... LoadCalibration should be called at least once !! " << dolog;
} 

BaseDanteWatcher::~BaseDanteWatcher()
{       
}

// should add nmore test to be sure the file is consistant ...
Bool_t BaseDanteWatcher::LoadCalibration(const char *calib_file)
{
        std::ifstream in(calib_file); std::string tmp;
        
        if ( in.is_open() == false ) 
                return false;
        
        getline(in,tmp);
        tmp += " ";
        while ( in.good() && !in.eof() ) { // read input stream line by line
                
                if ( tmp[0] == '#' ) {
                        getline(in,tmp); 
                        continue; 
                } // this line is a comment 
                {
                        Float_t dz, rz1, ry, rz2, rx, tx1, tx2,ty1, ty2, ax, bx, cx, ay, by, cy ;
                        Int_t which; 
                        std::istringstream decode(tmp);
                        
                        decode >> which;
                        if ( which == -1 ) { // this is LUT
                                for (Short_t i = 0; i < fNumberOfModules; i++ )  
                                        for (Short_t j = 0; j < gNumberOfItems; j++ )  
                                                decode >> gLUT[j][i] ;
                                
                                /*                              cout << " LUT " << endl;
                                 for (Short_t i = 0; i < fNumberOfModules; i++ )  
                                 for (Short_t j = 0; j < gNumberOfItems; j++ )  
                                 cout << j << " " <<  i << " " <<  gLUT[j][i] << endl ; */
                        }
                        else { // this a module
                                if ( which < gMaxNumberOfModules ) {
                                        
                                        decode >> dz >> rz1 >> ry >> rz2 >> rx >> tx1 >> tx2 >> ty1 >> ty2 >> ax >> bx >> cx >> ay >> by >> cy ;
                                        
                                        fNumberOfModules++;
                                        
                                        fDante[which].Z  = dz;
                                        
                                        fDante[which].AX = ax;
                                        fDante[which].BX = bx;
                                        fDante[which].CX = cx;
                                        
                                        fDante[which].AY = ay;
                                        fDante[which].BY = by;
                                        fDante[which].CY = cy;
                                        
                                        /*
                                         fDante[which].AX = ay;
                                         fDante[which].BX = by;
                                         fDante[which].CX = cy;
                                         
                                         fDante[which].AY = ax;
                                         fDante[which].BY = bx;
                                         fDante[which].CY = cx; */
                                        
                                        
                                        fDante[which].threshold_x_left  = tx1;
                                        fDante[which].threshold_x_right = tx2;
                                        fDante[which].threshold_y_down  = ty1;
                                        fDante[which].threshold_y_up    = ty2;
                                        
                                        fRotation_tmp[0][which].RotateZ(rz1*TMath::Pi()/180.);
                                        fRotation_tmp[1][which].RotateY(ry*TMath::Pi()/180.);
                                        fRotation_tmp[2][which].RotateZ(rz2*TMath::Pi()/180.);
                                        fRotation_tmp[3][which].RotateX((360.-rx)*TMath::Pi()/180.);
                                        
                                        fRotation[which].RotateZ(rz1*TMath::Pi()/180.);
                                        fRotation[which].RotateY(ry*TMath::Pi()/180.);
                                        fRotation[which].RotateZ(rz2*TMath::Pi()/180.);
                                        fRotation[which].RotateX((360.-rx)*TMath::Pi()/180.);                           
                           }
                           else {;}
                        }
                }
                // add null character at the end since getline
                // does not store \n and it may give troubles when decoding
                getline(in,tmp);
                tmp += " ";
        }       
        in.close();
        
        fLog << " Dante has been configured with " << fNumberOfModules << dolog;
        
        return true;    
}

void BaseDanteWatcher::InitLUT()
{
        // module 1
        gLUT[0][0] = 96;
        gLUT[1][0] = 97;
        gLUT[2][0] = 0;
        // module 2
        gLUT[0][1] = 98;
        gLUT[1][1] = 99;
        gLUT[2][1] = 0;
        // module 3
        gLUT[0][2] = 100;
        gLUT[1][2] = 101;
        gLUT[2][2] = 0;
        // module 2
        gLUT[0][3] = 0;
        gLUT[1][3] = 0;
        gLUT[2][3] = 0;
        // module 3
        gLUT[0][4] = 0;
        gLUT[1][4] = 0;
        gLUT[2][4] = 0;
}

void BaseDanteWatcher::ComputeRealPosition(Short_t which)
{
        if ( !(which < fNumberOfModules) )
                return;
        
        TVector3 v(fDante[which].X,fDante[which].Y,fDante[which].Z);
        
        v = fRotation_tmp[0][which]*v; v = fRotation_tmp[1][which]*v; v = fRotation_tmp[2][which]*v; v = fRotation_tmp[3][which]*v;
        //      v = fRotation[which]*v;
        
        fDante[which].Xr = v.X();
        fDante[which].Yr = v.Y();
        fDante[which].Zr = v.Z();
}


void BaseDanteWatcher::Exec(Option_t *option) 
{       
        // be sure the frame has been set properly
        if ( fFrame == 0x0 ) 
        { SetLastExecStatus(1); return; }
        else {
                if ( !fFrame->IsValid() )
                { SetLastExecStatus(1); return; }
        }
        
        // fill the array of floats
        RancLegnaroWatcher::Exec(option);
        
        // fill up the different modules : calibrates and check the threshold
        // to know the number of detector fired in this event
        for (Int_t i = 0u; i < fNumberOfModules; i++) {
                
                fDante[i].has_fired = true;
                
                if (   rawBuf[gLUT[0][i]] < fDante[i].threshold_adc 
                    || rawBuf[gLUT[1][i]] < fDante[i].threshold_adc 
                    || rawBuf[gLUT[2][i]] < fDante[i].threshold_adc ) {
                        
                        //                      fDante[i].has_fired = false;
                }
                
                // fill the data part
                fDante[i].X = rawBuf[gLUT[0][i]];
                fDante[i].Y = rawBuf[gLUT[1][i]];
                fDante[i].T = rawBuf[gLUT[2][i]];
                fDante[i].A = rawBuf[gLUT[3][i]];
                
                fDante[i].Calibrate(); ComputeRealPosition(i);
        }
}

#ifndef Gw_AgataEventContainer
#include "AgataEventContainer.h"
#endif

#ifndef Gw_AgataEventDisplay
#include "AgataEventDisplay.h"
#endif

#ifndef Gw_AgataGeometryTransformer
#include "AgataGeometryTransformer.h"
#endif

ClassImp(DisplayDante);

Float_t DisplayDante::fgMetric = 10.;

DisplayDante::DisplayDante(const char *name, const char *title,TDirectory *sp_dir, TDirectory *tag_dir): 
BaseDanteWatcher(name,title,sp_dir,tag_dir),
fCounts(0x0),
fPositions(0x0),
fPositionsC(0x0),
fPos3D(0x0),
fR(0x0),
fTheta(0x0),
fPhi(0x0),
fThetaxPhi(0x0),
fMaxCycle(100),
fConfFileName(Gw::AgataEventDisplay::GetDefaultDanteConf()),
fPathFile(Gw::AgataEventDisplay::GetDefaultAgataPath())
{
        
        fTypeName =  AgataEventDisplay::Instance()->RegisterBranch();
        
        AgataEventDisplay::Instance()->SetTrackStyle("Rectangle", fTypeName);
        
        fCounts = MakeTH1<TH1F>("NbModules","To get the number of modules fired",gMaxNumberOfModules,0,gMaxNumberOfModules);
        
        // one histogram for each module
        fPositions = new TObjArray(gMaxNumberOfModules); 
        fPositions->SetOwner(false);
        fPositionsC = new TObjArray(gMaxNumberOfModules); 
        fPositionsC->SetOwner(false);
        
        for (Short_t i = 0u; i < gMaxNumberOfModules; i++) {
                
                TString tmp; TH2F *h;
                
                tmp = Form("Dante_X_Y_%d",i);
                h = MakeTH2<TH2F>(tmp.Data(),"Raw Position in Dante",2048,0,4096,2048,0,4096);
                fPositions->Add ( h ) ;
                
                tmp = Form("Dante_X_Y_%d_C",i);
                h = MakeTH2<TH2F>(tmp.Data(),"Calibrated Positions in Dante;X;Y",1000,-30,30,1000,-30,30);
                fPositionsC->Add ( h ) ;
        }
        fPos3D = MakeTH3<TH3F>("Pos3D","Positions of recoils in beam system;X;Y;Z",100,-150,150,100,-150,150,100,-150,150);
        
        fR = MakeTH1<TH1F>("R","R position",100,0,200);
        fTheta = MakeTH1<TH1F>("Theta","Theta angles",90,0,180);
        fPhi = MakeTH1<TH1F>("Phi","Phi angles",90,-180,180);
        
        fThetaxPhi = MakeTH2<TH2F>("ThetaxPhi","Raw Position in Dante",180,0,180,360,-180,180);
} 

DisplayDante::~DisplayDante()
{ 
        
}

void DisplayDante::SetActive(Bool_t active) 
{
        TTask::SetActive(active); 
        AgataEventDisplay::Instance()->Reset();
} 

void DisplayDante::DoCanvas(TCanvas *c, Option_t *o)
{
        TString option = o;
        if ( option == "" ) {
                
                UShort_t which_pad = 1u;
                c->Divide(fNumberOfModules, 2, 0.001, 0.001); // 
                
                for (Short_t i = 0u; i < fNumberOfModules; i++) {
                        c->cd(which_pad++);
                        fPositions->At(i)->Draw();
                }
                for (Short_t i = 0u; i < fNumberOfModules; i++) {
                        c->cd(which_pad++);
                        fPositionsC->At(i)->Draw();
                }
        }
        if ( option.Contains("3D") ) {
                printf(" ==> OPTION %s \n",o);
                if ( option.Contains("OGL") )
                        fPos3D->Draw("glbox1");
                else {
                        fPos3D->Draw("col");
                }
        }       
        if ( option == "2D" )  {
                c->Divide(2, 2, 0.001, 0.001); // 
                
                c->cd(1);
                fR->Draw();
                c->cd(2);
                fTheta->Draw();         
                c->cd(3);
                fPhi->Draw();
                c->cd(4);
                fThetaxPhi->Draw("col");
        }
}

void DisplayDante::ShowEve(const Char_t* danteConfFile, const Char_t* agataPathFile)
{
   AgataEventDisplay::SetDefaultDanteConf(danteConfFile);
   AgataEventDisplay::SetDefaultAgataPath(agataPathFile);   
   AgataEventDisplay::Instance()->BuildDefaultGeometry(true, false);
   AgataEventDisplay::Instance()->ShowDisplay();
}

void DisplayDante::Exec(Option_t *option) 
{       
        // fill the array of floats
        BaseDanteWatcher::Exec(option);
        
        // count the nummber of detectors that have fired
        UShort_t nb_fired = 0u;
        /*
         AgataEventContainer* agataContainer = AgataEventDisplay::Instance()->GetEventContainer(); 
         
         if (agataContainer->GetNofEvents() > fMaxCycle) {
         TTask::SetActive(false);
         }
         */ 
        for (Int_t i = 0u; i < fNumberOfModules; i++) {
                // TrackHit* hit = 0x0; 
                if ( fDante[i].has_fired ) {
                        //                hit = agataContainer->NewTrackHit(fTypeName);
                        nb_fired++;
                } 
                TH2F *h;
                // raw positions
                h = (TH2F *)fPositions->At(i);                  
                h->Fill(rawBuf[gLUT[0][i]],rawBuf[gLUT[1][i]]);
                
                // calibrated
                if ( fDante[i].has_fired ) {
                        h = (TH2F *)fPositionsC->At(i); h->Fill(fDante[i].X,fDante[i].Y);
                        /*                      StdHit* thit = hit->NewHit();
                         thit->SetX(-fDante[i].X/fgMetric);
                         thit->SetY(fDante[i].Y/fgMetric);
                         thit->SetZ(fDante[i].Z/fgMetric);
                         thit->SetE(fDante[i].A);
                         StdHit* shit = agataContainer->NewHit(fTypeName);
                         shit->SetX(-fDante[i].X/fgMetric);
                         shit->SetY(fDante[i].Y/fgMetric);
                         shit->SetZ(fDante[i].Z/fgMetric);
                         shit->SetE(fDante[i].A); */
                }
                //cout << fDante[i].Xr << " " << fDante[i].Yr <<  " " << fDante[i].Zr << endl;
                
                if ( fDante[i].has_fired ) {
                        fPos3D->Fill(fDante[i].Xr,fDante[i].Yr,fDante[i].Zr);
                        
                        TVector3 vtmp(fDante[i].Xr,fDante[i].Yr,fDante[i].Zr);
                        fR->Fill(vtmp.Mag());
                        fTheta->Fill(vtmp.Theta()*180./TMath::Pi());
                        fPhi->Fill(vtmp.Phi()*180./TMath::Pi());
                        fThetaxPhi->Fill(vtmp.Theta()*180./TMath::Pi(),vtmp.Phi()*180./TMath::Pi());
                }
        }
        fCounts->Fill(nb_fired);
        //      agataContainer->FillTracks(fTypeName);
        //      agataContainer->FillHits(fTypeName);
}


ClassImp(MyDanteTree);

MyDanteTree::MyDanteTree(const char *name, const char *title, TTree *tree) : 
BaseDanteWatcher(name,title), TTreeBuilder(name,title,tree)
{
        fWhichDante = new UShort_t [gMaxNumberOfModules]; 
        // hit in the dante module in global reference system
        fDanteXr = new Float_t[gMaxNumberOfModules]; //[fNbModule]   
        fDanteYr = new Float_t[gMaxNumberOfModules]; //[fNbModule]   
        fDanteZr = new Float_t[gMaxNumberOfModules]; //[fNbModule]   
        // amplitude of the signal
        fDanteA = new Float_t[gMaxNumberOfModules]; //[fNbModule]   
        // calibrated poistion in the module reference system
        fDanteXc = new Float_t[gMaxNumberOfModules]; //[fNbModule]   
        fDanteYc = new Float_t[gMaxNumberOfModules]; //[fNbModule]   
        
        if ( GetTree() ) { // just add the branch for that watcher
                SetBranches();
        }
}

void MyDanteTree::SetBranches()
{
        GetTree()->Branch("NbDante", &fNbDante, "NbDante/I");
        GetTree()->Branch("WhichDante", fWhichDante, "WhichDante[NbDante]/s");
        
        GetTree()->Branch("DanteXc", fDanteXc, "Xc[NbDante]/F");
        GetTree()->Branch("DanteYc", fDanteYc, "Yc[NbDante]/F");
        
        GetTree()->Branch("DanteA", fDanteA, "A[NbDante]/F");
        
        GetTree()->Branch("DanteXr", fDanteXr, "DanteXr[NbDante]/F");
        GetTree()->Branch("DanteYr", fDanteYr, "DanteYr[NbDante]/F");
        GetTree()->Branch("DanteZr", fDanteZr, "DanteZr[NbDante]/F");
        
        GetTree()->Branch("DanteBeta", &fDanteBeta, "DanteBeta/F");
}

MyDanteTree::~MyDanteTree()
{
        delete [] fWhichDante;
        delete [] fDanteXr;
        delete [] fDanteYr;
        delete [] fDanteZr;
        delete [] fDanteA;
        delete [] fDanteXc;
        delete [] fDanteYc;
}

void MyDanteTree::Exec(Option_t *option) 
{
        if ( gDebug > 3 ) 
                printf("MyDanteTree::Exec is called \n");
        
        // fill the array of floats and compute the physics
        BaseDanteWatcher::Exec(option);
        
        fNbDante = 0;
        for (Int_t i = 0u; i < fNumberOfModules; i++) {
                if ( !fDante[i].has_fired ) {
                        continue;
                } 
                
                fWhichDante[fNbDante] = i;
                fDanteXc[fNbDante] = fDante[i].X;
                fDanteYc[fNbDante] = fDante[i].Y;
                
                fDanteXr[fNbDante] = fDante[i].Xr;
                fDanteYr[fNbDante] = fDante[i].Yr;
                fDanteZr[fNbDante] = fDante[i].Zr;
                
                fDanteA[fNbDante] = fDante[i].A;
                fDanteBeta = 0.5;
                
                fNbDante++;
                // printf("MyDanteTree::Exec is called with %d %d\n", fDante[i].X, fDante[i].Y);
        }
        
        // fill the tree only if this is the owner of the tree ... otherwise by TreeMaster
        FillTree(); 
}

ClassImp(DoDanteDoppler);

DoDanteDoppler::DoDanteDoppler(const char *name, const char *title,TDirectory *sp_dir, TDirectory *tag_dir): 
        BaseDanteWatcher(name,title,sp_dir,tag_dir), 
        VertexBuilder()
{
        fVertexFrame = MainFrameFactory::theMainFactory().
        NewSharedFrame( FactoryItem("Agata","meta:vertex",Version(0,0)), FactoryItem("Agata","meta:vertex",Version(1,1)) );
        
        fIsToBeDeleted = true; 
        
        fMeanBeta = 0.06;
} 

void DoDanteDoppler::ChangeMeanBeta(Double_t mean)
{
        fMeanBeta = mean;
} 

void DoDanteDoppler::Exec(Option_t *option) 
{       
        // fill the array of floats and compute the dante physics
        BaseDanteWatcher::Exec(option);
        
        // look to get the one that fired
        Short_t nb_fired = 0u, which = -1;
        for (Int_t i = 0u; i < fNumberOfModules; i++) {
                
                if ( !fDante[i].has_fired )
                        continue;
                
                nb_fired++; which = i;
        }
        if ( nb_fired == 1 && which == 1 )  {
                
                // fill the current vertex
                
                VertexInterface *gdata = GetVertex();;
                gdata->SetPosition(0,0,0);
                gdata->SetDirection(fDante[which].Xr,fDante[which].Yr,fDante[which].Zr);        
                gdata->SetBeta(fMeanBeta); 
        }
        else {
                VertexInterface *gdata = GetVertex();
                gdata->SetPosition(0,0,0);
                gdata->SetDirection(fDante[which].Xr,fDante[which].Yr,fDante[which].Zr);        
                gdata->SetBeta(0); 
        }
        
}