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