25d4cd214e2b49f5429a3ba27d11f0e7748f4e30
[u/mrichter/AliRoot.git] / PWG4 / macros / AddTaskJets.C
1 AliJetReader *CreateJetReader(Char_t *jr,UInt_t filterMask); // Common config\r
2 AliJetFinder *CreateJetFinder(Char_t *jf,Float_t radius = -1);\r
3 \r
4 AliAnalysisTaskJets *AddTaskJets(Char_t *jr, Char_t *jf,Float_t radius = -1,UInt_t filterMask = 0); // for the new AF\r
5 Int_t AddTaskJetsDelta(char *nonStdFile = "",UInt_t filterMask = 0,Bool_t kUseAODMC = kTRUE,UInt_t runFlag = 1|4|32|128|256);     \r
6 AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask = 0);\r
7 \r
8 AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask ){\r
9   // fills the standard "jets" branch in the AOD\r
10   // need the ESDFilter to run before, to access the AODtracks\r
11   // Tracks selected by the first Filter (1<<0)\r
12   // needs to be adapted for different cuts\r
13   \r
14   // UA1 as standard chosen, since it is the most robust and simple JF\r
15   // R = 0.4 suffficient to provide accurate jet axis for correlation studies\r
16   // energy resolution suffers a little\r
17   // Acceptance of jets not limited by the Jet Finder but should be done\r
18   // by user to abs(eta) < 0.5 \r
19 \r
20   return AddTaskJets("AOD","UA1",0.4,filterMask);\r
21 \r
22 }\r
23 \r
24 \r
25 \r
26 Int_t AddTaskJetsDelta(char *nonStdFile,UInt_t filterMask,Bool_t kUseAODMC,UInt_t runFlag){\r
27 \r
28   // Adds a whole set of jet finders  all to be written\r
29   // to a delta AOD\r
30   \r
31   // this can in principle be done also with on the fly \r
32   // if we have an ESD input jet fidner task does automatically fetch the ouput aod\r
33 \r
34    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
35    if (!mgr) {\r
36       ::Error("AddTaskJetsDelta", "No analysis manager to connect to.");\r
37       return 0;\r
38    }  \r
39    \r
40    // Check the analysis type using the event handlers connected to the analysis manager.\r
41    //==============================================================================\r
42    if (!mgr->GetInputEventHandler()) {\r
43       ::Error("AddTaskJetsDelta", "This task requires an input event handler");\r
44       return 0;\r
45    }\r
46 \r
47   AliAODHandler *aodh = (AliAODHandler*)mgr->GetOutputEventHandler();\r
48   if (!aodh) {\r
49     ::Error("AddTaskJetsDelta", "This task needs an output event handler");\r
50     return 0;\r
51   }   \r
52 \r
53   TString outFile(nonStdFile);\r
54 \r
55   AliAnalysisTaskJets *jetana = 0;\r
56   Int_t iCount = 0;\r
57 \r
58   // Jet Fidners Selected by run flag first bit 2^0 second by 2^1 etc\r
59   const char *cJF[14]        = {"UA1","UA1","UA1", "CDF", "CDF","DA","DA","SISCONE","SISCONE","FASTJET","FASTJET","FASTKT","FASTKT","UA1LO"};\r
60   const Float_t radius[14]   = {  0.4,  0.7,  1.0,  0.4,    0.7, 0.4, 0.7,      0.4,     0.7,      0.4,       0.7,     0.4,     0.7,    0.4};\r
61   UInt_t  flag[14]           = {    6,    7,    7,    7,      7,   7,   7,        7,       7,        7,         7,       7,       7,      7};\r
62   // flag[5] = 0; // set siscone to 0 for proof mode...\r
63   // flag first bit AOD, second bit AODMC2 third bit AODMC2 third (8) bit AOODMC2b (limited acceptance)\r
64   // i.e. 7 all, 6 only MC2 and MC\r
65   // this stay at three\r
66   const char *cReader[4] = {"AOD","AODMC","AODMC2","AODMC2b"};  \r
67 \r
68   \r
69 \r
70   for(int i = 0; i< 14;i++){\r
71     if(!(runFlag&(1<<i)))continue;\r
72     if(!kUseAODMC)flag[i]&=1; // switch OFF MC if we do not have it\r
73     for(int ib = 0;ib<4;ib++){      \r
74       if(flag[i]&(1<<ib)){\r
75         jetana = AddTaskJets(cReader[ib],cJF[i],radius[i],filterMask);\r
76         if(jetana){\r
77           char *cRadius = "";\r
78           if(radius[i]>0)cRadius = Form("%02d",(int)((radius[i]+0.01)*10.)); // add an offset beacuse of precision\r
79           //      jetana->SetNonStdBranch(Form("jets%s_%s%s",cReader[ib],cJF[i],cRadius)); // done in addtask jets\r
80           if(outFile.Length()>0)jetana->SetNonStdOutputFile(outFile.Data());\r
81           iCount++;\r
82         }\r
83       }\r
84     }\r
85   }\r
86     \r
87   Printf("Added %d JetFinders",iCount);\r
88 }\r
89 \r
90 \r
91 \r
92 \r
93 \r
94 \r
95 AliAnalysisTaskJets *AddTaskJets(Char_t *jr, Char_t *jf, Float_t radius,UInt_t filterMask)\r
96 {\r
97   // Creates a jet finder task, configures it and adds it to the analysis manager.\r
98 \r
99    // Get the pointer to the existing analysis manager via the static access method.\r
100    //==============================================================================\r
101    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
102    if (!mgr) {\r
103       ::Error("AddTaskJets", "No analysis manager to connect to.");\r
104       return NULL;\r
105    }  \r
106    \r
107    // Check the analysis type using the event handlers connected to the analysis manager.\r
108    //==============================================================================\r
109    if (!mgr->GetInputEventHandler()) {\r
110       ::Error("AddTaskJets", "This task requires an input event handler");\r
111       return NULL;\r
112    }\r
113 \r
114   AliAODHandler *aodh = (AliAODHandler*)mgr->GetOutputEventHandler();\r
115   if (!aodh) {\r
116     ::Error("AddTaskJets", "This task needs an output event handler");\r
117     return NULL;\r
118   }   \r
119 \r
120 \r
121    // Create the task and configure it.\r
122    //===========================================================================\r
123    AliAnalysisTaskJets *jetana;\r
124    AliJetReader *er = CreateJetReader(jr,filterMask);\r
125     // Define jet header and jet finder\r
126    AliJetFinder *jetFinder = CreateJetFinder(jf,radius);\r
127 \r
128    if (jetFinder){\r
129        if (er) jetFinder->SetJetReader(er);\r
130    }\r
131 \r
132    char *cRadius = "";\r
133    if(radius>0)cRadius = Form("%02d",(int)((radius+0.01)*10.));\r
134 \r
135    jetana = new AliAnalysisTaskJets(Form("JetAnalysis%s_%s%s",jr,jf,cRadius));\r
136 \r
137    TString cAdd = "";\r
138    cAdd += Form("_Filter%05d",filterMask);\r
139 \r
140    TString c_jr(jr);\r
141    c_jr.ToLower();\r
142    TString c_jf(jf);\r
143    c_jf.ToLower();\r
144 \r
145    if(c_jf.CompareTo("ua1")==0 && TMath::Abs(radius-0.4) < 0.01 && c_jr.CompareTo("aod") == 0){\r
146      // do nothing, this is the standard jet finder R = 0.4, UA1 on AOD\r
147    }\r
148    else{\r
149      TString bName =  Form("jets%s_%s%s%s",jr,jf,cRadius,cAdd.Data());\r
150      jetana->SetNonStdBranch(bName.Data());\r
151      Printf("Set jet branchname \"%s\"",bName.Data());\r
152      \r
153  }\r
154 \r
155 \r
156    AliAnalysisDataContainer *cout_jet = mgr->CreateContainer(\r
157                                                              Form("jethist_%s_%s%s%s",\r
158                                                                   c_jr.Data(),\r
159                                                                   c_jf.Data(),\r
160                                                                   cRadius,\r
161                                                                   cAdd.Data()), \r
162                                                              TList::Class(),\r
163                                                              AliAnalysisManager::kOutputContainer, \r
164                                                              Form("%s:PWG4_jethist_%s_%s%s",AliAnalysisManager::GetCommonFileName(),\r
165                                                                   c_jr.Data(),\r
166                                                                   c_jf.Data(),\r
167                                                                   cRadius,\r
168                                                                   cAdd.Data()));\r
169 \r
170    // Connect jet finder to task.\r
171    jetana->SetJetFinder(jetFinder);\r
172    jetana->SetConfigFile("");\r
173    jetana->SetDebugLevel(2);\r
174    if(TMath::Abs((radius-0.4))< 0.02&&c_jf.Contains("fastjet")){\r
175      //     jetana->SetFilterPt(20.);\r
176    }\r
177 \r
178 \r
179    mgr->AddTask(jetana);\r
180 \r
181 \r
182    // Create ONLY the output containers for the data produced by the task.\r
183    // Get and connect other common input/output containers via the manager as below\r
184    //==============================================================================\r
185    mgr->ConnectInput  (jetana, 0, mgr->GetCommonInputContainer());\r
186 // AOD output slot will be used in a different way in future\r
187    mgr->ConnectOutput (jetana, 0, mgr->GetCommonOutputContainer());\r
188    mgr->ConnectOutput (jetana, 1, cout_jet);\r
189    \r
190    return jetana;\r
191 }\r
192 \r
193 AliJetFinder *CreateJetFinder(Char_t *jf,Float_t radius){\r
194   AliJetFinder *jetFinder = 0;\r
195 \r
196   switch (jf) {\r
197   case "CDF":\r
198     AliCdfJetHeader *jh = new AliCdfJetHeader();\r
199     jh->SetRadius(0.7);\r
200     if(radius>0)jh->SetRadius(radius);    \r
201     jh->SetAODwrite(kTRUE);\r
202     jh->SetAODtracksWrite(kTRUE);\r
203     //    jh->SetDebugCDF(1);\r
204     jetFinder = new AliCdfJetFinder();\r
205     if (jh) jetFinder->SetJetHeader(jh);\r
206     break;\r
207 \r
208   case "DA":\r
209     AliDAJetHeader *jh=new AliDAJetHeader();\r
210     jh->SetComment("DA jet code with default parameters");\r
211     jh->SelectJets(kTRUE);\r
212 //    jh->SetNeff(200);\r
213 //    jh->SetEtaEff(2.2);\r
214     if(radius>0)jh->SetRadius(radius);\r
215     jh->SetEtMin(5.);\r
216     jetFinder = new AliDAJetFinder();\r
217     if (jh) jetFinder->SetJetHeader(jh);\r
218     break;\r
219 \r
220   case "FASTJET":\r
221     // DEFAULT is ANTI KT\r
222     AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();\r
223     jh->SetRparam(0.4); // setup parameters                                  \r
224     Double_t rBkg = 0.6;\r
225     if(radius>0)jh->SetRparam(radius);\r
226     jh->SetAlgorithm(2); // antikt from fastjet/JetDefinition.hh\r
227     jetFinder = new AliFastJetFinder();\r
228     jh->SetPtMin(1);\r
229     jh->SetRparamBkg(rBkg);\r
230     jh->SetGhostEtaMax(0.9+rBkg);\r
231     if (jh) jetFinder->SetJetHeader(jh);\r
232     break;\r
233 \r
234   case "FASTKT":\r
235     AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();\r
236     jh->SetRparam(0.4); // setup parameters                                  \r
237     Double_t rBkg = 0.6;\r
238     if(radius>0)jh->SetRparam(radius);\r
239     jh->SetAlgorithm(0); // kt from fastjet/JetDefinition.hh\r
240     jh->SetPtMin(1);\r
241     jh->SetRparamBkg(rBkg);\r
242     jh->SetGhostEtaMax(0.9+rBkg);\r
243     jetFinder = new AliFastJetFinder();\r
244     if (jh) jetFinder->SetJetHeader(jh);\r
245     break;\r
246 \r
247   case "SISCONE":\r
248     AliSISConeJetHeader * jh = new AliSISConeJetHeader();\r
249 \r
250     jh->SetJetEtaMax(1.5);\r
251     jh->SetJetEtaMin(-1.5);\r
252 \r
253     //siscone parameters\r
254     jh->SetConeRadius(0.4);                   // default cone radius\r
255     if(radius>0)jh->SetConeRadius(radius);   // cone radius\r
256 \r
257     jh->SetOverlapThreshold(0.75);            // overlap parameter, between 0 and 1 excluded!! 0.75 value is advised\r
258     jh->SetPtProtojetMin(0);                  // pt min of protojets\r
259     jh->SetMinJetPt(5);                      // Ptmin of jets (GeV)\r
260 \r
261     //do you want to subtract BG (0 = no, 1 = yes)\r
262     jh->SetBGMode(0); // if 1 set also the radius for the background determination..\r
263 \r
264     //for background\r
265     jh->SetRapRange( -0.9, 0.9);              // rapidity range for subtracting background must be < ghostmaxrap-0.95*R\r
266     jh->SetPhiRange(0 , 2*TMath::Pi());       // phi range for subtracting background\r
267     \r
268     //to determine jets area\r
269     jh->SetBGAlgorithm(1);                    // algorithm for rho calculation : 0 = kT, 1 = Cambridge\r
270     jh->SetGhostEtaMax(4);                    // eta max where ghosts are generated \r
271     jh->SetGhostArea(0.05);                   // area of a ghost \r
272     jh->SetMeanGhostKt(1e-100);               // average transverse momentum of the ghosts.\r
273     jh->SetAreaTypeNumber(4);                 // from 1 to 5 : 1 = active_area, 2 = active_area_explicit_ghosts, 3 = one_ghost_passive_area, 4 = passive_area, 5 = voronoi_area\r
274     jetFinder = new AliSISConeJetFinder();\r
275     if (jh) jetFinder->SetJetHeader(jh);\r
276     break;\r
277 \r
278   case "UA1":\r
279     AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();\r
280     jh->SetComment("UA1 jet code with default parameters");\r
281     jh->BackgMode(0);\r
282     jh->SetRadius(0.4);\r
283     if(radius>0)jh->SetRadius(radius);\r
284     jh->SetEtSeed(2.);\r
285     jh->SetNAcceptJets(6);\r
286     jh->SetLegoNbinPhi(432);\r
287     jh->SetLegoNbinEta(274);\r
288     jh->SetLegoEtaMin(-2);\r
289     jh->SetLegoEtaMax(+2);\r
290     jh->SetMinJetEt(5.);\r
291     jh->SetJetEtaMax(1.5);\r
292     jh->SetJetEtaMin(-1.5);\r
293 \r
294     jetFinder = new AliUA1JetFinderV1();\r
295     if (jh) jetFinder->SetJetHeader(jh);\r
296     break;\r
297   case "UA1LO":\r
298     AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();\r
299     jh->SetComment("UA1 jet code with Lo Pt settings parameters");\r
300     jh->BackgMode(0);\r
301     jh->SetRadius(0.4);\r
302     if(radius>0)jh->SetRadius(radius);\r
303     jh->SetEtSeed(1.);\r
304     jh->SetNAcceptJets(6);\r
305     jh->SetLegoNbinPhi(432);\r
306     jh->SetLegoNbinEta(274);\r
307     jh->SetLegoEtaMin(-2);\r
308     jh->SetLegoEtaMax(+2);\r
309     jh->SetMinJetEt(1.);\r
310     jh->SetJetEtaMax(1.5);\r
311     jh->SetJetEtaMin(-1.5);\r
312 \r
313     jetFinder = new AliUA1JetFinderV1();\r
314     if (jh) jetFinder->SetJetHeader(jh);\r
315     break;\r
316 \r
317 \r
318 \r
319   case "UA1MC":\r
320     AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();\r
321     jh->SetComment("UA1 jet code with default MC parameters");\r
322     jh->BackgMode(0);\r
323     jh->SetRadius(1.0);\r
324     if(radius>0)jh->SetRadius(radius);\r
325     jh->SetEtSeed(4.);\r
326     jh->SetNAcceptJets(6);\r
327     jh->SetLegoNbinPhi(432);\r
328     jh->SetLegoNbinEta(274);\r
329     jh->SetLegoEtaMin(-2);\r
330     jh->SetLegoEtaMax(+2);\r
331     jh->SetMinJetEt(5.);\r
332     jh->SetJetEtaMax(1.5);\r
333     jh->SetJetEtaMin(-1.5);\r
334     jetFinder = new AliUA1JetFinderV1();\r
335     if (jh) jetFinder->SetJetHeader(jh);\r
336     break;\r
337   default:\r
338     Printf("\n >>>>>>> AddTaskJets Error Wrong jet finder selected\n");\r
339     break;\r
340   }\r
341 \r
342   return jetFinder;\r
343 \r
344 }\r
345 \r
346 AliJetReader *CreateJetReader(Char_t *jr,UInt_t filterMask){\r
347   AliJetReader *er = 0;\r
348   const Float_t ptCut  = 0.15 ; // cut on track p_T\r
349   switch (jr) {\r
350   case "MC":\r
351     AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();\r
352     jrh->SetComment("MC full Kinematics");\r
353     jrh->SetFastSimTPC(kFALSE);\r
354     jrh->SetFastSimEMCAL(kFALSE);\r
355     jrh->SetPtCut(ptCut);\r
356     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0 .9                                                                             \r
357     // Define reader and set its header                                     \r
358     er = new AliJetKineReader();\r
359     er->SetReaderHeader(jrh);\r
360     break;\r
361   case "MC2":\r
362     AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();\r
363     jrh->SetComment("MC full Kinematics spearate config charged only");\r
364     jrh->SetFastSimTPC(kFALSE);\r
365     jrh->SetFastSimEMCAL(kFALSE);\r
366     jrh->SetChargedOnly(kTRUE);\r
367     jrh->SetPtCut(ptCut);\r
368     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0 .9                                                                             \r
369     // Define reader and set its header                                     \r
370     er = new AliJetKineReader();\r
371     er->SetReaderHeader(jrh);\r
372     break;\r
373   case "ESD":\r
374     AliJetESDReaderHeader *jrh = new AliJetESDReaderHeader();\r
375     jrh->SetComment("Testing");\r
376     jrh->SetFirstEvent(0.);\r
377     jrh->SetLastEvent(1000);\r
378     jrh->SetPtCut(ptCut);\r
379     jrh->SetReadSignalOnly(kFALSE);\r
380     // Define reader and set its header                                     \r
381     er = new AliJetESDReader();\r
382     er->SetReaderHeader(jrh);\r
383     break;\r
384 \r
385   case "AOD":\r
386     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
387     jrh->SetComment("AOD Reader");\r
388     jrh->SetPtCut(0.15); // set low p_T cut of to 150 MeV\r
389     jrh->SetTestFilterMask(32); // Change this one for a different set of cuts\r
390     if(filterMask>0)jrh->SetTestFilterMask(filterMask); \r
391     // Define reader and set its header\r
392     er = new AliJetAODReader();\r
393     er->SetReaderHeader(jrh);\r
394     break;\r
395   case "AODMC":\r
396     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
397     jrh->SetComment("AOD MC Reader");\r
398     jrh->SetPtCut(ptCut);\r
399     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9\r
400     jrh->SetReadAODMC(1);// 1 all primary MC , 2 all primary charged\r
401     // Define reader and set its header\r
402     er = new AliJetAODReader();\r
403     er->SetReaderHeader(jrh);\r
404     break;\r
405   case "AODMC2":\r
406     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
407     jrh->SetComment("AOD MC Reader");\r
408     jrh->SetPtCut(ptCut);\r
409     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9\r
410     jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged\r
411     // Define reader and set its header\r
412     er = new AliJetAODReader();\r
413     er->SetReaderHeader(jrh);\r
414     break;\r
415   case "AODMC2b":\r
416     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
417     jrh->SetComment("AOD MC Reader");\r
418     jrh->SetPtCut(ptCut);\r
419     jrh->SetFiducialEta(-0.9,0.9); // to take all MC particles default is 0.9\r
420     jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged\r
421     // Define reader and set its header\r
422     er = new AliJetAODReader();\r
423     er->SetReaderHeader(jrh);\r
424     break;\r
425 \r
426   default:\r
427     ::Error("AddTaskJets", "Wrong jet reader selected\n");\r
428     return 0;\r
429   }\r
430 \r
431   return er;\r
432 \r
433 }\r