updates for differnt jet configs (including background), reduced number of UA1 Jets
[u/mrichter/AliRoot.git] / PWG4 / macros / AddTaskJets.C
CommitLineData
5bf1593b 1AliJetReader *CreateJetReader(Char_t *jr,UInt_t filterMask); // Common config\r
76f9015f 2AliJetFinder *CreateJetFinder(Char_t *jf,Float_t radius = -1);\r
0651dd18 3\r
070f06d2 4AliAnalysisTaskJets *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
625fb1d6 5Int_t AddTaskJetsDelta(char *nonStdFile = "",UInt_t filterMask = 0,Bool_t kUseAODMC = kTRUE,UInt_t runFlag = 1|4|32|128|256); \r
5bf1593b 6AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask = 0);\r
0651dd18 7\r
070f06d2 8Float_t kPtTrackMin = 0.15;\r
9Int_t kBackgroundMode = 0;\r
10\r
11\r
5bf1593b 12AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask ){\r
9902ce1e 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
5bf1593b 24 return AddTaskJets("AOD","UA1",0.4,filterMask);\r
9902ce1e 25\r
26}\r
27\r
28\r
29\r
121e4bd5 30Int_t AddTaskJetsDelta(char *nonStdFile,UInt_t filterMask,Bool_t kUseAODMC,UInt_t runFlag){\r
b486b6c5 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
f6d6265b 57 TString outFile(nonStdFile);\r
b486b6c5 58\r
b486b6c5 59 AliAnalysisTaskJets *jetana = 0;\r
60 Int_t iCount = 0;\r
61\r
121e4bd5 62 // Jet Fidners Selected by run flag first bit 2^0 second by 2^1 etc\r
784b191d 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
0ff308fd 66 // flag[5] = 0; // set siscone to 0 for proof mode...\r
3a34b9a9 67 // flag first bit AOD, second bit AODMC2 third bit AODMC2 third (8) bit AOODMC2b (limited acceptance)\r
b486b6c5 68 // i.e. 7 all, 6 only MC2 and MC\r
69 // this stay at three\r
3a34b9a9 70 const char *cReader[4] = {"AOD","AODMC","AODMC2","AODMC2b"}; \r
b486b6c5 71\r
121e4bd5 72 \r
73\r
3a34b9a9 74 for(int i = 0; i< 14;i++){\r
121e4bd5 75 if(!(runFlag&(1<<i)))continue;\r
76 if(!kUseAODMC)flag[i]&=1; // switch OFF MC if we do not have it\r
3a34b9a9 77 for(int ib = 0;ib<4;ib++){ \r
b486b6c5 78 if(flag[i]&(1<<ib)){\r
5bf1593b 79 jetana = AddTaskJets(cReader[ib],cJF[i],radius[i],filterMask);\r
b486b6c5 80 if(jetana){\r
81 char *cRadius = "";\r
8a4430e1 82 if(radius[i]>0)cRadius = Form("%02d",(int)((radius[i]+0.01)*10.)); // add an offset beacuse of precision\r
f4279ef1 83 // jetana->SetNonStdBranch(Form("jets%s_%s%s",cReader[ib],cJF[i],cRadius)); // done in addtask jets\r
f6d6265b 84 if(outFile.Length()>0)jetana->SetNonStdOutputFile(outFile.Data());\r
b486b6c5 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
070f06d2 99AliAnalysisTaskJets *AddTaskJets(Char_t *jr, Char_t *jf, Float_t radius,UInt_t filterMask,float ptTrackMin,int iBack)\r
672f1183 100{\r
0651dd18 101 // Creates a jet finder task, configures it and adds it to the analysis manager.\r
672f1183 102\r
070f06d2 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
672f1183 118 }\r
070f06d2 119 \r
0511798f 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
070f06d2 125 \r
0511798f 126\r
070f06d2 127 // Create the task and configure it.\r
128 //===========================================================================\r
129 AliAnalysisTaskJets *jetana;\r
5bf1593b 130 AliJetReader *er = CreateJetReader(jr,filterMask);\r
070f06d2 131 // Define jet header and jet finder\r
76f9015f 132 AliJetFinder *jetFinder = CreateJetFinder(jf,radius);\r
0651dd18 133\r
672f1183 134 if (jetFinder){\r
672f1183 135 if (er) jetFinder->SetJetReader(er);\r
136 }\r
137\r
76f9015f 138\r
58ef946c 139\r
5325cd71 140 TString cAdd = "";\r
070f06d2 141 cAdd += Form("%02d",(int)((radius+0.01)*10.));\r
142 cAdd += Form("_B%d",(int)((kBackgroundMode)));\r
39e7e8ab 143 cAdd += Form("_Filter%05d",filterMask);\r
070f06d2 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
f4279ef1 148\r
be60d38b 149 TString c_jr(jr);\r
150 c_jr.ToLower();\r
151 TString c_jf(jf);\r
152 c_jf.ToLower();\r
153\r
070f06d2 154 /*\r
334e8d90 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
070f06d2 159\r
27cdf9c4 160 \r
161 }\r
070f06d2 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
f4279ef1 166\r
334e8d90 167\r
5325cd71 168 AliAnalysisDataContainer *cout_jet = mgr->CreateContainer(\r
070f06d2 169 Form("jethist_%s_%s%s",\r
5325cd71 170 c_jr.Data(),\r
171 c_jf.Data(),\r
5325cd71 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
5325cd71 178 cAdd.Data()));\r
179\r
672f1183 180 // Connect jet finder to task.\r
181 jetana->SetJetFinder(jetFinder);\r
182 jetana->SetConfigFile("");\r
23adcac8 183 jetana->SetDebugLevel(2);\r
0b64955a 184 if(TMath::Abs((radius-0.4))< 0.02&&c_jf.Contains("fastjet")){\r
80a790fd 185 // jetana->SetFilterPt(20.);\r
0b64955a 186 }\r
187\r
188\r
672f1183 189 mgr->AddTask(jetana);\r
190\r
0b64955a 191\r
672f1183 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
0651dd18 202\r
76f9015f 203AliJetFinder *CreateJetFinder(Char_t *jf,Float_t radius){\r
70144f07 204 AliJetFinder *jetFinder = 0;\r
0651dd18 205\r
206 switch (jf) {\r
207 case "CDF":\r
208 AliCdfJetHeader *jh = new AliCdfJetHeader();\r
209 jh->SetRadius(0.7);\r
8a4430e1 210 if(radius>0)jh->SetRadius(radius); \r
5bf1593b 211 jh->SetAODwrite(kTRUE);\r
58ef946c 212 jh->SetAODtracksWrite(kTRUE);\r
5bf1593b 213 // jh->SetDebugCDF(1);\r
0651dd18 214 jetFinder = new AliCdfJetFinder();\r
0651dd18 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
f25f827d 222// jh->SetNeff(200);\r
223// jh->SetEtaEff(2.2);\r
8a4430e1 224 if(radius>0)jh->SetRadius(radius);\r
198864bf 225 jh->SetEtMin(5.);\r
0651dd18 226 jetFinder = new AliDAJetFinder();\r
227 if (jh) jetFinder->SetJetHeader(jh);\r
228 break;\r
229\r
70144f07 230 case "FASTJET":\r
b486b6c5 231 // DEFAULT is ANTI KT\r
70144f07 232 AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();\r
233 jh->SetRparam(0.4); // setup parameters \r
23adcac8 234 Double_t rBkg = 0.6;\r
70144f07 235 if(radius>0)jh->SetRparam(radius);\r
236 jh->SetAlgorithm(2); // antikt from fastjet/JetDefinition.hh\r
0651dd18 237 jetFinder = new AliFastJetFinder();\r
0ff308fd 238 jh->SetPtMin(1);\r
23adcac8 239 jh->SetRparamBkg(rBkg);\r
240 jh->SetGhostEtaMax(0.9+rBkg);\r
0651dd18 241 if (jh) jetFinder->SetJetHeader(jh);\r
242 break;\r
243\r
8a4430e1 244 case "FASTKT":\r
245 AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();\r
246 jh->SetRparam(0.4); // setup parameters \r
23adcac8 247 Double_t rBkg = 0.6;\r
8a4430e1 248 if(radius>0)jh->SetRparam(radius);\r
249 jh->SetAlgorithm(0); // kt from fastjet/JetDefinition.hh\r
0ff308fd 250 jh->SetPtMin(1);\r
23adcac8 251 jh->SetRparamBkg(rBkg);\r
252 jh->SetGhostEtaMax(0.9+rBkg);\r
8a4430e1 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
53c63663 269 jh->SetMinJetPt(5); // Ptmin of jets (GeV)\r
8a4430e1 270\r
271 //do you want to subtract BG (0 = no, 1 = yes)\r
80a790fd 272 jh->SetBGMode(0); // if 1 set also the radius for the background determination..\r
8a4430e1 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
0651dd18 288 case "UA1":\r
289 AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();\r
290 jh->SetComment("UA1 jet code with default parameters");\r
ba146e16 291 jh->BackgMode(0);\r
292 jh->SetRadius(0.4);\r
76f9015f 293 if(radius>0)jh->SetRadius(radius);\r
070f06d2 294 jh->SetEtSeed(4.);\r
c40ab9fa 295 jh->SetNAcceptJets(6);\r
0651dd18 296 jh->SetLegoNbinPhi(432);\r
297 jh->SetLegoNbinEta(274);\r
298 jh->SetLegoEtaMin(-2);\r
299 jh->SetLegoEtaMax(+2);\r
9902ce1e 300 jh->SetMinJetEt(5.);\r
ba146e16 301 jh->SetJetEtaMax(1.5);\r
302 jh->SetJetEtaMin(-1.5);\r
070f06d2 303 jh->BackgMode(kBackgroundMode);\r
0651dd18 304\r
305 jetFinder = new AliUA1JetFinderV1();\r
306 if (jh) jetFinder->SetJetHeader(jh);\r
307 break;\r
0ff308fd 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
070f06d2 323 jh->BackgMode(kBackgroundMode);\r
0ff308fd 324\r
325 jetFinder = new AliUA1JetFinderV1();\r
326 if (jh) jetFinder->SetJetHeader(jh);\r
327 break;\r
328\r
329\r
0651dd18 330\r
ba146e16 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
76f9015f 336 if(radius>0)jh->SetRadius(radius);\r
ba146e16 337 jh->SetEtSeed(4.);\r
c40ab9fa 338 jh->SetNAcceptJets(6);\r
ba146e16 339 jh->SetLegoNbinPhi(432);\r
340 jh->SetLegoNbinEta(274);\r
341 jh->SetLegoEtaMin(-2);\r
342 jh->SetLegoEtaMax(+2);\r
9902ce1e 343 jh->SetMinJetEt(5.);\r
ba146e16 344 jh->SetJetEtaMax(1.5);\r
345 jh->SetJetEtaMin(-1.5);\r
070f06d2 346 jh->BackgMode(kBackgroundMode);\r
ba146e16 347 jetFinder = new AliUA1JetFinderV1();\r
348 if (jh) jetFinder->SetJetHeader(jh);\r
349 break;\r
70144f07 350 default:\r
351 Printf("\n >>>>>>> AddTaskJets Error Wrong jet finder selected\n");\r
352 break;\r
0651dd18 353 }\r
354\r
355 return jetFinder;\r
356\r
357}\r
358\r
5bf1593b 359AliJetReader *CreateJetReader(Char_t *jr,UInt_t filterMask){\r
0651dd18 360 AliJetReader *er = 0;\r
070f06d2 361 Float_t ptCut = kPtTrackMin ; // cut on track p_T\r
0651dd18 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
e0c120d9 368 jrh->SetPtCut(ptCut);\r
0651dd18 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
ba146e16 374 case "MC2":\r
375 AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();\r
5dbee319 376 jrh->SetComment("MC full Kinematics spearate config charged only");\r
ba146e16 377 jrh->SetFastSimTPC(kFALSE);\r
378 jrh->SetFastSimEMCAL(kFALSE);\r
5dbee319 379 jrh->SetChargedOnly(kTRUE);\r
e0c120d9 380 jrh->SetPtCut(ptCut);\r
ba146e16 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
0651dd18 386 case "ESD":\r
387 AliJetESDReaderHeader *jrh = new AliJetESDReaderHeader();\r
388 jrh->SetComment("Testing");\r
e0c120d9 389 jrh->SetFirstEvent(0.);\r
0651dd18 390 jrh->SetLastEvent(1000);\r
e0c120d9 391 jrh->SetPtCut(ptCut);\r
0651dd18 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
070f06d2 401 jrh->SetPtCut(ptCut); // set low p_T cut of to 150 MeV\r
3456f527 402 jrh->SetTestFilterMask(32); // Change this one for a different set of cuts\r
5bf1593b 403 if(filterMask>0)jrh->SetTestFilterMask(filterMask); \r
0651dd18 404 // Define reader and set its header\r
405 er = new AliJetAODReader();\r
406 er->SetReaderHeader(jrh);\r
407 break;\r
0c43aaa9 408 case "AODMC":\r
409 AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
410 jrh->SetComment("AOD MC Reader");\r
e0c120d9 411 jrh->SetPtCut(ptCut);\r
0c43aaa9 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
5dbee319 418 case "AODMC2":\r
419 AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
420 jrh->SetComment("AOD MC Reader");\r
e0c120d9 421 jrh->SetPtCut(ptCut);\r
5dbee319 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
3a34b9a9 428 case "AODMC2b":\r
429 AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
430 jrh->SetComment("AOD MC Reader");\r
e0c120d9 431 jrh->SetPtCut(ptCut);\r
3a34b9a9 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
5dbee319 438\r
0651dd18 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