testing patch for kaon flow
authorrbertens <rbertens@cern.ch>
Tue, 13 Jan 2015 14:30:31 +0000 (15:30 +0100)
committerrbertens <rbertens@cern.ch>
Tue, 13 Jan 2015 14:30:54 +0000 (15:30 +0100)
PWG/FLOW/Tasks/AliFlowTrackCuts.cxx
PWG/FLOW/Tasks/AliFlowTrackCuts.h
PWGCF/FLOW/macros/AddTaskPIDFlowSP.C

index 853b6da..3a26302 100644 (file)
@@ -35,6 +35,7 @@
 // AliFlowTrack is made using MakeFlowTrack() method, its an 'object factory'
 // so caller needs to take care of the freshly created object.
 
+#include "TFile.h"
 #include <limits.h>
 #include <float.h>
 #include <TMatrix.h>
@@ -116,6 +117,8 @@ AliFlowTrackCuts::AliFlowTrackCuts():
   fCutMinimalTPCdedx(kFALSE),
   fMinimalTPCdedx(0.),
   fLinearizeVZEROresponse(kFALSE),
+  fCentralityPercentileMin(0.),
+  fCentralityPercentileMax(5.),
   fCutPmdDet(kFALSE),
   fPmdDet(0),
   fCutPmdAdc(kFALSE),
@@ -156,6 +159,7 @@ AliFlowTrackCuts::AliFlowTrackCuts():
   fAllowTOFmismatchFlag(kFALSE),
   fRequireStrictTOFTPCagreement(kFALSE),
   fCutRejectElectronsWithTPCpid(kFALSE),
+  fUseTPCTOFNsigmaCutContours(kFALSE),
   fProbBayes(0.0),
   fCurrCentr(0.0),
   fVZEROgainEqualization(NULL),
@@ -166,11 +170,14 @@ AliFlowTrackCuts::AliFlowTrackCuts():
   fChi3A(0x0),
   fChi3C(0x0),
   fPIDResponse(NULL),
-  fNsigmaCut2(9)
+  fNsigmaCut2(9),
+  fContoursFile(0),
+  fCutContourList(0),
+  fMaxITSclusterShared(0),
+  fMaxITSChi2(0)
 {
   //io constructor 
   SetPriors(); //init arrays
-  
   // New PID procedure (Bayesian Combined PID)
   // allocating here is necessary because we don't 
   // stream this member
@@ -181,7 +188,12 @@ AliFlowTrackCuts::AliFlowTrackCuts():
       fVZEROApol[i] = 0;
       fVZEROCpol[i] = 0;
   }
-  for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = kTRUE;
+  for(Int_t i(0); i < 8; i++){ fUseVZERORing[i] = kTRUE;}
+    
+  for(int i=0;i<50;i++){
+    fCutContour[i]= NULL;
+    fCutGraph[i]=NULL;
+  }
 
 }
 
@@ -227,6 +239,8 @@ AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
   fCutMinimalTPCdedx(kFALSE),
   fMinimalTPCdedx(0.),
   fLinearizeVZEROresponse(kFALSE),
+  fCentralityPercentileMin(0.),
+  fCentralityPercentileMax(5.),
   fCutPmdDet(kFALSE),
   fPmdDet(0),
   fCutPmdAdc(kFALSE),
@@ -267,6 +281,7 @@ AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
   fAllowTOFmismatchFlag(kFALSE),
   fRequireStrictTOFTPCagreement(kFALSE),
   fCutRejectElectronsWithTPCpid(kFALSE),
+  fUseTPCTOFNsigmaCutContours(kFALSE),
   fProbBayes(0.0),
   fCurrCentr(0.0),
   fVZEROgainEqualization(NULL),
@@ -277,9 +292,13 @@ AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
   fChi3A(0x0),
   fChi3C(0x0),
   fPIDResponse(NULL),
-  fNsigmaCut2(9)
+  fNsigmaCut2(9),
+  fContoursFile(0),
+  fCutContourList(0),
+  fMaxITSclusterShared(0),
+  fMaxITSChi2(0)
 {
-  //constructor 
+  //constructor
   SetTitle("AliFlowTrackCuts");
   fESDpid.GetTPCResponse().SetBetheBlochParameters( 0.0283086,
                                                     2.63394e+01,
@@ -287,7 +306,6 @@ AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
                                                     2.12543e+00,
                                                     4.88663e+00 );
   SetPriors(); //init arrays
-
   // New PID procedure (Bayesian Combined PID)
   fBayesianResponse = new AliFlowBayesianPID();
   fBayesianResponse->SetNewTrackParam();
@@ -296,7 +314,7 @@ AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
       fVZEROCpol[i] = 0;
   }
   for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = kTRUE;
-}
+  }
 
 //-----------------------------------------------------------------------
 AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
@@ -340,6 +358,8 @@ AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
   fCutMinimalTPCdedx(that.fCutMinimalTPCdedx),
   fMinimalTPCdedx(that.fMinimalTPCdedx),
   fLinearizeVZEROresponse(that.fLinearizeVZEROresponse),
+  fCentralityPercentileMin(that.fCentralityPercentileMin),
+  fCentralityPercentileMax(that.fCentralityPercentileMax),
   fCutPmdDet(that.fCutPmdDet),
   fPmdDet(that.fPmdDet),
   fCutPmdAdc(that.fCutPmdAdc),
@@ -380,6 +400,7 @@ AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
   fAllowTOFmismatchFlag(that.fAllowTOFmismatchFlag),
   fRequireStrictTOFTPCagreement(that.fRequireStrictTOFTPCagreement),
   fCutRejectElectronsWithTPCpid(that.fCutRejectElectronsWithTPCpid),
+  fUseTPCTOFNsigmaCutContours(that.fUseTPCTOFNsigmaCutContours),
   fProbBayes(0.0),
   fCurrCentr(0.0),
   fVZEROgainEqualization(NULL),
@@ -390,7 +411,11 @@ AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
   fChi3A(0x0),
   fChi3C(0x0),
   fPIDResponse(that.fPIDResponse),
-  fNsigmaCut2(that.fNsigmaCut2)
+  fNsigmaCut2(that.fNsigmaCut2),
+  fContoursFile(0),
+  fCutContourList(0),
+  fMaxITSclusterShared(0),
+  fMaxITSChi2(0)
 {
   //copy constructor
   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
@@ -398,6 +423,7 @@ AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
   if (that.fAliESDtrackCuts) fAliESDtrackCuts = new AliESDtrackCuts(*(that.fAliESDtrackCuts));
   if (that.fMuonTrackCuts)   fMuonTrackCuts   = new AliMuonTrackCuts(*(that.fMuonTrackCuts));  // XZhang 20120604
   SetPriors(); //init arrays
+  if(fUseTPCTOFNsigmaCutContours) GetTPCTOFPIDContours();
   if (that.fQA) DefineHistograms();
   // would be neater via copy ctor of TArrayD
   fChi2A = new TArrayD(that.fChi2A->GetSize(), that.fChi2A->GetArray());
@@ -407,7 +433,7 @@ AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
   // New PID procedure (Bayesian Combined PID)
   fBayesianResponse = new AliFlowBayesianPID();
   fBayesianResponse->SetNewTrackParam();
-
+    
   // VZERO gain calibration
   // no reason to init fVZEROgainEqualizationPerRing, will be initialized on node if necessary
   // pointer is set to NULL in initialization list of this constructor
@@ -417,6 +443,12 @@ AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
       fVZEROCpol[i] = 0.;
   }
   for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = that.fUseVZERORing[i];
+  
+    for(int i=0;i<50;i++){
+        fCutContour[i]= that.fCutContour[i];
+        fCutGraph[i]=that.fCutGraph[i];
+    }
+
 
 }
 
@@ -475,6 +507,8 @@ AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
   fCutMinimalTPCdedx=that.fCutMinimalTPCdedx;
   fMinimalTPCdedx=that.fMinimalTPCdedx;
   fLinearizeVZEROresponse=that.fLinearizeVZEROresponse;
+  fCentralityPercentileMin=that.fCentralityPercentileMin;
+  fCentralityPercentileMax=that.fCentralityPercentileMax;
   fCutPmdDet=that.fCutPmdDet;
   fPmdDet=that.fPmdDet;
   fCutPmdAdc=that.fCutPmdAdc;
@@ -520,9 +554,10 @@ AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
   fAllowTOFmismatchFlag=that.fAllowTOFmismatchFlag;
   fRequireStrictTOFTPCagreement=that.fRequireStrictTOFTPCagreement;
   fCutRejectElectronsWithTPCpid=that.fCutRejectElectronsWithTPCpid;
+  fUseTPCTOFNsigmaCutContours=that.fUseTPCTOFNsigmaCutContours;
   fProbBayes = that.fProbBayes;
   fCurrCentr = that.fCurrCentr;
-  
+    
   fApplyRecentering = that.fApplyRecentering;
   fVZEROgainEqualizationPerRing = that.fVZEROgainEqualizationPerRing;
 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)           
@@ -531,6 +566,7 @@ AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
   //PH Lets try Clone, however the result might be wrong
   if (that.fVZEROgainEqualization) fVZEROgainEqualization = (TH1*)that.fVZEROgainEqualization->Clone();
 #endif
+    
   for(Int_t i(0); i < 4; i++) { // no use to copy these guys since they're only initialized on worked node
       fVZEROApol[i] = that.fVZEROApol[i];
       fVZEROCpol[i] = that.fVZEROCpol[i];
@@ -577,6 +613,11 @@ AliFlowTrackCuts::~AliFlowTrackCuts()
       delete fChi3C;
       fChi3C = 0x0;
   }
+  fContoursFile->Close();
+  for(int i=0;i<50;i++){
+     delete fCutContour[i];
+     delete fCutGraph[i];
+  }
 }
 
 //-----------------------------------------------------------------------
@@ -620,6 +661,8 @@ void AliFlowTrackCuts::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
   
   if(fPIDsource==kTOFbayesian) fBayesianResponse->SetDetAND(1);
   else if(fPIDsource==kTPCbayesian) fBayesianResponse->ResetDetOR(1);
+    
+  if(fUseTPCTOFNsigmaCutContours) GetTPCTOFPIDContours();
 }
 
 //-----------------------------------------------------------------------
@@ -1240,7 +1283,25 @@ Bool_t AliFlowTrackCuts::PassesAODcuts(const AliAODTrack* track, Bool_t passedFi
     {
       if (!PassesAODpidCut(track)) pass=kFALSE;
     }
-
+/*
+  if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
+     Double_t c = TMath::C()*1.E-9;
+     Double_t length = -999., beta =-999, tofTime = -999., tof = -999.;
+     tofTime = track->GetTOFsignal();//in ps
+     length = track->GetIntegratedLength();
+        
+     tof = tofTime*1E-3; // ns
+     if (tof > 0 && length > 0){
+        length = length*0.01; // in meters
+         tof = tof*c;
+         beta = length/tof;
+         if (fQA){
+             QAbefore(18)->Fill(track->P()*track->Charge(),beta);
+             if(pass) QAafter(18)->Fill(track->P()*track->Charge(),beta);
+         }
+     }
+  }*/
   if (fQA) {
     // changed 04062014 used to be filled before possible PID cut
     Double_t momTPC = track->GetTPCmomentum();
@@ -1262,6 +1323,10 @@ Bool_t AliFlowTrackCuts::PassesAODcuts(const AliAODTrack* track, Bool_t passedFi
     if (pass) QAafter( 11)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kKaon]));
     QAbefore(12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
     if (pass) QAafter( 12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
+      
+    //QAbefore(19)->Fill(track->P()*track->Charge(),track->GetDetPid()->GetTPCsignal());
+    //if(pass) QAafter(19)->Fill(track->P()*track->Charge(),track->GetDetPid()->GetTPCsignal());
+      
   }
 
 
@@ -1278,6 +1343,7 @@ Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
   track->GetImpactParameters(dcaxy,dcaz);
   const AliExternalTrackParam* pout = track->GetOuterParam();
   const AliExternalTrackParam* pin = track->GetInnerParam();
+  
   if (fIgnoreTPCzRange)
   {
     if (pin&&pout)
@@ -1343,6 +1409,13 @@ Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
     track->GetTPCpid(pidTPC);
     if (pidTPC[AliPID::kElectron]<fParticleProbability) pass=kFALSE;
   }
+    
+  Int_t counterForSharedCluster=MaxSharedITSClusterCuts(fMaxITSclusterShared, track);
+  if(counterForSharedCluster >= fMaxITSclusterShared) pass=kFALSE;
+    
+  Double_t chi2perClusterITS = MaxChi2perITSClusterCuts(fMaxITSChi2, track);
+  if(chi2perClusterITS >= fMaxITSChi2) pass=kFALSE;
+    
   if (fQA)
   {
     if (pass) QAafter(0)->Fill(track->GetP(),beta);
@@ -1427,6 +1500,7 @@ Int_t AliFlowTrackCuts::Count(AliVEvent* event)
   //calculate the number of track in given event.
   //if argument is NULL(default) take the event attached
   //by SetEvent()
+    
   Int_t multiplicity = 0;
   if (!event)
   {
@@ -2475,6 +2549,13 @@ void AliFlowTrackCuts::DefineHistograms()
 
   before->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
   after->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
+   /*
+  before->Add(new TH2F("Pionbeta",";p_{t}[GeV/c];#beta",kNbinsP,binsP,1000,0,1000 ));//18
+  after->Add(new TH2F("Pionbeta",";p_{t}[GeV/c];#beta",kNbinsP,binsP,1000,0,1000 ));//18
+    
+  before->Add(new TH2F("PiondEdX",";p_{t}[GeV/c];dE/dX",kNbinsP,binsP,1000,0,1000 ));//19
+  after->Add(new TH2F("PiondEdX",";p_{t}[GeV/c];dE/dX",kNbinsP,binsP,1000,0,1000 ));//19
+*/
   TH1::AddDirectory(adddirstatus);
 }
 
@@ -2632,6 +2713,9 @@ Bool_t AliFlowTrackCuts::PassesAODpidCut(const AliAODTrack* track )
   case kTPCTOFNsigma:
       if (!PassesTPCTOFNsigmaCut(track)) pass = kFALSE;
       break;
+  case kTPCTOFNsigmaPuritybased:
+      if(!PassesTPCTOFNsigmaCutPuritybased(track)) pass = kFALSE;
+      break;
   default:
     return kTRUE;
     break;
@@ -2792,6 +2876,7 @@ Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliAODTrack* track )
 
   //construct the pid index because it's not AliPID::EParticleType
   Int_t pid = 0;
+    cout<<"TOFbeta: fParticleID = "<<fParticleID<<endl;
   switch (fParticleID)
   {
     case AliPID::kPion:
@@ -3579,7 +3664,76 @@ Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaCut(const AliESDtrack* track)
     return (nsigma2 < fNsigmaCut2);
 
 }
+//-----------------------------------------------------------------------------
+Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaCutPuritybased(const AliAODTrack* track)
+{
+    // do a combined cut on the n sigma from tpc and tof based on purity of the identified particle. (Particle is selected if Purity>0.8)
+    // with information of the pid response object (needs pid response task)
 
+    if(!fPIDResponse) return kFALSE;
+    if(!track) return kFALSE;
+    if(!fUseTPCTOFNsigmaCutContours) return kFALSE;
+
+    
+    Bool_t pWithinRange = kFALSE;
+    Int_t p_bin = -999;
+    Double_t pBins[50];
+    for(int b=0;b<50;b++){pBins[b] = 0.1*b;}
+    for(int i=0;i<50;i++){
+        if(track->P()>pBins[i] && track->P()<(pBins[i+1])){
+            p_bin = i;
+        }
+    }
+    Int_t ispecie = 0;
+    switch (fParticleID)
+    {
+        case AliPID::kPion:
+            ispecie=0;
+            break;
+        case AliPID::kKaon:
+            ispecie=1;
+            break;
+        case AliPID::kProton:
+            ispecie=2;
+            break;
+        default:
+            return kFALSE;
+    }
+    
+    Double_t LowPtPIDTPCnsigLow_Pion[6] = {-3,-3,-3,-3,-3,-3};
+    Double_t LowPtPIDTPCnsigLow_Kaon[6] = {-3,-2,0,-1.8,-1.2,-0.8}; //for0.4<Pt<0.5 the purity is lower than 0.7
+    Double_t LowPtPIDTPCnsigHigh_Pion[6] ={2.4,3,3,3,2,1.4};
+    Double_t LowPtPIDTPCnsigHigh_Kaon[6] ={3,2.2,0,-0.2,1,1.8}; //for 0.4<Pt<0.5 he purity is lower than 0.7
+    
+    Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
+    Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
+
+    
+    if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid))){
+        if(p_bin<0) return kFALSE;
+        if(!fCutContour[p_bin]) cout<<"fCutContour[p_bin] does not exist"<<endl;
+        
+        if(p_bin==4 || p_bin>7){
+            if(fCutContour[p_bin]->IsInside(nSigmaTOF,nSigmaTPC)){
+                return kTRUE;
+            }
+            else{
+                return kFALSE;
+            }
+        }
+        
+        if(p_bin<8 && p_bin!=4){
+            if(fParticleID==AliPID::kPion && nSigmaTPC>LowPtPIDTPCnsigLow_Pion[p_bin-2] && nSigmaTPC<LowPtPIDTPCnsigHigh_Pion[p_bin-2]){
+                return kTRUE;
+            }
+            if(fParticleID==AliPID::kKaon && nSigmaTPC>LowPtPIDTPCnsigLow_Kaon[p_bin-2] && nSigmaTPC<LowPtPIDTPCnsigHigh_Kaon[p_bin-2]){return kTRUE;}
+            
+            if(fParticleID==AliPID::kProton && nSigmaTPC>-3 && nSigmaTPC<3){return kTRUE;}
+        }
+    }
+    
+    return kFALSE;
+}
 //-----------------------------------------------------------------------
 void AliFlowTrackCuts::SetPriors(Float_t centrCur){
  //set priors for the bayesian pid selection
@@ -4718,6 +4872,8 @@ const char* AliFlowTrackCuts::PIDsourceName(PIDsource s)
       return "TPCnuclei";
     case kTPCTOFNsigma:
       return "TPCTOFNsigma";
+    case kTPCTOFNsigmaPuritybased:
+      return "TPCTOFNsigmaPuritybased";
     default:
       return "NOPID";
   }
@@ -4923,6 +5079,85 @@ Long64_t AliFlowTrackCuts::Merge(TCollection* list)
   }
   return fQA->Merge(&tmplist);
 }
+//________________________________________________________________//
+void AliFlowTrackCuts::GetTPCTOFPIDContours()
+{
+    fContoursFile = TFile::Open(Form("$ALICE_ROOT/PWGCF/FLOW/database/PIDCutContours_%i-%i.root",fCentralityPercentileMin,fCentralityPercentileMax));
+    fCutContourList=(TDirectory*)fContoursFile->Get("Filterbit1");
+    if(!fCutContourList){printf("The contour file is empty"); return;}
+
+    Double_t pBinning[50];
+    for(int b=0;b<50;b++){pBinning[b]=b;}
+    TString species[3] = {"pion","kaon","proton"};
+    TList *Species_contours[3];
+    Int_t ispecie = 0;
+    for(ispecie = 0; ispecie < 3; ispecie++) {
+       Species_contours[ispecie] = (TList*)fCutContourList->Get(species[ispecie]);
+        if(!Species_contours[ispecie]) {
+            cout<<"Contours for species: "<<species[ispecie]<<" not found!!!"<<endl;
+            return;
+        }
+    }
+    Int_t iParticle = 0;
+
+    switch (fParticleID)
+    {
+        case AliPID::kPion:
+            iParticle=0;
+            break;
+        case AliPID::kKaon:
+            iParticle=1;
+            break;
+        case AliPID::kProton:
+            iParticle=2;
+            break;
+        default:
+            return;
+    }
+
+    
+    for(int i=0;i<50;i++){
+        TString Graph_Name = "contourlines_";
+        Graph_Name += species[iParticle];
+        Graph_Name += Form("%.f%.f-%i%icent",pBinning[i],pBinning[i]+1,fCentralityPercentileMin,fCentralityPercentileMax);
+        fCutGraph[i] = (TGraph*)Species_contours[iParticle]->FindObject(Graph_Name);
+        
+        if(!fCutGraph[i]){cout<<"Contour Graph does not exist"<<endl; continue;}
+        
+        fCutContour[i] = new TCutG(Graph_Name.Data(),fCutGraph[i]->GetN(),fCutGraph[i]->GetX(),fCutGraph[i]->GetY());
+        
+    }
+    
+}
+//__________________
+Int_t AliFlowTrackCuts::MaxSharedITSClusterCuts(Int_t maxITSclusterShared, AliESDtrack* track){
+    
+    Int_t counterForSharedCluster = 0;
+    for(int i =0;i<6;i++){
+        Bool_t sharedITSCluster = track->HasSharedPointOnITSLayer(i);
+        Bool_t PointOnITSLayer = track->HasPointOnITSLayer(i);
+        //   cout << "has a point???    " << PointOnITSLayer << "     is it shared or not??? "  << sharedITSCluster << endl;
+        if(sharedITSCluster == 1) counterForSharedCluster++;
+    }
+    //if(counterForSharedCluster >= maxITSclusterShared) pass=kFALSE;
+    
+    return counterForSharedCluster;
+
+}
+//___________________
+Double_t AliFlowTrackCuts::MaxChi2perITSClusterCuts(Double_t maxITSChi2, AliESDtrack* track){
+    
+    // cout << "  Total Shared Cluster ==   "  <<counterForSharedCluster<< endl;
+    Double_t chi2perClusterITS = track->GetITSchi2()/track->GetNcls(0);
+    //  cout << "  chi2perClusterITS ==   "  <<chi2perClusterITS <<endl;
+    
+    
+    //if(chi2perClusterITS >= maxITSChi2) pass=kFALSE;
+
+    return chi2perClusterITS;
+}
+
 
 
 
index 29b210a..b7c4285 100644 (file)
@@ -20,6 +20,8 @@
 #include "AliMuonTrackCuts.h"  // XZhang 20120604
 #include "AliPID.h"
 #include "AliESDpid.h"
+#include "TCutG.h"
+
 
 class TBrowser;
 class TArrayD;
@@ -91,10 +93,17 @@ class AliFlowTrackCuts : public AliFlowTrackSimpleCuts {
                    kTOFbetaSimple, //simple TOF only cut
                    kTPCbayesian, //bayesian cutTPC
                   kTPCNuclei,   // added by Natasha for Nuclei
-                   kTPCTOFNsigma // simple cut on combined tpc tof nsigma
+                   kTPCTOFNsigma, // simple cut on combined tpc tof nsigma
+                   kTPCTOFNsigmaPuritybased // purity>0.7 cut on combined tpc tof nsigma
                    };
 
   //setters (interface to AliESDtrackCuts)
+  Int_t MaxSharedITSClusterCuts(Int_t maxITSclusterShared, AliESDtrack* track);
+  Double_t MaxChi2perITSClusterCuts(Double_t maxITSChi2, AliESDtrack* track);
+
+  void SetMaxSharedITSCluster(Int_t b){fMaxITSclusterShared = b;}
+  void SetMaxChi2perITSCluster(Double_t b){fMaxITSChi2 = b;}
+    
   void SetMinNClustersTPC( Int_t a ) {fCutNClustersTPC=kTRUE; fNClustersTPCMin=a;}
   void SetMinNClustersITS( Int_t a ) {fCutNClustersITS=kTRUE; fNClustersITSMin=a;}
   void SetClusterRequirementITS( AliESDtrackCuts::Detector det,
@@ -239,6 +248,7 @@ class AliFlowTrackCuts : public AliFlowTrackSimpleCuts {
   //PID
   void SetPID(AliPID::EParticleType pid, PIDsource s=kTOFpid, Double_t prob=0.9)
              {fParticleID=pid; fPIDsource=s; fParticleProbability=prob; fCutPID=kTRUE; InitPIDcuts();}
+  void SetTPCTOFNsigmaPIDCutContours(Bool_t b,Double_t cMin,Double_t cMax){fUseTPCTOFNsigmaCutContours = b; fCentralityPercentileMin=cMin; fCentralityPercentileMax=cMax;}
   AliPID::EParticleType GetParticleID() const {return fParticleID;}
   Bool_t GetCutPID() const {return fCutPID;}
   void SetTPCpidCuts(const TMatrixF* mat) {fTPCpidCuts=new TMatrixF(*mat);}
@@ -280,7 +290,7 @@ class AliFlowTrackCuts : public AliFlowTrackSimpleCuts {
 
   void Browse(TBrowser* b);
   Long64_t Merge(TCollection* list);
-  
+  void GetTPCTOFPIDContours();
   //gain equalization and recentering
   void SetVZEROgainEqualisation(TH1* g) {fVZEROgainEqualization=g;}
   void SetVZEROApol(Int_t ring, Float_t f) {fVZEROApol[ring]=f;}
@@ -333,9 +343,10 @@ class AliFlowTrackCuts : public AliFlowTrackSimpleCuts {
   Bool_t PassesNucleiSelection(const AliESDtrack* track);   // added by Natasha
   Bool_t PassesTPCTOFNsigmaCut(const AliAODTrack* track); 
   Bool_t PassesTPCTOFNsigmaCut(const AliESDtrack* track);
+  Bool_t PassesTPCTOFNsigmaCutPuritybased(const AliAODTrack* track);
   Bool_t TPCTOFagree(const AliVTrack *track);
   // end part added by F. Noferini
-
+    
   //the cuts
   AliESDtrackCuts* fAliESDtrackCuts; //alianalysis cuts
   AliMuonTrackCuts* fMuonTrackCuts;  // muon selection cuts // XZhang 20120604
@@ -377,6 +388,9 @@ class AliFlowTrackCuts : public AliFlowTrackSimpleCuts {
   Double_t fMinimalTPCdedx;       //value for minimal TPC dedx
   Bool_t fLinearizeVZEROresponse; //linearize VZERO response using AliESDUtil
   
+  Int_t fCentralityPercentileMin; //centrality min
+  Int_t fCentralityPercentileMax; //centrality max
+    
   Bool_t  fCutPmdDet;   //cut on PMD detector plane 
   Int_t   fPmdDet;      // value of PMD detector plane
   Bool_t  fCutPmdAdc;   //cut on cluster ADC
@@ -422,6 +436,7 @@ class AliFlowTrackCuts : public AliFlowTrackSimpleCuts {
   Bool_t fAllowTOFmismatchFlag; //allow TOFmismatch flag=1 in ESD
   Bool_t fRequireStrictTOFTPCagreement; //require stricter than TOFmismatch flag TOF-TPC agreement
   Bool_t fCutRejectElectronsWithTPCpid; //reject electrons with TPC pid
+  Bool_t fUseTPCTOFNsigmaCutContours;  // Flag to use purity based TPC-TOF nsigma cuts
 
   // part added by F. Noferini
   static const Int_t fgkPIDptBin = 20; // pT bins for priors
@@ -445,6 +460,24 @@ class AliFlowTrackCuts : public AliFlowTrackSimpleCuts {
 
   AliPIDResponse *fPIDResponse;            //! Pid reponse to manage Nsigma cuts
   Float_t fNsigmaCut2;                     // Number of sigma^2 (cut value) for TPC+TOF nsigma cut
+    
+    
+  //TPC TOF nsigma Purity based cut contours
+  TFile                 *fContoursFile;       //! contours file
+  TDirectory            *fCutContourList;     //! contour list
+  
+  Int_t  fMaxITSclusterShared;
+  Double_t  fMaxITSChi2;
+  
+    
+  TCutG                 *fCutContour[50];    //! TCutG contours
+  //TCutG                 *fCutContour_kaon[50];    //! TCutG contours
+  //TCutG                 *fCutContour_proton[50];    //! TCutG contours
+  TGraph                *fCutGraph[50];      //! graphs
+  //TGraph                *fCutGraph_kaon[50];      //! graphs
+  //TGraph                *fCutGraph_proton[50];      //! graphs
+    
+    
 
   ClassDef(AliFlowTrackCuts,15)
 };
index 259e382..a716a3b 100644 (file)
@@ -17,10 +17,14 @@ void AddTaskPIDFlowSP(Int_t triggerSelectionString=AliVEvent::kMB,
                                    Int_t charge=0,
                                    Int_t MinTPCdedx = 10,
                                    Int_t ncentrality = 6,
+                                   Int_t maxITSCls = 7,
+                                   Double_t maxChi2ITSCls = 37,
                                    Bool_t doQA=kTRUE,
                                    Bool_t isPID = kTRUE,
                                    Bool_t VZERO = kFALSE, // use vzero sp method
                                    Bool_t is2011 = kFALSE,
+                                   Bool_t isAOD = kTRUE,
+                                   Bool_t UsePIDParContours = kFALSE,
                                    AliPID::EParticleType particleType=AliPID::kPion,
                                    AliFlowTrackCuts::PIDsource sourcePID=AliFlowTrackCuts::kTOFbayesian) {
 
@@ -115,12 +119,25 @@ for(int icentr=0;icentr<ncentr;icentr++){
     SP_POI[icentr] = DefinePOIcuts(icentr);
 
     SP_POI[icentr]->GetBayesianResponse()->ForceOldDedx(); // for 2010 data to use old TPC PID Response instead of the official one
-    //SP_POI[icentr]->SetParamType(poitype);
+    if(!isAOD){
+        SP_POI[icentr]->SetMaxSharedITSCluster(maxITSCls);
+        SP_POI[icentr]->SetMaxChi2perITSCluster(maxChi2ITSCls);
+        SP_POI[icentr]->SetMinNClustersITS(2);
+        SP_POI[icentr]->SetRequireITSRefit(kTRUE);
+        SP_POI[icentr]->SetRequireTPCRefit(kTRUE);
+        SP_POI[icentr]->SetMaxDCAToVertexXY(0.3);
+        SP_POI[icentr]->SetMaxDCAToVertexZ(0.3);
+        SP_POI[icentr]->SetAcceptKinkDaughters(kFALSE);
+        SP_POI[icentr]->SetMinimalTPCdedx(10.);
+    }
     //SP_POI[icentr]->SetParamMix(poimix);
-    SP_POI[icentr]->SetPtRange(0.2,5.);//
+    SP_POI[icentr]->SetPtRange(0.2,6.);//
     SP_POI[icentr]->SetMinNClustersTPC(70);
     SP_POI[icentr]->SetMinChi2PerClusterTPC(0.1);
     SP_POI[icentr]->SetMaxChi2PerClusterTPC(4.0);
+    
+    
     if(!VZERO && Qvector=="Qa"){
         SP_POI[icentr]->SetEtaRange( +0.5*EtaGap, etamax );
         printf(" > NOTE: Using half TPC (Qa) as POI selection u < \n");
@@ -145,12 +162,16 @@ for(int icentr=0;icentr<ncentr;icentr++){
     //SP_POI->SetMaxNsigmaToVertex(1.e+10);
     //SP_POI->SetRequireSigmaToVertex(kFALSE);
     SP_POI[icentr]->SetAcceptKinkDaughters(kFALSE);
-    if(isPID) SP_POI[icentr]->SetPID(particleType, sourcePID);//particleType, sourcePID
+    if(isPID){
+        SP_POI[icentr]->SetPID(particleType, sourcePID);//particleType, sourcePID
+        SP_POI[icentr]->SetTPCTOFNsigmaPIDCutContours(UsePIDParContours,centrMin[icentr],centrMax[icentr]);
+    }
+    
     if (charge!=0) SP_POI[icentr]->SetCharge(charge);
     //SP_POI->SetAllowTOFmismatch(kFALSE);
     SP_POI[icentr]->SetRequireStrictTOFTPCagreement(kTRUE);
     SP_POI[icentr]->SetMinimalTPCdedx(MinTPCdedx);
-    SP_POI[icentr]->SetAODfilterBit(AODfilterBit);
+    if(isAOD) SP_POI[icentr]->SetAODfilterBit(AODfilterBit);
     SP_POI[icentr]->SetQA(doQA);
     SP_POI[icentr]->SetPriors((centrMin[icentr]+centrMax[icentr])*0.5);
 
@@ -170,7 +191,9 @@ for(int icentr=0;icentr<ncentr;icentr++){
     if(isPID){
         suffixName[icentr]+=AliFlowTrackCuts::PIDsourceName(sourcePID);
         suffixName[icentr]+="_";
-        suffixName[icentr]+=AliPID::ParticleName(particleType);//particleType
+        suffixName[icentr]+=uniqueStr;
+        suffixName[icentr]+="_";
+        //suffixName[icentr]+=AliPID::ParticleName(particleType);//particleType
     }
     else{
         suffixName[icentr]+="AllCharged";
@@ -342,11 +365,11 @@ AliFlowTrackCuts* DefineRPcuts(Int_t icentr){
     AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts(Form("RP_%d",icentr));
     return cutsRP;
 }
-AliFlowTrackCuts* DefinePOIcuts(Int_t icentr){
+AliFlowTrackCuts* DefinePOIcuts(Int_t    icentr){
     AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts(Form("POI_%d",icentr));
     return cutsPOI;
 }
 
-
+