PWGLF/RESONANCES/macros/train_old/RsnConfigPhi.C

   1 //
   2 // *** Configuration script for phi->KK analysis with 2010 runs ***
   3 //
   4 // A configuration script for RSN package needs to define the followings:
   5 //
   6 // (1) decay tree of each resonance to be studied, which is needed to select
   7 //     true pairs and to assign the right mass to all candidate daughters
   8 // (2) cuts at all levels: single daughters, tracks, events
   9 // (3) output objects: histograms or trees
  10 //
  11 Bool_t RsnConfigPhi
  12 (
  13    Bool_t      isMC,
  14    const char *options,
  15    const char *path     = "$(ALICE_ROOT)/PWG2/RESONANCES/macros/train",
  16    const char *taskName = "RSNtask"
  17 )
  18 {
  19    // ==================================================================================================================
  20    // == OPTIONS =======================================================================================================
  21    // ==================================================================================================================
  22
  23    // Instead or getting confused with plenty of arguments in the macro (with default values),
  24    // we use a unique string of options with a set of conventional strings to set up the job:
  25    // -- "MC"/"DATA" --> what kind of sample
  26    // -- "ITS"/"TPC" --> what tracks to use (ITS standalone and/or TPC+ITS)
  27    // -- "xxxPID"    --> add the PID cuts for the detector xxx.
  28    //
  29    // In this point, these options are converted into boolean variables.
  30
  31    TString opt(options);
  32    opt.ToUpper();
  33    opt.ReplaceAll(" ", "");
  34
  35    Bool_t addITS = opt.Contains("ITS");
  36    Bool_t addTPC = opt.Contains("TPC");
  37    Bool_t useITS = opt.Contains("ITSPID");
  38    Bool_t useTPC = opt.Contains("TPCPID");
  39    Bool_t useTOF = opt.Contains("TOFPID");
  40
  41    // correct options when needed
  42    if (!addITS) useITS = kFALSE;
  43    if (!addTPC) useTPC = useTOF = kFALSE;
  44
  45    // ==================================================================================================================
  46    // == DEFINITIONS ===================================================================================================
  47    // ==================================================================================================================
  48
  49    // We put here the definitions of all objects which are needed in the following, in order to have then
  50    // a more readable code in the points where these objects are added to the analysis manager.
  51
  52    // pair definitions --> decay channels:
  53    // in our case, unlike-charged KK pairs for the signal, and like-charged ones for background
  54    AliRsnPairDef *pairDef[3];
  55    pairDef[0] = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '-', 333, 1.019455); // unlike
  56    pairDef[1] = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '+', 333, 1.019455); // like ++
  57    pairDef[2] = new AliRsnPairDef(AliRsnDaughter::kKaon, '-', AliRsnDaughter::kKaon, '-', 333, 1.019455); // like --
  58
  59    // computation objects:
  60    // these are the objects which drive the computations, whatever it is (histos or tree filling)
  61    // and all tracks passed to them will be given the mass according to the reference pair definition
  62    // we create two unlike-sign pair computators, one for all pairs and another for true pairs (useful in MC)
  63    AliRsnPairFunctions *pair[4];
  64    pair[0] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_phi", opt.Data()), pairDef[0]); // unlike - true
  65    pair[1] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_all", opt.Data()), pairDef[0]); // unlike - all
  66    pair[2] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_all", opt.Data()), pairDef[1]); // like ++
  67    pair[3] = new AliRsnPairFunctions(Form("%s_kaonP_kaonM_all", opt.Data()), pairDef[2]); // like --
  68
  69    // set additional option for true pairs (slot [0])
  70    pair[0]->SetOnlyTrue(kTRUE);
  71    pair[0]->SetCheckDecay(kTRUE);
  72
  73    // ==================================================================================================================
  74    // == CUTS AND AXES =================================================================================================
  75    // ==================================================================================================================
  76
  77    //
  78    // Track quality for ITS standalone:
  79    // this cut is used to select tracks of good quality, irrespective of the PID.
  80    // When adding status flags, the second argument tells if each considered flag
  81    // must be active or not in the track status, since the ITS-SA tracks need that
  82    // some of them are OFF (e.g.: kTPCin)
  83    //
  84    AliRsnCutTrackQuality *cutQualityITS = new AliRsnCutTrackQuality("cutQualityITS");
  85    cutQualityITS->AddStatusFlag(AliESDtrack::kITSin    , kTRUE);
  86    cutQualityITS->AddStatusFlag(AliESDtrack::kTPCin    , kFALSE);
  87    cutQualityITS->AddStatusFlag(AliESDtrack::kITSrefit , kTRUE);
  88    cutQualityITS->AddStatusFlag(AliESDtrack::kTPCrefit , kFALSE);
  89    cutQualityITS->AddStatusFlag(AliESDtrack::kITSpureSA, kFALSE);
  90    cutQualityITS->AddStatusFlag(AliESDtrack::kITSpid   , kTRUE);
  91    cutQualityITS->SetPtRange(0.15, 1E+20);
  92    cutQualityITS->SetEtaRange(-0.8, 0.8);
  93    cutQualityITS->SetDCARPtFormula("0.0595+0.0182/pt^1.55");
  94    cutQualityITS->SetDCAZmax(2.0);
  95    cutQualityITS->SetSPDminNClusters(1);
  96    cutQualityITS->SetITSminNClusters(4);
  97    cutQualityITS->SetITSmaxChi2(2.0);
  98    cutQualityITS->SetTPCminNClusters(0);
  99    cutQualityITS->SetTPCmaxChi2(1E+10);
 100    cutQualityITS->SetRejectKinkDaughters();
 101
 102    //
 103    // Track quality for TPC+ITS:
 104    // works exactly like the one above, but has settings for selecting TPC+ITS tracks
 105    // in this case, the flags required are all necessary, so here the procedure is simpler
 106    //
 107    AliRsnCutTrackQuality *cutQualityTPC = new AliRsnCutTrackQuality("cutQualityTPC");
 108    cutQualityTPC->AddStatusFlag(AliESDtrack::kTPCin   , kTRUE);
 109    cutQualityTPC->AddStatusFlag(AliESDtrack::kTPCrefit, kTRUE);
 110    cutQualityTPC->AddStatusFlag(AliESDtrack::kITSrefit, kTRUE);
 111    cutQualityTPC->SetPtRange(0.15, 1E+20);
 112    cutQualityTPC->SetEtaRange(-0.8, 0.8);
 113    cutQualityTPC->SetDCARPtFormula("0.0182+0.0350/pt^1.01");
 114    cutQualityTPC->SetDCAZmax(2.0);
 115    cutQualityTPC->SetSPDminNClusters(1);
 116    cutQualityTPC->SetITSminNClusters(0);
 117    cutQualityTPC->SetITSmaxChi2(1E+20);
 118    cutQualityTPC->SetTPCminNClusters(70);
 119    cutQualityTPC->SetTPCmaxChi2(4.0);
 120    cutQualityTPC->SetRejectKinkDaughters();
 121
 122    //
 123    // ITS PID
 124    // In this implementation, it is a 3sigma cut around the Bethe-Bloch value.
 125    //
 126    // NOTE:
 127    // --> The initialization of the BB is different between data and MC.
 128    // --> The cut is the same for all momenta.
 129    //
 130    AliRsnCutPIDITS *cutPIDITSkaon = new AliRsnCutPIDITS("cutPIDITSkaon", AliPID::kKaon, -3.0, 3.0);
 131    cutPIDITSkaon->SetMC(isMC);
 132    cutPIDITSkaon->SetMomentumRange(0.0, 1E20);
 133
 134    //
 135    // TPC PID
 136    // In this implementation, there are two instances:
 137    // - below 350 MeV --> 5 sigma cut
 138    // - above 350 MeV --> 3 sigma cut
 139    //
 140    // NOTE:
 141    // --> The initialization of the BB is different between data and MC.
 142    //
 143    AliRsnCutPIDTPC *cutPIDTPCkaonLow  = new AliRsnCutPIDTPC("cutPIDTPCkaonLow" , AliPID::kKaon, -5.0, 5.0);
 144    AliRsnCutPIDTPC *cutPIDTPCkaonHigh = new AliRsnCutPIDTPC("cutPIDTPCkaonHigh", AliPID::kKaon, -3.0, 3.0);
 145
 146    // assign the momentum range and tell to reject tracks outside it
 147    cutPIDTPCkaonLow ->SetMomentumRange(0.00, 0.35);
 148    cutPIDTPCkaonHigh->SetMomentumRange(0.35, 1E20);
 149    cutPIDTPCkaonLow ->SetRejectOutside(kTRUE);
 150    cutPIDTPCkaonHigh->SetRejectOutside(kTRUE);
 151
 152    // BB parameterization depends on data sample (MC, data)
 153    // the momentum range is passed and tracks outside it are rejected
 154    Double_t bbPar[5];
 155    if (isMC) {
 156       bbPar[0] = 2.15898 / 50.0;
 157       bbPar[1] = 1.75295E1;
 158       bbPar[2] = 3.40030E-9;
 159       bbPar[3] = 1.96178;
 160       bbPar[4] = 3.91720;
 161    } else {
 162       bbPar[0] = 1.41543 / 50.0;
 163       bbPar[1] = 2.63394E1;
 164       bbPar[2] = 5.0411E-11;
 165       bbPar[3] = 2.12543;
 166       bbPar[4] = 4.88663;
 167    }
 168    cutPIDTPCkaonLow ->SetBBParam(bbPar);
 169    cutPIDTPCkaonHigh->SetBBParam(bbPar);
 170
 171    //
 172    // TOF PID
 173    // In this implementation it is a 3sigma cout aroung expected kaon time.
 174    //
 175    // NOTE:
 176    // --> It is important to choose if this cut must reject tracks not matched in TOF.
 177    //     Usually, if TPC pid is used, we can accept them, otherwise we must reject.
 178    //     (here we assume TPC is used)
 179    //
 180    AliRsnCutPIDTOF *cutPIDTOFkaon = new AliRsnCutPIDTOF("cutPIDTOFkaon", AliPID::kKaon, -3.0, 3.0);
 181    cutPIDTOFkaon->SetRejectUnmatched(!useTPC);
 182
 183    //
 184    // Rapidity cut
 185    // Only thing to consider is that it needs a support object to define mass
 186    //
 187    AliRsnCutValue *cutRapidity = new AliRsnCutValue("cutY", AliRsnValue::kPairY, -0.5, 0.5);
 188    cutRapidity->GetValueObj()->SetSupportObject(pairDef[0]);
 189
 190    //
 191    // Axes
 192    //
 193    // NOTE:
 194    // --> multiplicity has variable bins, defined by array below
 195    //
 196
 197    Double_t mult[] = { 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13.,  14.,  15.,  16.,  17.,  18.,  19.,
 198                       20., 21., 22., 23., 24., 25., 30., 35., 40., 50., 60., 70., 80., 90., 100., 120., 140., 160., 180., 200., 500.};
 199    Int_t    nmult  = sizeof(mult) / sizeof(mult[0]);
 200
 201    AliRsnValue *axisIM      = new AliRsnValue("IM"  , AliRsnValue::kPairInvMass   ,  0.9, 1.4, 0.001);
 202    AliRsnValue *axisRes     = new AliRsnValue("Res" , AliRsnValue::kPairInvMassRes, -0.5, 0.5, 0.001);
 203    AliRsnValue *axisPt      = new AliRsnValue("PT"  , AliRsnValue::kPairPt        ,  0.0, 5.0, 0.1  );
 204    AliRsnValue *axisY       = new AliRsnValue("Y"   , AliRsnValue::kPairY         , -1.1, 1.1, 0.1  );
 205    AliRsnValue *axisMultESD = new AliRsnValue("MESD", AliRsnValue::kEventMultESDCuts, nmult, mult);
 206    AliRsnValue *axisMultSPD = new AliRsnValue("MSPD", AliRsnValue::kEventMultSPD    , nmult, mult);
 207    AliRsnValue *axisMultMC  = new AliRsnValue("MMC" , AliRsnValue::kEventMultMC     , nmult, mult);
 208
 209    // ==================================================================================================================
 210    // == PRELIMINARY OPERATIONS ========================================================================================
 211    // ==================================================================================================================
 212
 213    // retrieve task from manager, using its name
 214    AliAnalysisManager *mgr  = AliAnalysisManager::GetAnalysisManager();
 215    AliRsnAnalysisTask *task = (AliRsnAnalysisTask*)mgr->GetTask(taskName);
 216    if (!task) {
 217       Error("RsnConfigMonitor", "Task not found");
 218       return kFALSE;
 219    }
 220
 221    // ==================================================================================================================
 222    // == SINGLE DAUGHTER CUTS ==========================================================================================
 223    // ==================================================================================================================
 224
 225    // we need to define cuts for track quality and for PID
 226    // which one will be used, depend on the analysis type
 227
 228    Char_t qualityITS[255], qualityTPC[255];
 229    Char_t pidITS[255], pidTPC[255], pidTOF[255];
 230    Char_t schemeITS[255], schemeTPC[255], scheme[255];
 231
 232    sprintf(qualityITS, "%s"     , cutQualityITS->GetName());
 233    sprintf(qualityTPC, "%s"     , cutQualityTPC->GetName());
 234    sprintf(pidITS    , "%s"     , cutPIDITSkaon->GetName());
 235    sprintf(pidTPC    , "(%s|%s)", cutPIDTPCkaonHigh->GetName(), cutPIDTPCkaonLow->GetName());
 236    sprintf(pidTOF    , "%s"     , cutPIDTOFkaon->GetName());
 237    sprintf(schemeITS , "");
 238    sprintf(schemeTPC , "");
 239
 240    // choose cuts to add depending on used tracks
 241    for (Int_t i = 0; i < 4; i++) {
 242       // assign cuts to common manager for daughters
 243       AliRsnCutSet *cutSet = pair[i]->GetCutManager()->GetCommonDaughterCuts();
 244       // add cuts and define schemes depending on options
 245       if (addITS) {
 246          // if adding ITS standalone,
 247          // by default we use the quality check
 248          // and if required we add the PID (option "ITSPID")
 249          cutSet->AddCut(cutQualityITS);
 250          sprintf(schemeITS, "%s", qualityITS);
 251          if (useITS) {
 252             sprintf(schemeITS, "%s & %s", qualityITS, pidITS);
 253             cutSet->AddCut(cutPIDITSkaon);
 254          }
 255       }
 256       if (addTPC) {
 257          // if adding TPC+ITS tracks,
 258          // by default we use the quality check
 259          // and then we can use no PID, or both,
 260          // or only one among TPC and TOF
 261          // the TPC PID cut is always an 'OR' of
 262          // the two ones, in order to cover the full momentum range
 263          cutSet->AddCut(cutQualityTPC);
 264          sprintf(schemeTPC, "%s", qualityTPC);
 265          if (useTPC && useTOF) {
 266             cutSet->AddCut(cutPIDTPCkaonLow);
 267             cutSet->AddCut(cutPIDTPCkaonHigh);
 268             cutSet->AddCut(cutPIDTOFkaon);
 269             sprintf(schemeTPC, "%s & %s & %s", qualityTPC, pidTPC, pidTOF);
 270          } else if (useTPC) {
 271             cutSet->AddCut(cutPIDTPCkaonLow);
 272             cutSet->AddCut(cutPIDTPCkaonHigh);
 273             sprintf(schemeTPC, "%s & %s", qualityTPC, pidTPC);
 274          } else if (useTOF) {
 275             cutSet->AddCut(cutPIDTOFkaon);
 276             sprintf(schemeTPC, "%s & %s", qualityTPC, pidTOF);
 277          }
 278       }
 279
 280       // final scheme depends on what of the above were added
 281       // in case both ITS-SA and TPC tracks are added, the scheme
 282       // is the OR of the cuts for the first and those for the second
 283       // category of tracks
 284       if (strlen(schemeITS) > 0 && strlen(schemeTPC) > 0) {
 285          sprintf(scheme, "(%s) | (%s)", schemeITS, schemeTPC);
 286       } else if (strlen(schemeITS) > 0) {
 287          sprintf(scheme, "%s", schemeITS);
 288       } else if (strlen(schemeTPC) > 0) {
 289          sprintf(scheme, "%s", schemeTPC);
 290       } else {
 291          ::Error("Scheme is empty!");
 292          return kFALSE;
 293       }
 294       cutSet->SetCutScheme(scheme);
 295       ::Info("RsnConfigPhi", "Scheme for daughter cuts: %s", cutSet->GetCutScheme().Data());
 296    }
 297
 298    // ==================================================================================================================
 299    // == PAIR CUTS =====================================================================================================
 300    // ==================================================================================================================
 301
 302    // in this case, we add the cut to the specific cut sets of all pairs
 303    // and we must then loop over all pairs, to add cuts to the related sets
 304    for (Int_t ipair = 0; ipair < 4; ipair++) {
 305       pair[ipair]->GetCutManager()->GetMotherCuts()->AddCut(cutRapidity);
 306       pair[ipair]->GetCutManager()->GetMotherCuts()->SetCutScheme(cutRapidity->GetName());
 307       ::Info("RsnConfigPhi", "Scheme for pair cuts: %s", pair[ipair]->GetCutManager()->GetMotherCuts()->GetCutScheme().Data());
 308    }
 309
 310    // ==================================================================================================================
 311    // == FUNCTIONS FOR HISTOGRAMS ======================================================================================
 312    // ==================================================================================================================
 313
 314    // we define a histogram with inv mass and other support binning values
 315    // and, when possible, a resolution on inv. mass
 316
 317    // create function for inv. mass and add axes
 318    AliRsnFunction *fcnIM = new AliRsnFunction;
 319    if (!fcnIM->AddAxis(axisIM)     ) return kFALSE;
 320    if (!fcnIM->AddAxis(axisPt)     ) return kFALSE;
 321    if (!fcnIM->AddAxis(axisMultSPD)) return kFALSE;
 322    fcnIM->UseSparse();
 323
 324    // create function for inv. mass and add axes
 325    AliRsnFunction *fcnIMRes = new AliRsnFunction;
 326    if (!fcnIMRes->AddAxis(axisIM)     ) return kFALSE;
 327    if (!fcnIMRes->AddAxis(axisPt)     ) return kFALSE;
 328    if (!fcnIMRes->AddAxis(axisMultSPD)) return kFALSE;
 329    if (!fcnIMRes->AddAxis(axisRes)    ) return kFALSE;
 330
 331    // add functions to pairs
 332    pair[0]->AddFunction(fcnIMRes);
 333    for (Int_t ipair = 1; ipair < 4; ipair++) pair[ipair]->AddFunction(fcnIM);
 334
 335    // ==================================================================================================================
 336    // == CONCLUSION ====================================================================================================
 337    // ==================================================================================================================
 338
 339    // add all created AliRsnPair objects to the AliRsnAnalysisManager in the task
 340    task->GetAnalysisManager()->Add(pair[1]);
 341    task->GetAnalysisManager()->Add(pair[2]);
 342    task->GetAnalysisManager()->Add(pair[3]);
 343    if (isMC) task->GetAnalysisManager()->Add(pair[0]);
 344
 345    return kTRUE;
 346 }