ITS/AliITSCompareHitsRecPoints.C

   1 Bool_t AliITSCompareHitsRecPoints(Char_t *rfn="galice.root"){
   2     // This macro compares the average location of the hits associated
   3     // with a specific track and the best match RecPoint location.
   4     // Inputs:
   5     //   none.
   6     // Outputs:
   7     //   none.
   8     // Return:
   9     //   kTRUE if no errors were encountered, otherwise kFALSE.
  10     TProfile *pxg,*pyg,*pzg,*pxl,*pyl,*pzl;
  11     Double_t hg[3],hl[3],tof,rg[3],rl[3];
  12     Float_t rgf[3];
  13
  14     // Dynamically link some shared libs
  15     if (gClassTable->GetID("AliRun") < 0) {
  16         gROOT->LoadMacro("loadlibs.C");
  17         loadlibs();
  18     } else {
  19         if(gAlice){
  20             delete gAlice->GetRunLoader();
  21             delete gAlice;
  22             gAlice=0;
  23         } // end if gAlice
  24     }// end if gClassTable->GetID()
  25
  26     gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSstandard.C");
  27
  28     AliRunLoader *runl = AccessFile(rfn); // Set up to read in Data
  29     Int_t retval = runl->LoadHeader();
  30     if (retval){
  31         cerr<<"AliITSPrintRecPoints.C : LoadHeader returned error"<<endl;
  32         return kFALSE;
  33     } // end if retval
  34
  35     AliITSLoader* ITSloader =  (AliITSLoader*) runl->GetLoader("ITSLoader");
  36     if(!ITSloader){
  37         cerr<<"AliITSPrintRecPoints.C :  ITS loader not found"<<endl;
  38         return kFALSE;
  39     } // end if !ITSloader
  40
  41     if(!gGeoManager){
  42         gGeoManger = new TGeoManager();
  43         gGeoManager->Import("geometry.root","");
  44     } // end if !gGeoManager
  45     if(ITSloader->LoadHits("read")!=0){
  46         cerr<<"Error Loading Hits"<<endl;
  47         return kFALSE;
  48     }// end if ITSloader->LoadHits
  49     if(ITSloader->LoadRecPoints("read")!=0){
  50         cerr<<"Error Loading RecPoints"<<endl;
  51         return kFALSE;
  52     }// end if ITSloader->LoadRecPoints
  53     AliITS *ITS = 0;
  54     ITS = (AliITS*)(gAlice->GetDetector("ITS"));
  55     if(!ITS){
  56         cout << "Error: no ITS found. Aborting"<<endl;
  57         return kFALSE;
  58     } // end if !ITS
  59     if(!(ITS->GetDetTypeSim())){
  60         cout <<"No AliITSDetTypeSim object found in ITS"<<endl;
  61         return kFALSE;
  62     } // end if
  63     AliITSmodule     *m   = 0;
  64     AliITShit        *h   = 0;
  65     AliITSDetTypeSim *sim = ITS->GetDetTypeSim();
  66     AliITSgeom *gm=0;
  67     gm = ITS->GetITSgeom();
  68     if(!gm){
  69         cout <<"No AliITSgeom object found in ITS"<<endl;
  70         if(!gGeoManager){
  71             cout <<"No gGeoManger. Aborting"<<endl;
  72             return kFALSE;
  73         }else{
  74             ITS->UpdateInternalGeometry(gGeoManager);
  75         } // end if
  76     } // end if !AliITSgeom
  77     //
  78     Int_t nMods= gm->GetIndexMax(),nEvents=gAlice->GetEventsPerRun();
  79     Int_t mod=0,evnt=0,size=-1,irp=0,ih=0,trkindexOld=-1;
  80     Double_t xmod,nHitPerTrack;
  81     TTree            *rpt = 0;
  82     TClonesArray     *rpa = 0;
  83     AliITSRecPoint   *rp  = 0;
  84     AliITSDetTypeRec *rec = new AliITSDetTypeRec(ITSloader);
  85     rec->SetDefaults();
  86     // We are going to use the following
  87     // <hit_x[i]> == { Sum_j hit_x[i][j] }/N_hits per track
  88     // d_x[i] == <hit_x[i]>-recpoint_x[i]
  89     // <d_x>  == {sum_i d_x[i]}/N_RecPoints   (excluding noise and no merging)
  90     //         = {sum_i (<hit_x[i]>-recpoint_x[i])}/N_RecPoints
  91     //         = {sum_i ((sum_j hit_x[i][j])/N_hits_per_track -recpoinnt_x[i])}/N_Recpoints
  92     //         = sum_i sum_j {hit_x[i][j]/N_hits_per_track}/N_Recpoints -
  93     //           sum_i {recpoint_x[i]}/N_Recpoints
  94     pxg = new TProfile("XdiffGlobal","Mean displacement in gX direction",nMods,
  95                        -0.5,(Double_t)(nMods)+0.5," ");
  96     pxg->SetXTitle("module");
  97     pxg->SetYTitle("Global X [cm]");
  98     pyg = new TProfile("YdiffGlobal","Mean displacement in gY direction",nMods,
  99                        -0.5,(Double_t)(nMods)+0.5," ");
 100     pyg->SetXTitle("module");
 101     pyg->SetYTitle("Global Y [cm]");
 102     pzg = new TProfile("ZdiffGlobal","Mean displacement in gZ direction",nMods,
 103                        -0.5,(Double_t)(nMods)+0.5," ");
 104     pzg->SetXTitle("module");
 105     pzg->SetYTitle("Global Z [cm]");
 106     pxl = new TProfile("XdiffLocal","Mean displacement in lX direction",nMods,
 107                        -0.5,(Double_t)(nMods)+0.5," ");
 108     pxl->SetXTitle("module");
 109     pxl->SetYTitle("local X [cm]");
 110     pyl = new TProfile("YdiffLocal","Mean displacement in lY direction",nMods,
 111                        -0.5,(Double_t)(nMods)+0.5," ");
 112     pyl->SetXTitle("module");
 113     pyl->SetYTitle("local Y [cm]");
 114     pzl = new TProfile("ZdiffLocal","Mean displacement in lZ direction",nMods,
 115                        -0.5,(Double_t)(nMods)+0.5," ");
 116     pzl->SetXTitle("module");
 117     pzl->SetYTitle("local Z [cm]");
 118     //
 119     for(evnt=0;evnt<nEvents;evnt++){
 120         runl->GetEvent(evnt);
 121         ITS->InitModules(size,nMods);
 122         ITS->FillModules(evnt,0,-1," "," ");
 123         rec->SetTreeAddress();
 124         for(mod=0;mod<nMods;mod++){
 125             xmod = (Double_t) mod;
 126             m = ITS->GetModule(mod);
 127             rec->ResetRecPoints();
 128             rpt = ITSloader->TreeR();
 129             rpt->GetEvent(mod);
 130             rpa = rec->RecPoints();
 131             for(irp=0;irp<rpa->GetEntriesFast();irp++){
 132                 rp = (AliITSRecPoint*)(rpa->At(irp));
 133                 rp->GetGlobalXYZ(rgf);rg[0]=rgf[0];rg[1]=rgf[1];rg[2]=rgf[2];
 134                 rl[0] = rp->GetDetLocalX();
 135                 rl[1] = 0.0;
 136                 rl[2] = rp->GetDetLocalZ();
 137                 pxg->Fill(xmod,-rg[0],0.5);
 138                 pyg->Fill(xmod,-rg[1],0.5);
 139                 pzg->Fill(xmod,-rg[2],0.5);
 140                 pxl->Fill(xmod,-rl[0],0.5);
 141                 //pyl->Fill(xmod,-rl[1],0.5); // assumed to be zero always.
 142                 pzl->Fill(xmod,-rl[2],0.5);
 143             }// rnf got itp
 144             trkindexOld=-1;
 145             nHitPerTrack = 0.0;
 146             for(ih=0;ih<m->GetNhits();ih++){ // We want the median hit location
 147                 // for each track.
 148                 h = m->GetHit(ih);
 149                 if(m->GetHitTrackIndex(ih)!=trkindexOld){// Enterence location
 150                     trkindexOld = m->GetHitTrackIndex(ih);
 151                     nHitPerTrack = 1.0;
 152                     do{
 153                         nHitPerTrack += 1.0;
 154                     }while (m->GetHitTrackIndex(ih+nHitPerTrack-1)==trkindexOld);
 155                     h->GetPositionG0(hg[0],hg[1],hg[2],tof);
 156                     h->GetPositionL0(hl[0],hl[1],hl[2],tof);
 157                     pxg->Fill(xmod,hg[0],0.5/nHitPerTrack);
 158                     pyg->Fill(xmod,hg[1],0.5/nHitPerTrack);
 159                     pzg->Fill(xmod,hg[2],0.5/nHitPerTrack);
 160                     pxl->Fill(xmod,hl[0],0.5/nHitPerTrack);
 161                     pyl->Fill(xmod,hl[1],1.0/nHitPerTrack); // rl[1]=0 always.
 162                     pzl->Fill(xmod,hl[2],0.5/nHitPerTrack);
 163                 } // end if
 164                 if(nHitPerTrack==0.0) continue;
 165                 h->GetPositionG(hg[0],hg[1],hg[2],tof);
 166                 h->GetPositionL(hl[0],hl[1],hl[2],tof);
 167                 pxg->Fill(xmod,hg[0],0.5/nHitPerTrack);
 168                 pyg->Fill(xmod,hg[1],0.5/nHitPerTrack);
 169                 pzg->Fill(xmod,hg[2],0.5/nHitPerTrack);
 170                 pxl->Fill(xmod,hl[0],0.5/nHitPerTrack);
 171                 pyl->Fill(xmod,hl[1],1.0/nHitPerTrack); // rl[1]=0 always.
 172                 pzl->Fill(xmod,hl[2],0.5/nHitPerTrack);
 173             } // end for ih
 174         } // end for mod
 175         ITS->ClearModules();
 176     } // end for evnt
 177     //
 178     Int_t wh=800;
 179     TCanvas *c0 = new TCanvas("c0","Average displacements between hits and RecPoints"
 180                               ,1.5*wh,wh);
 181     c0->Divide(3,2);
 182     c0->cd(1);
 183     pxg->Draw();
 184     c0->cd(2);
 185     pyg->Draw();
 186     c0->cd(3);
 187     pzg->Draw();
 188     c0->cd(4);
 189     pxl->Draw();
 190     c0->cd(5);
 191     pyl->Draw();
 192     c0->cd(6);
 193     pzl->Draw();
 194     c0->Update();
 195     return kTRUE;
 196 }