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