adding proper destructors and some small changes in analysis macro
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliSpectraBothTrackCuts.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 //-----------------------------------------------------------------
18 //         AliSpectraBothTrackCuts class
19 //-----------------------------------------------------------------
20
21 #include "TChain.h"
22 #include "TTree.h"
23 #include "TLegend.h"
24 #include "TH1F.h"
25 #include "TH1I.h"
26 #include "TH2F.h"
27 #include "TCanvas.h"
28 #include "AliAnalysisTask.h"
29 #include "AliAnalysisManager.h"
30 #include "AliAODTrack.h"
31 #include "AliVTrack.h"
32 #include "AliExternalTrackParam.h"
33 #include "AliAODMCParticle.h"
34 #include "AliAODEvent.h"
35 #include "AliAODInputHandler.h"
36 #include "AliAnalysisTaskESDfilter.h"
37 #include "AliAnalysisDataContainer.h"
38 #include "AliSpectraBothTrackCuts.h"
39 //#include "AliSpectraBothHistoManager.h"
40 #include <iostream>
41
42 using namespace std;
43
44 const char * AliSpectraBothTrackCuts::kBinLabel[] ={"TrkBit",
45                                                    "TrkCuts",
46                                                    "TrkEta",
47                                                    "TrkDCA",
48                                                    "TrkP",
49                                                    "TrkPt",
50                                                    "TrkPtTOF",
51                                                    "TOFMatching",
52                                                    "kTOFout",
53                                                    "kTIME",
54                                                    "kTOFpid",
55                                                    "Accepted"};
56
57
58 ClassImp(AliSpectraBothTrackCuts)
59
60
61 AliSpectraBothTrackCuts::AliSpectraBothTrackCuts(const char *name) : TNamed(name, "AOD Track Cuts"), fIsSelected(0), fTrackBits(0), fMinTPCcls(0), fEtaCutMin(0), fEtaCutMax(0), fDCACut(0), fPCut(0), fPtCut(0), fYCutMax(0),fYCutMin(0),
62   fPtCutTOFMatching(0),fAODtrack(kotherobject), fHashitinSPD1(0),fusedadditionalcuts(kTRUE),
63 fHistoCuts(0), fHistoNSelectedPos(0), fHistoNSelectedNeg(0), fHistoNMatchedPos(0), fHistoNMatchedNeg(0), fHistoEtaPhiHighPt(0), fHistoNclustersITS(0),fTrack(0),fCuts(0)
64   
65 {
66   Bool_t oldStatus = TH1::AddDirectoryStatus();
67   TH1::AddDirectory(kFALSE);    
68   // Constructor
69   fHistoCuts = new TH1I("fTrkCuts", "Track Cuts", kNTrkCuts, -0.5, kNTrkCuts - 0.5);
70   for(Int_t ibin=1;ibin<=kNTrkCuts;ibin++)fHistoCuts->GetXaxis()->SetBinLabel(ibin,kBinLabel[ibin-1]);
71   //standard histo
72   const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
73   Int_t nbinsTempl=52;
74   
75   fHistoNSelectedPos=new TH1F("fHistoNSelectedPos","fHistoNSelectedPos",nbinsTempl,templBins);
76   fHistoNSelectedPos->GetXaxis()->SetTitle("P_{T} (GeV / c)");
77   fHistoNSelectedNeg=new TH1F("fHistoNSelectedNeg","fHistoNSelectedNeg",nbinsTempl,templBins);
78   fHistoNSelectedNeg->GetXaxis()->SetTitle("P_{T} (GeV / c)");
79   fHistoNMatchedPos=new TH1F("fHistoNMatchedPos","fHistoNMatchedPos",nbinsTempl,templBins);
80   fHistoNMatchedPos->GetXaxis()->SetTitle("P_{T} (GeV / c)");
81   fHistoNMatchedNeg=new TH1F("fHistoNMatchedNeg","fHistoNMatchedNeg",nbinsTempl,templBins);
82   fHistoNMatchedNeg->GetXaxis()->SetTitle("P_{T} (GeV / c)");
83   fHistoEtaPhiHighPt=new TH2F("fHistoEtaPhiHighPt","fHistoEtaPhiHighPt",200,-1,1,400,0,7);
84   fHistoEtaPhiHighPt->SetXTitle("eta");
85   fHistoEtaPhiHighPt->SetYTitle("phi");
86   fHistoNclustersITS=new TH1F("fHistoNclustersITS","fHistoNclustersITS;N;ITSLayer",6,-0.5,5.5);
87   
88   fEtaCutMin = -100000.0; // default value of eta cut ~ no cut
89   fEtaCutMax = 100000.0; // default value of eta cut ~ no cut
90   fDCACut = 100000.0; // default value of dca cut ~ no cut
91   fPCut = 100000.0; // default value of p cut ~ no cut
92   fPtCut = 100000.0; // default value of pt cut ~ no cut 
93   fPtCutTOFMatching=0.6; //default value fot matching with TOF
94   fYCutMax       = 100000.0; // default value of y cut ~ no cut 
95   fYCutMin       = -100000.0; // default value of y cut ~ no cut 
96   fMinTPCcls=70; // ncls in TPC
97    TH1::AddDirectory(oldStatus);
98         
99 }
100 //_______________________________________________________
101 AliSpectraBothTrackCuts::~AliSpectraBothTrackCuts()
102 {
103         if(fHistoCuts)
104                 delete fHistoCuts;
105         if(fHistoNSelectedPos)
106                 delete fHistoNSelectedPos;
107         if(fHistoNSelectedNeg)
108                 delete fHistoNSelectedNeg;
109         if(fHistoNMatchedPos)
110                 delete fHistoNMatchedPos;
111         if(fHistoNMatchedNeg)
112                 delete fHistoNMatchedNeg;
113         if(fHistoEtaPhiHighPt)
114                 delete fHistoEtaPhiHighPt;
115         if(fHistoNclustersITS)
116                 delete fHistoNclustersITS;
117
118
119
120 }
121 //_______________________________________________________
122 Bool_t AliSpectraBothTrackCuts::IsSelected(AliVTrack * track,Bool_t FillHistStat)
123 {
124 // Returns true if Track Cuts are selected and applied
125   if (!track)
126     {
127       printf("ERROR: Could not receive track");
128       return kFALSE;
129     }
130     fTrack = track;
131    TString nameoftrack(track->ClassName());  
132     if(!nameoftrack.CompareTo("AliESDtrack"))
133                 fAODtrack=kESDobject;
134         else if(!nameoftrack.CompareTo("AliAODTrack"))
135                 fAODtrack=kAODobject;
136         else
137                 fAODtrack=kotherobject;
138   if(!CheckTrackType()){
139     return kFALSE;
140   }
141   if(FillHistStat)fHistoCuts->Fill(kTrkBit);
142   if(!CheckTrackCuts()){
143     return kFALSE;
144   }
145   if(FillHistStat)fHistoCuts->Fill(kTrkCuts);
146   if(!CheckEtaCut()){
147     return kFALSE;
148   }
149   if(FillHistStat)fHistoCuts->Fill(kTrkEta);
150   if(!CheckDCACut()){
151     return kFALSE;
152   }
153   if(FillHistStat)fHistoCuts->Fill(kTrkDCA);
154   if(!CheckPCut()){
155     return kFALSE;
156   }
157   if(FillHistStat)fHistoCuts->Fill(kTrkP);
158   if(!CheckPtCut()){
159     return kFALSE;
160   }
161   if(FillHistStat)fHistoCuts->Fill(kTrkPt);
162   if(!CheckTOFMatching(FillHistStat)){
163     return kFALSE;
164   }
165   if(FillHistStat)fHistoCuts->Fill(kAccepted);
166   //Printf("-------- %d,%d",kTOFMatching,kAccepted);
167   
168   return kTRUE;
169 }
170 //_________________________________________________________
171
172 Bool_t AliSpectraBothTrackCuts::CheckTrackType()
173 {
174   // Check track Type
175   if(fAODtrack==kESDobject)
176   {
177         AliESDtrack* esdtrack=dynamic_cast<AliESDtrack*>(fTrack);
178         if(!esdtrack)
179                 return kFALSE;  
180         if(fCuts->AcceptTrack(esdtrack)) return kTRUE;
181                 return kFALSE;
182  }
183   else if(fAODtrack==kAODobject)
184   {
185         AliAODTrack* aodtrack=dynamic_cast<AliAODTrack*>(fTrack);
186         if(!aodtrack)
187                 return kFALSE;
188         if (aodtrack->TestFilterBit(fTrackBits)) return kTRUE;
189                 return kFALSE;
190   }
191
192   else
193         return kFALSE;
194   
195 }
196 //_________________________________________________________
197
198 Bool_t AliSpectraBothTrackCuts::CheckTrackCuts()
199 {
200   // Check additional track Cuts
201   Bool_t PassTrackCuts=kTRUE;
202   if(!fusedadditionalcuts)
203         return PassTrackCuts;
204   AliAODTrack* aodtrack=0;
205   AliESDtrack* esdtrack=0;
206   if(fAODtrack==kESDobject)
207   {
208         esdtrack=dynamic_cast<AliESDtrack*>(fTrack);
209         if(!esdtrack)
210                 return kFALSE;
211         if (!esdtrack->HasPointOnITSLayer(0) && !esdtrack->HasPointOnITSLayer(1))PassTrackCuts=kFALSE; //FIXME 1 SPD for the moment
212         if (fHashitinSPD1&&!esdtrack->HasPointOnITSLayer(0)) PassTrackCuts=kFALSE;              
213         if (esdtrack->GetTPCNcls()<fMinTPCcls)PassTrackCuts=kFALSE;
214         if(!esdtrack->IsOn(AliESDtrack::kTPCrefit))PassTrackCuts=kFALSE;
215         if(!esdtrack->IsOn(AliESDtrack::kITSrefit))PassTrackCuts=kFALSE;
216         if(PassTrackCuts)
217         {
218                 for(int i=0;i<6;i++)
219                         if(esdtrack->HasPointOnITSLayer(i))
220                                 fHistoNclustersITS->Fill(i);
221         }       
222   }
223   else if (fAODtrack==kAODobject)
224         {
225         aodtrack=dynamic_cast<AliAODTrack*>(fTrack);
226         if(!aodtrack)
227                 return kFALSE;
228         if (!aodtrack->HasPointOnITSLayer(0) && !aodtrack->HasPointOnITSLayer(1))PassTrackCuts=kFALSE; //FIXME 1 SPD for the moment
229         if (fHashitinSPD1&&!aodtrack->HasPointOnITSLayer(0)) PassTrackCuts=kFALSE;      
230         if (aodtrack->GetTPCNcls()<fMinTPCcls)PassTrackCuts=kFALSE;
231         if(!aodtrack->IsOn(AliAODTrack::kTPCrefit))PassTrackCuts=kFALSE;
232         if(!aodtrack->IsOn(AliAODTrack::kITSrefit))PassTrackCuts=kFALSE;
233         if(PassTrackCuts)
234         {
235                 for(int i=0;i<6;i++)
236                         if(aodtrack->HasPointOnITSLayer(i))
237                                 fHistoNclustersITS->Fill(i);
238         }       
239   }
240   else
241         return kFALSE;
242     
243   
244   
245   return PassTrackCuts;
246 }
247 //________________________________________________________
248 Bool_t AliSpectraBothTrackCuts::CheckEtaCut()
249 {
250    // Check eta cut
251    if (fTrack->Eta() < fEtaCutMax && fTrack->Eta() > fEtaCutMin) return kTRUE;
252     return kFALSE;
253 }
254
255 Bool_t AliSpectraBothTrackCuts::CheckYCut(BothParticleSpecies_t species) 
256 {
257   // check if the rapidity is within the set range
258   Double_t y;
259   
260   Double_t pz=fTrack->Pz();
261   Double_t p=fTrack->P();
262   Double_t mass=-1.0;
263   /*
264   if (species == kSpProton) { y = fTrack->Y(9.38271999999999995e-01); }
265   if ( species == kSpKaon ) { y = fTrack->Y(4.93676999999999977e-01); }
266   if ( species == kSpPion)  { y = fTrack->Y(1.39570000000000000e-01); }
267   
268   */
269     if (species == kSpProton) { mass=9.38271999999999995e-01; }
270   if ( species == kSpKaon ) { mass=4.93676999999999977e-01; }
271   if ( species == kSpPion)  { mass=1.39570000000000000e-01; }
272   if(mass<0.0)
273         y =-999.0 ;
274   else
275         y=0.5*TMath::Log((TMath::Sqrt(mass*mass+p*p)+pz)/(TMath::Sqrt(mass*mass+p*p)-pz));
276   if (y > fYCutMax || y<fYCutMin||y < -998.) return kFALSE;
277         return kTRUE;
278 }
279 //_______________________________________________________
280 Bool_t AliSpectraBothTrackCuts::CheckDCACut()
281 {
282    // Check DCA cut
283  // if (TMath::Abs(fTrack->DCA()) < fDCACut) return kTRUE; //FIXME for newest AOD fTrack->DCA() always gives -999
284    
285     AliAODTrack* aodtrack=0;
286   AliESDtrack* esdtrack=0;
287   if(fAODtrack==kESDobject)
288   {
289         esdtrack=dynamic_cast<AliESDtrack*>(fTrack);
290         if(!esdtrack)
291                 return kFALSE;
292         Float_t dcaxy=0.0; 
293         Float_t dcaz=0.0;
294         esdtrack->GetImpactParameters(dcaxy,dcaz);      
295         if (TMath::Abs(dcaxy) < fDCACut) 
296                 return kTRUE;
297         else 
298                 return kFALSE;
299    }
300   else if (fAODtrack==kAODobject)
301   {
302         aodtrack=dynamic_cast<AliAODTrack*>(fTrack);
303         if(!aodtrack)
304                 return kFALSE;
305         if (TMath::Abs(aodtrack->DCA()) < fDCACut) return kTRUE;
306         else 
307                 return kFALSE;
308                 
309    }
310         else
311         return kFALSE;
312 }
313 //________________________________________________________
314 Bool_t AliSpectraBothTrackCuts::CheckPCut()
315 {
316    // Check P cut
317    if (fTrack->P() < fPCut) return kTRUE;
318    return kFALSE;
319 }
320 //_______________________________________________________
321 Bool_t AliSpectraBothTrackCuts::CheckPtCut()
322 {
323     // check Pt cut
324 //    if ((fTrack->Pt() < fPtCut) && (fTrack->Pt() > 0.3 )) return kTRUE;
325    if (fTrack->Pt() < fPtCut) return kTRUE;
326     return kFALSE;
327 }
328
329 //_______________________________________________________
330 Bool_t AliSpectraBothTrackCuts::CheckTOFMatching(Bool_t FillHistStat)
331 {
332   // check Pt cut
333   //    if ((fTrack->Pt() < fPtCut) && (fTrack->Pt() > 0.3 )) return kTRUE;
334   
335   if (fTrack->Pt() < fPtCutTOFMatching) return kTRUE;
336   else{
337     if(FillHistStat)fHistoCuts->Fill(kTrkPtTOF);
338     if(fTrack->Charge()>0)fHistoNSelectedPos->Fill(fTrack->Pt());
339     else fHistoNSelectedNeg->Fill(fTrack->Pt());
340     UInt_t status; 
341     status=fTrack->GetStatus();
342     if((status&AliAODTrack::kTOFout)&&FillHistStat)fHistoCuts->Fill(kTrTOFout);
343     if((status&AliAODTrack::kTIME)&&FillHistStat)fHistoCuts->Fill(kTrTIME);
344     if((status&AliAODTrack::kTOFpid)&&FillHistStat)fHistoCuts->Fill(kTrTOFpid);
345     
346     if((status&AliAODTrack::kTOFout)==0 || (status&AliAODTrack::kTIME)==0){//kTOFout and kTIME
347       return kFALSE; 
348     } 
349     if(FillHistStat)fHistoCuts->Fill(kTOFMatching);
350     if(fTrack->Charge()>0)fHistoNMatchedPos->Fill(fTrack->Pt());
351     else fHistoNMatchedNeg->Fill(fTrack->Pt());
352     if(fTrack->Pt()>1.5){
353       //fHistoEtaPhiHighPt->Fill(fTrack->GetOuterParam()->Eta(),fTrack->GetOuterParam()->Phi());
354       //Printf("AliExternalTrackParam * extpar=(AliExternalTrackParam*)fTrack->GetOuterParam();");
355       //AliExternalTrackParam * extpar=(AliExternalTrackParam*)fTrack->GetOuterParam();
356       fHistoEtaPhiHighPt->Fill(fTrack->Eta(),fTrack->Phi());
357       //Printf("fHistoEtaPhiHighPt->Fill(extpar->Eta(),extpar->Phi());");
358       //fHistoEtaPhiHighPt->Fill(extpar->Eta(),extpar->Phi());
359       //delete extpar;
360     }
361     return kTRUE;
362   }
363 }
364 //_______________________________________________________
365 void AliSpectraBothTrackCuts::PrintCuts() const
366 {
367   // Print cuts
368     cout << "Track Cuts" << endl;
369     cout << " > TrackBit\t" << fTrackBits << endl;
370     cout << " > Eta cut\t" << fEtaCutMin <<","<< fEtaCutMax << endl;
371     cout << " > DCA cut\t" << fDCACut << endl;
372     cout << " > P cut\t" << fPCut << endl;
373     cout << " > Pt cut \t" << fPtCut << endl;
374     cout << " > TPC cls \t" << fMinTPCcls << endl;
375 }
376 //_______________________________________________________
377 void AliSpectraBothTrackCuts::SetTrackType(UInt_t bit)
378 {
379    // Set the type of track to be used. The argument should be the bit number. The mask is produced automatically.
380    fTrackBits = (0x1 << (bit - 1));
381 }
382 //_______________________________________________________
383
384 Long64_t AliSpectraBothTrackCuts::Merge(TCollection* list)
385 {
386   // Merge a list of AliSpectraBothTrackCuts objects with this.
387   // Returns the number of merged objects (including this).
388
389   //  AliInfo("Merging");
390
391
392   if (!list)
393     return 0;
394
395   if (list->IsEmpty())
396     return 1;
397
398   TIterator* iter = list->MakeIterator();
399   TObject* obj;
400
401   // collections of all histograms
402   TList collections;//FIXME we should only 1 collection
403   TList collections_histoNSelectedPos;
404   TList collections_histoNSelectedNeg;
405   TList collections_histoNMatchedPos;
406   TList collections_histoNMatchedNeg;
407   TList collections_histoEtaPhiHighPt;
408
409   Int_t count = 0;
410
411   while ((obj = iter->Next())) {
412     AliSpectraBothTrackCuts* entry = dynamic_cast<AliSpectraBothTrackCuts*> (obj);
413     if (entry == 0) 
414       continue;
415     
416     TH1I * histo = entry->GetHistoCuts();      
417     collections.Add(histo);
418     TH1F * histoNSelectedPos = entry->GetHistoNSelectedPos();      
419     collections_histoNSelectedPos.Add(histoNSelectedPos);
420     TH1F * histoNSelectedNeg = entry->GetHistoNSelectedNeg();      
421     collections_histoNSelectedNeg.Add(histoNSelectedNeg);
422     TH1F * histoNMatchedPos = entry->GetHistoNMatchedPos();      
423     collections_histoNMatchedPos.Add(histoNMatchedPos);
424     TH1F * histoNMatchedNeg = entry->GetHistoNMatchedNeg();      
425     collections_histoNMatchedNeg.Add(histoNMatchedNeg);
426     TH2F * histoEtaPhiHighPt = entry->GetHistoEtaPhiHighPt();      
427     collections_histoEtaPhiHighPt.Add(histoEtaPhiHighPt);
428     count++;
429   }
430   
431   fHistoCuts->Merge(&collections);
432   fHistoNSelectedPos->Merge(&collections_histoNSelectedPos);
433   fHistoNSelectedNeg->Merge(&collections_histoNSelectedNeg);
434   fHistoNMatchedPos->Merge(&collections_histoNMatchedPos);
435   fHistoNMatchedNeg->Merge(&collections_histoNMatchedNeg);
436   fHistoEtaPhiHighPt->Merge(&collections_histoEtaPhiHighPt);
437   
438   delete iter;
439
440   return count+1;
441 }
442