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