Add new version of macros for RSN analysis
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / macros / train / RsnConfigPhiTPC.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 RsnConfigPhiTPC
12 (
13    AliRsnAnalysisTask *task,
14    Bool_t              isMC,
15    Bool_t              useCentrality,
16    AliRsnCutSet       *eventCuts
17 )
18 {
19    if (!task) ::Error("RsnConfigPhiTPC", "NULL task");
20    
21    // we define here a suffix to differentiate names of different setups for the same resonance
22    // and we define also the name of the list of tracks we want to select for the analysis
23    // (if will fail if no lists with this name were added to the RsnInputHandler)
24    const char *suffix     = "tpc";
25    const char *listName   = "kaonTPC";
26    Bool_t      useCharged = kTRUE;
27    Int_t       listID     = -1;
28    
29    // find the index of the corresponding list in the RsnInputHandler
30    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
31    AliMultiInputEventHandler *multi = dynamic_cast<AliMultiInputEventHandler*>(mgr->GetInputEventHandler());
32    if (multi) {
33       TObjArray *array = multi->InputEventHandlers();
34       AliRsnInputHandler *rsn = (AliRsnInputHandler*)array->FindObject("rsnInputHandler");
35       if (rsn) {
36          AliRsnDaughterSelector *sel = rsn->GetSelector();
37          listID = sel->GetID(listName, useCharged);
38       }
39    }
40    if (listID >= 0)
41       ::Info("RsnConfigPhiTPC.C", "Required list '%s' stays in position %d", listName, listID);
42    else {
43       ::Error("RsnConfigPhiTPC.C", "Required list '%s' absent in handler!", listName);
44       return kFALSE;
45    }
46             
47    // ----------------------------------------------------------------------------------------------
48    // -- DEFINITIONS -------------------------------------------------------------------------------
49    // ----------------------------------------------------------------------------------------------
50    
51    // PAIR DEFINITIONS:
52    // this contains the definition of particle species and charge for both daughters of a resonance,
53    // which are used for the following purposes:
54    // --> species is used to assign the mass to the daughter (e.g. for building invariant mass)
55    // --> charge is used to select what tracks to use when doing the computation loops
56    // When a user wants to compute a like-sign background, he must define also a pair definition
57    // for each like-sign: in case of charged track decays, we need one for ++ and one for --
58    // Last two arguments are necessary only in some cases (but it is not bad to well initialize them):
59    // --> PDG code of resonance, which is used for selecting true pairs, when needed
60    // --> nominal resonance mass, which is used for computing quantities like Y or Mt
61    AliRsnPairDef *phi_kaonP_kaonM = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '-', 333, 1.019455);
62    AliRsnPairDef *phi_kaonP_kaonP = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '+', 333, 1.019455);
63    AliRsnPairDef *phi_kaonM_kaonM = new AliRsnPairDef(AliRsnDaughter::kKaon, '-', AliRsnDaughter::kKaon, '-', 333, 1.019455);
64
65    // PAIR LOOPS:
66    // these are the objects which drive the computations and fill the output histograms
67    // each one requires to be initialized with an AliRsnPairDef object, which provided masses,
68    // last argument tells if the pair is for mixing or not (this can be also set afterwards, anyway)
69    const Int_t     nPairs = 5;
70    Bool_t          addPair[nPairs] = {1, 1, 1, 1, 1};
71    AliRsnLoopPair *phiLoop[nPairs];
72    phiLoop[0] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonM"     , suffix), phi_kaonP_kaonM, kFALSE);
73    phiLoop[1] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonM_true", suffix), phi_kaonP_kaonM, kFALSE);
74    phiLoop[2] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonM_mix" , suffix), phi_kaonP_kaonM, kTRUE );
75    phiLoop[3] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonP"     , suffix), phi_kaonP_kaonP, kFALSE);
76    phiLoop[4] = new AliRsnLoopPair(Form("%s_phi_kaonM_kaonM"     , suffix), phi_kaonM_kaonM, kFALSE);
77
78    // set additional option for true pairs
79    // 1) we select only pairs coming from the same mother, which must have the right PDG code (from pairDef)
80    // 2) we select only pairs decaying according to the right channel (from pairDef species+charge definitions)
81    phiLoop[1]->SetOnlyTrue(kTRUE);
82    phiLoop[1]->SetCheckDecay(kTRUE);
83    
84    // don't add true pairs if not MC
85    if (!isMC) addPair[1] = 0;
86    
87    // ----------------------------------------------------------------------------------------------
88    // -- PAIR CUTS ---------------------------------------------------------------------------------
89    // ----------------------------------------------------------------------------------------------
90    
91    // for pairs we define a rapidity windows, defined through a cut
92    // --> NOTE: it needs a support AliRsnPairDef from which it takes the mass
93    AliRsnValueStd *valRapidity = new AliRsnValueStd("valY", AliRsnValueStd::kPairY);
94    AliRsnCutValue *cutRapidity = new AliRsnCutValue("phi_cutY", -0.5, 0.5, isMC);
95    valRapidity->SetSupportObject(phi_kaonP_kaonM);
96    cutRapidity->SetValueObj(valRapidity);
97    
98    // add the cut to a cut set (will be simple, there's only one)
99    AliRsnCutSet *pairCuts = new AliRsnCutSet("phi_pairCuts", AliRsnTarget::kMother);
100    pairCuts->AddCut(cutRapidity);
101    pairCuts->SetCutScheme(cutRapidity->GetName());
102    
103    // ----------------------------------------------------------------------------------------------
104    // -- COMPUTED VALUES & OUTPUTS -----------------------------------------------------------------
105    // ----------------------------------------------------------------------------------------------
106    
107    // All values which should be computed are defined here and passed to the computation objects,
108    // since they define all that is computed bye each one, and, in case one output is a histogram
109    // they define the binning and range for that value
110    //
111    // NOTE:
112    // --> multiplicity bins have variable size
113    
114    Double_t mult[] = { 0.,  1.,  2.,  3.,  4.,  5.,   6.,   7.,   8.,   9.,  10., 11., 12., 13., 
115                       14., 15., 16., 17., 18., 19.,  20.,  21.,  22.,  23.,  24.,  25., 30., 35., 
116                       40., 50., 60., 70., 80., 90., 100., 120., 140., 160., 180., 200., 500.};
117    Int_t    nmult  = sizeof(mult) / sizeof(mult[0]);
118    
119    AliRsnValueStd *axisIM      = new AliRsnValueStd("pIM"  , AliRsnValueStd::kPairInvMass       ,   0.9,   1.4, 0.001);
120    AliRsnValueStd *axisRes     = new AliRsnValueStd("RES"  , AliRsnValueStd::kPairInvMassRes    ,  -0.5,   0.5, 0.001);
121    AliRsnValueStd *axisPt      = new AliRsnValueStd("PT"   , AliRsnValueStd::kPairPt            ,   0.0,   5.0, 0.1  );
122    AliRsnValueStd *axisCentV0  = new AliRsnValueStd("CNT"  , AliRsnValueStd::kEventCentralityV0 ,   0.0, 100.0, 5.0  );
123    AliRsnValueStd *axisMultESD = new AliRsnValueStd("MESD" , AliRsnValueStd::kEventMultESDCuts  , nmult, mult);
124    AliRsnValueStd *axisMultSPD = new AliRsnValueStd("MSPD" , AliRsnValueStd::kEventMultSPD      , nmult, mult);
125    AliRsnValueStd *axisMultTRK = new AliRsnValueStd("MTRK" , AliRsnValueStd::kEventMult         , nmult, mult);
126    AliRsnValueStd *axisMultMC  = new AliRsnValueStd("MMC"  , AliRsnValueStd::kEventMultMC       , nmult, mult);
127
128    // create outputs:
129    // we define one for true pairs, where we add resolution, and another without it, for all others
130    // it seems that it is much advantageous to use sparse histograms when adding more than 2 axes
131    AliRsnListOutput *out[2];
132    out[0] = new AliRsnListOutput("res"  , AliRsnListOutput::kHistoSparse);
133    out[1] = new AliRsnListOutput("nores", AliRsnListOutput::kHistoSparse);
134    
135    // add values to outputs:
136    // if centrality is required, we add it only, otherwise we add all multiplicities
137    // other axes (invmass, pt) are always added
138    for (Int_t i = 0; i < 2; i++) {
139       out[i]->AddValue(axisIM);
140       out[i]->AddValue(axisPt);
141       if (useCentrality) {
142          ::Info("RsnConfigPhiTPC.C", "Adding centrality axis");
143          out[i]->AddValue(axisCentV0);
144       } else {
145          ::Info("RsnConfigPhiTPC.C", "Adding multiplicity axes");
146          //out[i]->AddValue(axisMultESD);
147          //out[i]->AddValue(axisMultSPD);
148          out[i]->AddValue(axisMultTRK);
149          if (isMC) out[i]->AddValue(axisMultMC);
150       }
151    }
152    // resolution only in the first
153    out[0]->AddValue(axisRes);
154    
155    // ----------------------------------------------------------------------------------------------
156    // -- ADD SETTINGS TO LOOPS AND LOOPS TO TASK ---------------------------------------------------
157    // ----------------------------------------------------------------------------------------------
158    
159    for (Int_t ip = 0; ip < nPairs; ip++) {
160       // skip pairs not to be added
161       if (!addPair[ip]) continue;
162       // assign list IDs
163       phiLoop[ip]->SetListID(0, listID);
164       phiLoop[ip]->SetListID(1, listID);
165       // assign event cuts
166       phiLoop[ip]->SetEventCuts(eventCuts);
167       // assign pair cuts
168       phiLoop[ip]->SetPairCuts(pairCuts);
169       // assign outputs
170       if (ip != 1)
171          phiLoop[ip]->AddOutput(out[1]);
172       else
173          phiLoop[ip]->AddOutput(out[0]);
174       // add to task
175       task->Add(phiLoop[ip]);
176    }
177    
178    return kTRUE;
179 }