]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/macros/AddTaskJets.C
Allow for differerent selection flag than kMB, take into account filterbit in the...
[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
5bf1593b 4AliAnalysisTaskJets *AddTaskJets(Char_t *jr, Char_t *jf,Float_t radius = -1,UInt_t filterMask = 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
5bf1593b 8AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask ){\r
9902ce1e 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
5bf1593b 20 return AddTaskJets("AOD","UA1",0.4,filterMask);\r
9902ce1e 21\r
22}\r
23\r
24\r
25\r
121e4bd5 26Int_t AddTaskJetsDelta(char *nonStdFile,UInt_t filterMask,Bool_t kUseAODMC,UInt_t runFlag){\r
b486b6c5 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
f6d6265b 53 TString outFile(nonStdFile);\r
b486b6c5 54\r
b486b6c5 55 AliAnalysisTaskJets *jetana = 0;\r
56 Int_t iCount = 0;\r
57\r
121e4bd5 58 // Jet Fidners Selected by run flag first bit 2^0 second by 2^1 etc\r
784b191d 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
0ff308fd 62 // flag[5] = 0; // set siscone to 0 for proof mode...\r
3a34b9a9 63 // flag first bit AOD, second bit AODMC2 third bit AODMC2 third (8) bit AOODMC2b (limited acceptance)\r
b486b6c5 64 // i.e. 7 all, 6 only MC2 and MC\r
65 // this stay at three\r
3a34b9a9 66 const char *cReader[4] = {"AOD","AODMC","AODMC2","AODMC2b"}; \r
b486b6c5 67\r
121e4bd5 68 \r
69\r
3a34b9a9 70 for(int i = 0; i< 14;i++){\r
121e4bd5 71 if(!(runFlag&(1<<i)))continue;\r
72 if(!kUseAODMC)flag[i]&=1; // switch OFF MC if we do not have it\r
3a34b9a9 73 for(int ib = 0;ib<4;ib++){ \r
b486b6c5 74 if(flag[i]&(1<<ib)){\r
5bf1593b 75 jetana = AddTaskJets(cReader[ib],cJF[i],radius[i],filterMask);\r
b486b6c5 76 if(jetana){\r
77 char *cRadius = "";\r
8a4430e1 78 if(radius[i]>0)cRadius = Form("%02d",(int)((radius[i]+0.01)*10.)); // add an offset beacuse of precision\r
f4279ef1 79 // jetana->SetNonStdBranch(Form("jets%s_%s%s",cReader[ib],cJF[i],cRadius)); // done in addtask jets\r
f6d6265b 80 if(outFile.Length()>0)jetana->SetNonStdOutputFile(outFile.Data());\r
b486b6c5 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
5bf1593b 95AliAnalysisTaskJets *AddTaskJets(Char_t *jr, Char_t *jf, Float_t radius,UInt_t filterMask)\r
672f1183 96{\r
0651dd18 97 // Creates a jet finder task, configures it and adds it to the analysis manager.\r
672f1183 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
0511798f 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
672f1183 121 // Create the task and configure it.\r
122 //===========================================================================\r
123 AliAnalysisTaskJets *jetana;\r
5bf1593b 124 AliJetReader *er = CreateJetReader(jr,filterMask);\r
672f1183 125 // Define jet header and jet finder\r
76f9015f 126 AliJetFinder *jetFinder = CreateJetFinder(jf,radius);\r
0651dd18 127\r
672f1183 128 if (jetFinder){\r
672f1183 129 if (er) jetFinder->SetJetReader(er);\r
130 }\r
131\r
76f9015f 132 char *cRadius = "";\r
8a4430e1 133 if(radius>0)cRadius = Form("%02d",(int)((radius+0.01)*10.));\r
76f9015f 134\r
b486b6c5 135 jetana = new AliAnalysisTaskJets(Form("JetAnalysis%s_%s%s",jr,jf,cRadius));\r
58ef946c 136\r
5325cd71 137 TString cAdd = "";\r
39e7e8ab 138 cAdd += Form("_Filter%05d",filterMask);\r
f4279ef1 139\r
be60d38b 140 TString c_jr(jr);\r
141 c_jr.ToLower();\r
142 TString c_jf(jf);\r
143 c_jf.ToLower();\r
144\r
334e8d90 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
27cdf9c4 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
f4279ef1 154\r
334e8d90 155\r
5325cd71 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
672f1183 170 // Connect jet finder to task.\r
171 jetana->SetJetFinder(jetFinder);\r
172 jetana->SetConfigFile("");\r
23adcac8 173 jetana->SetDebugLevel(2);\r
0b64955a 174 if(TMath::Abs((radius-0.4))< 0.02&&c_jf.Contains("fastjet")){\r
39e7e8ab 175 jetana->SetFilterPt(100.);\r
0b64955a 176 }\r
177\r
178\r
672f1183 179 mgr->AddTask(jetana);\r
180\r
0b64955a 181\r
672f1183 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
0651dd18 192\r
76f9015f 193AliJetFinder *CreateJetFinder(Char_t *jf,Float_t radius){\r
70144f07 194 AliJetFinder *jetFinder = 0;\r
0651dd18 195\r
196 switch (jf) {\r
197 case "CDF":\r
198 AliCdfJetHeader *jh = new AliCdfJetHeader();\r
199 jh->SetRadius(0.7);\r
8a4430e1 200 if(radius>0)jh->SetRadius(radius); \r
5bf1593b 201 jh->SetAODwrite(kTRUE);\r
58ef946c 202 jh->SetAODtracksWrite(kTRUE);\r
5bf1593b 203 // jh->SetDebugCDF(1);\r
0651dd18 204 jetFinder = new AliCdfJetFinder();\r
0651dd18 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
f25f827d 212// jh->SetNeff(200);\r
213// jh->SetEtaEff(2.2);\r
8a4430e1 214 if(radius>0)jh->SetRadius(radius);\r
198864bf 215 jh->SetEtMin(5.);\r
0651dd18 216 jetFinder = new AliDAJetFinder();\r
217 if (jh) jetFinder->SetJetHeader(jh);\r
218 break;\r
219\r
70144f07 220 case "FASTJET":\r
b486b6c5 221 // DEFAULT is ANTI KT\r
70144f07 222 AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();\r
223 jh->SetRparam(0.4); // setup parameters \r
23adcac8 224 Double_t rBkg = 0.6;\r
70144f07 225 if(radius>0)jh->SetRparam(radius);\r
226 jh->SetAlgorithm(2); // antikt from fastjet/JetDefinition.hh\r
0651dd18 227 jetFinder = new AliFastJetFinder();\r
0ff308fd 228 jh->SetPtMin(1);\r
23adcac8 229 jh->SetRparamBkg(rBkg);\r
230 jh->SetGhostEtaMax(0.9+rBkg);\r
0651dd18 231 if (jh) jetFinder->SetJetHeader(jh);\r
232 break;\r
233\r
8a4430e1 234 case "FASTKT":\r
235 AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();\r
236 jh->SetRparam(0.4); // setup parameters \r
23adcac8 237 Double_t rBkg = 0.6;\r
8a4430e1 238 if(radius>0)jh->SetRparam(radius);\r
239 jh->SetAlgorithm(0); // kt from fastjet/JetDefinition.hh\r
0ff308fd 240 jh->SetPtMin(1);\r
23adcac8 241 jh->SetRparamBkg(rBkg);\r
242 jh->SetGhostEtaMax(0.9+rBkg);\r
8a4430e1 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
53c63663 259 jh->SetMinJetPt(5); // Ptmin of jets (GeV)\r
8a4430e1 260\r
261 //do you want to subtract BG (0 = no, 1 = yes)\r
262 jh->SetBGMode(0);\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
0651dd18 278 case "UA1":\r
279 AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();\r
280 jh->SetComment("UA1 jet code with default parameters");\r
ba146e16 281 jh->BackgMode(0);\r
282 jh->SetRadius(0.4);\r
76f9015f 283 if(radius>0)jh->SetRadius(radius);\r
cfd02305 284 jh->SetEtSeed(2.);\r
c40ab9fa 285 jh->SetNAcceptJets(6);\r
0651dd18 286 jh->SetLegoNbinPhi(432);\r
287 jh->SetLegoNbinEta(274);\r
288 jh->SetLegoEtaMin(-2);\r
289 jh->SetLegoEtaMax(+2);\r
9902ce1e 290 jh->SetMinJetEt(5.);\r
ba146e16 291 jh->SetJetEtaMax(1.5);\r
292 jh->SetJetEtaMin(-1.5);\r
0651dd18 293\r
294 jetFinder = new AliUA1JetFinderV1();\r
295 if (jh) jetFinder->SetJetHeader(jh);\r
296 break;\r
0ff308fd 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
0651dd18 318\r
ba146e16 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
76f9015f 324 if(radius>0)jh->SetRadius(radius);\r
ba146e16 325 jh->SetEtSeed(4.);\r
c40ab9fa 326 jh->SetNAcceptJets(6);\r
ba146e16 327 jh->SetLegoNbinPhi(432);\r
328 jh->SetLegoNbinEta(274);\r
329 jh->SetLegoEtaMin(-2);\r
330 jh->SetLegoEtaMax(+2);\r
9902ce1e 331 jh->SetMinJetEt(5.);\r
ba146e16 332 jh->SetJetEtaMax(1.5);\r
333 jh->SetJetEtaMin(-1.5);\r
ba146e16 334 jetFinder = new AliUA1JetFinderV1();\r
335 if (jh) jetFinder->SetJetHeader(jh);\r
336 break;\r
70144f07 337 default:\r
338 Printf("\n >>>>>>> AddTaskJets Error Wrong jet finder selected\n");\r
339 break;\r
0651dd18 340 }\r
341\r
342 return jetFinder;\r
343\r
344}\r
345\r
5bf1593b 346AliJetReader *CreateJetReader(Char_t *jr,UInt_t filterMask){\r
0651dd18 347 AliJetReader *er = 0;\r
e0c120d9 348 const Float_t ptCut = 0.15 ; // cut on track p_T\r
0651dd18 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
e0c120d9 355 jrh->SetPtCut(ptCut);\r
0651dd18 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
ba146e16 361 case "MC2":\r
362 AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();\r
5dbee319 363 jrh->SetComment("MC full Kinematics spearate config charged only");\r
ba146e16 364 jrh->SetFastSimTPC(kFALSE);\r
365 jrh->SetFastSimEMCAL(kFALSE);\r
5dbee319 366 jrh->SetChargedOnly(kTRUE);\r
e0c120d9 367 jrh->SetPtCut(ptCut);\r
ba146e16 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
0651dd18 373 case "ESD":\r
374 AliJetESDReaderHeader *jrh = new AliJetESDReaderHeader();\r
375 jrh->SetComment("Testing");\r
e0c120d9 376 jrh->SetFirstEvent(0.);\r
0651dd18 377 jrh->SetLastEvent(1000);\r
e0c120d9 378 jrh->SetPtCut(ptCut);\r
0651dd18 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
e0c120d9 388 jrh->SetPtCut(0.15); // set low p_T cut of to 150 MeV\r
3456f527 389 jrh->SetTestFilterMask(32); // Change this one for a different set of cuts\r
5bf1593b 390 if(filterMask>0)jrh->SetTestFilterMask(filterMask); \r
0651dd18 391 // Define reader and set its header\r
392 er = new AliJetAODReader();\r
393 er->SetReaderHeader(jrh);\r
394 break;\r
0c43aaa9 395 case "AODMC":\r
396 AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
397 jrh->SetComment("AOD MC Reader");\r
e0c120d9 398 jrh->SetPtCut(ptCut);\r
0c43aaa9 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
5dbee319 405 case "AODMC2":\r
406 AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
407 jrh->SetComment("AOD MC Reader");\r
e0c120d9 408 jrh->SetPtCut(ptCut);\r
5dbee319 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
3a34b9a9 415 case "AODMC2b":\r
416 AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();\r
417 jrh->SetComment("AOD MC Reader");\r
e0c120d9 418 jrh->SetPtCut(ptCut);\r
3a34b9a9 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
5dbee319 425\r
0651dd18 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