PWG2/RESONANCES/macros/AddAnalysisTaskRsn.C

   1 //
   2 // Macro to create the full analysis manager for Resonances
   3 //
   4
   5 static Double_t  cov11 = 2;
   6 static Double_t  cov22 = 2;
   7 static Double_t  cov33 = 0.5;
   8 static Double_t  cov44 = 0.5;
   9 static Double_t  cov55 = 2;
  10 static Double_t  nSigmaToVertex = 4;
  11 static Double_t  dcaToVertex = 3.0;
  12 static Double_t  maxChi2PerClusterTPC = 3.5;
  13 static Bool_t    requireTPCRefit = kTRUE;
  14 static Bool_t    requireSigmaToVertex = kTRUE;
  15 static Bool_t    acceptKinkDaughters = kFALSE;
  16 static Int_t     minNClustersTPC = 50;
  17
  18 Bool_t AddAnalysisTaskRsn
  19 (
  20   AliLog::EType_t  debugType  = AliLog::kInfo, // debug depth for some classes
  21   Bool_t           useTPCOnly = kFALSE,        // use TPC only PID
  22   const char      *outFile    = "rsn.root",    // output file name
  23   Bool_t           sourceESD  = kTRUE          // if true, the source of data is ESD, otherwise is AOD from filter task
  24 )
  25 {
  26   // retrieve analysis manager
  27   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  28
  29   // initialize task
  30   AliRsnAnalysisSE *task = new AliRsnAnalysisSE("AliRsnAnalysisSE");
  31   //task->SetLogType(debugType, "AliRsnAnalysisManager:AliRsnPairManager:AliRsnPairManager:AliRsnPair");
  32   task->SetLogType(debugType, "AliRsnCut:AliRsnCutPrimaryVertex");
  33
  34   // set prior probabilities for PID
  35   task->SetPriorProbability(AliPID::kElectron, 0.02);
  36   task->SetPriorProbability(AliPID::kMuon,     0.02);
  37   task->SetPriorProbability(AliPID::kPion,     0.83);
  38   task->SetPriorProbability(AliPID::kKaon,     0.07);
  39   task->SetPriorProbability(AliPID::kProton,   0.06);
  40   task->DumpPriors();
  41
  42   // initialize analysis manager with pairs from config
  43   AliRsnAnalysisManager *anaMgr = task->GetAnalysisManager("MyAnalysisSE");
  44   anaMgr->Add(CreatePairsNoPID("PHI_NoPID_0"   , 333, AliPID::kKaon, AliPID::kKaon, 10000.0));
  45   anaMgr->Add(CreatePairsNoPID("PHI_NoPID_10"  , 333, AliPID::kKaon, AliPID::kKaon,     0.1));
  46   anaMgr->Add(CreatePairsNoPID("PHI_NoPID_20"  , 333, AliPID::kKaon, AliPID::kKaon,     0.2));
  47   anaMgr->Add(CreatePairsPID  ("PHI_PID"       , 333, AliPID::kKaon, AliPID::kKaon));
  48   //anaMgr->Add(CreatePairsNoPID("KSTAR_NoPID_0" , 313, AliPID::kKaon, AliPID::kPion, 10000.0));
  49   //anaMgr->Add(CreatePairsNoPID("KSTAR_NoPID_10", 313, AliPID::kKaon, AliPID::kPion,     0.1));
  50   //anaMgr->Add(CreatePairsNoPID("KSTAR_NoPID_20", 313, AliPID::kKaon, AliPID::kPion,     0.2));
  51   //anaMgr->Add(CreatePairsPID  ("KSTAR_PID"     , 313, AliPID::kKaon, AliPID::kPion));
  52   cout << "CREATED PAIRS" << endl;
  53
  54   // setup cuts for ESD tracks
  55   if (sourceESD)
  56   {
  57     AliESDtrackCuts *esdCuts = new AliESDtrackCuts;
  58     esdCuts->SetMaxCovDiagonalElements(cov11, cov22, cov33, cov44, cov55);
  59     esdCuts->SetRequireSigmaToVertex(requireSigmaToVertex);
  60     if (requireSigmaToVertex) esdCuts->SetMaxNsigmaToVertex(nSigmaToVertex);
  61     else
  62     {
  63       esdCuts->SetDCAToVertexZ(dcaToVertex);
  64       esdCuts->SetDCAToVertexXY(dcaToVertex);
  65     }
  66     esdCuts->SetRequireTPCRefit(requireTPCRefit);
  67     esdCuts->SetAcceptKinkDaughters(acceptKinkDaughters);
  68     esdCuts->SetMinNClustersTPC(minNClustersTPC);
  69     esdCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
  70     task->SetESDtrackCuts(esdCuts);
  71   }
  72   cout << "SET ESD CUTS" << endl;
  73
  74   // set PID customization if necessary
  75   if (useTPCOnly) {
  76     Info("Using TPC only PID");
  77     task->GetPIDDef()->SetScheme(AliRsnPIDDefESD::kSchemeTPC);
  78   }
  79   cout << "SET PID SCHEME" << endl;
  80
  81   // setup cuts for events (good primary vertex)
  82   AliRsnCutPrimaryVertex *cutVertex = new AliRsnCutPrimaryVertex("cutVertex", 3);
  83   AliRsnCutSet *cutSetEvent = new AliRsnCutSet("eventCuts");
  84   cutSetEvent->AddCut(cutVertex);
  85   cutSetEvent->SetCutScheme("cutVertex");
  86   task->SetEventCuts(cutSetEvent);
  87   cout << "SET EVENT CUT SCHEME" << endl;
  88
  89   // add the task to manager
  90   mgr->AddTask(task);
  91   cout << "ADD TASK" << endl;
  92
  93   // connect input container according to source choice
  94   if (sourceESD) mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
  95   else mgr->ConnectInput(task, 0, mgr->GetCommonOutputContainer());
  96   cout << "CONNECT INPUT" << endl;
  97
  98   // initialize and connect container for the output
  99   AliAnalysisDataContainer *out = mgr->CreateContainer("RSN", TList::Class(), AliAnalysisManager::kOutputContainer, outFile);
 100   mgr->ConnectOutput(task, 1, out);
 101   cout << "CONNECT OUTPUT" << endl;
 102 }
 103
 104 AliRsnFunction* DefineFunctionIM()
 105 {
 106 //
 107 // In general, for all processed pairs in one analysis the same functions are computed.
 108 // Then, they are defined separately here and added in the same way to all pairs.
 109 //
 110
 111   // define all binnings
 112   AliRsnFunctionAxis *axisIM   = new AliRsnFunctionAxis(AliRsnFunctionAxis::kPairInvMass,    1000,  0.0,   2.0);
 113   AliRsnFunctionAxis *axisPt   = new AliRsnFunctionAxis(AliRsnFunctionAxis::kPairPt,           10,  0.0,  10.0);
 114   AliRsnFunctionAxis *axisEta  = new AliRsnFunctionAxis(AliRsnFunctionAxis::kPairEta,          10, -1.0,   1.0);
 115   AliRsnFunctionAxis *axisMult = new AliRsnFunctionAxis(AliRsnFunctionAxis::kEventMult,         8,  0.0, 200.0);
 116
 117   // define function
 118   AliRsnFunction *fcn = new AliRsnFunction;
 119   fcn->AddAxis(axisIM);
 120   fcn->AddAxis(axisPt);
 121   fcn->AddAxis(axisEta);
 122   fcn->AddAxis(axisMult);
 123
 124   return fcn;
 125 }
 126
 127 AliRsnFunction* DefineFunctionP1P2()
 128 {
 129 //
 130 // In general, for all processed pairs in one analysis the same functions are computed.
 131 // Then, they are defined separately here and added in the same way to all pairs.
 132 //
 133
 134   // define all binnings
 135   AliRsnFunctionAxis *axisP1   = new AliRsnFunctionAxis(AliRsnFunctionAxis::kTrack1P,          50,  0.0,   5.0);
 136   AliRsnFunctionAxis *axisP2   = new AliRsnFunctionAxis(AliRsnFunctionAxis::kTrack2P,          50,  0.0,   5.0);
 137   AliRsnFunctionAxis *axisPt   = new AliRsnFunctionAxis(AliRsnFunctionAxis::kPairPt,           10,  0.0,  10.0);
 138   AliRsnFunctionAxis *axisEta  = new AliRsnFunctionAxis(AliRsnFunctionAxis::kPairEta,          10, -1.0,   1.0);
 139   AliRsnFunctionAxis *axisMult = new AliRsnFunctionAxis(AliRsnFunctionAxis::kEventMult,         8,  0.0, 200.0);
 140
 141   // define function
 142   AliRsnFunction *fcn = new AliRsnFunction;
 143   fcn->AddAxis(axisP1);
 144   fcn->AddAxis(axisP2);
 145   fcn->AddAxis(axisPt);
 146   fcn->AddAxis(axisEta);
 147   fcn->AddAxis(axisMult);
 148
 149   return fcn;
 150 }
 151
 152 AliRsnPairManager* CreatePairsNoPID
 153 (
 154   const char            *pairMgrName,    // name for the pair manager
 155   Int_t                  resonancePDG,   // PDG code of resonance (for true pairs)
 156   AliPID::EParticleType  type1,          // PID of one member of decay (+)
 157   AliPID::EParticleType  type2,          // PID of other member of decay (-)
 158   Double_t               bbCut = 10000.0 // cut on Bethe-Bloch
 159 )
 160 {
 161 //
 162 // Creates an AliRsnPairMgr for a specified resonance, which contains:
 163 // - signal (inv. mass)
 164 // - event mixing (inv. mass)
 165 // - like-sign (inv. mass)
 166 // - true pairs (inv. mass, resolution)
 167 //
 168 // For all pairs, a binning in Pt and Eta is provided, and a cut in multiplicity
 169 // which defines a multiplicity bin where the analysis is computed.
 170 //
 171 // Arguments define how the pair manager must be built, and are explained above
 172 //
 173
 174   AliRsnPairManager  *pairMgr  = new AliRsnPairManager(pairMgrName);
 175
 176   // === PAIR DEFINITIONS =========================================================================
 177
 178   // if particle #1 and #2 are different, two histograms must be built
 179   // for each scheme (signal, true, mixed, like-sign) exchanging both particles and signs
 180   Int_t i, j, nArray = 1;
 181   if (type1 != type2) nArray = 2;
 182
 183   AliRsnPairDef *defUnlike[2] = {0, 0};
 184   AliRsnPairDef *defLikePP[2] = {0, 0};
 185   AliRsnPairDef *defLikeMM[2] = {0, 0};
 186
 187   defUnlike[0] = new AliRsnPairDef(type1, '+', type2, '-', resonancePDG);
 188   defLikePP[0] = new AliRsnPairDef(type1, '+', type2, '+', resonancePDG);
 189   defLikeMM[0] = new AliRsnPairDef(type1, '-', type2, '-', resonancePDG);
 190
 191   defUnlike[1] = new AliRsnPairDef(type2, '+', type1, '-', resonancePDG);
 192   defLikePP[1] = new AliRsnPairDef(type2, '+', type1, '+', resonancePDG);
 193   defLikeMM[1] = new AliRsnPairDef(type2, '-', type1, '-', resonancePDG);
 194
 195   // === PAIR ANALYSIS ENGINES ====================================================================
 196
 197   // define null (dummy) objects and initialize only the ones which are needed,
 198   // depending again on particle types;
 199   // array is organized as follows:
 200   // [0] - true pairs
 201   // [1] - signal
 202   // [2] - like PP
 203   // [3] - like MM
 204   AliRsnPair *noPID[2][4] = {0,0,0,0,0,0,0,0};
 205
 206   for (i = 0; i < nArray; i++) {
 207     noPID[i][0] = new AliRsnPair(AliRsnPair::kNoPID, defUnlike[i]);
 208     noPID[i][1] = new AliRsnPair(AliRsnPair::kNoPID, defUnlike[i]);
 209     noPID[i][2] = new AliRsnPair(AliRsnPair::kNoPID, defLikePP[i]);
 210     noPID[i][3] = new AliRsnPair(AliRsnPair::kNoPID, defLikeMM[i]);
 211   }
 212
 213   // === CUTS =====================================================================================
 214
 215   // cuts for tracks:
 216   // -- Bethe-Bloch & required kaon PID
 217   AliRsnCutBetheBloch *cutKaonBB = new AliRsnCutBetheBloch("cutKaon", bbCut, AliPID::kKaon);
 218   cutKaonBB->SetCalibConstant(0, 0.76176e-1);
 219   cutKaonBB->SetCalibConstant(1, 10.632);
 220   cutKaonBB->SetCalibConstant(2, 0.13279e-4);
 221   cutKaonBB->SetCalibConstant(3, 1.8631);
 222   cutKaonBB->SetCalibConstant(4, 1.9479);
 223
 224   // cuts on pairs:
 225   // -- true daughters of a phi resonance (only for true pairs histogram)cutSetPairTrue->AddCut(cutTrue);
 226   AliRsnCutStd *cutTruePair = new AliRsnCutStd("cutTrue", AliRsnCutStd::kTruePair, resonancePDG);
 227
 228   // cuts on event:
 229   // -- none (specified for whole task)
 230
 231   // cut set definition for all pairs
 232   AliRsnCutSet *cutSetParticle = new AliRsnCutSet("trackCuts");
 233   cutSetParticle->AddCut(cutKaonBB);
 234   cutSetParticle->SetCutScheme("cutKaonBB");
 235
 236   // cut set definition for true pairs
 237   AliRsnCutSet *cutSetPairTrue = new AliRsnCutSet("truePairs");
 238   cutSetPairTrue->AddCut(cutTruePair);
 239   cutSetPairTrue->SetCutScheme("cutTrue");
 240
 241   // cut manager for all pairs
 242   AliRsnCutMgr *cutMgrAll = new AliRsnCutMgr("std", "All");
 243   cutMgrAll->SetCutSet(AliRsnCut::kParticle, cutSetParticle);
 244
 245   // cut manager for all pairs
 246   AliRsnCutMgr *cutMgrTrue = new AliRsnCutMgr("true", "True");
 247   cutMgrTrue->SetCutSet(AliRsnCut::kParticle, cutSetParticle);
 248   cutMgrTrue->SetCutSet(AliRsnCut::kPair, cutSetPairTrue);
 249
 250   for (i = 0; i < nArray; i++) {
 251     noPID[i][0]->SetCutMgr(cutMgrTrue);
 252     for (j = 1; j < 4; j++) {
 253       noPID[i][j]->SetCutMgr(cutMgrAll);
 254     }
 255   }
 256
 257   // === FUNCTIONS ================================================================================
 258
 259   AliRsnFunction *fcn   = DefineFunctionIM();
 260   AliRsnFunction *fcnPP = DefineFunctionP1P2();
 261
 262   for (i = 0; i < nArray; i++) {
 263     for (j = 0; j < 4; j++) {
 264       noPID[i][j]->AddFunction(fcn);
 265       if (j < 2) noPID[i][j]->AddFunction(fcnPP);
 266     }
 267   }
 268
 269   // === ADD TO PAIR MANAGER ======================================================================
 270
 271   for (i = 0; i < nArray; i++) {
 272     for (j = 0; j < 4; j++) {
 273       pairMgr->AddPair(noPID[i][j]);
 274     }
 275   }
 276
 277   return pairMgr;
 278 }
 279
 280 AliRsnPairManager* CreatePairsPID
 281 (
 282   const char            *pairMgrName,    // name for the pair manager
 283   Int_t                  resonancePDG,   // PDG code of resonance (for true pairs)
 284   AliPID::EParticleType  type1,          // PID of one member of decay (+)
 285   AliPID::EParticleType  type2           // PID of other member of decay (-)
 286 )
 287 {
 288 //
 289 // Creates an AliRsnPairMgr for a specified resonance, which contains:
 290 // - signal (inv. mass)
 291 // - event mixing (inv. mass)
 292 // - like-sign (inv. mass)
 293 // - true pairs (inv. mass, resolution)
 294 //
 295 // For all pairs, a binning in Pt and Eta is provided, and a cut in multiplicity
 296 // which defines a multiplicity bin where the analysis is computed.
 297 //
 298 // Arguments define how the pair manager must be built, and are explained above
 299 //
 300
 301   AliRsnPairManager  *pairMgr  = new AliRsnPairManager(pairMgrName);
 302
 303   // === PAIR DEFINITIONS =========================================================================
 304
 305   // if particle #1 and #2 are different, two histograms must be built
 306   // for each scheme (signal, true, mixed, like-sign) exchanging both particles and signs
 307   Int_t i, j, nArray = 1;
 308   if (type1 != type2) nArray = 2;
 309
 310   AliRsnPairDef *defUnlike[2] = {0, 0};
 311   AliRsnPairDef *defLikePP[2] = {0, 0};
 312   AliRsnPairDef *defLikeMM[2] = {0, 0};
 313
 314   defUnlike[0] = new AliRsnPairDef(type1, '+', type2, '-', resonancePDG);
 315   defLikePP[0] = new AliRsnPairDef(type1, '+', type2, '+', resonancePDG);
 316   defLikeMM[0] = new AliRsnPairDef(type1, '-', type2, '-', resonancePDG);
 317
 318   defUnlike[1] = new AliRsnPairDef(type2, '+', type1, '-', resonancePDG);
 319   defLikePP[1] = new AliRsnPairDef(type2, '+', type1, '+', resonancePDG);
 320   defLikeMM[1] = new AliRsnPairDef(type2, '-', type1, '-', resonancePDG);
 321
 322   // === PAIR ANALYSIS ENGINES ====================================================================
 323
 324   // define null (dummy) objects and initialize only the ones which are needed,
 325   // depending again on particle types;
 326   // array is organized as follows:
 327   // [0] - true pairs
 328   // [1] - signal
 329   // [2] - like PP
 330   // [3] - like MM
 331   AliRsnPair *perfectPID[2][4]   = {0,0,0,0,0,0,0,0};
 332   AliRsnPair *realisticPID[2][4] = {0,0,0,0,0,0,0,0};
 333
 334   for (i = 0; i < nArray; i++) {
 335     perfectPID[i][0] = new AliRsnPair(AliRsnPair::kPerfectPID, defUnlike[i]);
 336     perfectPID[i][1] = new AliRsnPair(AliRsnPair::kPerfectPID, defUnlike[i]);
 337     perfectPID[i][2] = new AliRsnPair(AliRsnPair::kPerfectPID, defLikePP[i]);
 338     perfectPID[i][3] = new AliRsnPair(AliRsnPair::kPerfectPID, defLikeMM[i]);
 339
 340     realisticPID[i][0] = new AliRsnPair(AliRsnPair::kRealisticPID, defUnlike[i]);
 341     realisticPID[i][1] = new AliRsnPair(AliRsnPair::kRealisticPID, defUnlike[i]);
 342     realisticPID[i][2] = new AliRsnPair(AliRsnPair::kRealisticPID, defLikePP[i]);
 343     realisticPID[i][3] = new AliRsnPair(AliRsnPair::kRealisticPID, defLikeMM[i]);
 344   }
 345
 346   // === CUTS =====================================================================================
 347
 348   // cuts for tracks:
 349   // -- nothing
 350
 351   // cuts on pairs:
 352   // -- true daughters of a phi resonance (only for true pairs histogram)cutSetPairTrue->AddCut(cutTrue);
 353   AliRsnCutStd *cutTruePair = new AliRsnCutStd("cutTrue", AliRsnCutStd::kTruePair, resonancePDG);
 354
 355   // cut set definition for true pairs
 356   AliRsnCutSet *cutSetPairTrue = new AliRsnCutSet("truePairs");
 357   cutSetPairTrue->AddCut(cutTruePair);
 358   cutSetPairTrue->SetCutScheme("cutTrue");
 359
 360   // cut manager for all pairs
 361   AliRsnCutMgr *cutMgrAll = new AliRsnCutMgr("std", "All");
 362
 363   // cut manager for all pairs
 364   AliRsnCutMgr *cutMgrTrue = new AliRsnCutMgr("true", "True");
 365   cutMgrTrue->SetCutSet(AliRsnCut::kPair, cutSetPairTrue);
 366
 367   for (i = 0; i < nArray; i++) {
 368     perfectPID[i][0]->SetCutMgr(cutMgrTrue);
 369     realisticPID[i][0]->SetCutMgr(cutMgrTrue);
 370     for (j = 1; j < 4; j++) {
 371       perfectPID[i][j]->SetCutMgr(cutMgrAll);
 372       realisticPID[i][j]->SetCutMgr(cutMgrAll);
 373     }
 374   }
 375
 376   // === FUNCTIONS ================================================================================
 377
 378   AliRsnFunction *fcn   = DefineFunctionIM();
 379   AliRsnFunction *fcnPP = DefineFunctionP1P2();
 380
 381   for (i = 0; i < nArray; i++) {
 382     for (j = 0; j < 4; j++) {
 383       perfectPID[i][j]->AddFunction(fcn);
 384       realisticPID[i][j]->AddFunction(fcn);
 385       if (j < 2) {
 386         perfectPID[i][j]->AddFunction(fcnPP);
 387         realisticPID[i][j]->AddFunction(fcnPP);
 388       }
 389     }
 390   }
 391
 392   // === ADD TO PAIR MANAGER ======================================================================
 393
 394   for (i = 0; i < nArray; i++) {
 395     for (j = 0; j < 4; j++) {
 396       pairMgr->AddPair(perfectPID[i][j]);
 397       pairMgr->AddPair(realisticPID[i][j]);
 398     }
 399   }
 400
 401   return pairMgr;
 402 }
 403