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(const char *outFile = "rsn.root", Bool_t useKine = kTRUE)
  19 {
  20   // Retrieve analysis manager.
  21   // Since it is usually created in the steering macro,
  22   // then we don't ever need to initialize a new one
  23   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  24
  25   // create handlers for input and MonteCarlo
  26   AliMCEventHandler  *mcHandler  = new AliMCEventHandler();
  27   AliESDInputHandler *esdHandler = new AliESDInputHandler();
  28   //AliAODHandler      *aodHandler = new AliAODHandler();
  29   //aodHandler->SetOutputFileName("AliAOD.root");
  30   if (useKine) mgr->SetMCtruthEventHandler(mcHandler);
  31   mgr->SetInputEventHandler(esdHandler);
  32   //mgr->SetOutputEventHandler(aodHandler);
  33
  34   // initialize task
  35   AliRsnAnalysisSE *task = new AliRsnAnalysisSE("AliRsnAnalysisSE");
  36
  37   // set prior probabilities for PID
  38   task->SetPriorProbability(AliPID::kElectron, 0.02);
  39   task->SetPriorProbability(AliPID::kMuon,     0.02);
  40   task->SetPriorProbability(AliPID::kPion,     0.83);
  41   task->SetPriorProbability(AliPID::kKaon,     0.07);
  42   task->SetPriorProbability(AliPID::kProton,   0.06);
  43   task->DumpPriors();
  44
  45   // initialize analysis manager with pairs from config
  46   AliRsnAnalysisManager *anaMgr = task->GetAnalysisManager("MyAnalysisSE");
  47   anaMgr->Add(CreatePairsNoPID("PHI_NoPID_0"   , 333, AliPID::kKaon, AliPID::kKaon, 10000.0));
  48   anaMgr->Add(CreatePairsNoPID("PHI_NoPID_10"  , 333, AliPID::kKaon, AliPID::kKaon,     0.1));
  49   anaMgr->Add(CreatePairsNoPID("PHI_NoPID_20"  , 333, AliPID::kKaon, AliPID::kKaon,     0.2));
  50   anaMgr->Add(CreatePairsPID  ("PHI_PID"       , 333, AliPID::kKaon, AliPID::kKaon));
  51   anaMgr->Add(CreatePairsNoPID("KSTAR_NoPID_0" , 313, AliPID::kKaon, AliPID::kPion, 10000.0));
  52   anaMgr->Add(CreatePairsNoPID("KSTAR_NoPID_10", 313, AliPID::kKaon, AliPID::kPion,     0.1));
  53   anaMgr->Add(CreatePairsNoPID("KSTAR_NoPID_20", 313, AliPID::kKaon, AliPID::kPion,     0.2));
  54   anaMgr->Add(CreatePairsPID  ("KSTAR_PID"     , 313, AliPID::kKaon, AliPID::kPion));
  55
  56   // setup cuts
  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   // add the task to manager
  73   mgr->AddTask(task);
  74
  75   // initialize container for the output
  76   // (input container is common)
  77   AliAnalysisDataContainer *out = mgr->CreateContainer("RSN", TList::Class(), AliAnalysisManager::kOutputContainer, outFile);
  78
  79   // connect input and output slots for this task
  80   mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
  81   mgr->ConnectOutput(task, 1, out);
  82 }
  83
  84 void AddRsnFunctionsToPair(AliRsnPair *pair)
  85 {
  86 //
  87 // In general, for all processed pairs in one analysis the same functions are computed.
  88 // Then, they are defined separately here and added in the same way to all pairs.
  89 //
  90
  91   // define histogram templates
  92   AliRsnHistoDef *hdIM   = new AliRsnHistoDef(1000,  0.0,   2.0);
  93   AliRsnHistoDef *hdPt   = new AliRsnHistoDef(  10,  0.0,  10.0);
  94   AliRsnHistoDef *hdMult = new AliRsnHistoDef(   8,  0.0, 200.0);
  95
  96   // define functions:
  97   // -- one computed w.r. to Pt and Mult
  98   // -- one computed w.r. to Eta and Mult
  99   AliRsnFunction *fcnIM1 = new AliRsnFunction(AliRsnFunction::kInvMass, AliRsnFunction::kBinPairPt, AliRsnFunction::kBinEventMult, hdIM, hdPt, hdMult);
 100
 101   // add to pair
 102   pair->AddFunction(fcnIM1);
 103 }
 104
 105 AliRsnPairManager* CreatePairsNoPID
 106 (
 107   const char            *pairMgrName,    // name for the pair manager
 108   Int_t                  resonancePDG,   // PDG code of resonance (for true pairs)
 109   AliPID::EParticleType  type1,          // PID of one member of decay (+)
 110   AliPID::EParticleType  type2,          // PID of other member of decay (-)
 111   Double_t               bbCut = 10000.0 // cut on Bethe-Bloch
 112 )
 113 {
 114 //
 115 // Creates an AliRsnPairMgr for a specified resonance, which contains:
 116 // - signal (inv. mass)
 117 // - event mixing (inv. mass)
 118 // - like-sign (inv. mass)
 119 // - true pairs (inv. mass, resolution)
 120 //
 121 // For all pairs, a binning in Pt and Eta is provided, and a cut in multiplicity
 122 // which defines a multiplicity bin where the analysis is computed.
 123 //
 124 // Arguments define how the pair manager must be built, and are explained above
 125 //
 126
 127   AliRsnPairManager  *pairMgr  = new AliRsnPairManager(pairMgrName);
 128
 129   // === PAIR DEFINITIONS =========================================================================
 130
 131   // if particle #1 and #2 are different, two histograms must be built
 132   // for each scheme (signal, true, mixed, like-sign) exchanging both particles and signs
 133   Int_t i, j, nArray = 1;
 134   if (type1 != type2) nArray = 2;
 135
 136   AliRsnPairDef *defUnlike[2] = {0, 0};
 137   AliRsnPairDef *defLikePP[2] = {0, 0};
 138   AliRsnPairDef *defLikeMM[2] = {0, 0};
 139
 140   defUnlike[0] = new AliRsnPairDef(type1, '+', type2, '-', resonancePDG);
 141   defLikePP[0] = new AliRsnPairDef(type1, '+', type2, '+', resonancePDG);
 142   defLikeMM[0] = new AliRsnPairDef(type1, '-', type2, '-', resonancePDG);
 143
 144   defUnlike[1] = new AliRsnPairDef(type2, '+', type1, '-', resonancePDG);
 145   defLikePP[1] = new AliRsnPairDef(type2, '+', type1, '+', resonancePDG);
 146   defLikeMM[1] = new AliRsnPairDef(type2, '-', type1, '-', resonancePDG);
 147
 148   // === PAIR ANALYSIS ENGINES ====================================================================
 149
 150   // define null (dummy) objects and initialize only the ones which are needed,
 151   // depending again on particle types;
 152   // array is organized as follows:
 153   // [0] - true pairs
 154   // [1] - signal
 155   // [2] - like PP
 156   // [3] - like MM
 157   AliRsnPair *noPID[2][4] = {0,0,0,0,0,0,0,0};
 158
 159   for (i = 0; i < nArray; i++) {
 160     noPID[i][0] = new AliRsnPair(AliRsnPair::kNoPID, defUnlike[i]);
 161     noPID[i][1] = new AliRsnPair(AliRsnPair::kNoPID, defUnlike[i]);
 162     noPID[i][2] = new AliRsnPair(AliRsnPair::kNoPID, defLikePP[i]);
 163     noPID[i][3] = new AliRsnPair(AliRsnPair::kNoPID, defLikeMM[i]);
 164   }
 165
 166   // === CUTS =====================================================================================
 167
 168   // cuts for tracks:
 169   // -- Bethe-Bloch & required kaon PID
 170   AliRsnCutStd *cutReqKaon = new AliRsnCutStd("cutReqKaon", AliRsnCutStd::kRequiredPID, (Int_t)AliPID::kKaon);
 171   AliRsnCutBetheBloch *cutKaonBB = new AliRsnCutBetheBloch("cutKaon", bbCut, AliPID::kKaon);
 172   cutKaonBB->SetCalibConstant(0, 0.76176e-1);
 173   cutKaonBB->SetCalibConstant(1, 10.632);
 174   cutKaonBB->SetCalibConstant(2, 0.13279e-4);
 175   cutKaonBB->SetCalibConstant(3, 1.8631);
 176   cutKaonBB->SetCalibConstant(4, 1.9479);
 177
 178   // cuts on pairs:
 179   // -- true daughters of a phi resonance (only for true pairs histogram)cutSetPairTrue->AddCut(cutTrue);
 180   AliRsnCutStd *cutTruePair = new AliRsnCutStd("cutTrue", AliRsnCutStd::kTruePair, resonancePDG);
 181
 182   // cut set definition for all pairs
 183   AliRsnCutSet *cutSetParticle = new AliRsnCutSet("trackCuts");
 184   cutSetParticle->AddCut(cutKaonBB);
 185   cutSetParticle->AddCut(cutReqKaon);
 186   cutSetParticle->SetCutScheme("(cutReqKaon&cutKaonBB)|(!cutReqKaon)");
 187
 188   // cut set definition for true pairs
 189   AliRsnCutSet *cutSetPairTrue = new AliRsnCutSet("truePairs");
 190   cutSetPairTrue->AddCut(cutTruePair);
 191   cutSetPairTrue->SetCutScheme("cutTrue");
 192
 193   // cut manager for all pairs
 194   AliRsnCutMgr *cutMgrAll = new AliRsnCutMgr("std", "All");
 195   cutMgrAll->SetCutSet(AliRsnCut::kParticle, cutSetParticle);
 196
 197   // cut manager for all pairs
 198   AliRsnCutMgr *cutMgrTrue = new AliRsnCutMgr("true", "True");
 199   cutMgrTrue->SetCutSet(AliRsnCut::kParticle, cutSetParticle);
 200   cutMgrTrue->SetCutSet(AliRsnCut::kPair, cutSetPairTrue);
 201
 202   for (i = 0; i < nArray; i++) {
 203     noPID[i][0]->SetCutMgr(cutMgrTrue);
 204     for (j = 1; j < 4; j++) {
 205       noPID[i][j]->SetCutMgr(cutMgrAll);
 206     }
 207   }
 208
 209   // === FUNCTIONS ================================================================================
 210
 211   for (i = 0; i < nArray; i++) {
 212     for (j = 0; j < 4; j++) {
 213       AddRsnFunctionsToPair(noPID[i][j]);
 214     }
 215   }
 216
 217   // === ADD TO PAIR MANAGER ======================================================================
 218
 219   for (i = 0; i < nArray; i++) {
 220     for (j = 0; j < 4; j++) {
 221       pairMgr->AddPair(noPID[i][j]);
 222     }
 223   }
 224
 225   return pairMgr;
 226 }
 227
 228 AliRsnPairManager* CreatePairsPID
 229 (
 230   const char            *pairMgrName,    // name for the pair manager
 231   Int_t                  resonancePDG,   // PDG code of resonance (for true pairs)
 232   AliPID::EParticleType  type1,          // PID of one member of decay (+)
 233   AliPID::EParticleType  type2           // PID of other member of decay (-)
 234 )
 235 {
 236 //
 237 // Creates an AliRsnPairMgr for a specified resonance, which contains:
 238 // - signal (inv. mass)
 239 // - event mixing (inv. mass)
 240 // - like-sign (inv. mass)
 241 // - true pairs (inv. mass, resolution)
 242 //
 243 // For all pairs, a binning in Pt and Eta is provided, and a cut in multiplicity
 244 // which defines a multiplicity bin where the analysis is computed.
 245 //
 246 // Arguments define how the pair manager must be built, and are explained above
 247 //
 248
 249   AliRsnPairManager  *pairMgr  = new AliRsnPairManager(pairMgrName);
 250
 251   // === PAIR DEFINITIONS =========================================================================
 252
 253   // if particle #1 and #2 are different, two histograms must be built
 254   // for each scheme (signal, true, mixed, like-sign) exchanging both particles and signs
 255   Int_t i, j, nArray = 1;
 256   if (type1 != type2) nArray = 2;
 257
 258   AliRsnPairDef *defUnlike[2] = {0, 0};
 259   AliRsnPairDef *defLikePP[2] = {0, 0};
 260   AliRsnPairDef *defLikeMM[2] = {0, 0};
 261
 262   defUnlike[0] = new AliRsnPairDef(type1, '+', type2, '-', resonancePDG);
 263   defLikePP[0] = new AliRsnPairDef(type1, '+', type2, '+', resonancePDG);
 264   defLikeMM[0] = new AliRsnPairDef(type1, '-', type2, '-', resonancePDG);
 265
 266   defUnlike[1] = new AliRsnPairDef(type2, '+', type1, '-', resonancePDG);
 267   defLikePP[1] = new AliRsnPairDef(type2, '+', type1, '+', resonancePDG);
 268   defLikeMM[1] = new AliRsnPairDef(type2, '-', type1, '-', resonancePDG);
 269
 270   // === PAIR ANALYSIS ENGINES ====================================================================
 271
 272   // define null (dummy) objects and initialize only the ones which are needed,
 273   // depending again on particle types;
 274   // array is organized as follows:
 275   // [0] - true pairs
 276   // [1] - signal
 277   // [2] - like PP
 278   // [3] - like MM
 279   AliRsnPair *perfectPID[2][4]   = {0,0,0,0,0,0,0,0};
 280   AliRsnPair *realisticPID[2][4] = {0,0,0,0,0,0,0,0};
 281
 282   for (i = 0; i < nArray; i++) {
 283     perfectPID[i][0] = new AliRsnPair(AliRsnPair::kPerfectPID, defUnlike[i]);
 284     perfectPID[i][1] = new AliRsnPair(AliRsnPair::kPerfectPID, defUnlike[i]);
 285     perfectPID[i][2] = new AliRsnPair(AliRsnPair::kPerfectPID, defLikePP[i]);
 286     perfectPID[i][3] = new AliRsnPair(AliRsnPair::kPerfectPID, defLikeMM[i]);
 287
 288     realisticPID[i][0] = new AliRsnPair(AliRsnPair::kRealisticPID, defUnlike[i]);
 289     realisticPID[i][1] = new AliRsnPair(AliRsnPair::kRealisticPID, defUnlike[i]);
 290     realisticPID[i][2] = new AliRsnPair(AliRsnPair::kRealisticPID, defLikePP[i]);
 291     realisticPID[i][3] = new AliRsnPair(AliRsnPair::kRealisticPID, defLikeMM[i]);
 292   }
 293
 294   // === CUTS =====================================================================================
 295
 296   // cuts for tracks:
 297   // -- nothing
 298
 299   // cuts on pairs:
 300   // -- true daughters of a phi resonance (only for true pairs histogram)cutSetPairTrue->AddCut(cutTrue);
 301   AliRsnCutStd *cutTruePair = new AliRsnCutStd("cutTrue", AliRsnCutStd::kTruePair, resonancePDG);
 302
 303   // cut set definition for true pairs
 304   AliRsnCutSet *cutSetPairTrue = new AliRsnCutSet("truePairs");
 305   cutSetPairTrue->AddCut(cutTruePair);
 306   cutSetPairTrue->SetCutScheme("cutTrue");
 307
 308   // cut manager for all pairs
 309   AliRsnCutMgr *cutMgrAll = new AliRsnCutMgr("std", "All");
 310
 311   // cut manager for all pairs
 312   AliRsnCutMgr *cutMgrTrue = new AliRsnCutMgr("true", "True");
 313   cutMgrTrue->SetCutSet(AliRsnCut::kPair, cutSetPairTrue);
 314
 315   for (i = 0; i < nArray; i++) {
 316     perfectPID[i][0]->SetCutMgr(cutMgrTrue);
 317     realisticPID[i][0]->SetCutMgr(cutMgrTrue);
 318     for (j = 1; j < 4; j++) {
 319       perfectPID[i][j]->SetCutMgr(cutMgrAll);
 320       realisticPID[i][j]->SetCutMgr(cutMgrAll);
 321     }
 322   }
 323
 324   // === FUNCTIONS ================================================================================
 325
 326   for (i = 0; i < nArray; i++) {
 327     for (j = 0; j < 4; j++) {
 328       AddRsnFunctionsToPair(perfectPID[i][j]);
 329       AddRsnFunctionsToPair(realisticPID[i][j]);
 330     }
 331   }
 332
 333   // === ADD TO PAIR MANAGER ======================================================================
 334
 335   for (i = 0; i < nArray; i++) {
 336     for (j = 0; j < 4; j++) {
 337       pairMgr->AddPair(perfectPID[i][j]);
 338       pairMgr->AddPair(realisticPID[i][j]);
 339     }
 340   }
 341
 342   return pairMgr;
 343 }
 344