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{
 116     ::Error("AddTaskBF", "analysis type NOT known.");
 117     return NULL;
 118   }
 119
 120   // Create the task, add it to manager and configure it.
 121   //===========================================================================
 122   AliAnalysisTaskBFPsi *taskBF = new AliAnalysisTaskBFPsi(Form("TaskBFPsi_%.0f-%.0f_Bit%d_%s%s",centrMin,centrMax,AODfilterBit,centralityEstimator.Data(),dirNameExtra.Data()));
 123
 124   //Event characteristics scheme
 125   taskBF->SetEventClass(fArgEventClass);
 126   //taskBF->SetCustomBinning("centralityVertex:0,80");
 127   //taskBF->SetCustomBinning("multiplicity:0,260");
 128
 129   if(fArgEventClass == "Multiplicity") {
 130     taskBF->SetMultiplicityRange(centrMin,centrMax);
 131     taskBF->SetMultiplicityEstimator(centralityEstimator);
 132     cout<<"Multiplicity estimator "<<centralityEstimator.Data()<<endl;
 133   }
 134   else if(fArgEventClass == "Centrality") {
 135     if(analysisType == "MC")
 136       taskBF->SetImpactParameterRange(centrMin,centrMax);
 137     else {
 138       taskBF->SetCentralityPercentileRange(centrMin,centrMax);
 139       // centrality estimator (default = V0M)
 140       taskBF->SetCentralityEstimator(centralityEstimator);
 141       cout<<"Centrality estimator "<<centralityEstimator.Data()<<endl;
 142     }
 143   }
 144
 145   //++++++++++++++++++++++
 146   // Efficiency + Contamination corrections
 147   // If correctionFileName = "", do not use corrections
 148   if(correctionFileName != "")
 149     taskBF->SetInputCorrection(Form("$ALICE_ROOT/PWGCF/EBYE/BalanceFunctions/Corrections/%s",correctionFileName.Data()),nCentralityArrayBinsForCorrection,gCentralityArrayForCorrections);
 150
 151   //+++++++++++++++++++++
 152
 153   taskBF->SetAnalysisObject(bf);
 154   if(gRunShuffling) taskBF->SetShufflingObject(bfs);
 155   if(gRunMixing){
 156     taskBF->SetMixingObject(bfm);
 157     taskBF->SetMixingTracks(50000);
 158     if(gRunMixingWithEventPlane){
 159       taskBF->SetMixingWithEventPlane(gRunMixingWithEventPlane);
 160     }
 161   }
 162
 163   if(analysisType == "ESD") {
 164     AliESDtrackCuts *trackCuts = GetTrackCutsObject(ptMin,ptMax,etaMin,etaMax,maxTPCchi2,DCAxy,DCAz,minNClustersTPC);
 165     taskBF->SetAnalysisCutObject(trackCuts);
 166     if(kUsePID) {
 167       if(kUseBayesianPID)
 168         taskBF->SetUseBayesianPID(gMinAcceptedProbability);
 169       else if(kUseNSigmaPID)
 170         taskBF->SetUseNSigmaPID(nSigmaMax);
 171       taskBF->SetParticleOfInterest(AliAnalysisTaskBFPsi::kPion);
 172       taskBF->SetDetectorUsedForPID(AliAnalysisTaskBFPsi::kTOFpid);
 173     }
 174   }
 175   else if(analysisType == "AOD") {
 176     // pt and eta cut (pt_min, pt_max, eta_min, eta_max)
 177     taskBF->SetAODtrackCutBit(AODfilterBit);
 178     taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);
 179
 180     // set extra DCA cuts (-1 no extra cut)
 181     taskBF->SetExtraDCACutsAOD(DCAxy,DCAz);
 182
 183     // set extra TPC chi2 / nr of clusters cut
 184     taskBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC);
 185
 186     // electron rejection (so far only for AOD), <0 --> no rejection
 187     if(sigmaElectronRejection > 0){
 188       if(electronExclusiveRejection) taskBF->SetElectronOnlyRejection(sigmaElectronRejection); // no other particle in nsigma
 189       else                           taskBF->SetElectronRejection(sigmaElectronRejection); // check only if electrons in nsigma
 190     }
 191
 192     //++++++++++++++++//
 193     if(kUsePID) {
 194       if(kUseBayesianPID)
 195         taskBF->SetUseBayesianPID(gMinAcceptedProbability);
 196       else if(kUseNSigmaPID)
 197         taskBF->SetUseNSigmaPID(nSigmaMax);
 198       taskBF->SetParticleOfInterest(AliAnalysisTaskBFPsi::kPion);
 199       taskBF->SetDetectorUsedForPID(AliAnalysisTaskBFPsi::kTPCTOF); //TOFpid,TPCpid
 200     }
 201     //++++++++++++++++//
 202
 203   }
 204   else if(analysisType == "MC") {
 205     taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);
 206   }
 207   else if(analysisType == "MCAOD") {
 208     // pt and eta cut (pt_min, pt_max, eta_min, eta_max)
 209     taskBF->SetAODtrackCutBit(AODfilterBit);
 210     taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);
 211   }
 212   else if(analysisType == "MCAODrec") {     //++++++++++++++++
 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     // set extra DCA cuts (-1 no extra cut)
 218     taskBF->SetExtraDCACutsAOD(DCAxy,DCAz);
 219
 220     // set extra TPC chi2 / nr of clusters cut
 221     taskBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC);
 222
 223     // electron rejection (so far only for AOD), <0 --> no rejection
 224     if(sigmaElectronRejection > 0){
 225       if(electronExclusiveRejection) taskBF->SetElectronOnlyRejection(sigmaElectronRejection); // no other particle in nsigma
 226       else                           taskBF->SetElectronRejection(sigmaElectronRejection); // check only if electrons in nsigma
 227     }
 228   }//++++++++++++++++
 229
 230   // offline trigger selection (AliVEvent.h)
 231   // taskBF->UseOfflineTrigger(); // NOT used (selection is done with the AliAnalysisTaskSE::SelectCollisionCandidates())
 232   // with this only selected events are analyzed (first 2 bins in event QA histogram are the same))
 233   // documentation in https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PWG1EvSelDocumentation
 234   if(bCentralTrigger) taskBF->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral);
 235   else                taskBF->SelectCollisionCandidates(AliVEvent::kMB);
 236
 237   // centrality estimator (default = V0M)
 238   taskBF->SetCentralityEstimator(centralityEstimator);
 239
 240   // vertex cut (x,y,z)
 241   taskBF->SetVertexDiamond(3.,3.,vertexZ);
 242
 243
 244
 245   //bf->PrintAnalysisSettings();
 246   mgr->AddTask(taskBF);
 247
 248   // Create ONLY the output containers for the data produced by the task.
 249   // Get and connect other common input/output containers via the manager as below
 250   //==============================================================================
 251   TString outputFileName = AliAnalysisManager::GetCommonFileName();
 252   outputFileName += ":PWGCFEbyE.outputBalanceFunctionPsiAnalysis";
 253   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());
 254   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());
 255   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());
 256   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());
 257   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());
 258
 259   mgr->ConnectInput(taskBF, 0, mgr->GetCommonInputContainer());
 260   mgr->ConnectOutput(taskBF, 1, coutQA);
 261   mgr->ConnectOutput(taskBF, 2, coutBF);
 262   if(gRunShuffling) mgr->ConnectOutput(taskBF, 3, coutBFS);
 263   if(gRunMixing) mgr->ConnectOutput(taskBF, 4, coutBFM);
 264   if(kUsePID||sigmaElectronRejection > 0) mgr->ConnectOutput(taskBF, 5, coutQAPID);
 265   //if((kUsePID && analysisType == "AOD")||sigmaElectronRejection > 0) mgr->ConnectOutput(taskBF, 5, coutQAPID);
 266   //if((kUsePID && analysisType == "ESD")||sigmaElectronRejection > 0) mgr->ConnectOutput(taskBF, 5, coutQAPID);
 267
 268   return taskBF;
 269 }