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