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