PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C

   1 // now in options
   2 //=============================================//
   3 //const char* centralityEstimator = "V0M";
   4 //const char* centralityEstimator = "CL1";
   5 //const char* centralityEstimator = "TRK";
   6 //=============================================//
   7 //Bool_t gRunShuffling = kFALSE;
   8 //Bool_t gRunShuffling = kTRUE;
   9 //=============================================//
  10
  11 //PID config
  12 Bool_t kUseNSigmaPID = kTRUE;
  13 Double_t nSigmaMax = 3.0;
  14 Bool_t kUseBayesianPID = kFALSE;
  15 Double_t gMinAcceptedProbability = 0.7;
  16
  17 //_________________________________________________________//
  18 AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,
  19                                                        Double_t centrMax=100.,
  20                                                        Bool_t gRunShuffling=kFALSE,
  21                                                        Bool_t gRunMixing=kTRUE,
  22                                                        Bool_t gRunMixingWithEventPlane=kFALSE,
  23                                                        TString centralityEstimator="V0M",
  24                                                        Double_t vertexZ=10.,
  25                                                        Double_t DCAxy=-1,
  26                                                        Double_t DCAz=-1,
  27                                                        Double_t ptMin=0.3,
  28                                                        Double_t ptMax=1.5,
  29                                                        Double_t etaMin=-0.8,
  30                                                        Double_t etaMax=0.8,
  31                                                        Double_t maxTPCchi2 = -1,
  32                                                        Int_t minNClustersTPC = -1,
  33                                                        Bool_t kUsePID = kTRUE,
  34                                                        Bool_t bResonancesCut = kTRUE,
  35                                                        Bool_t bHBTcut = kTRUE,
  36                                                        Double_t HBTCutValue = 0.02,
  37                                                        Bool_t bConversionCut = kTRUE,
  38                                                        Double_t invMassForConversionCut = 0.04,
  39                                                        Bool_t bMomentumDifferenceCut = kTRUE,
  40                                                        Double_t fQCutMin = 0.0,
  41                                                        Int_t AODfilterBit = 128,
  42                                                        Bool_t bCentralTrigger = kFALSE,
  43                                                        TString fileNameBase="AnalysisResults",
  44                                                        TString dirNameExtra="",
  45                                                        TString fArgEventClass="Centrality",
  46                                                        TString analysisTypeUser="AOD",
  47                                                        Bool_t bVertexBinning=kTRUE,
  48                                                        Double_t sigmaElectronRejection=3,
  49                                                        Bool_t electronExclusiveRejection=kFALSE,
  50                                                        TString correctionFileName = "",
  51                                                        Int_t nCentralityArrayBinsForCorrection,
  52                                                        Double_t *gCentralityArrayForCorrections) {
  53   // Creates a balance function analysis task and adds it to the analysis manager.
  54   // Get the pointer to the existing analysis manager via the static access method.
  55   TString outputFileName(fileNameBase);
  56   outputFileName.Append(".root");
  57
  58   //===========================================================================
  59   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  60   if (!mgr) {
  61     ::Error("AddTaskBF", "No analysis manager to connect to.");
  62     return NULL;
  63   }
  64
  65   // Check the analysis type using the event handlers connected to the analysis manager.
  66   //===========================================================================
  67   if (!mgr->GetInputEventHandler()) {
  68     ::Error("AddTaskBF", "This task requires an input event handler");
  69     return NULL;
  70   }
  71   TString analysisType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
  72   if(dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())) analysisType = "MC";
  73
  74   // to set the analysis type manually
  75   if(analysisTypeUser != ""){
  76     analysisType = analysisTypeUser;
  77     ::Info("AddTaskBF",Form("Analysis Type manually set to %s",analysisType.Data()));
  78   }
  79
  80   // for local changed BF configuration
  81   //gROOT->LoadMacro("./configBalanceFunctionPsiAnalysis.C");
  82   gROOT->LoadMacro("$ALICE_ROOT/PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C");
  83   AliBalancePsi *bf  = 0;  // Balance Function object
  84   AliBalancePsi *bfs = 0;  // shuffled Balance function object
  85   AliBalancePsi *bfm = 0;  // mixing Balance function object
  86
  87   //maximum Delta eta range
  88   Double_t deltaEtaMax=TMath::Abs(etaMax-etaMin);
  89
  90   if (analysisType=="ESD"){
  91     bf  = GetBalanceFunctionObject("ESD",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
  92     if(gRunShuffling) bfs = GetBalanceFunctionObject("ESD",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
  93     if(gRunMixing)    bfm = GetBalanceFunctionObject("ESD",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
  94   }
  95   else if (analysisType=="AOD"){
  96     bf  = GetBalanceFunctionObject("AOD",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
  97     if(gRunShuffling) bfs = GetBalanceFunctionObject("AOD",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
  98     if(gRunMixing)    bfm = GetBalanceFunctionObject("AOD",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
  99   }
 100   else if (analysisType=="MC"){
 101     bf  = GetBalanceFunctionObject("MC",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 102     if(gRunShuffling) bfs = GetBalanceFunctionObject("MC",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 103     if(gRunMixing)    bfm = GetBalanceFunctionObject("MC",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 104   }
 105   else if (analysisType=="MCAOD"){
 106     bf  = GetBalanceFunctionObject("MCAOD",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 107     if(gRunShuffling) bfs = GetBalanceFunctionObject("MCAOD",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 108     if(gRunMixing)    bfm = GetBalanceFunctionObject("MCAOD",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 109   }
 110   else if (analysisType=="MCAODrec"){
 111     bf  = GetBalanceFunctionObject("MCAODrec",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 112     if(gRunShuffling) bfs = GetBalanceFunctionObject("MCAODrec",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 113     if(gRunMixing)    bfm = GetBalanceFunctionObject("MCAODrec",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 114   }
 115   else if (analysisType=="AODnano"){
 116     bf  = GetBalanceFunctionObject("AODnano",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 117     if(gRunShuffling) bfs = GetBalanceFunctionObject("AODnano",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 118     if(gRunMixing)    bfm = GetBalanceFunctionObject("AODnano",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,HBTCutValue,bConversionCut,invMassForConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);
 119   }
 120   else{
 121     ::Error("AddTaskBF", "analysis type NOT known.");
 122     return NULL;
 123   }
 124
 125   // Create the task, add it to manager and configure it.
 126   //===========================================================================
 127   AliAnalysisTaskBFPsi *taskBF = new AliAnalysisTaskBFPsi(Form("TaskBFPsi_%.0f-%.0f_Bit%d_%s%s",centrMin,centrMax,AODfilterBit,centralityEstimator.Data(),dirNameExtra.Data()));
 128
 129   //Event characteristics scheme
 130   taskBF->SetEventClass(fArgEventClass);
 131   //taskBF->SetCustomBinning("centralityVertex:0,80");
 132   //taskBF->SetCustomBinning("multiplicity:0,260");
 133
 134   if(fArgEventClass == "Multiplicity") {
 135     taskBF->SetMultiplicityRange(centrMin,centrMax);
 136     taskBF->SetMultiplicityEstimator(centralityEstimator);
 137     cout<<"Multiplicity estimator "<<centralityEstimator.Data()<<endl;
 138   }
 139   else if(fArgEventClass == "Centrality") {
 140     if(analysisType == "MC")
 141       taskBF->SetImpactParameterRange(centrMin,centrMax);
 142     else {
 143       taskBF->SetCentralityPercentileRange(centrMin,centrMax);
 144       // centrality estimator (default = V0M)
 145       taskBF->SetCentralityEstimator(centralityEstimator);
 146       cout<<"Centrality estimator "<<centralityEstimator.Data()<<endl;
 147     }
 148   }
 149
 150   //++++++++++++++++++++++
 151   // Efficiency + Contamination corrections
 152   // If correctionFileName = "", do not use corrections
 153   if(correctionFileName != "")
 154     taskBF->SetInputCorrection(Form("$ALICE_ROOT/PWGCF/EBYE/BalanceFunctions/Corrections/%s",correctionFileName.Data()),nCentralityArrayBinsForCorrection,gCentralityArrayForCorrections);
 155
 156   //+++++++++++++++++++++
 157
 158   taskBF->SetAnalysisObject(bf);
 159   if(gRunShuffling) taskBF->SetShufflingObject(bfs);
 160   if(gRunMixing){
 161     taskBF->SetMixingObject(bfm);
 162     taskBF->SetMixingTracks(50000);
 163     if(gRunMixingWithEventPlane){
 164       taskBF->SetMixingWithEventPlane(gRunMixingWithEventPlane);
 165     }
 166   }
 167
 168   if(analysisType == "ESD") {
 169     AliESDtrackCuts *trackCuts = GetTrackCutsObject(ptMin,ptMax,etaMin,etaMax,maxTPCchi2,DCAxy,DCAz,minNClustersTPC);
 170     taskBF->SetAnalysisCutObject(trackCuts);
 171     if(kUsePID) {
 172       if(kUseBayesianPID)
 173         taskBF->SetUseBayesianPID(gMinAcceptedProbability);
 174       else if(kUseNSigmaPID)
 175         taskBF->SetUseNSigmaPID(nSigmaMax);
 176       taskBF->SetParticleOfInterest(AliAnalysisTaskBFPsi::kPion);
 177       taskBF->SetDetectorUsedForPID(AliAnalysisTaskBFPsi::kTOFpid);
 178     }
 179   }
 180   else if(analysisType == "AOD" || analysisType == "AODnano") {
 181     // pt and eta cut (pt_min, pt_max, eta_min, eta_max)
 182     taskBF->SetAODtrackCutBit(AODfilterBit);
 183     taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);
 184
 185     // set extra DCA cuts (-1 no extra cut)
 186     taskBF->SetExtraDCACutsAOD(DCAxy,DCAz);
 187
 188     // set extra TPC chi2 / nr of clusters cut
 189     taskBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC);
 190
 191     // electron rejection (so far only for AOD), <0 --> no rejection
 192     if(sigmaElectronRejection > 0){
 193       if(electronExclusiveRejection) taskBF->SetElectronOnlyRejection(sigmaElectronRejection); // no other particle in nsigma
 194       else                           taskBF->SetElectronRejection(sigmaElectronRejection); // check only if electrons in nsigma
 195     }
 196
 197     //++++++++++++++++//
 198     if(kUsePID) {
 199       if(kUseBayesianPID)
 200         taskBF->SetUseBayesianPID(gMinAcceptedProbability);
 201       else if(kUseNSigmaPID)
 202         taskBF->SetUseNSigmaPID(nSigmaMax);
 203       taskBF->SetParticleOfInterest(AliAnalysisTaskBFPsi::kKaon);
 204       taskBF->SetDetectorUsedForPID(AliAnalysisTaskBFPsi::kTPCTOF); //TOFpid,TPCpid
 205     }
 206     //++++++++++++++++//
 207
 208   }
 209   else if(analysisType == "MC") {
 210     taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);
 211   }
 212   else if(analysisType == "MCAOD") {
 213     // pt and eta cut (pt_min, pt_max, eta_min, eta_max)
 214     taskBF->SetAODtrackCutBit(AODfilterBit);
 215     taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);
 216   }
 217   else if(analysisType == "MCAODrec") {     //++++++++++++++++
 218     // pt and eta cut (pt_min, pt_max, eta_min, eta_max)
 219     taskBF->SetAODtrackCutBit(AODfilterBit);
 220     taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);
 221
 222     // set extra DCA cuts (-1 no extra cut)
 223     taskBF->SetExtraDCACutsAOD(DCAxy,DCAz);
 224
 225     // set extra TPC chi2 / nr of clusters cut
 226     taskBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC);
 227
 228     // electron rejection (so far only for AOD), <0 --> no rejection
 229     if(sigmaElectronRejection > 0){
 230       if(electronExclusiveRejection) taskBF->SetElectronOnlyRejection(sigmaElectronRejection); // no other particle in nsigma
 231       else                           taskBF->SetElectronRejection(sigmaElectronRejection); // check only if electrons in nsigma
 232     }
 233   }//++++++++++++++++
 234
 235   // offline trigger selection (AliVEvent.h)
 236   // taskBF->UseOfflineTrigger(); // NOT used (selection is done with the AliAnalysisTaskSE::SelectCollisionCandidates())
 237   // with this only selected events are analyzed (first 2 bins in event QA histogram are the same))
 238   // documentation in https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PWG1EvSelDocumentation
 239   if(bCentralTrigger) taskBF->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral);
 240   else                taskBF->SelectCollisionCandidates(AliVEvent::kMB);
 241
 242   // centrality estimator (default = V0M)
 243   taskBF->SetCentralityEstimator(centralityEstimator);
 244
 245   // vertex cut (x,y,z)
 246   taskBF->SetVertexDiamond(3.,3.,vertexZ);
 247
 248
 249
 250   //bf->PrintAnalysisSettings();
 251   mgr->AddTask(taskBF);
 252
 253   // Create ONLY the output containers for the data produced by the task.
 254   // Get and connect other common input/output containers via the manager as below
 255   //==============================================================================
 256   TString outputFileName = AliAnalysisManager::GetCommonFileName();
 257   outputFileName += ":PWGCFEbyE.outputBalanceFunctionPsiAnalysis";
 258   AliAnalysisDataContainer *coutQA = mgr->CreateContainer(Form("listQAPsi_%.0f-%.0f_Bit%d_%s%s",centrMin,centrMax,AODfilterBit,centralityEstimator.Data(),dirNameExtra.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());
 259   AliAnalysisDataContainer *coutBF = mgr->CreateContainer(Form("listBFPsi_%.0f-%.0f_Bit%d_%s%s",centrMin,centrMax,AODfilterBit,centralityEstimator.Data(),dirNameExtra.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());
 260   if(gRunShuffling) AliAnalysisDataContainer *coutBFS = mgr->CreateContainer(Form("listBFPsiShuffled_%.0f-%.0f_Bit%d_%s%s",centrMin,centrMax,AODfilterBit,centralityEstimator.Data(),dirNameExtra.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());
 261   if(gRunMixing) AliAnalysisDataContainer *coutBFM = mgr->CreateContainer(Form("listBFPsiMixed_%.0f-%.0f_Bit%d_%s%s",centrMin,centrMax,AODfilterBit,centralityEstimator.Data(),dirNameExtra.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());
 262   if(kUsePID || sigmaElectronRejection > 0) AliAnalysisDataContainer *coutQAPID = mgr->CreateContainer(Form("listQAPIDPsi_%.0f-%.0f_Bit%d_%s%s",centrMin,centrMax,AODfilterBit,centralityEstimator.Data(),dirNameExtra.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());
 263
 264   mgr->ConnectInput(taskBF, 0, mgr->GetCommonInputContainer());
 265   mgr->ConnectOutput(taskBF, 1, coutQA);
 266   mgr->ConnectOutput(taskBF, 2, coutBF);
 267   if(gRunShuffling) mgr->ConnectOutput(taskBF, 3, coutBFS);
 268   if(gRunMixing) mgr->ConnectOutput(taskBF, 4, coutBFM);
 269   if(kUsePID||sigmaElectronRejection > 0) mgr->ConnectOutput(taskBF, 5, coutQAPID);
 270   //if((kUsePID && analysisType == "AOD")||sigmaElectronRejection > 0) mgr->ConnectOutput(taskBF, 5, coutQAPID);
 271   //if((kUsePID && analysisType == "ESD")||sigmaElectronRejection > 0) mgr->ConnectOutput(taskBF, 5, coutQAPID);
 272
 273   return taskBF;
 274 }