Overlaps corrected, new shape of sectors
[u/mrichter/AliRoot.git] / 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 AliRunLoader::Instance();
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=AliRunLoader::GetNumberOfEvents();
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 }