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",
  51         "kClsITSCut",
  52         "kClsTPCCut",
  53         "kDCAcut",
  54         "kPtCut",
  55         "kPiDedxSigmaITSCut",
  56         "kPiDedxSigmaTPCCut",
  57         "kPiTOFSigmaCut",
  58         "kMassCut"
  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         if (fHistTrackDCAxyPtbefore) fHistTrackDCAxyPtbefore->Fill(dcaToVertexXY,lTrack->Pt());
 305         if (fHistTrackDCAzPtbefore) fHistTrackDCAzPtbefore->Fill( dcaToVertexZ, lTrack->Pt());
 306         if (fHistTrackNFindClsPtTPCbefore) fHistTrackNFindClsPtTPCbefore->Fill( clsToF, lTrack->Pt());
 307
 308         if (fHistCutIndex) fHistCutIndex->Fill(kPionIn);
 309
 310         if (lTrack == NULL){
 311                 if (fHistCutIndex) fHistCutIndex->Fill(kNoTracks);
 312                 return kFALSE;
 313         }
 314
 315         if ( ! lTrack->GetConstrainedParam() ){
 316                 return kFALSE;
 317         }
 318         AliVTrack * track = dynamic_cast<AliVTrack*>(lTrack);
 319
 320         // Track Cuts
 321         if( !TrackIsSelected(lTrack) ){
 322                 if (fHistCutIndex) fHistCutIndex->Fill(kTrackCuts);
 323                 return kFALSE;
 324         }
 325
 326         // dEdx Cuts
 327         if( ! dEdxCuts( track ) ) {
 328                 if(fHistCutIndex)fHistCutIndex->Fill(kdEdxCuts);
 329                 return kFALSE;
 330
 331         }
 332
 333         //Pion passed the cuts
 334         if (fHistCutIndex) fHistCutIndex->Fill(kPionOut);
 335         if (fHistTrackDCAxyPtafter) fHistTrackDCAxyPtafter->Fill(dcaToVertexXY,lTrack->Pt());
 336         if (fHistTrackDCAzPtafter) fHistTrackDCAzPtafter->Fill(dcaToVertexZ,lTrack->Pt());
 337         if (fHistTrackNFindClsPtTPCafter) fHistTrackNFindClsPtTPCafter->Fill( clsToF, lTrack->Pt());
 338
 339         return kTRUE;
 340 }
 341
 342 ///________________________________________________________________________
 343 Bool_t AliPrimaryPionCuts::TrackIsSelected(AliESDtrack* lTrack) {
 344         // Track Selection for Photon Reconstruction
 345
 346         Double_t clsToF = GetNFindableClustersTPC(lTrack);
 347
 348         if( ! fEsdTrackCuts->AcceptTrack(lTrack) ){
 349                 return kFALSE;
 350         }
 351
 352         if( fDoEtaCut ) {
 353                 if(  lTrack->Eta() > (fEtaCut + fEtaShift) || lTrack->Eta() < (-fEtaCut + fEtaShift) ) {
 354                         return kFALSE;
 355                 }
 356         }
 357
 358         if( lTrack->Pt() < fPtCut ) {
 359                 return kFALSE;
 360         }
 361
 362         if( clsToF < fMinClsTPCToF){
 363                 return kFALSE;
 364         }
 365
 366         return kTRUE;
 367 }
 368
 369 ///________________________________________________________________________
 370 Bool_t AliPrimaryPionCuts::dEdxCuts(AliVTrack *fCurrentTrack){
 371
 372         // Pion Identification Cuts for Photon reconstruction
 373
 374         if(!fPIDResponse){  InitPIDResponse();  }// Try to reinitialize PID Response
 375         if(!fPIDResponse){  AliError("No PID Response"); return kFALSE;}// if still missing fatal error
 376
 377         Int_t cutIndex=0;
 378
 379         if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
 380         if(fHistITSdEdxbefore)fHistITSdEdxbefore->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasITS(fCurrentTrack, AliPID::kPion));
 381         if(fHistTPCdEdxbefore)fHistTPCdEdxbefore->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTPC(fCurrentTrack, AliPID::kPion));
 382         if(fHistTPCdEdxSignalbefore)fHistTPCdEdxSignalbefore->Fill(fCurrentTrack->P(),TMath::Abs(fCurrentTrack->GetTPCsignal()));
 383
 384         cutIndex++;
 385
 386         if( fDodEdxSigmaITSCut == kTRUE ){
 387                 if( fPIDResponse->NumberOfSigmasITS(fCurrentTrack,AliPID::kPion)<fPIDnSigmaBelowPionLineITS ||
 388                                 fPIDResponse->NumberOfSigmasITS(fCurrentTrack,AliPID::kPion)> fPIDnSigmaAbovePionLineITS ){
 389                         if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
 390                         return kFALSE;
 391                 }
 392         }
 393
 394         if(fHistITSdEdxafter)fHistITSdEdxafter->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasITS(fCurrentTrack, AliPID::kPion));
 395
 396
 397         cutIndex++;
 398
 399
 400         if(fDodEdxSigmaTPCCut == kTRUE){
 401                 // TPC Pion Line
 402                 if( fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)<fPIDnSigmaBelowPionLineTPC ||
 403                         fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)>fPIDnSigmaAbovePionLineTPC){
 404                         if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
 405                         return kFALSE;
 406                 }
 407                 cutIndex++;
 408         } else { cutIndex+=1; }
 409
 410         if( ( fCurrentTrack->GetStatus() & AliESDtrack::kTOFpid ) && ( !( fCurrentTrack->GetStatus() & AliESDtrack::kTOFmismatch) ) ){
 411                 if(fHistTOFbefore) fHistTOFbefore->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kPion));
 412                 if(fUseTOFpid){
 413                         if( fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kPion)>fPIDnSigmaAbovePionLineTOF ||
 414                                 fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kPion)<fPIDnSigmaBelowPionLineTOF ){
 415                                 if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
 416                                 return kFALSE;
 417                         }
 418                 }
 419                 if(fHistTOFafter)fHistTOFafter->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTOF(fCurrentTrack, AliPID::kPion));
 420         } else if ( fRequireTOF == kTRUE ) {
 421                 if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
 422                 return kFALSE;
 423         }
 424         cutIndex++;
 425
 426         if(fHistdEdxCuts)fHistdEdxCuts->Fill(cutIndex);
 427         if(fHistTPCdEdxafter)fHistTPCdEdxafter->Fill(fCurrentTrack->P(),fPIDResponse->NumberOfSigmasTPC(fCurrentTrack, AliPID::kPion));
 428         if(fHistTPCdEdxSignalafter)fHistTPCdEdxSignalafter->Fill(fCurrentTrack->P(),TMath::Abs(fCurrentTrack->GetTPCsignal()));
 429
 430         return kTRUE;
 431 }
 432
 433 ///________________________________________________________________________
 434 AliVTrack *AliPrimaryPionCuts::GetTrack(AliVEvent * event, Int_t label){
 435     //Returns pointer to the track with given ESD label
 436     //(Important for AOD implementation, since Track array in AOD data is different
 437     //from ESD array, but ESD tracklabels are stored in AOD Tracks)
 438
 439         AliESDEvent * esdEvent = dynamic_cast<AliESDEvent*>(event);
 440         if(esdEvent) {
 441                 if(label > event->GetNumberOfTracks() ) return NULL;
 442                 AliESDtrack * track = esdEvent->GetTrack(label);
 443                 return track;
 444         } else {
 445                 for(Int_t ii=0; ii<event->GetNumberOfTracks(); ii++) {
 446                         AliVTrack * track = dynamic_cast<AliVTrack*>(event->GetTrack(ii));
 447                         if(track) {
 448                                 if(track->GetID() == label) {
 449                                         return track;
 450                                 }
 451                         }
 452                 }
 453         }
 454         cout << "track not found " << label << " " << event->GetNumberOfTracks() << endl;
 455         return NULL;
 456 }
 457
 458 ///________________________________________________________________________
 459 Double_t AliPrimaryPionCuts::GetNFindableClustersTPC(AliESDtrack* lTrack){
 460
 461         Double_t clsToF=0;
 462         if ( !fUseCorrectedTPCClsInfo ){
 463                 if(lTrack->GetTPCNclsF()!=0){
 464                         clsToF = (Double_t)lTrack->GetNcls(1)/(Double_t)lTrack->GetTPCNclsF();
 465                 }
 466         } else {
 467                 clsToF = lTrack->GetTPCClusterInfo(2,0); //NOTE ask friederike
 468         }
 469         return clsToF;
 470
 471 }
 472
 473 ///________________________________________________________________________
 474 Bool_t AliPrimaryPionCuts::UpdateCutString() {
 475 ///Update the cut string (if it has been created yet)
 476
 477         if(fCutString && fCutString->GetString().Length() == kNCuts) {
 478                 fCutString->SetString(GetCutNumber());
 479         } else {
 480                 return kFALSE;
 481         }
 482         return kTRUE;
 483
 484 }
 485
 486 ///________________________________________________________________________
 487 Bool_t AliPrimaryPionCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
 488         // Initialize Cuts from a given Cut string
 489
 490         AliInfo(Form("Set PionCuts Number: %s",analysisCutSelection.Data()));
 491
 492         if(analysisCutSelection.Length()!=kNCuts) {
 493                 AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
 494                 return kFALSE;
 495         }
 496         if(!analysisCutSelection.IsDigit()){
 497                 AliError("Cut selection contains characters");
 498                 return kFALSE;
 499         }
 500
 501         const char *cutSelection = analysisCutSelection.Data();
 502         #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
 503         for(Int_t ii=0;ii<kNCuts;ii++){
 504                 ASSIGNARRAY(ii);
 505         }
 506
 507         // Set Individual Cuts
 508         for(Int_t ii=0;ii<kNCuts;ii++){
 509                 if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
 510         }
 511
 512         PrintCuts();
 513         return kTRUE;
 514 }
 515 ///________________________________________________________________________
 516 Bool_t AliPrimaryPionCuts::SetCut(cutIds cutID, const Int_t value) {
 517         ///Set individual cut ID
 518
 519         //cout << "Updating cut  " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
 520
 521         switch (cutID) {
 522                 case kEtaCut:
 523                         if( SetEtaCut(value)) {
 524                                 fCuts[kEtaCut] = value;
 525                                 UpdateCutString();
 526                                 return kTRUE;
 527                         } else return kFALSE;
 528                 case kClsITSCut:
 529                         if( SetITSClusterCut(value) ) {
 530                                 fCuts[kClsITSCut] = value;
 531                                 UpdateCutString();
 532                                 return kTRUE;
 533                         } else return kFALSE;
 534                 case kClsTPCCut:
 535                         if( SetTPCClusterCut(value)) {
 536                                 fCuts[kClsTPCCut] = value;
 537                                 UpdateCutString();
 538                                 return kTRUE;
 539                         } else return kFALSE;
 540                 case kDCACut:
 541                         if( SetDCACut(value)) {
 542                                 fCuts[kDCACut] = value;
 543                                 UpdateCutString();
 544                                 return kTRUE;
 545                         } else return kFALSE;
 546                 case kPtCut:
 547                         if( SetPtCut(value)) {
 548                                 fCuts[kPtCut] = value;
 549                                 UpdateCutString();
 550                                 return kTRUE;
 551                         } else return kFALSE;
 552                 case kPidedxSigmaITSCut:
 553                         if( SetITSdEdxCutPionLine(value)) {
 554                                 fCuts[kPidedxSigmaITSCut] = value;
 555                                 UpdateCutString();
 556                                 return kTRUE;
 557                         } else return kFALSE;
 558                 case kPidedxSigmaTPCCut:
 559                         if( SetTPCdEdxCutPionLine(value)) {
 560                                 fCuts[kPidedxSigmaTPCCut] = value;
 561                                 UpdateCutString();
 562                                 return kTRUE;
 563                         } else return kFALSE;
 564                 case kPiTOFSigmaPID:
 565                         if( SetTOFPionPIDCut(value)) {
 566                                 fCuts[kPiTOFSigmaPID] = value;
 567                                 UpdateCutString();
 568                                 return kTRUE;
 569                         } else return kFALSE;
 570                 case kMassCut:
 571                         if( SetMassCut(value)) {
 572                                 fCuts[kMassCut] = value;
 573                                 UpdateCutString();
 574                                 return kTRUE;
 575                         } else return kFALSE;
 576                 case kNCuts:
 577                         cout << "Error:: Cut id out of range"<< endl;
 578                         return kFALSE;
 579         }
 580
 581         cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
 582         return kFALSE;
 583 }
 584
 585 ///________________________________________________________________________
 586
 587 void AliPrimaryPionCuts::PrintCuts() {
 588     // Print out current Cut Selection
 589         for(Int_t ic = 0; ic < kNCuts; ic++) {
 590                 printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
 591         }
 592 }
 593
 594 ///________________________________________________________________________
 595 Bool_t AliPrimaryPionCuts::SetITSdEdxCutPionLine(Int_t ededxSigmaCut){
 596         switch(ededxSigmaCut){
 597                 case 0:
 598                         fDodEdxSigmaITSCut = kFALSE;
 599                         fPIDnSigmaBelowPionLineITS=-100;
 600                         fPIDnSigmaAbovePionLineITS= 100;
 601                         break;
 602                 case 1: // -10,10
 603                         fDodEdxSigmaITSCut = kTRUE;
 604                         fPIDnSigmaBelowPionLineITS=-10;
 605                         fPIDnSigmaAbovePionLineITS=10;
 606                         break;
 607                 case 2: // -6,7
 608                         fDodEdxSigmaITSCut = kTRUE;
 609                         fPIDnSigmaBelowPionLineITS=-6;
 610                         fPIDnSigmaAbovePionLineITS=7;
 611                         break;
 612                 case 3: // -5,5
 613                         fDodEdxSigmaITSCut = kTRUE;
 614                         fPIDnSigmaBelowPionLineITS=-5;
 615                         fPIDnSigmaAbovePionLineITS=5;
 616                         break;
 617                 case 4: // -4,5
 618                         fDodEdxSigmaITSCut = kTRUE;
 619                         fPIDnSigmaBelowPionLineITS=-4;
 620                         fPIDnSigmaAbovePionLineITS=5;
 621                         break;
 622                 case 5: // -3,5
 623                         fDodEdxSigmaITSCut = kTRUE;
 624                         fPIDnSigmaBelowPionLineITS=-3;
 625                         fPIDnSigmaAbovePionLineITS=5;
 626                         break;
 627                 case 6: // -4,4
 628                         fDodEdxSigmaITSCut = kTRUE;
 629                         fPIDnSigmaBelowPionLineITS=-4;
 630                         fPIDnSigmaAbovePionLineITS=4;
 631                         break;
 632                 case 7: // -2.5,4
 633                         fDodEdxSigmaITSCut = kTRUE;
 634                         fPIDnSigmaBelowPionLineITS=-2.5;
 635                         fPIDnSigmaAbovePionLineITS=4;
 636                         break;
 637                 case 8: // -2,3.5
 638                         fDodEdxSigmaITSCut = kTRUE;
 639                         fPIDnSigmaBelowPionLineITS=-2;
 640                         fPIDnSigmaAbovePionLineITS=3.5;
 641                         break;
 642                 default:
 643                         cout<<"Warning: ITSdEdxCutPionLine not defined"<<ededxSigmaCut<<endl;
 644                         return kFALSE;
 645
 646         }
 647         return kTRUE;
 648 }
 649
 650 ///________________________________________________________________________
 651 Bool_t AliPrimaryPionCuts::SetTPCdEdxCutPionLine(Int_t ededxSigmaCut){
 652         switch(ededxSigmaCut){
 653                 case 0:
 654                         fDodEdxSigmaTPCCut = kFALSE;
 655                         fPIDnSigmaBelowPionLineTPC=-10;
 656                         fPIDnSigmaAbovePionLineTPC=10;
 657                         break;
 658                 case 1: // -10,10
 659                         fDodEdxSigmaTPCCut = kTRUE;
 660                         fPIDnSigmaBelowPionLineTPC=-10;
 661                         fPIDnSigmaAbovePionLineTPC=10;
 662                         break;
 663                 case 2: // -6,7
 664                         fDodEdxSigmaTPCCut = kTRUE;
 665                         fPIDnSigmaBelowPionLineTPC=-6;
 666                         fPIDnSigmaAbovePionLineTPC=7;
 667                         break;
 668                 case 3: // -5,5
 669                         fDodEdxSigmaTPCCut = kTRUE;
 670                         fPIDnSigmaBelowPionLineTPC=-5;
 671                         fPIDnSigmaAbovePionLineTPC=5;
 672                         break;
 673                 case 4: // -4,5
 674                         fDodEdxSigmaTPCCut = kTRUE;
 675                         fPIDnSigmaBelowPionLineTPC=-4;
 676                         fPIDnSigmaAbovePionLineTPC=5;
 677                         break;
 678                 case 5: // -3,5
 679                         fDodEdxSigmaTPCCut = kTRUE;
 680                         fPIDnSigmaBelowPionLineTPC=-3;
 681                         fPIDnSigmaAbovePionLineTPC=5;
 682                         break;
 683                 case 6: // -4,4
 684                         fDodEdxSigmaTPCCut = kTRUE;
 685                         fPIDnSigmaBelowPionLineTPC=-4;
 686                         fPIDnSigmaAbovePionLineTPC=4;
 687                         break;
 688                 case 7: // -2.5,4
 689                         fDodEdxSigmaTPCCut = kTRUE;
 690                         fPIDnSigmaBelowPionLineTPC=-2.5;
 691                         fPIDnSigmaAbovePionLineTPC=4;
 692                         break;
 693                 case 8: // -2,3.5
 694                         fDodEdxSigmaTPCCut = kTRUE;
 695                         fPIDnSigmaBelowPionLineTPC=-2;
 696                         fPIDnSigmaAbovePionLineTPC=3.5;
 697                         break;
 698                 default:
 699                         cout<<"Warning: TPCdEdxCutPionLine not defined"<<ededxSigmaCut<<endl;
 700                         return kFALSE;
 701         }
 702         return kTRUE;
 703 }
 704
 705 ///________________________________________________________________________
 706 Bool_t AliPrimaryPionCuts::SetITSClusterCut(Int_t clsITSCut){
 707         if( !fEsdTrackCuts ) {
 708                 cout<<"Warning: AliESDtrackCut is not initialized "<<endl;
 709                 return kFALSE;
 710         }
 711
 712         switch(clsITSCut){
 713                 case 0:
 714                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
 715                         break;
 716                 case 1:
 717                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
 718                         break;  //1 hit first layer of SPD
 719                 case 2:
 720                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 721                         break; //1 hit in any layer of SPD
 722                 case 3:
 723                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
 724                         fEsdTrackCuts->SetMinNClustersITS(4);
 725                         // 4 hits in total in the ITS. At least 1 hit in the first layer of SPD
 726                         break;
 727                 case 4:
 728                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 729                         fEsdTrackCuts->SetMinNClustersITS(3);
 730                         // 3 hits in total in the ITS. At least 1 hit in any layer of SPD
 731                         break;
 732                 case 5:
 733                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 734                         fEsdTrackCuts->SetMinNClustersITS(4);
 735                         // 4 hits in total in the ITS. At least 1 hit in any layer of SPD
 736                         break;
 737                 case 6:
 738                         fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 739                         fEsdTrackCuts->SetMinNClustersITS(5);
 740                         // 5 hits in total in the ITS. At least 1 hit in any layer of SPD
 741                         break;
 742                 default:
 743                         cout<<"Warning: clsITSCut not defined "<<clsITSCut<<endl;
 744                         return kFALSE;
 745         }
 746         return kTRUE;
 747 }
 748
 749 ///________________________________________________________________________
 750 Bool_t AliPrimaryPionCuts::SetTPCClusterCut(Int_t clsTPCCut){
 751         switch(clsTPCCut){
 752                 case 0: // 0
 753                         fMinClsTPC= 0.;
 754                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 755                         break;
 756                 case 1:  // 70
 757                         fMinClsTPC= 70.;
 758                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 759                         break;
 760                 case 2:  // 80
 761                         fMinClsTPC= 80.;
 762                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 763                         break;
 764                 case 3:  // 100
 765                         fMinClsTPC= 100.;
 766                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 767                         break;
 768                 case 4:  // 0% of findable clusters
 769                                 fMinClsTPC= 70.;
 770                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 771                         fMinClsTPCToF= 0.0;
 772                         fUseCorrectedTPCClsInfo=0;
 773                         break;
 774                 case 5:  // 35% of findable clusters
 775                         fMinClsTPC = 70.;
 776                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 777                         fMinClsTPCToF= 0.35;
 778                         fUseCorrectedTPCClsInfo=0;
 779                         break;
 780                 case 6:  // 60% of findable clusters
 781                         fMinClsTPC= 70.;
 782                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 783                         fMinClsTPCToF= 0.6;
 784                         fUseCorrectedTPCClsInfo=0;
 785                         break;
 786                 case 7:  // 70% of findable clusters
 787                         fMinClsTPC= 70.;
 788                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 789                         fMinClsTPCToF= 0.7;
 790                         fUseCorrectedTPCClsInfo=0;
 791                         break;
 792                 case 8: fMinClsTPC = 0.;
 793                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 794                         fMinClsTPCToF= 0.35;
 795                         fUseCorrectedTPCClsInfo=0;
 796                         break;
 797                 case 9:  // 35% of findable clusters
 798                         fMinClsTPC = 70.;
 799                         fEsdTrackCuts->SetMinNClustersTPC(fMinClsTPC);
 800                         fMinClsTPCToF= 0.35;
 801                         fUseCorrectedTPCClsInfo=1;
 802                         break;
 803
 804                 default:
 805                         cout<<"Warning: clsTPCCut not defined "<<clsTPCCut<<endl;
 806                         return kFALSE;
 807         }
 808         return kTRUE;
 809 }
 810
 811 ///________________________________________________________________________
 812 Bool_t AliPrimaryPionCuts::SetEtaCut(Int_t etaCut){
 813         // Set eta Cut
 814         switch(etaCut){
 815                 case 0:
 816                         fEtaCut = 100.;
 817                         fDoEtaCut = kFALSE;
 818                         break;
 819                 case 1: // 1.4
 820                         fEtaCut = 1.4;
 821                         fDoEtaCut = kTRUE;
 822                         break;
 823                 case 2: // 1.2
 824                         fEtaCut = 1.2;
 825                         fDoEtaCut = kTRUE;
 826                         break;
 827                 case 3: // 0.9
 828                         fEtaCut = 0.9;
 829                         fDoEtaCut = kTRUE;
 830                         break;
 831                 case 4: // 0.8
 832                         fEtaCut = 0.8;
 833                         fDoEtaCut = kTRUE;
 834                         break;
 835                 case 5: // 0.75
 836                         fEtaCut = 0.75;
 837                         fDoEtaCut = kTRUE;
 838                         break;
 839                 case 6: //0.6
 840                         fEtaCut = 0.6; //changed from 0.4 to 0.6 2013.06.10
 841                         fDoEtaCut = kTRUE;
 842                         break;
 843                 case 7: //0.5
 844                         fEtaCut = 0.5; //changed from 0.3 to 0.5 2013.06.10
 845                         fDoEtaCut = kTRUE;
 846                         break;
 847                 case 8:
 848                         fEtaCut = 0.4;
 849                         fDoEtaCut = kTRUE;
 850                         break;
 851                 default:
 852                         cout<<"Warning: EtaCut not defined "<<etaCut<<endl;
 853                         return kFALSE;
 854         }
 855         return kTRUE;
 856 }
 857
 858 ///________________________________________________________________________
 859 Bool_t AliPrimaryPionCuts::SetPtCut(Int_t ptCut){
 860         switch(ptCut){
 861                 case 0:
 862                         fPtCut = 0.075;
 863                         break;
 864                 case 1:  // 0.1
 865                         fPtCut  = 0.1;
 866                         break;
 867                 case 2:  // 0.125 GeV
 868                         fPtCut  = 0.125;
 869                         break;
 870                 case 3: // 0.15 GeV
 871                         fPtCut  = 0.15;
 872                         break;
 873                 default:
 874                         cout<<"Warning: PtCut not defined "<<ptCut<<endl;
 875                         return kFALSE;
 876         }
 877         return kTRUE;
 878 }
 879
 880
 881 ///________________________________________________________________________
 882 Bool_t AliPrimaryPionCuts::SetDCACut(Int_t dcaCut)
 883 {
 884         // Set DCA Cut
 885         if( !fEsdTrackCuts ) {
 886                 cout<<"Warning: AliESDtrackCut is not initialized "<<endl;
 887                 return kFALSE;
 888         }
 889
 890         switch(dcaCut){
 891                 case 0:
 892                         //Open cuts//
 893                         fEsdTrackCuts->SetMaxDCAToVertexZ(1000);
 894                         fEsdTrackCuts->SetMaxDCAToVertexXY(1000);
 895                         fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 896                         break;
 897                 case 1:
 898                         fEsdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 899                         fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 900                         break;
 901                 case 2: fEsdTrackCuts->SetMaxDCAToVertexZ(2);
 902                         fEsdTrackCuts->SetMaxDCAToVertexXY(1);
 903                         fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
 904                         break;
 905                 default:
 906                         cout<<"Warning: dcaCut not defined "<<dcaCut<<endl;
 907                         return kFALSE;
 908         }
 909         return kTRUE;
 910 }
 911
 912 ///________________________________________________________________________
 913 Bool_t AliPrimaryPionCuts::SetTOFPionPIDCut(Int_t TOFelectronPID){
 914     // Set Cut
 915         switch(TOFelectronPID){
 916                 case 0: // no cut
 917                         fRequireTOF = kFALSE;
 918                         fUseTOFpid = kFALSE;
 919                         fPIDnSigmaBelowPionLineTOF=-100;
 920                         fPIDnSigmaAbovePionLineTOF=100;
 921                         break;
 922                 case 1: // -7,7
 923                         fRequireTOF = kFALSE;
 924                         fUseTOFpid = kTRUE;
 925                         fPIDnSigmaBelowPionLineTOF=-7;
 926                         fPIDnSigmaAbovePionLineTOF=7;
 927                         break;
 928                 case 2: // -5,5
 929                         fRequireTOF = kFALSE;
 930                         fUseTOFpid = kTRUE;
 931                         fPIDnSigmaBelowPionLineTOF=-5;
 932                         fPIDnSigmaAbovePionLineTOF=5;
 933                         break;
 934                 case 3: // -3,5
 935                         fRequireTOF = kFALSE;
 936                         fUseTOFpid = kTRUE;
 937                         fPIDnSigmaBelowPionLineTOF=-3;
 938                         fPIDnSigmaAbovePionLineTOF=5;
 939                         break;
 940                 case 4: // -2,3
 941                         fRequireTOF = kFALSE;
 942                         fUseTOFpid = kTRUE;
 943                         fPIDnSigmaBelowPionLineTOF=-2;
 944                         fPIDnSigmaAbovePionLineTOF=3;
 945                         break;
 946                 case 5: // -3, 3 TOF mandatory
 947                         fRequireTOF = kTRUE;
 948                         fUseTOFpid  = kTRUE;
 949                         fPIDnSigmaBelowPionLineTOF= -3;
 950                         fPIDnSigmaAbovePionLineTOF=  3;
 951                         break;
 952                 default:
 953                         cout<<"Warning: TOFPionCut not defined "<<TOFelectronPID<<endl;
 954                 return kFALSE;
 955     }
 956     return kTRUE;
 957 }
 958
 959 ///________________________________________________________________________
 960 Bool_t AliPrimaryPionCuts::SetMassCut(Int_t massCut){
 961     // Set Cut
 962         switch(massCut){
 963                 case 0: // no cut
 964                         fDoMassCut = kFALSE;
 965                         fMassCut = 10;
 966                         break;
 967                 case 1: // cut at 1 GeV/c^2
 968                         fDoMassCut = kTRUE;
 969                         fMassCut = 1;
 970                         break;
 971                 case 2: // cut at 0.7 GeV/c^2
 972                         fDoMassCut = kTRUE;
 973                         fMassCut = 0.7;
 974                         break;
 975                 case 3: // cut at 0.6 GeV/c^2
 976                         fDoMassCut = kTRUE;
 977                         fMassCut = 0.6;
 978                         break;
 979                 case 4: // cut at eta mass
 980                         fDoMassCut = kTRUE;
 981                         fMassCut = 0.547853;
 982                         break;
 983                 case 5: // cut at 0.5 GeV/c^2
 984                         fDoMassCut = kTRUE;
 985                         fMassCut = 0.5;
 986                         break;
 987                 default:
 988                         cout<<"Warning: MassCut not defined "<<massCut<<endl;
 989                 return kFALSE;
 990     }
 991     return kTRUE;
 992 }
 993
 994
 995 ///________________________________________________________________________
 996 TString AliPrimaryPionCuts::GetCutNumber(){
 997         // returns TString with current cut number
 998         TString a(kNCuts);
 999         for(Int_t ii=0;ii<kNCuts;ii++){
1000                 a.Append(Form("%d",fCuts[ii]));
1001         }
1002         return a;
1003 }
1004
1005
1006 ///________________________________________________________________________
1007 AliPrimaryPionCuts* AliPrimaryPionCuts::GetStandardCuts2010PbPb(){
1008     //Create and return standard 2010 PbPb cuts
1009     AliPrimaryPionCuts *cuts=new AliPrimaryPionCuts("StandardCuts2010PbPb","StandardCuts2010PbPb");
1010     if(!cuts->InitializeCutsFromCutString("000000400")){
1011                 cout<<"Warning: Initialization of Standardcuts2010PbPb failed"<<endl;
1012         }
1013     return cuts;
1014 }
1015
1016 ///________________________________________________________________________
1017 AliPrimaryPionCuts* AliPrimaryPionCuts::GetStandardCuts2010pp(){
1018     //Create and return standard 2010 PbPb cuts
1019     AliPrimaryPionCuts *cuts=new AliPrimaryPionCuts("StandardCuts2010pp","StandardCuts2010pp");
1020
1021     if(!cuts->InitializeCutsFromCutString("000000400")){
1022                 cout<<"Warning: Initialization of Standardcuts2010pp failed"<<endl;
1023         }
1024     return cuts;
1025 }
1026