64bbed85092f685cab59761434541244ea04fb4f
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliPrimaryPionCuts.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved.       *
4  *                                                                                                                                              *
5  * Authors: Friederike Bock                                                                                                     *
6  * Version 1.0                                                                                                                          *
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 //--------------------------------------------- 
19 // Class handling all kinds of selection cuts for
20 // Gamma Conversion analysis
21 //---------------------------------------------
22 ////////////////////////////////////////////////
23
24
25 #include "AliPrimaryPionCuts.h"
26 #include "AliAODConversionPhoton.h"
27 #include "AliKFVertex.h"
28 #include "AliAODTrack.h"
29 #include "AliESDtrack.h"
30 #include "AliAnalysisManager.h"
31 #include "AliInputEventHandler.h"
32 #include "AliMCEventHandler.h"
33 #include "AliAODHandler.h"
34 #include "AliPIDResponse.h"
35 #include "TH1.h"
36 #include "TH2.h"
37 #include "AliStack.h"
38 #include "TObjString.h"
39 #include "AliAODEvent.h"
40 #include "AliESDEvent.h"
41 #include "TList.h"
42 class iostream;
43
44 using namespace std;
45
46 ClassImp(AliPrimaryPionCuts)
47
48
49 const char* AliPrimaryPionCuts::fgkCutNames[AliPrimaryPionCuts::kNCuts] = {
50         "kEtaCut",                              // 0
51         "kClsITSCut",                   // 1
52         "kClsTPCCut",                   // 2
53         "kDCAcut",                              // 3
54         "kPtCut",                               // 4
55         "kPiDedxSigmaITSCut",   // 5
56         "kPiDedxSigmaTPCCut",   // 6
57         "kPiTOFSigmaCut",               // 7 
58         "kMassCut"                              // 8
59 };
60
61 //________________________________________________________________________
62 AliPrimaryPionCuts::AliPrimaryPionCuts(const char *name,const char *title) : AliAnalysisCuts(name,title),
63         fHistograms(NULL),
64         fPIDResponse(NULL),
65         fEsdTrackCuts(NULL),
66         fEtaCut(0.9),
67         fEtaShift(0.0),
68         fDoEtaCut(kFALSE),
69         fPtCut(0.0),
70         fMinClsTPC(0), // minimum clusters in the TPC
71         fMinClsTPCToF(0), // minimum clusters to findable clusters
72         fDodEdxSigmaITSCut(kFALSE),
73         fDodEdxSigmaTPCCut(kTRUE),
74         fDoTOFsigmaCut(kFALSE), // RRnewTOF
75         fPIDnSigmaAbovePionLineITS(100),
76         fPIDnSigmaBelowPionLineITS(-100),
77         fPIDnSigmaAbovePionLineTPC(100),
78         fPIDnSigmaBelowPionLineTPC(-100),
79         fPIDnSigmaAbovePionLineTOF(100), // RRnewTOF
80         fPIDnSigmaBelowPionLineTOF(-100), // RRnewTOF
81         fUseCorrectedTPCClsInfo(kFALSE),
82         fUseTOFpid(kFALSE),
83         fRequireTOF(kFALSE),
84         fDoMassCut(kFALSE),
85         fMassCut(10),
86         fDoWeights(kFALSE),
87         fCutString(NULL),
88         fHistCutIndex(NULL),
89         fHistdEdxCuts(NULL),
90         fHistITSdEdxbefore(NULL),
91         fHistITSdEdxafter(NULL),
92         fHistTPCdEdxbefore(NULL),
93         fHistTPCdEdxafter(NULL),
94         fHistTPCdEdxSignalbefore(NULL),
95         fHistTPCdEdxSignalafter(NULL),
96         fHistTOFbefore(NULL),
97         fHistTOFafter(NULL),
98         fHistTrackDCAxyPtbefore(NULL),
99         fHistTrackDCAxyPtafter(NULL),
100         fHistTrackDCAzPtbefore(NULL),
101         fHistTrackDCAzPtafter(NULL),
102         fHistTrackNFindClsPtTPCbefore(NULL),
103         fHistTrackNFindClsPtTPCafter(NULL)
104 {
105         InitPIDResponse();
106         for(Int_t jj=0;jj<kNCuts;jj++){ fCuts[jj]=0; }
107         fCutString=new TObjString((GetCutNumber()).Data());
108
109         // Using standard function for setting Cuts
110         Bool_t selectPrimaries=kFALSE;
111         fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
112 }
113
114 //________________________________________________________________________
115 AliPrimaryPionCuts::~AliPrimaryPionCuts() {
116                 // Destructor
117         //Deleting fHistograms leads to seg fault it it's added to output collection of a task
118         // if(fHistograms)
119         //      delete fHistograms;
120         // fHistograms = NULL;
121
122         if(fCutString != NULL){
123                 delete fCutString;
124                 fCutString = NULL;
125         }
126 }
127
128 //________________________________________________________________________
129 void AliPrimaryPionCuts::InitCutHistograms(TString name, Bool_t preCut,TString cutNumber){
130
131         // Initialize Cut Histograms for QA (only initialized and filled if function is called)
132
133         TString cutName = "";
134         
135         if( cutNumber==""){
136                 cutName = GetCutNumber().Data();
137         }
138         else {
139                 cutName = cutNumber.Data();
140         } 
141
142         if(fHistograms != NULL){
143                 delete fHistograms;
144                 fHistograms=NULL;
145         }
146         if(fHistograms==NULL){
147                 fHistograms=new TList();
148                 if(name=="")fHistograms->SetName(Form("PionCuts_%s",cutName.Data()));
149                 else fHistograms->SetName(Form("%s_%s",name.Data(),cutName.Data()));
150         }
151
152
153         fHistCutIndex=new TH1F(Form("IsPionSelected %s",cutName.Data()),"IsPionSelected",10,-0.5,9.5);
154         fHistCutIndex->GetXaxis()->SetBinLabel(kPionIn+1,"in");
155         fHistCutIndex->GetXaxis()->SetBinLabel(kNoTracks+1,"no tracks");
156         fHistCutIndex->GetXaxis()->SetBinLabel(kTrackCuts+1,"Track cuts");
157         fHistCutIndex->GetXaxis()->SetBinLabel(kdEdxCuts+1,"dEdx");
158         fHistCutIndex->GetXaxis()->SetBinLabel(kPionOut+1,"out");
159         fHistograms->Add(fHistCutIndex);
160
161         // dEdx Cuts
162         fHistdEdxCuts=new TH1F(Form("dEdxCuts %s",cutName.Data()),"dEdxCuts",5,-0.5,4.5);
163         fHistdEdxCuts->GetXaxis()->SetBinLabel(1,"in");
164         fHistdEdxCuts->GetXaxis()->SetBinLabel(2,"ITSpion");
165         fHistdEdxCuts->GetXaxis()->SetBinLabel(3,"TPCpion");
166         fHistdEdxCuts->GetXaxis()->SetBinLabel(4,"TOFpion");
167         fHistdEdxCuts->GetXaxis()->SetBinLabel(5,"out");
168         fHistograms->Add(fHistdEdxCuts);
169         
170         TAxis *axisBeforeITS  = NULL;
171         TAxis *axisBeforedEdx = NULL;
172         TAxis *axisBeforeTOF  = NULL;
173         TAxis *axisBeforedEdxSignal = NULL;
174
175         if(preCut){
176                 fHistITSdEdxbefore=new TH2F(Form("Pion_ITS_before %s",cutName.Data()),"ITS dEdx pion before" ,150,0.05,20,400,-10,10);
177                 fHistograms->Add(fHistITSdEdxbefore);
178                 axisBeforeITS = fHistITSdEdxbefore->GetXaxis();
179
180                 fHistTPCdEdxbefore=new TH2F(Form("Pion_dEdx_before %s",cutName.Data()),"dEdx pion before" ,150,0.05,20,400,-10,10);
181                 fHistograms->Add(fHistTPCdEdxbefore);
182                 axisBeforedEdx = fHistTPCdEdxbefore->GetXaxis();
183
184                 fHistTPCdEdxSignalbefore=new TH2F(Form("Pion_dEdxSignal_before %s",cutName.Data()),"dEdx pion signal before" ,150,0.05,20.0,800,0.0,200);
185                 fHistograms->Add(fHistTPCdEdxSignalbefore);
186                 axisBeforedEdxSignal = fHistTPCdEdxSignalbefore->GetXaxis();
187
188                 fHistTOFbefore=new TH2F(Form("Pion_TOF_before %s",cutName.Data()),"TOF pion before" ,150,0.05,20,400,-6,10);
189                 fHistograms->Add(fHistTOFbefore);
190                 axisBeforeTOF = fHistTOFbefore->GetXaxis();
191                 
192                 fHistTrackDCAxyPtbefore = new TH2F(Form("hTrack_DCAxy_Pt_before %s",cutName.Data()),"DCAxy Vs Pt of tracks before",800,-4.0,4.0,400,0.,10.);
193                 fHistograms->Add(fHistTrackDCAxyPtbefore);      
194                 
195                 fHistTrackDCAzPtbefore  = new TH2F(Form("hTrack_DCAz_Pt_before %s",cutName.Data()), "DCAz  Vs Pt of tracks before",800,-4.0,4.0,400,0.,10.);
196                 fHistograms->Add(fHistTrackDCAzPtbefore); 
197                 
198                 fHistTrackNFindClsPtTPCbefore = new TH2F(Form("hTrack_NFindCls_Pt_TPC_before %s",cutName.Data()),"Track: N Findable Cls TPC Vs Pt before",100,0,1,400,0.,10.);
199                 fHistograms->Add(fHistTrackNFindClsPtTPCbefore); 
200         }
201
202         fHistITSdEdxafter=new TH2F(Form("Pion_ITS_after %s",cutName.Data()),"ITS dEdx pion after" ,150,0.05,20,400, -10,10);
203         fHistograms->Add(fHistITSdEdxafter);
204
205         fHistTPCdEdxafter=new TH2F(Form("Pion_dEdx_after %s",cutName.Data()),"dEdx pion after" ,150,0.05,20,400, -10,10);
206         fHistograms->Add(fHistTPCdEdxafter);
207
208         fHistTPCdEdxSignalafter=new TH2F(Form("Pion_dEdxSignal_after %s",cutName.Data()),"dEdx pion signal after" ,150,0.05,20.0,800,0.0,200);
209         fHistograms->Add(fHistTPCdEdxSignalafter);
210
211         fHistTOFafter=new TH2F(Form("Pion_TOF_after %s",cutName.Data()),"TOF pion after" ,150,0.05,20,400,-6,10);
212         fHistograms->Add(fHistTOFafter);
213         
214         fHistTrackDCAxyPtafter  = new TH2F(Form("hTrack_DCAxy_Pt_after %s",cutName.Data()),"DCAxy Vs Pt of tracks after",800,-4.0,4.0,400,0.,10.);
215         fHistograms->Add(fHistTrackDCAxyPtafter); 
216         
217         fHistTrackDCAzPtafter  = new TH2F(Form("hTrack_DCAz_Pt_after %s",cutName.Data()), "DCAz Vs Pt of tracks  after",800,-4.0,4.0,400,0.,10.);
218         fHistograms->Add(fHistTrackDCAzPtafter); 
219         
220         fHistTrackNFindClsPtTPCafter = new TH2F(Form("hTrack_NFindCls_Pt_TPC_after %s",cutName.Data()),"Track: N Findable Cls TPC Vs Pt after",100,0,1,400,0.,10.);
221         fHistograms->Add(fHistTrackNFindClsPtTPCafter); 
222         
223         TAxis *AxisAfter = fHistTPCdEdxafter->GetXaxis(); 
224         Int_t bins = AxisAfter->GetNbins();
225         Double_t from = AxisAfter->GetXmin();
226         Double_t to = AxisAfter->GetXmax();
227         Double_t *newBins = new Double_t[bins+1];
228         newBins[0] = from;
229         Double_t factor = TMath::Power(to/from, 1./bins);
230         for(Int_t i=1; i<=bins; ++i) newBins[i] = factor * newBins[i-1];
231         AxisAfter->Set(bins, newBins);
232         AxisAfter = fHistTOFafter->GetXaxis(); 
233         AxisAfter->Set(bins, newBins);
234         AxisAfter = fHistITSdEdxafter->GetXaxis();
235         AxisAfter->Set(bins,newBins); 
236         AxisAfter = fHistTPCdEdxSignalafter->GetXaxis();
237         AxisAfter->Set(bins,newBins);
238         
239         if(preCut){
240                 axisBeforeITS->Set(bins, newBins);
241                 axisBeforedEdx->Set(bins, newBins);
242                 axisBeforedEdxSignal->Set(bins,newBins);
243                 axisBeforeTOF->Set(bins, newBins);
244         
245         }
246         delete [] newBins;
247         
248         // Event Cuts and Info
249 }
250
251
252 //________________________________________________________________________
253 Bool_t AliPrimaryPionCuts::InitPIDResponse(){
254
255 // Set Pointer to AliPIDResponse
256
257   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
258
259   if(man) {
260     AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
261     fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
262     if(fPIDResponse)return kTRUE;
263
264   }
265
266   return kFALSE;
267 }
268 ///________________________________________________________________________
269 Bool_t AliPrimaryPionCuts::PionIsSelectedMC(Int_t labelParticle,AliStack *fMCStack){
270         
271         if( labelParticle < 0 || labelParticle >= fMCStack->GetNtrack() ) return kFALSE;
272         if( fMCStack->IsPhysicalPrimary(labelParticle) == kFALSE ) return kFALSE; 
273
274         TParticle* particle = fMCStack->Particle(labelParticle);
275
276         if( TMath::Abs( particle->GetPdgCode() ) != 211 )  return kFALSE;
277         
278         if( fDoEtaCut ){
279         if( particle->Eta() > (fEtaCut + fEtaShift) || particle->Eta() < (-fEtaCut + fEtaShift) )
280                 return kFALSE;
281         }
282         
283         return kTRUE;
284 }
285
286
287 ///________________________________________________________________________
288 Bool_t AliPrimaryPionCuts::PionIsSelected(AliESDtrack* lTrack){
289     //Selection of Reconstructed electrons
290         
291         Float_t b[2];
292         Float_t bCov[3];
293         lTrack->GetImpactParameters(b,bCov);
294
295         if (bCov[0]<=0 || bCov[2]<=0) {
296                 AliDebug(1, "Estimated b resolution lower or equal zero!");
297                 bCov[0]=0; bCov[2]=0;
298         }
299
300         Float_t dcaToVertexXY = b[0];
301         Float_t dcaToVertexZ  = b[1];
302         Double_t clsToF = GetNFindableClustersTPC(lTrack);
303
304
305         if (fHistCutIndex) fHistCutIndex->Fill(kPionIn);
306
307         if (lTrack == NULL){
308                 if (fHistCutIndex) fHistCutIndex->Fill(kNoTracks);
309                 return kFALSE;  
310         }   
311         if (fHistTrackDCAxyPtbefore) fHistTrackDCAxyPtbefore->Fill(dcaToVertexXY,lTrack->Pt());
312         if (fHistTrackDCAzPtbefore) fHistTrackDCAzPtbefore->Fill( dcaToVertexZ, lTrack->Pt());
313         if (fHistTrackNFindClsPtTPCbefore) fHistTrackNFindClsPtTPCbefore->Fill( clsToF, lTrack->Pt());
314
315         
316         if ( ! lTrack->GetConstrainedParam() ){
317                 return kFALSE;
318         }
319         AliVTrack * track = dynamic_cast<AliVTrack*>(lTrack);
320
321         // Track Cuts
322         if( !TrackIsSelected(lTrack) ){
323                 if (fHistCutIndex) fHistCutIndex->Fill(kTrackCuts);
324                 return kFALSE;
325         }
326
327         // dEdx Cuts
328         if( ! dEdxCuts( track ) ) {
329                 if(fHistCutIndex)fHistCutIndex->Fill(kdEdxCuts);
330                 return kFALSE;
331
332         }
333
334         //Pion passed the cuts
335         if (fHistCutIndex) fHistCutIndex->Fill(kPionOut);
336         if (fHistTrackDCAxyPtafter) fHistTrackDCAxyPtafter->Fill(dcaToVertexXY,lTrack->Pt());
337         if (fHistTrackDCAzPtafter) fHistTrackDCAzPtafter->Fill(dcaToVertexZ,lTrack->Pt());
338         if (fHistTrackNFindClsPtTPCafter) fHistTrackNFindClsPtTPCafter->Fill( clsToF, lTrack->Pt());
339
340         return kTRUE;
341 }
342
343 ///________________________________________________________________________
344 Bool_t AliPrimaryPionCuts::TrackIsSelected(AliESDtrack* lTrack) {
345         // Track Selection for Photon Reconstruction
346         
347         Double_t clsToF = GetNFindableClustersTPC(lTrack);
348         
349         if( ! fEsdTrackCuts->AcceptTrack(lTrack) ){
350                 return kFALSE;
351         }
352                 
353         if( fDoEtaCut ) {
354                 if(  lTrack->Eta() > (fEtaCut + fEtaShift) || lTrack->Eta() < (-fEtaCut + fEtaShift) ) {
355                         return kFALSE;
356                 }
357         }
358         
359         if( lTrack->Pt() < fPtCut ) {
360                 return kFALSE;          
361         }
362
363         if( clsToF < fMinClsTPCToF){
364                 return kFALSE;
365         }
366
367         return kTRUE;
368 }
369
370 ///________________________________________________________________________
371 Bool_t AliPrimaryPionCuts::dEdxCuts(AliVTrack *fCurrentTrack){
372
373         // Pion Identification Cuts for Photon reconstruction
374
375         if(!fPIDResponse){  InitPIDResponse();  }// Try to reinitialize PID Response
376         if(!fPIDResponse){  AliError("No PID Response"); return kFALSE;}// if still missing fatal error
377
378         Int_t cutIndex=0;
379
380         if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
381         if(fHistITSdEdxbefore)fHistITSdEdxbefore->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasITS(fCurrentTrack, AliPID::kPion));
382         if(fHistTPCdEdxbefore)fHistTPCdEdxbefore->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTPC(fCurrentTrack, AliPID::kPion));
383         if(fHistTPCdEdxSignalbefore)fHistTPCdEdxSignalbefore->Fill(fCurrentTrack->P(),TMath::Abs(fCurrentTrack->GetTPCsignal()));
384
385         cutIndex++;
386
387         if( fDodEdxSigmaITSCut == kTRUE ){
388                 if( fPIDResponse->NumberOfSigmasITS(fCurrentTrack,AliPID::kPion)<fPIDnSigmaBelowPionLineITS ||
389                                 fPIDResponse->NumberOfSigmasITS(fCurrentTrack,AliPID::kPion)> fPIDnSigmaAbovePionLineITS ){     
390                         if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
391                         return kFALSE;
392                 }
393         }
394                 
395         if(fHistITSdEdxafter)fHistITSdEdxafter->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasITS(fCurrentTrack, AliPID::kPion));
396         
397         
398         cutIndex++;
399                 
400                 
401         if(fDodEdxSigmaTPCCut == kTRUE){
402                 // TPC Pion Line
403                 if( fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)<fPIDnSigmaBelowPionLineTPC ||
404                         fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)>fPIDnSigmaAbovePionLineTPC){
405                         if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
406                         return kFALSE;
407                 }
408                 cutIndex++;
409         } else { cutIndex+=1; }
410         
411         if( ( fCurrentTrack->GetStatus() & AliESDtrack::kTOFpid ) && ( !( fCurrentTrack->GetStatus() & AliESDtrack::kTOFmismatch) ) ){
412                 if(fHistTOFbefore) fHistTOFbefore->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kPion));
413                 if(fUseTOFpid){
414                         if( fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kPion)>fPIDnSigmaAbovePionLineTOF ||
415                                 fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kPion)<fPIDnSigmaBelowPionLineTOF ){
416                                 if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
417                                 return kFALSE;
418                         }
419                 }
420                 if(fHistTOFafter)fHistTOFafter->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kPion));
421         } else if ( fRequireTOF == kTRUE ) {
422                 if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
423                 return kFALSE;
424         }
425         cutIndex++;
426                 
427         if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
428         if(fHistTPCdEdxafter)fHistTPCdEdxafter->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTPC(fCurrentTrack, AliPID::kPion));
429         if(fHistTPCdEdxSignalafter)fHistTPCdEdxSignalafter->Fill(fCurrentTrack->P(),TMath::Abs(fCurrentTrack->GetTPCsignal()));
430         
431         return kTRUE;
432 }
433
434 ///________________________________________________________________________
435 AliVTrack *AliPrimaryPionCuts::GetTrack(AliVEvent * event, Int_t label){
436     //Returns pointer to the track with given ESD label
437     //(Important for AOD implementation, since Track array in AOD data is different
438     //from ESD array, but ESD tracklabels are stored in AOD Tracks)
439
440         AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(event);
441         if(esdEvent) {
442                 if(label > event->GetNumberOfTracks() ) return NULL;
443                 AliESDtrack * track = esdEvent->GetTrack(label);
444                 return track;           
445         } else { 
446                 for(Int_t ii=0; ii<event->GetNumberOfTracks(); ii++) {
447                         AliVTrack * track = dynamic_cast<AliVTrack*>(event->GetTrack(ii));              
448                         if(track) { 
449                                 if(track->GetID() == label) {
450                                         return track;
451                                 }
452                         }
453                 }
454         }
455         cout << "track not found " << label << " " << event->GetNumberOfTracks() << endl;
456         return NULL;
457 }
458
459 ///________________________________________________________________________
460 Double_t AliPrimaryPionCuts::GetNFindableClustersTPC(AliESDtrack* lTrack){
461
462         Double_t clsToF=0;
463         if ( !fUseCorrectedTPCClsInfo ){
464                 if(lTrack->GetTPCNclsF()!=0){
465                         clsToF = (Double_t)lTrack->GetNcls(1)/(Double_t)lTrack->GetTPCNclsF();
466                 }
467         } else {
468                 clsToF = lTrack->GetTPCClusterInfo(2,0); //NOTE ask friederike  
469         }
470         return clsToF;
471
472 }
473
474 ///________________________________________________________________________
475 Bool_t AliPrimaryPionCuts::UpdateCutString() {
476 ///Update the cut string (if it has been created yet)
477
478         if(fCutString && fCutString->GetString().Length() == kNCuts) {
479                 fCutString->SetString(GetCutNumber());
480         } else {
481                 return kFALSE;
482         }
483         return kTRUE;
484
485 }
486
487 ///________________________________________________________________________
488 Bool_t AliPrimaryPionCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
489         // Initialize Cuts from a given Cut string
490
491         AliInfo(Form("Set PionCuts Number: %s",analysisCutSelection.Data()));
492         
493         if(analysisCutSelection.Length()!=kNCuts) {
494                 AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
495                 return kFALSE;
496         }
497         if(!analysisCutSelection.IsDigit()){
498                 AliError("Cut selection contains characters");
499                 return kFALSE;
500         }
501         
502         const char *cutSelection = analysisCutSelection.Data();
503         #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
504         for(Int_t ii=0;ii<kNCuts;ii++){
505                 ASSIGNARRAY(ii);
506         }
507
508         // Set Individual Cuts
509         for(Int_t ii=0;ii<kNCuts;ii++){
510                 if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
511         }
512
513         PrintCuts();
514         return kTRUE;
515 }
516 ///________________________________________________________________________
517 Bool_t AliPrimaryPionCuts::SetCut(cutIds cutID, const Int_t value) {
518         ///Set individual cut ID
519
520         //cout << "Updating cut  " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
521
522         switch (cutID) {
523                 case kEtaCut:
524                         if( SetEtaCut(value)) {
525                                 fCuts[kEtaCut] = value;
526                                 UpdateCutString();
527                                 return kTRUE;
528                         } else return kFALSE;
529                 case kClsITSCut:
530                         if( SetITSClusterCut(value) ) {
531                                 fCuts[kClsITSCut] = value;
532                                 UpdateCutString();
533                                 return kTRUE;                   
534                         } else return kFALSE;
535                 case kClsTPCCut:
536                         if( SetTPCClusterCut(value)) {
537                                 fCuts[kClsTPCCut] = value;
538                                 UpdateCutString();
539                                 return kTRUE;
540                         } else return kFALSE;
541                 case kDCACut:
542                         if( SetDCACut(value)) {
543                                 fCuts[kDCACut] = value;
544                                 UpdateCutString();
545                                 return kTRUE;
546                         } else return kFALSE;
547                 case kPtCut:    
548                         if( SetPtCut(value)) {
549                                 fCuts[kPtCut] = value;
550                                 UpdateCutString();
551                                 return kTRUE;
552                         } else return kFALSE;
553                 case kPidedxSigmaITSCut:
554                         if( SetITSdEdxCutPionLine(value)) { 
555                                 fCuts[kPidedxSigmaITSCut] = value;
556                                 UpdateCutString();
557                                 return kTRUE;
558                         } else return kFALSE;
559                 case kPidedxSigmaTPCCut:
560                         if( SetTPCdEdxCutPionLine(value)) { 
561                                 fCuts[kPidedxSigmaTPCCut] = value;
562                                 UpdateCutString();
563                                 return kTRUE;
564                         } else return kFALSE;
565                 case kPiTOFSigmaPID:
566                         if( SetTOFPionPIDCut(value)) {
567                                 fCuts[kPiTOFSigmaPID] = value;
568                                 UpdateCutString();
569                                 return kTRUE;
570                         } else return kFALSE;
571                 case kMassCut:
572                         if( SetMassCut(value)) {
573                                 fCuts[kMassCut] = value;
574                                 UpdateCutString();
575                                 return kTRUE;
576                         } else return kFALSE;
577                 case kNCuts:
578                         cout << "Error:: Cut id out of range"<< endl;
579                         return kFALSE;
580         }
581                 
582         cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
583         return kFALSE;
584 }
585
586 ///________________________________________________________________________
587
588 void AliPrimaryPionCuts::PrintCuts() {
589     // Print out current Cut Selection
590         for(Int_t ic = 0; ic < kNCuts; ic++) {
591                 printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
592         }
593 }
594
595 ///________________________________________________________________________
596 Bool_t AliPrimaryPionCuts::SetITSdEdxCutPionLine(Int_t ededxSigmaCut){ 
597         switch(ededxSigmaCut){
598                 case 0: 
599                         fDodEdxSigmaITSCut = kFALSE;
600                         fPIDnSigmaBelowPionLineITS=-100;
601                         fPIDnSigmaAbovePionLineITS= 100;
602                         break;
603                 case 1: // -10,10
604                         fDodEdxSigmaITSCut = kTRUE;
605                         fPIDnSigmaBelowPionLineITS=-10;
606                         fPIDnSigmaAbovePionLineITS=10;
607                         break;
608                 case 2: // -6,7
609                         fDodEdxSigmaITSCut = kTRUE;
610                         fPIDnSigmaBelowPionLineITS=-6;
611                         fPIDnSigmaAbovePionLineITS=7;
612                         break;
613                 case 3: // -5,5
614                         fDodEdxSigmaITSCut = kTRUE;
615                         fPIDnSigmaBelowPionLineITS=-5;
616                         fPIDnSigmaAbovePionLineITS=5;
617                         break;
618                 case 4: // -4,5
619                         fDodEdxSigmaITSCut = kTRUE;
620                         fPIDnSigmaBelowPionLineITS=-4;
621                         fPIDnSigmaAbovePionLineITS=5;
622                         break;
623                 case 5: // -3,5
624                         fDodEdxSigmaITSCut = kTRUE;
625                         fPIDnSigmaBelowPionLineITS=-3;
626                         fPIDnSigmaAbovePionLineITS=5;
627                         break;
628                 case 6: // -4,4
629                         fDodEdxSigmaITSCut = kTRUE;
630                         fPIDnSigmaBelowPionLineITS=-4;
631                         fPIDnSigmaAbovePionLineITS=4;
632                         break;
633                 case 7: // -2.5,4
634                         fDodEdxSigmaITSCut = kTRUE;
635                         fPIDnSigmaBelowPionLineITS=-2.5;
636                         fPIDnSigmaAbovePionLineITS=4;
637                         break;
638                 case 8: // -2,3.5
639                         fDodEdxSigmaITSCut = kTRUE;
640                         fPIDnSigmaBelowPionLineITS=-2;
641                         fPIDnSigmaAbovePionLineITS=3.5;
642                         break;
643                 default:
644                         cout<<"Warning: ITSdEdxCutPionLine not defined"<<ededxSigmaCut<<endl;
645                         return kFALSE;
646                 
647         }
648         return kTRUE;
649 }
650
651 ///________________________________________________________________________
652 Bool_t AliPrimaryPionCuts::SetTPCdEdxCutPionLine(Int_t ededxSigmaCut){
653         switch(ededxSigmaCut){
654                 case 0: 
655                         fDodEdxSigmaTPCCut = kFALSE;
656                         fPIDnSigmaBelowPionLineTPC=-10;
657                         fPIDnSigmaAbovePionLineTPC=10;
658                         break;
659                 case 1: // -10,10
660                         fDodEdxSigmaTPCCut = kTRUE;
661                         fPIDnSigmaBelowPionLineTPC=-10;
662                         fPIDnSigmaAbovePionLineTPC=10;
663                         break;
664                 case 2: // -6,7
665                         fDodEdxSigmaTPCCut = kTRUE;
666                         fPIDnSigmaBelowPionLineTPC=-6;
667                         fPIDnSigmaAbovePionLineTPC=7;
668                         break;
669                 case 3: // -5,5
670                         fDodEdxSigmaTPCCut = kTRUE;
671                         fPIDnSigmaBelowPionLineTPC=-5;
672                         fPIDnSigmaAbovePionLineTPC=5;
673                         break;
674                 case 4: // -4,5
675                         fDodEdxSigmaTPCCut = kTRUE;
676                         fPIDnSigmaBelowPionLineTPC=-4;
677                         fPIDnSigmaAbovePionLineTPC=5;
678                         break;  
679                 case 5: // -4,4
680                         fDodEdxSigmaTPCCut = kTRUE;
681                         fPIDnSigmaBelowPionLineTPC=-4;
682                         fPIDnSigmaAbovePionLineTPC=4;
683                         break;
684                 case 6: // -3,4
685                         fDodEdxSigmaTPCCut = kTRUE;
686                         fPIDnSigmaBelowPionLineTPC=-3;
687                         fPIDnSigmaAbovePionLineTPC=4;
688                         break;
689                 case 7: // -3,3
690                         fDodEdxSigmaTPCCut = kTRUE;
691                         fPIDnSigmaBelowPionLineTPC=-3;
692                         fPIDnSigmaAbovePionLineTPC=3;
693                         break;
694                 case 8: // -2,3.
695                         fDodEdxSigmaTPCCut = kTRUE;
696                         fPIDnSigmaBelowPionLineTPC=-2;
697                         fPIDnSigmaAbovePionLineTPC=3.;
698                         break;
699                 default:
700                         cout<<"Warning: TPCdEdxCutPionLine not defined"<<ededxSigmaCut<<endl;
701                         return kFALSE;
702         }
703         return kTRUE;
704 }
705
706 ///________________________________________________________________________
707 Bool_t AliPrimaryPionCuts::SetITSClusterCut(Int_t clsITSCut){
708         if( !fEsdTrackCuts ) {
709                 cout<<"Warning: AliESDtrackCut is not initialized "<<endl;
710                 return kFALSE;
711         }
712
713         switch(clsITSCut){
714                 case 0: 
715                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
716                         break;
717                 case 1: 
718                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
719                         break;  //1 hit first layer of SPD
720                 case 2: 
721                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
722                         break; //1 hit in any layer of SPD
723                 case 3: 
724                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
725                         fEsdTrackCuts->SetMinNClustersITS(4);
726                         // 4 hits in total in the ITS. At least 1 hit in the first layer of SPD  
727                         break;
728                 case 4: 
729                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
730                         fEsdTrackCuts->SetMinNClustersITS(3);
731                         // 3 hits in total in the ITS. At least 1 hit in any layer of SPD
732                         break;
733                 case 5: 
734                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
735                         fEsdTrackCuts->SetMinNClustersITS(4);
736                         // 4 hits in total in the ITS. At least 1 hit in any layer of SPD
737                         break;
738                 case 6: 
739                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
740                         fEsdTrackCuts->SetMinNClustersITS(5);
741                         // 5 hits in total in the ITS. At least 1 hit in any layer of SPD
742                         break;
743                 default:
744                         cout<<"Warning: clsITSCut not defined "<<clsITSCut<<endl;
745                         return kFALSE;
746         }
747         return kTRUE;
748 }
749
750 ///________________________________________________________________________
751 Bool_t AliPrimaryPionCuts::SetTPCClusterCut(Int_t clsTPCCut){  
752         switch(clsTPCCut){
753                 case 0: // 0
754                         fMinClsTPC= 0.;
755                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
756                         break;
757                 case 1:  // 70
758                         fMinClsTPC= 70.;
759                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
760                         break;
761                 case 2:  // 80
762                         fMinClsTPC= 80.;
763                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
764                         break;
765                 case 3:  // 100
766                         fMinClsTPC= 100.;
767                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
768                         break;
769                 case 4:  // 0% of findable clusters
770                                 fMinClsTPC= 70.;  
771                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
772                         fMinClsTPCToF= 0.0;
773                         fUseCorrectedTPCClsInfo=0;
774                         break;
775                 case 5:  // 35% of findable clusters
776                         fMinClsTPC = 70.;  
777                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
778                         fMinClsTPCToF= 0.35;
779                         fUseCorrectedTPCClsInfo=0;
780                         break;
781                 case 6:  // 60% of findable clusters
782                         fMinClsTPC= 70.;  
783                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
784                         fMinClsTPCToF= 0.6;
785                         fUseCorrectedTPCClsInfo=0;
786                         break;
787                 case 7:  // 70% of findable clusters
788                         fMinClsTPC= 70.;  
789                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
790                         fMinClsTPCToF= 0.7;
791                         fUseCorrectedTPCClsInfo=0;
792                         break;
793                 case 8: fMinClsTPC = 0.;  
794                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
795                         fMinClsTPCToF= 0.35;
796                         fUseCorrectedTPCClsInfo=0;
797                         break;
798                 case 9:  // 35% of findable clusters
799                         fMinClsTPC = 70.;  
800                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
801                         fMinClsTPCToF= 0.35;
802                         fUseCorrectedTPCClsInfo=1;
803                         break;
804                 
805                 default:
806                         cout<<"Warning: clsTPCCut not defined "<<clsTPCCut<<endl;
807                         return kFALSE;
808         }
809         return kTRUE;
810 }
811
812 ///________________________________________________________________________
813 Bool_t AliPrimaryPionCuts::SetEtaCut(Int_t etaCut){ 
814         // Set eta Cut
815         switch(etaCut){
816                 case 0: 
817                         fEtaCut = 100.;
818                         fDoEtaCut = kFALSE;
819                         break;
820                 case 1: // 1.4
821                         fEtaCut = 1.4;
822                         fDoEtaCut = kTRUE;
823                         break;
824                 case 2: // 1.2
825                         fEtaCut = 1.2;
826                         fDoEtaCut = kTRUE;
827                         break;
828                 case 3: // 0.9
829                         fEtaCut = 0.9;
830                         fDoEtaCut = kTRUE;
831                         break;
832                 case 4: // 0.8
833                         fEtaCut = 0.8;
834                         fDoEtaCut = kTRUE;
835                         break;
836                 case 5: // 0.75
837                         fEtaCut = 0.75;
838                         fDoEtaCut = kTRUE;
839                         break;
840                 case 6: //0.6
841                         fEtaCut = 0.6; //changed from 0.4 to 0.6 2013.06.10
842                         fDoEtaCut = kTRUE;
843                         break;
844                 case 7: //0.5
845                         fEtaCut = 0.5; //changed from 0.3 to 0.5 2013.06.10
846                         fDoEtaCut = kTRUE;
847                         break;
848                 case 8: 
849                         fEtaCut = 0.4;
850                         fDoEtaCut = kTRUE;
851                         break;
852                 default:
853                         cout<<"Warning: EtaCut not defined "<<etaCut<<endl;
854                         return kFALSE;
855         }
856         return kTRUE;
857 }
858
859 ///________________________________________________________________________
860 Bool_t AliPrimaryPionCuts::SetPtCut(Int_t ptCut){ 
861         switch(ptCut){    
862                 case 0: 
863                         fPtCut = 0.075;         
864                         break;
865                 case 1:  // 0.1
866                         fPtCut  = 0.1;  
867                         break;
868                 case 2:  // 0.125 GeV
869                         fPtCut  = 0.125;                
870                         break;
871                 case 3: // 0.15 GeV
872                         fPtCut  = 0.15;
873                         break;
874                 default:
875                         cout<<"Warning: PtCut not defined "<<ptCut<<endl;
876                         return kFALSE;
877         }
878         return kTRUE;
879 }
880
881
882 ///________________________________________________________________________
883 Bool_t AliPrimaryPionCuts::SetDCACut(Int_t dcaCut)
884
885         // Set DCA Cut
886         if( !fEsdTrackCuts ) {
887                 cout<<"Warning: AliESDtrackCut is not initialized "<<endl;
888                 return kFALSE;
889         }
890   
891         switch(dcaCut){   
892                 case 0: 
893                         //Open cuts//
894                         fEsdTrackCuts->SetMaxDCAToVertexZ(1000);
895                         fEsdTrackCuts->SetMaxDCAToVertexXY(1000);
896                         fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
897                         break;
898                 case 1: 
899                         fEsdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
900                         fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
901                         break;
902                 case 2: fEsdTrackCuts->SetMaxDCAToVertexZ(2);
903                         fEsdTrackCuts->SetMaxDCAToVertexXY(1);
904                         fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
905                         break; 
906                 default:
907                         cout<<"Warning: dcaCut not defined "<<dcaCut<<endl;
908                         return kFALSE;
909         }
910         return kTRUE;
911 }
912
913 ///________________________________________________________________________
914 Bool_t AliPrimaryPionCuts::SetTOFPionPIDCut(Int_t TOFelectronPID){
915     // Set Cut
916         switch(TOFelectronPID){ 
917                 case 0: // no cut
918                         fRequireTOF = kFALSE;
919                         fUseTOFpid = kFALSE;
920                         fPIDnSigmaBelowPionLineTOF=-100;
921                         fPIDnSigmaAbovePionLineTOF=100;
922                         break;
923                 case 1: // -7,7
924                         fRequireTOF = kFALSE;
925                         fUseTOFpid = kTRUE;
926                         fPIDnSigmaBelowPionLineTOF=-7;
927                         fPIDnSigmaAbovePionLineTOF=7;
928                         break;
929                 case 2: // -5,5
930                         fRequireTOF = kFALSE;
931                         fUseTOFpid = kTRUE;
932                         fPIDnSigmaBelowPionLineTOF=-5;
933                         fPIDnSigmaAbovePionLineTOF=5;
934                         break;
935                 case 3: // -3,5
936                         fRequireTOF = kFALSE;
937                         fUseTOFpid = kTRUE;
938                         fPIDnSigmaBelowPionLineTOF=-3;
939                         fPIDnSigmaAbovePionLineTOF=5;
940                         break;
941                 case 4: // -2,3
942                         fRequireTOF = kFALSE;
943                         fUseTOFpid = kTRUE;
944                         fPIDnSigmaBelowPionLineTOF=-2;
945                         fPIDnSigmaAbovePionLineTOF=3;
946                         break;
947                 case 5: // -3, 3 TOF mandatory
948                         fRequireTOF = kTRUE;
949                         fUseTOFpid  = kTRUE;
950                         fPIDnSigmaBelowPionLineTOF= -3;
951                         fPIDnSigmaAbovePionLineTOF=  3;
952                         break;
953                 default:
954                         cout<<"Warning: TOFPionCut not defined "<<TOFelectronPID<<endl;
955                 return kFALSE;
956     } 
957     return kTRUE;
958 }
959
960 ///________________________________________________________________________
961 Bool_t AliPrimaryPionCuts::SetMassCut(Int_t massCut){
962     // Set Cut
963         switch(massCut){ 
964                 case 0: // no cut
965                         fDoMassCut = kFALSE;
966                         fMassCut = 10;
967                         break;
968                 case 1: // cut at 1 GeV/c^2
969                         fDoMassCut = kTRUE;
970                         fMassCut = 1;
971                         break;
972                 case 2: // cut at 0.7 GeV/c^2
973                         fDoMassCut = kTRUE;
974                         fMassCut = 0.7;
975                         break;
976                 case 3: // cut at 0.6 GeV/c^2
977                         fDoMassCut = kTRUE;
978                         fMassCut = 0.6;
979                         break;
980                 case 4: // cut at eta mass
981                         fDoMassCut = kTRUE;
982                         fMassCut = 0.547853;
983                         break;
984                 case 5: // cut at 0.5 GeV/c^2
985                         fDoMassCut = kTRUE;
986                         fMassCut = 0.5;
987                         break;
988                 default:
989                         cout<<"Warning: MassCut not defined "<<massCut<<endl;
990                 return kFALSE;
991     } 
992     return kTRUE;
993 }
994
995
996 ///________________________________________________________________________
997 TString AliPrimaryPionCuts::GetCutNumber(){
998         // returns TString with current cut number
999         TString a(kNCuts);
1000         for(Int_t ii=0;ii<kNCuts;ii++){
1001                 a.Append(Form("%d",fCuts[ii]));
1002         }
1003         return a;
1004 }
1005
1006
1007 ///________________________________________________________________________
1008 AliPrimaryPionCuts* AliPrimaryPionCuts::GetStandardCuts2010PbPb(){
1009     //Create and return standard 2010 PbPb cuts
1010     AliPrimaryPionCuts *cuts=new AliPrimaryPionCuts("StandardCuts2010PbPb","StandardCuts2010PbPb");
1011     if(!cuts->InitializeCutsFromCutString("000000400")){
1012                 cout<<"Warning: Initialization of Standardcuts2010PbPb failed"<<endl;   
1013         }
1014     return cuts;
1015 }
1016
1017 ///________________________________________________________________________
1018 AliPrimaryPionCuts* AliPrimaryPionCuts::GetStandardCuts2010pp(){
1019     //Create and return standard 2010 PbPb cuts
1020     AliPrimaryPionCuts *cuts=new AliPrimaryPionCuts("StandardCuts2010pp","StandardCuts2010pp");
1021                                           
1022     if(!cuts->InitializeCutsFromCutString("000000400")){
1023                 cout<<"Warning: Initialization of Standardcuts2010pp failed"<<endl;
1024         }
1025     return cuts;
1026 }
1027