PWGLF/RESONANCES/macros/train_old/AddRsnAnalysisTaskEffPhi.C

   1 //
   2 // This macro add an analysis task for computing efficiency.
   3 // It will have as output an AliCFContainer with several steps:
   4 //
   5 //  0) all resonances in MC which decay in the pair specified
   6 //  1) subset of (0) whose daughters are in acceptance
   7 //  2) subset of (1) whose daughters satisfy quality track cuts (covariance, chi square && nTPCclusters)
   8 //  3) subset of (2) whose daughters satisfy primary track cuts (nsigma to vertex, no kink daughters)
   9 //  4) subset of (3) whose daughters satisty the BB TPC compatibility cut at 3 sigma
  10 //
  11 Bool_t AddRsnAnalysisTaskEffPhi
  12 (
  13    const char *evtopts,
  14    const char *options,
  15    const char *path = "$(ALICE_ROOT)/PWG2/RESONANCES/macros/train"
  16 )
  17 {
  18    // ==================================================================================================================
  19    // == OPTIONS =======================================================================================================
  20    // ==================================================================================================================
  21
  22    // Instead or getting confused with plenty of arguments in the macro (with default values),
  23    // we use a unique string of options with a set of conventional strings to set up the job:
  24    // -- "MC"/"DATA" --> what kind of sample
  25    // -- "ITS"/"TPC" --> what tracks to use (ITS standalone and/or TPC+ITS)
  26    // -- "xxxPID"    --> add the PID cuts for the detector xxx.
  27    //
  28    // In this point, these options are converted into boolean variables.
  29
  30    TString opt(options);
  31    opt.ToUpper();
  32    opt.ReplaceAll(" ", "");
  33
  34    Bool_t addITS = opt.Contains("ITS");
  35    Bool_t addTPC = opt.Contains("TPC");
  36    Bool_t useITS = opt.Contains("ITSPID");
  37    Bool_t useTPC = opt.Contains("TPCPID");
  38    Bool_t useTOF = opt.Contains("TOFPID");
  39
  40    // correct options when needed
  41    if (!addITS) useITS = kFALSE;
  42    if (!addTPC) useTPC = useTOF = kFALSE;
  43
  44    // ==================================================================================================================
  45    // == DEFINITIONS ===================================================================================================
  46    // ==================================================================================================================
  47
  48    // We put here the definitions of all objects which are needed in the following, in order to have then
  49    // a more readable code in the points where these objects are added to the analysis manager.
  50
  51    // pair definition:
  52    // phi --> K+ K-
  53    AliRsnPairDef *pairPhi = new AliRsnPairDef(AliPID::kKaon, '+', AliPID::kKaon, '-', 333, 1.019455);
  54
  55    // ==================================================================================================================
  56    // == CUTS AND AXES =================================================================================================
  57    // ==================================================================================================================
  58
  59    //
  60    // Track quality for ITS standalone:
  61    // this cut is used to select tracks of good quality, irrespective of the PID.
  62    // When adding status flags, the second argument tells if each considered flag
  63    // must be active or not in the track status, since the ITS-SA tracks need that
  64    // some of them are OFF (e.g.: kTPCin)
  65    //
  66    AliRsnCutTrackQuality *cutQualityITS = new AliRsnCutTrackQuality("cutQualityITS");
  67    cutQualityITS->AddStatusFlag(AliESDtrack::kITSin    , kTRUE);
  68    cutQualityITS->AddStatusFlag(AliESDtrack::kTPCin    , kFALSE);
  69    cutQualityITS->AddStatusFlag(AliESDtrack::kITSrefit , kTRUE);
  70    cutQualityITS->AddStatusFlag(AliESDtrack::kTPCrefit , kFALSE);
  71    cutQualityITS->AddStatusFlag(AliESDtrack::kITSpureSA, kFALSE);
  72    cutQualityITS->AddStatusFlag(AliESDtrack::kITSpid   , kTRUE);
  73    cutQualityITS->SetPtRange(0.15, 1E+20);
  74    cutQualityITS->SetEtaRange(-0.8, 0.8);
  75    cutQualityITS->SetDCARPtFormula("0.0595+0.0182/pt^1.55");
  76    cutQualityITS->SetDCAZmax(2.0);
  77    cutQualityITS->SetSPDminNClusters(1);
  78    cutQualityITS->SetITSminNClusters(4);
  79    cutQualityITS->SetITSmaxChi2(2.0);
  80    cutQualityITS->SetTPCminNClusters(0);
  81    cutQualityITS->SetTPCmaxChi2(1E+10);
  82    cutQualityITS->SetRejectKinkDaughters();
  83
  84    //
  85    // Track quality for TPC+ITS:
  86    // works exactly like the one above, but has settings for selecting TPC+ITS tracks
  87    // in this case, the flags required are all necessary, so here the procedure is simpler
  88    //
  89    AliRsnCutTrackQuality *cutQualityTPC = new AliRsnCutTrackQuality("cutQualityTPC");
  90    cutQualityTPC->AddStatusFlag(AliESDtrack::kTPCin   , kTRUE);
  91    cutQualityTPC->AddStatusFlag(AliESDtrack::kTPCrefit, kTRUE);
  92    cutQualityTPC->AddStatusFlag(AliESDtrack::kITSrefit, kTRUE);
  93    cutQualityTPC->SetPtRange(0.15, 1E+20);
  94    cutQualityTPC->SetEtaRange(-0.8, 0.8);
  95    cutQualityTPC->SetDCARPtFormula("0.0182+0.0350/pt^1.01");
  96    cutQualityTPC->SetDCAZmax(2.0);
  97    cutQualityTPC->SetSPDminNClusters(1);
  98    cutQualityTPC->SetITSminNClusters(0);
  99    cutQualityTPC->SetITSmaxChi2(1E+20);
 100    cutQualityTPC->SetTPCminNClusters(70);
 101    cutQualityTPC->SetTPCmaxChi2(4.0);
 102    cutQualityTPC->SetRejectKinkDaughters();
 103
 104    //
 105    // ITS PID
 106    // In this implementation, it is a 3sigma cut around the Bethe-Bloch value.
 107    //
 108    // NOTE:
 109    // --> The initialization of the BB is different between data and MC.
 110    // --> The cut is the same for all momenta.
 111    //
 112    AliRsnCutPIDITS *cutPIDITSkaon = new AliRsnCutPIDITS("cutPIDITSkaon", AliPID::kKaon, -3.0, 3.0);
 113    cutPIDITSkaon->SetMC(kTRUE);
 114    cutPIDITSkaon->SetMomentumRange(0.0, 1E20);
 115
 116    //
 117    // TPC PID
 118    // In this implementation, there are two instances:
 119    // - below 350 MeV --> 5 sigma cut
 120    // - above 350 MeV --> 3 sigma cut
 121    //
 122    // NOTE:
 123    // --> The initialization of the BB is different between data and MC.
 124    //
 125    AliRsnCutPIDTPC *cutPIDTPCkaonLow  = new AliRsnCutPIDTPC("cutPIDTPCkaonLow" , AliPID::kKaon, -5.0, 5.0);
 126    AliRsnCutPIDTPC *cutPIDTPCkaonHigh = new AliRsnCutPIDTPC("cutPIDTPCkaonHigh", AliPID::kKaon, -3.0, 3.0);
 127
 128    // assign the momentum range and tell to reject tracks outside it
 129    cutPIDTPCkaonLow ->SetMomentumRange(0.00, 0.35);
 130    cutPIDTPCkaonHigh->SetMomentumRange(0.35, 1E20);
 131    cutPIDTPCkaonLow ->SetRejectOutside(kTRUE);
 132    cutPIDTPCkaonHigh->SetRejectOutside(kTRUE);
 133
 134    // BB parameterization depends on data sample (MC, data)
 135    // the momentum range is passed and tracks outside it are rejected
 136    Double_t bbPar[5];
 137    bbPar[0] = 2.15898 / 50.0;
 138    bbPar[1] = 1.75295E1;
 139    bbPar[2] = 3.40030E-9;
 140    bbPar[3] = 1.96178;
 141    bbPar[4] = 3.91720;
 142    cutPIDTPCkaonLow ->SetBBParam(bbPar);
 143    cutPIDTPCkaonHigh->SetBBParam(bbPar);
 144
 145    //
 146    // TOF PID
 147    // In this implementation it is a 3sigma cout aroung expected kaon time.
 148    //
 149    // NOTE:
 150    // --> It is important to choose if this cut must reject tracks not matched in TOF.
 151    //     Usually, if TPC pid is used, we can accept them, otherwise we must reject.
 152    //     (here we assume TPC is used)
 153    //
 154    AliRsnCutPIDTOF *cutPIDTOFkaon = new AliRsnCutPIDTOF("cutPIDTOFkaon", AliPID::kKaon, -3.0, 3.0);
 155    cutPIDTOFkaon->SetRejectUnmatched(!useTPC);
 156
 157    //
 158    // Rapidity cut
 159    // Only thing to consider is that it needs a support object to define mass
 160    //
 161    AliRsnCutValue *cutRapidity = new AliRsnCutValue("cutY", AliRsnValue::kPairY, -0.5, 0.5);
 162    cutRapidity->GetValueObj()->SetSupportObject(pairPhi);
 163
 164    //
 165    // Axes
 166    //
 167    // NOTE:
 168    // --> multiplicity has variable bins, defined by array below
 169    //
 170
 171    Double_t mult[] = { 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13.,  14.,  15.,  16.,  17.,  18.,  19.,
 172                       20., 21., 22., 23., 24., 25., 30., 35., 40., 50., 60., 70., 80., 90., 100., 120., 140., 160., 180., 200., 500.};
 173    Int_t    nmult  = sizeof(mult) / sizeof(mult[0]);
 174
 175    AliRsnValue *axisIM      = new AliRsnValue("IM"  , AliRsnValue::kPairInvMass   ,  0.9, 1.4, 0.001);
 176    AliRsnValue *axisRes     = new AliRsnValue("Res" , AliRsnValue::kPairInvMassRes, -0.5, 0.5, 0.001);
 177    AliRsnValue *axisPt      = new AliRsnValue("PT"  , AliRsnValue::kPairPt        ,  0.0, 5.0, 0.1  );
 178    AliRsnValue *axisY       = new AliRsnValue("Y"   , AliRsnValue::kPairY         , -1.1, 1.1, 0.1  );
 179    AliRsnValue *axisMultESD = new AliRsnValue("MESD", AliRsnValue::kEventMultESDCuts, nmult, mult);
 180    AliRsnValue *axisMultSPD = new AliRsnValue("MSPD", AliRsnValue::kEventMultSPD    , nmult, mult);
 181    AliRsnValue *axisMultMC  = new AliRsnValue("MMC" , AliRsnValue::kEventMultMC     , nmult, mult);
 182
 183    // ==================================================================================================================
 184    // == PRELIMINARY OPERATIONS ========================================================================================
 185    // ==================================================================================================================
 186
 187    // retrieve analysis manager
 188    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 189
 190    // create task
 191    AliRsnAnalysisTaskEffPair *task = new AliRsnAnalysisTaskEffPair(Form("RsnEff_%s", options));
 192    task->SelectCollisionCandidates();
 193
 194    // add pair definition, to choose the checked resonance
 195    task->AddDef(pairPhi);
 196
 197    // add the output histogram axis
 198    task->AddAxis(axisIM);
 199    task->AddAxis(axisRes);
 200    task->AddAxis(axisPt);
 201    task->AddAxis(axisY);
 202    task->AddAxis(axisMultSPD);
 203    task->AddAxis(axisMultMC);
 204
 205    // ==================================================================================================================
 206    // == EVENT CUTS ====================================================================================================
 207    // ==================================================================================================================
 208
 209    gROOT->LoadMacro(Form("%s/AddRsnEventComputations.C", path));
 210    AddRsnEventComputations(kTRUE, evtopts);
 211
 212    // ==================================================================================================================
 213    // == STEPS =========================================================================================================
 214    // ==================================================================================================================
 215
 216    Char_t qualityITS[255], qualityTPC[255];
 217    Char_t pidITS[255], pidTPC[255], pidTOF[255];
 218    Char_t schemeITS[255], schemeTPC[255], scheme[255];
 219
 220    sprintf(qualityITS, "%s"     , cutQualityITS->GetName());
 221    sprintf(qualityTPC, "%s"     , cutQualityTPC->GetName());
 222    sprintf(pidITS    , "%s"     , cutPIDITSkaon->GetName());
 223    sprintf(pidTPC    , "(%s|%s)", cutPIDTPCkaonHigh->GetName(), cutPIDTPCkaonLow->GetName());
 224    sprintf(pidTOF    , "%s"     , cutPIDTOFkaon->GetName());
 225    sprintf(schemeITS , "");
 226    sprintf(schemeTPC , "");
 227    sprintf(scheme    , "");
 228
 229    //
 230    // *** STEP 0 - All resonances which decay in the specified pair
 231    //
 232    // This step does not need any kind of definition, since
 233    // its requirement is automatically checked during execution,
 234    // but to avoid segfaults, it is needed to initialize a cut manager.
 235    //
 236    AliRsnCutManager *mgr_step0 = new AliRsnCutManager("mc_step0", "");
 237
 238    //
 239    // *** STEP 1 - Track quality
 240    //
 241    // All resonances whose daughters were reconstructed
 242    // and pass quality track cuts will enter this step
 243    //
 244    AliRsnCutManager *mgr_step1 = new AliRsnCutManager("rec_step1", "");
 245    AliRsnCutSet     *set_step1 = mgr_step1->GetCommonDaughterCuts();
 246
 247    if (addTPC && addITS) {
 248       set_step1->AddCut(cutQualityTPC);
 249       set_step1->AddCut(cutQualityITS);
 250       set_step1->SetCutScheme(Form("%s|%s", qualityTPC, qualityITS));
 251    } else if (addTPC) {
 252       set_step1->AddCut(cutQualityTPC);
 253       set_step1->SetCutScheme(qualityTPC);
 254    } else if (addITS) {
 255       set_step1->AddCut(cutQualityITS);
 256       set_step1->SetCutScheme(qualityITS);
 257    } else {
 258       ::Error("Need to ad at least one between ITS and TPC tracks");
 259       return kFALSE;
 260    }
 261    ::Info("AddRsnAnalysisTaskEffPhi", "Cut scheme for step #1: %s", set_step1->GetCutScheme().Data());
 262
 263    //
 264    // *** STEP 2 - PID
 265    //
 266    // Add all TPC cuts, according to options
 267    //
 268    AliRsnCutManager *mgr_step2 = new AliRsnCutManager("esd_step2", "");
 269    AliRsnCutSet     *set_step2 = mgr_step2->GetCommonDaughterCuts();
 270
 271    if (addITS && useITS) {
 272       sprintf(schemeITS, "%s & %s", qualityITS, pidITS);
 273       set_step2->AddCut(cutPIDITSkaon);
 274    }
 275    if (addTPC) {
 276       if (useTPC && useTOF) {
 277          set_step2->AddCut(cutPIDTPCkaonLow);
 278          set_step2->AddCut(cutPIDTPCkaonHigh);
 279          set_step2->AddCut(cutPIDTOFkaon);
 280          sprintf(schemeTPC, "%s & %s", pidTPC, pidTOF);
 281       } else if (useTPC) {
 282          set_step2->AddCut(cutPIDTPCkaonLow);
 283          set_step2->AddCut(cutPIDTPCkaonHigh);
 284          sprintf(schemeTPC, "%s & %s", pidTPC);
 285       } else if (useTOF) {
 286          set_step2->AddCut(cutPIDTOFkaon);
 287          sprintf(schemeTPC, "%s & %s", pidTOF);
 288       }
 289    }
 290
 291    // final scheme depends on what of the above were added
 292    // in case both ITS-SA and TPC tracks are added, the scheme
 293    // is the OR of the cuts for the first and those for the second
 294    // category of tracks
 295    if (strlen(schemeITS) > 0 && strlen(schemeTPC) > 0) {
 296       sprintf(scheme, "(%s) | (%s)", schemeITS, schemeTPC);
 297    } else if (strlen(schemeITS) > 0) {
 298       sprintf(scheme, "%s", schemeITS);
 299    } else if (strlen(schemeTPC) > 0) {
 300       sprintf(scheme, "%s", schemeTPC);
 301    } else {
 302       ::Error("Scheme is empty!");
 303       return kFALSE;
 304    }
 305    // check scheme
 306    set_step2->SetCutScheme(scheme);
 307    ::Info("AddRsnAnalysisTaskEffPhi", "Cut scheme for step #1: %s", set_step2->GetCutScheme().Data());
 308
 309    // add all steps to the task:
 310    // - first step computed on MC
 311    // - all other steps computed on reconstruction
 312    task->AddStepMC(mgr_step0);
 313    task->AddStepRec(mgr_step1);
 314    task->AddStepRec(mgr_step2);
 315
 316    // add the task to the manager and connect to input
 317    mgr->AddTask(task);
 318    mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
 319
 320    // create paths for the output in the common file
 321    TString infoname(task->GetName());
 322    TString histname(task->GetName());
 323    infoname.Append("_Info");
 324    histname.Append("_Hist");
 325    AliAnalysisDataContainer *outputInfo = mgr->CreateContainer(infoname.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName());
 326    AliAnalysisDataContainer *outputHist = mgr->CreateContainer(histname.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName());
 327    mgr->ConnectOutput(task, 1, outputInfo);
 328    mgr->ConnectOutput(task, 2, outputHist);
 329
 330    return kTRUE;
 331 }