2557f333e1588bba3f0823e540a04a27e28d16ad
[u/mrichter/AliRoot.git] / PWG1 / TOF / AliAnalysisTaskTOFqa.cxx
1 /*  created by fbellini@cern.ch on 14/09/2010 */
2 /*  last modified by fbellini   on 21/10/2010 */
3
4
5 #ifndef ALIANALYSISTASKTOFQA_CXX
6 #define ALIANALYSISTASKTOFQA_CXX
7
8 #include "TChain.h"
9 #include "TTree.h"
10 #include "TH1F.h"
11 #include "TH2F.h"
12 #include "TCanvas.h"
13 #include "AliAnalysisTaskSE.h"
14 #include "AliAnalysisManager.h"
15 #include "AliESDEvent.h"
16 #include "AliESDInputHandler.h"
17 #include "AliAnalysisTaskTOFqa.h"
18 #include "AliAnalysisFilter.h"
19 #include "AliESDtrackCuts.h"
20 #include "AliLog.h"
21 #include "AliTOFRawStream.h"
22 #include "AliTOFGeometry.h"
23
24 ClassImp(AliAnalysisTaskTOFqa)
25
26 //________________________________________________________________________
27 AliAnalysisTaskTOFqa::AliAnalysisTaskTOFqa() :
28   fRunNumber(0), 
29   fESD(0x0), 
30   fTrackFilter(0x0), 
31   fVertex(0x0),
32   fNTOFtracks(0), 
33   fNPrimaryTracks(0), 
34   fT0(0), 
35   fHlist(0x0),
36   fHlistExperts(0x0)
37  {
38   // Default constructor
39 }
40 //________________________________________________________________________
41 AliAnalysisTaskTOFqa::AliAnalysisTaskTOFqa(const char *name) : 
42   AliAnalysisTaskSE(name), 
43   fRunNumber(0), 
44   fESD(0x0), 
45   fTrackFilter(0x0),
46   fVertex(0x0),
47   fNTOFtracks(0), 
48   fNPrimaryTracks(0), 
49   fT0(0), 
50   fHlist(0x0),
51   fHlistExperts(0)
52  {
53   // Constructor
54   // Define input and output slots here
55    Info("AliAnalysisTaskTOFqa","Calling Constructor");
56    
57    // Input slot #0 works with a TChain
58    DefineInput(0, TChain::Class());
59    
60    // Output slot #0 writes into a TH1 container
61    // Output slot #1 writes into a user defined  container
62    DefineOutput(1, TList::Class());
63    DefineOutput(2, TList::Class());
64  }
65
66 //________________________________________________________________________
67 AliAnalysisTaskTOFqa::AliAnalysisTaskTOFqa(const AliAnalysisTaskTOFqa& copy) 
68 : AliAnalysisTaskSE(), 
69   fRunNumber(copy.fRunNumber), 
70   fESD(copy.fESD), 
71   fTrackFilter(copy.fTrackFilter), 
72   fVertex(copy.fVertex),
73   fNTOFtracks(copy.fNTOFtracks), 
74   fNPrimaryTracks(copy.fNPrimaryTracks), 
75   fT0(copy.fT0),
76   fHlist(copy.fHlist),
77   fHlistExperts(copy.fHlistExperts)
78 {
79   // Copy constructor
80 }
81
82 //___________________________________________________________________________
83 AliAnalysisTaskTOFqa& AliAnalysisTaskTOFqa::operator=(const AliAnalysisTaskTOFqa& copy) 
84 {
85   //
86   // Assignment operator
87   //
88   if (this!=&copy) {
89     AliAnalysisTaskSE::operator=(copy) ;
90     fRunNumber=copy.fRunNumber; 
91     fESD=copy.fESD;
92     fTrackFilter=copy.fTrackFilter;
93     fVertex=copy.fVertex;
94     fNTOFtracks=copy.fNTOFtracks; 
95     fNPrimaryTracks=copy.fNPrimaryTracks; 
96     fT0=copy.fT0;
97     fHlist=copy.fHlist;
98     fHlistExperts=copy.fHlistExperts;
99   }
100   return *this;
101 }
102
103 //___________________________________________________________________________
104 AliAnalysisTaskTOFqa::~AliAnalysisTaskTOFqa() {
105   //
106   //destructor
107   //
108
109   Info("~AliAnalysisTaskTOFqa","Calling Destructor");
110   if (fVertex) delete fVertex;
111   if (fTrackFilter) delete fTrackFilter;
112   if (fHlist) {
113     delete fHlist;
114     fHlist = 0;
115   }
116   if (fHlistExperts) {
117     delete fHlistExperts;
118     fHlistExperts = 0;
119   }
120 }
121
122 //________________________________________________________________________
123 void AliAnalysisTaskTOFqa::UserCreateOutputObjects()
124 {
125   //Defines output objects and histograms
126   Info("CreateOutputObjects","CreateOutputObjects (TList) of task %s", GetName());
127   OpenFile(1);
128   if (!fHlist) fHlist = new TList();    
129   fHlist->SetOwner(kTRUE);
130   if (!fHlistExperts) fHlistExperts = new TList();      
131   fHlistExperts->SetOwner(kTRUE);
132   //0
133   TH1I* hTOFmatchedESDperEvt = new TH1I("hTOFmatchedPerEvt", "Number of matched TOF tracks per event;Number of TOF matched ESD tracks;Counts", 100, 0, 100) ;  
134   hTOFmatchedESDperEvt->SetLineWidth(2);
135   hTOFmatchedESDperEvt->SetLineColor(kBlue);
136   fHlist->AddLast(hTOFmatchedESDperEvt) ;
137   //1
138   TH1F* hTOFmatchedESDtime = new TH1F("hTOFmatchedESDtime", "Matched  ESDs tracks: TOF Time spectrum; t [ns];Counts", 250, 0., 610. ) ; 
139   hTOFmatchedESDtime->SetLineWidth(2);
140   hTOFmatchedESDtime->SetLineColor(kBlue);
141   hTOFmatchedESDtime->SetFillColor(kBlue);
142   fHlist->AddLast(hTOFmatchedESDtime) ;
143   //2
144   TH1F* hTOFmatchedESDrawTime = new TH1F("hTOFmatchedESDrawTime", "Matched ESDs tracks: TOF raw Time spectrum;t_{raw} [ns];Counts", 250, 0., 610.) ; 
145   hTOFmatchedESDrawTime->SetLineWidth(2);
146   hTOFmatchedESDrawTime->SetLineColor(kGreen+2);
147   hTOFmatchedESDrawTime->SetFillColor(kGreen+2);
148   fHlist->AddLast(hTOFmatchedESDrawTime) ;
149   //3
150   TH1F* hTOFmatchedESDToT = new TH1F("hTOFmatchedESDToT", "Matched ESDs tracks: TOF ToT spectrum; ToT [ns];Counts",100, 0., 48.8) ; 
151   hTOFmatchedESDToT->SetLineWidth(2);
152   hTOFmatchedESDToT->SetLineColor(kBlue);
153   hTOFmatchedESDToT->SetFillColor(kBlue);
154   fHlist->AddLast(hTOFmatchedESDToT) ;
155   //4
156   TH1F* hTOFmatchedESDP  = new TH1F("hTOFmatchedESDP", "TPC-TOF matched tracks momentum distribution (GeV/c); P(GeV/c);tracks", 500,0.,5.) ;  
157   hTOFmatchedESDP->SetLineWidth(1);
158   hTOFmatchedESDP->SetLineColor(kBlue);
159   hTOFmatchedESDP->SetMarkerStyle(21);
160   hTOFmatchedESDP->SetMarkerSize(0.6);
161   hTOFmatchedESDP->SetMarkerColor(kBlue);
162   fHlist->AddLast(hTOFmatchedESDP) ; 
163   //5
164   TH1F* hTOFmatchedESDPt  = new TH1F("hTOFmatchedESDPt", "TPC-TOF matched tracks p_{T} distribution (GeV/c); p_{T}(GeV/c);tracks", 500,0.,5.) ;  
165   hTOFmatchedESDPt->SetLineWidth(1);
166   hTOFmatchedESDPt->SetLineColor(kBlue);
167   hTOFmatchedESDPt->SetMarkerStyle(4);
168   hTOFmatchedESDPt->SetMarkerSize(0.6);
169   hTOFmatchedESDPt->SetMarkerColor(kBlue);
170   fHlist->AddLast(hTOFmatchedESDPt) ; 
171   //6
172   TH1F* hTOFmatchedESDtrkLength  = new TH1F("hTOFmatchedESDtrkLength", "Matched ESDs tracks length; Track length [cm];Counts", 1600, -800., 800) ; 
173   hTOFmatchedESDtrkLength->SetLineWidth(1);
174   hTOFmatchedESDtrkLength->SetLineColor(kBlue);
175   fHlist->AddLast(hTOFmatchedESDtrkLength) ; 
176   //7
177   TH2F* hTOFmatchedESDpVsBeta  = new TH2F("hTOFmatchedESDpVsBeta", "Matched ESDs tracks p vs. beta; p(GeV/c); beta", 500, 0., 5.,500, 0., 5.) ; 
178   fHlist->AddLast(hTOFmatchedESDpVsBeta);
179   
180   //8
181   TH1F* hESDprimaryTrackP = new TH1F("hESDprimaryTrackP", "All ESDs tracks Pt distribution (GeV/c); p_{T}(GeV/c);tracks", 500, 0., 5.0) ;  
182   hESDprimaryTrackP->SetLineWidth(1);
183   hESDprimaryTrackP->SetMarkerStyle(8);
184   hESDprimaryTrackP->SetMarkerSize(0.6);
185   hESDprimaryTrackP->SetLineColor(kGray+1);
186   fHlist->AddLast(hESDprimaryTrackP);
187   //9
188   TH1F* hESDprimaryTrackPt = new TH1F("hESDprimaryTrackPt", "ESDs primary tracks Pt distribution (GeV/c); p_{T}(GeV/c);tracks", 500, 0.0, 5.0) ;  
189   hESDprimaryTrackPt->SetLineWidth(1);
190   hESDprimaryTrackPt->SetMarkerStyle(4);
191   hESDprimaryTrackPt->SetMarkerSize(0.6);
192   hESDprimaryTrackPt->SetLineColor(kBlack);
193   fHlist->AddLast(hESDprimaryTrackPt);
194   //10
195   TH1F* hTOFmatchedESDeta = new TH1F("hTOFmatchedESDeta", "Matched ESDtracks eta; eta;Counts", 500, -2.5, 2.5) ; 
196   fHlist->AddLast(hTOFmatchedESDeta) ; 
197   //11
198   TH1F* hTOFprimaryESDeta = new TH1F("hTOFprimaryESDeta", "Primary ESDtracks eta; eta;Counts", 500, -2.5, 2.5) ; 
199   fHlist->AddLast(hTOFprimaryESDeta) ; 
200   //12
201   TH1F* hTOFmatchedESDphi = new TH1F("hTOFmatchedESDphi", "Matched ESDtracks Phi; Phi (rad);Counts", 65, 0., 6.5) ; 
202   fHlist->AddLast(hTOFmatchedESDphi) ; 
203   //13
204   TH1F* hTOFprimaryESDphi = new TH1F("hTOFprimaryESDphi", "Primary ESDtracks Phi;Phi (rad);Counts", 65, 0., 6.5) ; 
205   fHlist->AddLast(hTOFprimaryESDphi) ; 
206   
207   //Experts 0
208   TH1F* hTOFESDsMatchingProb  = new TH1F("hTOFESDsMatchingProb", "TPC-TOF track-matching probability per event(|eta|<0.9 && pt>0.5GeV/c);TPC-TOF track-matching probability (%)  ;Counts",21, 0., 110.) ;  
209   hTOFESDsMatchingProb->SetLineColor(kRed);
210   fHlistExperts->AddLast(hTOFESDsMatchingProb) ; 
211   
212   //Experts 1
213   TH1F* hTOFmatchedESDdiffTime  = new TH1F("hTOFmatchedESDdiffTime", "ESDs t_{TOF}-t_{pi,exp} spectrum in TOF (ps); t_{TOF}-t_{pi,exp} [ps];Counts", 4000, -50000., 50000.) ; 
214   hTOFmatchedESDdiffTime->SetLineWidth(1);
215   hTOFmatchedESDdiffTime->SetLineColor(kBlack);
216   hTOFmatchedESDdiffTime->SetMarkerStyle(8);
217   hTOFmatchedESDdiffTime->SetMarkerSize(0.8);
218   hTOFmatchedESDdiffTime->SetMarkerColor(kAzure+7);
219   hTOFmatchedESDdiffTime->SetFillColor(kAzure-2);
220   fHlistExperts->AddLast(hTOFmatchedESDdiffTime) ; 
221   
222   //Experts 2
223   TH2F* hTOFmatchedESDdiffTimeVsStrip = new TH2F("hTOFmatchedESDdiffTimeVsStrip", "ESDs t_{TOF}-t_{pi,exp} vs strip number; strip (Eta); t_{TOF}-t_{pi,exp} [ps]", 92,0.,92,400, -5000., 5000.) ; 
224   fHlistExperts->AddLast(hTOFmatchedESDdiffTimeVsStrip) ; 
225   
226   //Experts 3
227   TH2F* hTOFmatchedESDdxVsEta = new TH2F("hTOFmatchedESDdxVsEta", "Matched ESD tracks Dx vs eta; strip(eta); Dx (cm)", 92,0.,92., 200,-10.,10.) ; 
228   fHlistExperts->AddLast(hTOFmatchedESDdxVsEta) ; 
229   
230   //Experts 4
231   TH2F* hTOFmatchedESDdzVsEta  = new TH2F("hTOFmatchedESDdzVsEta", "Matched ESDtracks Dz vs eta; strip(eta); Dz (cm)", 92,0.,92., 200,-10.,10.) ; 
232   fHlistExperts->AddLast(hTOFmatchedESDdzVsEta) ; 
233   
234   //Experts 5
235   TH1F* hTOFmatchedMass= new TH1F("hTOFmatchedMass","Matched ESD tracks mass distribution; M (GeV/c^{2}); entries", 500,0., 5. );
236   hTOFmatchedMass->SetLineWidth(2);
237   hTOFmatchedMass->SetLineColor(kBlue);
238   hTOFmatchedMass->SetLineColor(kBlue);
239   fHlistExperts->AddLast(hTOFmatchedMass);
240    
241   //Experts 6
242   TH1D* hEventT0DetAND = new TH1D("hEventT0DetAND", "Event T0 from T0 detector (A&C); t [ps];Counts", 4000, -50000., 50000. ) ; 
243   hEventT0DetAND->SetLineWidth(2);
244   hEventT0DetAND->SetLineColor(kRed);
245   hEventT0DetAND->SetFillColor(kRed);
246   fHlistExperts->AddLast(hEventT0DetAND) ;
247
248   //Experts 7
249   TH1D* hEventT0DetA = new TH1D("hEventT0DetA", "Event T0 from T0 detector (A side); t [ps];Counts", 4000, -50000., 50000. ) ; 
250   hEventT0DetA->SetLineWidth(2);
251   hEventT0DetA->SetLineColor(kBlue);
252   hEventT0DetA->SetFillColor(kBlue);
253   fHlistExperts->AddLast(hEventT0DetA) ;
254
255    //Experts 8
256   TH1D* hEventT0DetC = new TH1D("hEventT0DetC", "Event T0 from T0 detector (C side); t [ps];Counts", 4000, -50000., 50000. ) ; 
257   hEventT0DetC->SetLineWidth(2);
258   hEventT0DetC->SetLineColor(kGreen);
259   hEventT0DetC->SetFillColor(kGreen);
260   fHlistExperts->AddLast(hEventT0DetC);
261  
262   //Experts 9
263   TH1F* hTOFmatchedExpTime = new TH1F("hTOFmatchedExpTime", "Matched  ESDs tracks - pions expected  time; t [ns];Counts", 4000, -50000., 50000. ) ; 
264   hTOFmatchedExpTime->SetLineWidth(1);
265   hTOFmatchedExpTime->SetLineColor(kBlack);
266   hTOFmatchedExpTime->SetMarkerStyle(8);
267   hTOFmatchedExpTime->SetMarkerSize(0.8); 
268   hTOFmatchedExpTime->SetMarkerColor(kRed);
269   fHlistExperts->AddLast(hTOFmatchedExpTime) ;
270
271   
272   PostData(1, fHlist);
273   PostData(2, fHlistExperts);
274 }
275 //________________________________________________________________________
276 void AliAnalysisTaskTOFqa::UserExec(Option_t *) 
277
278   /* Main - executed for each event.
279     It extracts event information and track information after selecting 
280     primary tracks via standard cuts. */
281   const Double_t speedOfLight =  TMath::C()*1E2*1E-12; // cm/ps
282   
283   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
284   if (!esdH) {
285     Printf("ERROR: Could not get ESDInputHandler");
286     return;
287   } else {
288     fESD = (AliESDEvent*) esdH->GetEvent();
289   } 
290   
291   if (!fESD) {
292     Printf("ERROR: fESD not available");
293     return;
294   }
295   
296   // loop over ESD tracks
297   fNTOFtracks=0;
298   fNPrimaryTracks=0;
299   
300   //info from To detector
301   fT0= fESD->GetT0TOF(0);//ps
302   ((TH1D*)fHlistExperts->At(6))->Fill(fT0);//ps
303   ((TH1D*)fHlistExperts->At(7))->Fill((Double_t)fESD->GetT0TOF(1)); //ps
304   ((TH1D*)fHlistExperts->At(8))->Fill((Double_t)fESD->GetT0TOF(2));//ps
305       
306   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
307     AliESDtrack* track = fESD->GetTrack(iTracks);
308     if (!track) {
309       Printf("ERROR: Could not receive track %d", iTracks);
310       continue;
311     }
312
313     //primary tracks selection: kTPCrefit and std cuts
314     if(!fTrackFilter->IsSelected(track)) continue;
315    
316     Double_t tofTime=track->GetTOFsignal();//in ps
317     Double_t tofTimeRaw=track->GetTOFsignalRaw();//in ps
318     Double_t tofToT=track->GetTOFsignalToT(); //in ps
319     Double_t expTimes[5];
320     track->GetIntegratedTimes(expTimes);
321     Double_t length =track->GetIntegratedLength();
322     Double_t ptot[3];
323     track->GetConstrainedPxPyPz(ptot);    
324     Double_t pT = TMath::Sqrt(ptot[0]*ptot[0] + ptot[1]*ptot[1]);
325     Double_t P2 = pT*pT + ptot[2]*ptot[2];
326     Double_t eta=track->Eta();
327     Int_t channel=track->GetTOFCalChannel(); 
328     Int_t volId[5]; //(sector, plate,strip,padZ,padX)
329     AliTOFGeometry::GetVolumeIndices(channel,volId);
330     
331     if (TMath::Abs(eta)<0.9) { //cut for acceptance
332       if (P2>=0)
333         ((TH1F*)fHlist->At(8))->Fill(TMath::Sqrt(P2)); 
334       ((TH1F*)fHlist->At(9))->Fill(track->Pt()); //all esd tracks within acceptance
335       ((TH1F*)fHlist->At(11))->Fill(eta);
336       ((TH1F*)fHlist->At(13))->Fill(track->Phi());
337       if ((TMath::Abs(track->Eta())<0.9)&&(track->Pt()>0.5)) fNPrimaryTracks++;
338       
339       //matched tracks selection: kTOFout and kTIME
340       if ((track->GetStatus() & AliESDtrack::kTOFout) &&
341           (track->GetStatus() & AliESDtrack::kTIME)) {      
342         
343         if ((TMath::Abs(track->Eta())<0.9)&&(track->Pt()>0.5)) fNTOFtracks++; //matched counter
344         
345           ((TH1F*)fHlist->At(1))->Fill(tofTime*1E-3); //ns
346           ((TH1F*)fHlist->At(2))->Fill(tofTimeRaw*1E-3); //ns
347           ((TH1F*)fHlist->At(3))->Fill(tofToT);
348           if (P2>=0)
349             ((TH1F*)fHlist->At(4))->Fill(TMath::Sqrt(P2));
350           ((TH1F*)fHlist->At(5))->Fill(track->Pt());
351           ((TH1F*)fHlist->At(6))->Fill(length);  
352           ((TH1F*)fHlist->At(10))->Fill(eta);
353           ((TH1F*)fHlist->At(12))->Fill(track->Phi());
354           
355           ((TH1F*)fHlistExperts->At(1))->Fill(tofTime-expTimes[2]);//ps
356           ((TH1F*)fHlistExperts->At(2))->Fill((Int_t)GetStripIndex(volId),tofTime-expTimes[2]); //ps
357           ((TH1F*)fHlistExperts->At(3))->Fill((Int_t)GetStripIndex(volId),track->GetTOFsignalDx());
358           ((TH1F*)fHlistExperts->At(4))->Fill((Int_t)GetStripIndex(volId),track->GetTOFsignalDz());
359           ((TH1F*)fHlistExperts->At(9))->Fill(expTimes[2]);//ps
360
361           //basic PID performance check 
362           Double_t tof= tofTime; //average T0 fill subtracted, no info from T0detector 
363           if (length>350){
364             Double_t beta= length/(tof*speedOfLight);
365             //Double_t mass2=P2*((tof/length)*(tof/length)-(1/(speedOfLight*speedOfLight)));
366             ((TH1F*)fHlist->At(7))->Fill(TMath::Sqrt(P2),beta);
367             //if (mass2>=0)((TH1F*)fHlistExperts->At(5))->Fill(TMath::Sqrt(mass2));
368           }
369       }//matched
370     }//acceptance cut
371   }//end loop on tracks
372   ((TH1F*)fHlist->At(0))->Fill(fNTOFtracks) ;
373   //if (fNTOFtracks>fNPrimaryTracks) printf("Something strange!!!\n");
374   if(fNPrimaryTracks>0){
375     ((TH1F*)fHlistExperts->At(0))->Fill((fNTOFtracks/(Float_t)fNPrimaryTracks)*100) ;
376   }
377   
378  
379   PostData(1, fHlist);
380   PostData(2, fHlistExperts);
381   
382 }      
383
384 //________________________________________________________________________
385 void AliAnalysisTaskTOFqa::Terminate(Option_t *) 
386 {
387   //check on output validity
388   fHlist = dynamic_cast<TList*> (GetOutputData(1));
389   if (!fHlist || !fHlistExperts) {
390     Printf("ERROR: lists not available");
391     return;   
392   }   
393  
394 }
395
396 //---------------------------------------------------------------
397 Int_t AliAnalysisTaskTOFqa::GetStripIndex(const Int_t * const in)
398 {
399   /* return tof strip index between 0 and 91 */
400   
401   Int_t nStripA = AliTOFGeometry::NStripA();
402   Int_t nStripB = AliTOFGeometry::NStripB();
403   Int_t nStripC = AliTOFGeometry::NStripC();
404
405   Int_t iplate = in[1];
406   Int_t istrip = in[2];
407   
408   Int_t stripOffset = 0;
409   switch (iplate) {
410   case 0:
411     stripOffset = 0;
412       break;
413   case 1:
414     stripOffset = nStripC;
415     break;
416   case 2:
417     stripOffset = nStripC+nStripB;
418     break;
419   case 3:
420     stripOffset = nStripC+nStripB+nStripA;
421     break;
422   case 4:
423     stripOffset = nStripC+nStripB+nStripA+nStripB;
424     break;
425   default:
426     stripOffset=-1;
427     break;
428   };
429   
430   if (stripOffset<0 || stripOffset>92) return -1;
431   else 
432     return (stripOffset+istrip);
433 }
434
435 #endif