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