Changes required by Coding Conventions
[u/mrichter/AliRoot.git] / TOF / AliTOFanalyzeHits.C
CommitLineData
32bd3eb5 1Int_t AliTOFanalyzeHits(Int_t nevents = 1, Bool_t drawing = kFALSE)
42036c57 2{
3 //
4 // Analyzes the hits and fills QA-histograms
32bd3eb5 5 // Use case:
6 // start aliroot
7 // root [0] .L AliTOFanalyzeHits.C
8 // root [1] AliTOFanalyzeHits()
42036c57 9 //
32bd3eb5 10 // Updated to the new I/O by: A. De Caro, C. Zampolli
42036c57 11
12 Int_t rc = 0;
32bd3eb5 13
42036c57 14 // Define the histograms
15 // x,y,z, rho, tof, padx, padz, sector, plate, strip, (x vs y)
16
17 // hit-map in a plane
18 TH2F *h2hitMap = new TH2F("h2hitMap","Hit Map (projection on the plane)",2500,-12500.,12500.,800,-400.,400.);
19 // time of flight distribution for primaries and secondaries
20 TH1F *htofp = new TH1F("htofp","Time of Flight (primaries)",800,0.,80.);
21 TH1F *htofs = new TH1F("htofs","Time of Flight (secondaries)",800,0.,80.);
22 // momentum when striking the TOF for primaries and secondaries
23 TH1F *htofmomp = new TH1F("htofmomp","Momentum at TOF (primaries)",100,0.,10.);
24 TH1F *htofmoms = new TH1F("htofmoms","Momentum at TOF (secondaries)",100,0.,10.);
25 // TOF hit volumes
da3d3acd 26 TH1F *hsector = new TH1F("hsector","Sector",18,0.,18.);
27 TH1F *hplate = new TH1F("hplate","Plate ", 5,0., 5.);
28 TH1F *hstrip = new TH1F("hstrip","Strip ",20,0.,20.);
29 TH1F *hpadz = new TH1F("hpadz","Pad along z ",2,0.,2.);
30 TH1F *hpadx = new TH1F("hpadx","Pad along x",48,0.,48.);
42036c57 31 // track length when striking the TOF (used by AliTOFT0)
32 TH1F *htrackLenp= new TH1F("htrackLenp","Track Length on TOF for Primaries",800,0.,800.);
33
32bd3eb5 34 // Histograms added to control the right TOF element numbering:
35 // it should be increasing with the azimuthal and polar angles
36
da3d3acd 37 TH2F *hmoduleVStheta = new TH2F("hmoduleVStheta", "hmoduleVStheta", 180,0.,180.,5,0,5);
38 TH2F *hsectorVSphi = new TH2F("hsectorVSphi", "hsectorVSphi", 360,0.,360.,18,0,18);
39 TH2F *hstripVStheta = new TH2F("hstripVStheta", "hstripVStheta", 180,0.,180.,20,0,20);
40 //TH2F *hpadzVStheta = new TH2F("hpadzVStheta", "hpadzVStheta", 180,0.,180.,2,0,2);
41 TH2F *hpadxVSphi = new TH2F("hpadxVSphi", "hpadxVSphi", 360,0.,360.,48,0,48);
42 TH2F *hpadz2stripVStheta = new TH2F("hpadz2stripVStheta", "hpadz2stripVStheta", 180,0.,180.,40,0,40);
43 //TH2F *hdzVSpadz2strip = new TH2F("hdzVSpadz2strip", "hdzVSpadz2strip",40,0,40,70,-3.5,3.5);
32bd3eb5 44
45 // Dynamically link some shared libs
46 if (gClassTable->GetID("AliRun") < 0) {
47 gROOT->LoadMacro("loadlibs.C");
48 loadlibs();
49 }
50
51 if (gAlice)
52 {
53 delete gAlice->GetRunLoader();
54 delete gAlice;
55 gAlice = 0x0;
56 }
57
58 AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::fgkDefaultEventFolderName,"read");
59 if (!rl)
60 {
61 cerr<<"Can't load RunLoader from file"<<"!\n";
62 rc = 1;
63 return rc;
64 }
65
66 rl->LoadgAlice();
67 gAlice = rl->GetAliRun();
68
69 if (!gAlice)
70 {
71 cerr << "<AliTOFanalyzeHits> AliRun object not found on file\n ";
72 rc = 2;
73 return rc;
74 }
75
76 // Get the pointer to the TOF detector
77 AliTOF *tof = (AliTOF *) gAlice->GetDetector("TOF");
78 AliLoader *tofl = rl->GetLoader("TOFLoader");
79 if (tof == 0x0 || tofl == 0x0) {
80 cerr << "<AliTOFanalyzeHits> Can not find TOF or TOFLoader\n";
81 rc = 3;
42036c57 82 return rc;
83 }
84
85 Int_t countHits = 0;
32bd3eb5 86
87 rl->LoadHeader();
88
89 for (Int_t ievent=0; ievent<nevents; ievent++) {
90 printf ("Processing event %d \n", ievent);
91 rl->GetEvent(ievent);
92
93 // Get the pointer Hit tree
94 tofl->LoadHits();
95 TTree *hitTree = tofl->TreeH();
96 tof->SetTreeAddress();
97 if (!hitTree) {
98 cout << "<AliTOFanalyzeHits> No TreeH found" << endl;
99 rc = 4;
100 return rc;
101 }
102
103 rl->LoadKinematics();
104 //AliStack* stack = rl->Stack(); // it is not necessary to use the stack!
105
106 // Get the number of entries in the hit tree
107 // (Number of primary particles creating a hit somewhere)
108 Int_t nTrack = (Int_t) hitTree->GetEntries();
109 cout << "<AliTOFanalyzeHits> Found " << nTrack
110 << " primary particles with hits" << endl;
111
112 Int_t nPrimaryOnTof = 0;
113 Int_t nSecondaryOnTof = 0;
114 Int_t nelectron = 0;
115 Int_t npion = 0;
116 Int_t nkaon = 0;
117 Int_t nproton = 0;
118 Int_t nmuon = 0;
42036c57 119
32bd3eb5 120 // Loop through all entries in the tree
121 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
122
123 tof->ResetHits();
124 hitTree->GetEvent(iTrack);
125
126 // Loop through the TOF hits
127 AliTOFhitT0 *hit = (AliTOFhitT0 *) tof->FirstHit(-1);
128 while (hit) {
129
130 countHits++;
131
132 Float_t x = hit->X();
133 Float_t y = hit->Y();
134 Float_t z = hit->Z();
135 Float_t phiAngle=TMath::Pi() + TMath::ATan2(-y,-x);
136 Float_t rhoRadius=TMath::Sqrt(x*x+y*y);
137 Float_t thetaAngle=TMath::Pi() + TMath::ATan2(-rhoRadius,-z);
138 Float_t dummy=rhoRadius*phiAngle;
139 h2hitMap->Fill(dummy,z);
140
141 phiAngle*=180./TMath::Pi();
142 thetaAngle*=180./TMath::Pi();
143
144 Float_t flightTime = hit->GetTof(); // [s]
145 flightTime *= 1.e+09; // convert in [ns]
146 Float_t angle = hit->GetIncA();
147 Float_t tofmom = hit->GetMom(); // [GeV/c]
148 Float_t trackLen = hit->GetLen(); // [cm]
149
150 // TOF hit volumes
151 Int_t sector = hit->GetSector(); // range [1-18]
152 Int_t plate = hit->GetPlate(); // range [1- 5]
153 Int_t strip = hit->GetStrip(); // range [1-20]
154 Int_t padz = hit->GetPadz(); // range [1- 2]
155 Int_t padx = hit->GetPadx(); // range [1-48]
156 // it is QA, then I perform QA!
32bd3eb5 157
da3d3acd 158
159 Bool_t isHitBad = (sector<0 || sector>17 ||
160 plate<0 || plate>4 ||
161 padz<0 || padz>1 ||
162 padx<0 || padx>47 ||
163
164 ((strip<0 || strip>14) && plate == 2) ||
165 ((strip<0 || strip>18) && (plate == 1 || plate == 3)) ||
166 ((strip<0 || strip>19) && (plate == 0 || plate == 4)));
167
168 if (isHitBad) {
169 cout << "<AliTOFanalyzeHits> strange hit found \n";
170 cout << "sector = " << sector <<
171 " plate = " << plate <<
172 " strip = " << strip <<
173 " padx = " << padx <<
174 " padz = " << padz << endl;
32bd3eb5 175 rc = 5;
176 return rc;
177 }
da3d3acd 178
32bd3eb5 179 hmoduleVStheta->Fill(thetaAngle,plate);
180 hstripVStheta->Fill(thetaAngle,strip);
181 hsectorVSphi->Fill(phiAngle,sector);
da3d3acd 182 //hpadzVStheta->Fill(thetaAngle,padx);
32bd3eb5 183 hpadxVSphi->Fill(phiAngle,padx);
da3d3acd 184
185 Float_t dummy2 = 2*strip + padz;
186 hpadz2stripVStheta->Fill(thetaAngle,dummy2);
187
188 /*
189 Float_t dummy3;
190 if (hit->GetDz()<=0) dummy3 = hit->GetDz()-1.75;
191 else dummy3 = hit->GetDz()+1.75;
192 hdzVSpadz2strip->Fill(dummy2,dummy3);
193 */
32bd3eb5 194
195 // filling hit volume histos
196 hsector->Fill(sector);
197 hplate->Fill(plate);
198 hstrip->Fill(strip);
199 hpadx->Fill(padx);
200 hpadz->Fill(padz);
201
202 Int_t track = hit->Track();
203 TParticle *part = gAlice->Particle(track);
204
205 // getting MC info for the current track
206 if (part->GetFirstMother()<0){
207 Int_t icod = TMath::Abs(part->GetPdgCode());
208 switch (icod) {
209 case 211:
210 npion++;
211 break ;
212 case 321:
213 nkaon++;
214 break ;
215 case 2212:
216 nproton++;
217 break ;
218 case 11:
219 nelectron++;
220 break ;
221 case 13:
222 nmuon++;
223 break ;
224 }
225 htofp->Fill(flightTime);
226 htofmomp->Fill(tofmom);
227 htrackLenp->Fill(trackLen);
228 } else {
229 htofs->Fill(flightTime);
230 htofmoms->Fill(tofmom);
231 }
232
233 // go to next hit
234 hit = (AliTOFhitT0 *) tof->NextHit();
42036c57 235
32bd3eb5 236 }
42036c57 237
238 }
239
32bd3eb5 240 tofl->UnloadHits();
241 rl->UnloadKinematics();
42036c57 242
32bd3eb5 243 cout << "<AliTOFanalyzeHits> Found " << countHits << " hits in total" << endl;
244 cout << npion << " primary pions reached the TOF detector" << endl;
245 cout << nkaon << " primary kaons reached the TOF detector" << endl;
246 cout << nproton << " primary protons reached the TOF detector" << endl;
247 cout << nelectron << " primary electrons reached the TOF detector" << endl;
248 cout << nmuon << " primary muons reached the TOF detector" << endl;
42036c57 249
32bd3eb5 250 }
42036c57 251
da3d3acd 252 cout << "hpadx->GetEntries() = " << hpadx->GetEntries() << endl;
253
32bd3eb5 254 rl->UnloadHeader();
255 rl->UnloadgAlice();
256
257 if (drawing) {
258 TCanvas *cHits = new TCanvas("cHits","AliTOFanalyzeHits hit volumes",50,50,900,900);
259 cHits->Divide(3,2);
260 cHits->cd(1);
261 hsector->Draw();
262 cHits->cd(2);
263 hplate->Draw();
264 cHits->cd(3);
265 hstrip->Draw();
266 cHits->cd(4);
267 hpadz->Draw();
268 cHits->cd(5);
269 hpadx->Draw();
270
271 TCanvas *chitmap = new TCanvas("chitmap","AliTOFanalyzeHits Hit Map",50,50,600,600);
272 chitmap->cd();
273 h2hitMap->Draw();
274
275 TCanvas *ctrackLen = new TCanvas("ctrackLen","AliTOFanalyzeHits Track Length for primaries on TOF",50,50,400,400);
276 ctrackLen->cd();
277 htrackLenp->Draw();
278
279 TCanvas *ctofmom = new TCanvas("ctofmom","AliTOFanalyzeHits flight times",50,50,700,700);
280 ctofmom->Divide(2,2);
281 ctofmom->cd(1);
282 gPad->SetLogy();
283 htofp->Draw();
284 ctofmom->cd(2);
285 gPad->SetLogy();
286 htofs->Draw();
287 ctofmom->cd(3);
288 gPad->SetLogy();
289 htofmomp->Draw();
290 ctofmom->cd(4);
291 gPad->SetLogy();
292 htofmoms->Draw();
293 }
294
42036c57 295 // save histos into file TOF_hitsQA.root
296 TFile *fout = new TFile("TOF_hitsQA.root","RECREATE");
297 h2hitMap->Write();
298 htofp->Write();
299 htofs->Write();
300 htofmomp->Write();
301 htofmoms->Write();
302 hsector->Write();
303 hplate->Write();
304 hstrip->Write();
305 hpadz->Write();
306 hpadx->Write();
307 htrackLenp->Write();
32bd3eb5 308
309 hmoduleVStheta->Write();
310 hsectorVSphi->Write();
311 hstripVStheta->Write();
da3d3acd 312 //hpadzVStheta->Write();
32bd3eb5 313 hpadxVSphi->Write();
da3d3acd 314 hpadz2stripVStheta->Write();
315 //hdzVSpadz2strip->Write();
32bd3eb5 316
42036c57 317 fout->Close();
318
32bd3eb5 319 cout << " Finished AliTOFanalizeHits \n";
320
321 if (gAlice)
322 {
323 delete gAlice->GetRunLoader();
324 delete gAlice;
325 gAlice = 0x0;
326 }
327
42036c57 328 return rc;
329
330}