Compatibility with the Root trunk
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetHBOM.cxx
1 // **************************************
2 // Task used for the correction of detector effects for background fluctuations in jet spectra by the HBOM method
3 // *******************************************
4
5
6 /**************************************************************************
7  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
8  *                                                                        *
9  * Author: The ALICE Off-line Project.                                    *
10  * Contributors are mentioned in the code where appropriate.              *
11  *                                                                        *
12  * Permission to use, copy, modify and distribute this software and its   *
13  * documentation strictly for non-commercial purposes is hereby granted   *
14  * without fee, provided that the above copyright notice appears in all   *
15  * copies and that both the copyright notice and this permission notice   *
16  * appear in the supporting documentation. The authors make no claims     *
17  * about the suitability of this software for any purpose. It is          *
18  * provided "as is" without express or implied warranty.                  *
19  **************************************************************************/
20
21  
22 #include <TROOT.h>
23 #include <TRandom3.h>
24 #include <TSystem.h>
25 #include <TInterpreter.h>
26 #include <TChain.h>
27 #include <TRefArray.h>
28 #include <TFile.h>
29 #include <TKey.h>
30 #include <TH1F.h>
31 #include <TH2F.h>
32 #include <TH3F.h>
33 #include <TProfile.h>
34 #include <TF1.h>
35 #include <TList.h>
36 #include <TLorentzVector.h>
37 #include <TClonesArray.h>
38 #include  "TDatabasePDG.h"
39
40 #include "AliAnalysisTaskJetHBOM.h"
41 #include "AliAnalysisManager.h"
42 #include "AliESDEvent.h"
43 #include "AliAODEvent.h"
44 #include "AliAODHandler.h"
45 #include "AliAODExtension.h"
46 #include "AliAODTrack.h"
47 #include "AliAODJet.h"
48 #include "AliAODMCParticle.h"
49 #include "AliMCEventHandler.h"
50 #include "AliMCEvent.h"
51 #include "AliStack.h"
52 #include "AliGenPythiaEventHeader.h"
53 #include "AliGenCocktailEventHeader.h"
54 #include "AliInputEventHandler.h"
55 #include "AliAODJetEventBackground.h"
56
57
58 using std::cout;
59 using std::endl;
60 using std::vector;
61
62 ClassImp(AliAnalysisTaskJetHBOM)
63
64 AliAnalysisTaskJetHBOM::~AliAnalysisTaskJetHBOM(){
65   //
66   // Destructor
67   //
68
69   delete fRef;
70   delete fRandom;
71
72   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
73   delete fTCARandomConesOut;
74   
75   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
76   delete fTCARandomConesOutRan;
77
78 }
79
80 AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM(): 
81   AliAnalysisTaskSE(),
82   fAOD(0x0),
83   fAODExtension(0x0),
84   fRef(new TRefArray),
85   fUseAODTrackInput(kFALSE),
86   fUseAODMCInput(kFALSE),
87   fEventSelection(kFALSE),     
88   fFilterMask(0),
89   fFilterMaskBestPt(0),
90   fFilterType(0),
91   fJetTypes(kJet),
92   fTrackTypeRec(kTrackUndef),
93   fTrackTypeGen(kTrackUndef),  
94   fNSkipLeadingRan(0),
95   fNSkipLeadingCone(0),
96   fNRandomCones(0),
97   randCone_pos(0),
98   randCone_Eta(0),
99   randCone_Phi(0),
100   fNHBOM(0),
101   fTrackEtaWindow(0.9),    
102   //fRecEtaWindow(0.5),
103   fTrackPtCut(0.),                                                      
104   fJetOutputMinPt(0.150),
105   fMaxTrackPtInJet(100.),
106   fVtxZCut(8),
107   fVtxR2Cut(1),
108   fCentCutUp(0),
109   fCentCutLo(0),
110   fNonStdBranch(""),
111   fBackgroundBranch(""),
112   fNonStdFile(""),
113   fMomResH1(0x0),
114   fMomResH2(0x0),
115   fMomResH3(0x0),
116   fMomResH1Fit(0x0),
117   fMomResH2Fit(0x0),
118   fMomResH3Fit(0x0),
119   fhEffH1(0x0),
120   fhEffH2(0x0),
121   fhEffH3(0x0),
122   fUseTrMomentumSmearing(kFALSE),
123   fUseDiceEfficiency(kFALSE),
124   fRparam(1.0), 
125   fAlgorithm(fastjet::kt_algorithm),
126   fStrategy(fastjet::Best),
127   fRecombScheme(fastjet::BIpt_scheme),
128   fAreaType(fastjet::active_area), 
129   fGhostArea(0.01),
130   fActiveAreaRepeats(1),
131   fGhostEtamax(1.5),
132   background(0),
133   fTCARandomConesOut(0x0),
134   fTCARandomConesOutRan(0x0),
135   fRandom(0),
136   fh1Xsec(0x0),
137   fh1Trials(0x0),
138   fh1PtHard(0x0),
139   fh1PtHardNoW(0x0),  
140   fh1PtHardTrials(0x0),
141   fh1Nch(0x0),
142   fh1CentralityPhySel(0x0), 
143   fh1Centrality(0x0), 
144   fh1DeltapT(0x0),
145   fh1Rho(0x0),
146   fh1PtRandCone(0x0),
147   fh1Area(0x0),
148   fh1efficiencyPt(0x0),
149   fh2efficiencyPhi(0x0),
150   fh1ZPhySel(0x0), 
151   fh1Z(0x0), 
152   fHistList(0x0)  
153 {
154   
155 }
156
157 AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM(const char* name):
158   AliAnalysisTaskSE(name),
159   fAOD(0x0),
160   fAODExtension(0x0),
161   fRef(new TRefArray),
162   fUseAODTrackInput(kFALSE),
163   fUseAODMCInput(kFALSE),
164   fEventSelection(kFALSE),                                                        
165   fFilterMask(0),
166   fFilterMaskBestPt(0),
167   fFilterType(0),
168   fJetTypes(kJet),
169   fTrackTypeRec(kTrackUndef),
170   fTrackTypeGen(kTrackUndef),
171   fNSkipLeadingRan(0),
172   fNSkipLeadingCone(0),
173   fNRandomCones(0),
174   randCone_pos(0),
175   randCone_Eta(0),
176   randCone_Phi(0),
177   fNHBOM(0),
178   fTrackEtaWindow(0.9),    
179   //fRecEtaWindow(0.5),
180   fTrackPtCut(0.),                                                      
181   fJetOutputMinPt(0.150),
182   fMaxTrackPtInJet(100.),
183   fVtxZCut(8),
184   fVtxR2Cut(1),
185   fCentCutUp(0),
186   fCentCutLo(0),
187   fNonStdBranch(""),
188   fBackgroundBranch(""),
189   fNonStdFile(""),
190   fMomResH1(0x0),
191   fMomResH2(0x0),
192   fMomResH3(0x0),
193   fMomResH1Fit(0x0),
194   fMomResH2Fit(0x0),
195   fMomResH3Fit(0x0),
196   fhEffH1(0x0),
197   fhEffH2(0x0),
198   fhEffH3(0x0),
199   fUseTrMomentumSmearing(kFALSE),
200   fUseDiceEfficiency(kFALSE),
201   fRparam(1.0), 
202   fAlgorithm(fastjet::kt_algorithm),
203   fStrategy(fastjet::Best),
204   fRecombScheme(fastjet::BIpt_scheme),
205   fAreaType(fastjet::active_area), 
206   fGhostArea(0.01),
207   fActiveAreaRepeats(1),
208   fGhostEtamax(1.5),
209   background(0),
210   fTCARandomConesOut(0x0),
211   fTCARandomConesOutRan(0x0),
212   fRandom(0),
213   fh1Xsec(0x0),
214   fh1Trials(0x0),
215   fh1PtHard(0x0),
216   fh1PtHardNoW(0x0),  
217   fh1PtHardTrials(0x0),
218   fh1Nch(0x0),
219   fh1CentralityPhySel(0x0), 
220   fh1Centrality(0x0), 
221   fh1DeltapT(0x0),
222   fh1Rho(0x0),
223   fh1PtRandCone(0x0),
224   fh1Area(0x0),
225   fh1efficiencyPt(0x0),
226   fh2efficiencyPhi(0x0),
227   fh1ZPhySel(0x0), 
228   fh1Z(0x0), 
229   fHistList(0x0)
230 {
231  DefineOutput(1, TList::Class());  
232 }
233
234
235
236 Bool_t AliAnalysisTaskJetHBOM::Notify()
237 {
238   //
239   // Implemented Notify() to read the cross sections
240   // and number of trials from pyxsec.root
241   // 
242   return kTRUE;
243 }
244
245 void AliAnalysisTaskJetHBOM::UserCreateOutputObjects()
246 {
247
248   //
249   // Create the output container
250   //
251
252   fRandom = new TRandom3(0);
253
254
255   // Connect the AOD
256
257
258   if (fDebug > 1) printf("AnalysisTaskJetHBOM::UserCreateOutputObjects() \n");
259
260   
261
262   if(fNonStdBranch.Length()!=0)
263     {
264       // only create the output branch if we have a name
265       // Create a new branch for jets...
266       //  -> cleared in the UserExec....
267       // here we can also have the case that the brnaches are written to a separate file
268       
269       // create the branch for the random cones with the same R 
270       TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
271       if(fUseDiceEfficiency || fUseTrMomentumSmearing)
272         cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone);
273     
274       //create array for the random cones; Until now only one cone per event is used
275       if(!AODEvent()->FindListObject(cName.Data())){
276         fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
277         fTCARandomConesOut->SetName(cName.Data());
278         AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
279       }
280       // create the branch with the random for the random cones on the random event
281       cName = Form("%sRandomCone_random",fNonStdBranch.Data());
282       if(!AODEvent()->FindListObject(cName.Data())){
283         fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
284         fTCARandomConesOutRan->SetName(cName.Data());
285         AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
286       }
287       
288       if(fNonStdFile.Length()!=0){
289         // 
290         // case that we have an AOD extension we need to fetch the jets from the extended output
291         // we identify the extension aod event by looking for the branchname
292         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
293         // case that we have an AOD extension we need can fetch the background maybe from the extended output                                                                  
294         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
295       }
296     }
297
298   //  FitMomentumResolution();
299
300
301   if(!fHistList)fHistList = new TList();
302   fHistList->SetOwner();
303   PostData(1, fHistList); // post data in any case once
304
305   Bool_t oldStatus = TH1::AddDirectoryStatus();
306   TH1::AddDirectory(kFALSE);
307
308   //
309   //  Histogram
310     
311   const Int_t nBinPt = 100;
312   Double_t binLimitsPt[nBinPt+1];
313   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
314     if(iPt == 0){
315       binLimitsPt[iPt] = 0.0;
316     }
317     else {// 1.0
318       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
319     }
320   }
321   
322   const Int_t nBinPhi = 90;
323   Double_t binLimitsPhi[nBinPhi+1];
324   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
325     if(iPhi==0){
326       binLimitsPhi[iPhi] = -1.*TMath::Pi();
327     }
328     else{
329       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
330     }
331   }
332
333
334
335   const Int_t nBinEta = 40;
336   Double_t binLimitsEta[nBinEta+1];
337   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
338     if(iEta==0){
339       binLimitsEta[iEta] = -2.0;
340     }
341     else{
342       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
343     }
344   }
345
346   const int nChMax = 5000;
347
348   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
349   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
350
351   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
352   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
353
354
355   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
356   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
357   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
358
359   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
360
361   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
362   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
363
364   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
365   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
366   
367   fh1DeltapT = new TH1F("fh1DeltapT","DeltapT",100,-50,50);
368   fh1Rho = new TH1F("fh1Rho","Rho",100,0,200);
369   fh1PtRandCone = new TH1F("fh1PtRandCone","pt",100,0,200);
370   fh1Area = new TH1F("fh1Area","area",100,0,1);
371
372   const Int_t saveLevel = 3; // large save level more histos
373   if(saveLevel>0){
374     fHistList->Add(fh1Xsec);
375     fHistList->Add(fh1Trials);
376
377     fHistList->Add(fh1Nch);
378     fHistList->Add(fh1Centrality);
379     fHistList->Add(fh1CentralityPhySel);
380     fHistList->Add(fh1Z);
381     fHistList->Add(fh1ZPhySel);
382     fHistList->Add(fh1DeltapT);
383     fHistList->Add(fh1Rho);
384     fHistList->Add(fh1PtRandCone);
385     fHistList->Add(fh1Area);
386   }
387
388   // =========== Switch on Sumw2 for all histos ===========
389   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
390     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
391     if (h1){
392       h1->Sumw2();
393       continue;
394     }
395     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
396     if(hn)hn->Sumw2();
397   }
398   TH1::AddDirectory(oldStatus);
399 }
400
401 void AliAnalysisTaskJetHBOM::Init()
402 {
403   //
404   // Initialization
405   //
406
407   if (fDebug > 1) printf("AnalysisTaskJetHBOM::Init() \n");
408
409   FitMomentumResolution();
410
411 }
412
413 void AliAnalysisTaskJetHBOM::UserExec(Option_t */*option*/)
414 {
415
416   // handle and reset the output jet branch 
417
418   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
419   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
420
421   //
422   // Execute analysis for current event
423   //
424   AliESDEvent *fESD = 0;
425   if(fUseAODTrackInput){    
426     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
427     if(!fAOD){
428       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
429       return;
430     }
431     // fetch the header
432   }
433   else{
434     //  assume that the AOD is in the general output...
435     fAOD  = AODEvent();
436     if(!fAOD){
437       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
438       return;
439     }
440     if(fDebug>0){
441       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
442     }
443   }
444
445   //Check if information is provided detector level effects
446   if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
447   if(!fhEffH1 || !fhEffH2 || !fhEffH3)       fUseDiceEfficiency = kFALSE;
448   
449   Bool_t selectEvent =  false;
450   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
451
452   Float_t cent = 0;
453   Float_t zVtx  = 0;
454
455   if(fAOD){
456     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
457     TString vtxTitle(vtxAOD->GetTitle());
458     zVtx = vtxAOD->GetZ();
459
460     cent = fAOD->GetHeader()->GetCentrality();
461     if(physicsSelection){
462       fh1CentralityPhySel->Fill(cent);
463       fh1ZPhySel->Fill(zVtx);
464     }
465     // zVertex and centrality selection
466     if(fEventSelection){
467       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
468         Float_t yvtx = vtxAOD->GetY();
469         Float_t xvtx = vtxAOD->GetX();
470         Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
471         if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){//usual fVtxZCut=10 and fVtxR2Cut=1  // apply vertex cut later on 
472           if(physicsSelection){
473             selectEvent = true;
474           }
475         }
476       }
477       //centrality selection
478       if(fCentCutUp>0){
479         if(cent<fCentCutLo||cent>fCentCutUp){
480           selectEvent = false;
481         }
482       }
483     }else{
484       selectEvent = true;
485     }
486   }
487   
488
489   if(!selectEvent){
490     PostData(1, fHistList);
491     return;
492   }
493   fh1Centrality->Fill(cent);  
494   fh1Z->Fill(zVtx);
495   fh1Trials->Fill("#sum{ntrials}",1);
496   
497
498   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
499
500   // ==== General variables needed
501
502
503
504   // we simply fetch the tracks/mc particles as a list of AliVParticles
505
506   //reconstructed particles
507   TList recParticles;
508   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
509   Float_t nCh = recParticles.GetEntries(); 
510   fh1Nch->Fill(nCh);
511   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
512   //nT = GetListOfTracks(&genParticles,fTrackTypeGen);
513   //if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
514
515
516   //apply efficiency fNHBOM times
517   if(fNHBOM>0){
518     for(int particle=0;particle<recParticles.GetEntries();particle++){
519       // hier Effizienzen laden und ├╝berpr├╝fen ob das Teilchen nachgewiesen wird.
520       AliVParticle *vp = (AliVParticle*)recParticles.At(particle);
521       
522       Double_t pT = vp->Pt();
523       Double_t phi = vp->Phi();
524       
525       //load efficiency
526       Double_t efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
527       Double_t efficiencyPhi = fh2efficiencyPhi->GetBinContent(fh2efficiencyPhi->FindBin(phi,pT));
528       if(pT>10){
529         efficiencyPt = 0.857; //this is the result for the fit with pT>10; statistic is low for pT>10
530         //if the efficiency is not from fastMCInput_LHC10h_110719a.root this should be calculated new
531         if(fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(10))>0.9||fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(10))<0.8){
532           efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
533           Printf("%s:%d Wrong efficiency input. Check efficiency for pt>10GeV",(char*)__FILE__,__LINE__);
534         }
535       }
536       Double_t eff = efficiencyPt*efficiencyPhi; //over all efficiency
537       
538       // if ran<eff -> particle is detected eff^fNHBOM = efficiency to detect particle fNHBOM times
539       Double_t ran = fRandom->Rndm();
540       if(ran>TMath::Power(eff,fNHBOM)){
541         recParticles.Remove(vp);
542       }
543     }
544   }
545
546
547   // find the jets....
548
549   vector<fastjet::PseudoJet> inputParticlesRec;
550   vector<fastjet::PseudoJet> inputParticlesRecRan;
551   
552   //randomize particles
553   AliAODJet vTmpRan(1,0,0,1);
554   for(int i = 0; i < recParticles.GetEntries(); i++){
555     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
556
557     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
558     // we take total momentum here
559
560       //Add particles to fastjet in case we are not running toy model
561       fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
562       jInp.set_user_index(i);
563       inputParticlesRec.push_back(jInp);
564
565     // the randomized input changes eta and phi, but keeps the p_T
566       Double_t pT = vp->Pt();
567       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
568       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
569       
570       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));  
571       Double_t pZ = pT/TMath::Tan(theta);
572
573       Double_t pX = pT * TMath::Cos(phi);
574       Double_t pY = pT * TMath::Sin(phi);
575       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
576       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
577
578       jInpRan.set_user_index(i);
579       inputParticlesRecRan.push_back(jInpRan);
580       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
581
582   }// recparticles
583
584   if(inputParticlesRec.size()==0){
585     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
586     PostData(1, fHistList);
587     return;
588   }
589   
590   // run fast jet
591   // employ setters for these...
592
593   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
594   fastjet::AreaType areaType =   fastjet::active_area;
595   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
596   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
597   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
598   
599   //range where to compute background
600   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
601   phiMin = 0;
602   phiMax = 2*TMath::Pi();
603   rapMax = fGhostEtamax - fRparam;
604   rapMin = - fGhostEtamax + fRparam;
605   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
606  
607
608   const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
609   const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
610
611  
612  // loop over all jets and fill information, first one is the leading jet
613
614   if(inclusiveJets.size()>0){
615     
616     //background estimates:all bckg jets wo the 2 hardest
617     vector<fastjet::PseudoJet> jets2=sortedJets;
618     if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); //removes the two jets with the highest pT; +2 is correct ro remove 2 jets
619     Double_t bkg1=0;
620     Double_t sigma1=0.;
621     Double_t meanarea1=0.;
622     clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
623     background = bkg1;//sets background variable of the task to the correct value
624
625
626     // generate random cones
627     if(fTCARandomConesOut){
628       // create a random jet within the acceptance
629       Double_t etaMax = fTrackEtaWindow - fRparam;//0.9 - 0.4
630       Int_t nCone = 0;
631       Int_t nConeRan = 0;
632       Double_t pTC = 1; // small number
633       Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
634       Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
635       // use fixed position for random Cones
636       if(randCone_pos){
637         etaC = randCone_Eta;
638         phiC = randCone_Phi;
639       }
640       // massless jet
641       Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
642       Double_t pZC = pTC/TMath::Tan(thetaC);
643       Double_t pXC = pTC * TMath::Cos(phiC);
644       Double_t pYC = pTC * TMath::Sin(phiC);
645       Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
646       AliAODJet tmpRecC (pXC,pYC,pZC, pC); 
647       
648       tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
649       if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
650       if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
651   
652       // loop over the reconstructed particles and add up the pT in the random cones
653       // maybe better to loop over randomized particles not in the real jets...
654       // but this by definition brings dow average energy in the whole  event
655       AliAODJet vTmpRanR(1,0,0,1);
656       for(int i = 0; i < recParticles.GetEntries(); i++){
657         AliVParticle *vp = (AliVParticle*)recParticles.At(i);
658         //add up energy in cone
659         if(fTCARandomConesOut){
660           AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(0);  
661           if(jC&&jC->DeltaR(vp)<fRparam){
662             if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
663             jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
664           }
665         }// add up energy in cone
666         
667         // the randomized input changes eta and phi, but keeps the p_T
668         if(fTCARandomConesOutRan){
669           Double_t pTR = vp->Pt();
670           Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
671           Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
672           
673           Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));  
674           Double_t pZR = pTR/TMath::Tan(thetaR);
675         
676           Double_t pXR = pTR * TMath::Cos(phiR);
677           Double_t pYR = pTR * TMath::Sin(phiR);
678           Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR); 
679           vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
680           
681           AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(0);  
682           if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
683             if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
684             jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
685           }
686         }
687       }// loop over recparticles
688     }  //fTCARandomConesOut
689     Float_t jetArea = fRparam*fRparam*TMath::Pi();
690     if(fTCARandomConesOut){
691       // rescale the momentum vectors for the random cones
692       AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(0);
693       if(rC){
694         Double_t etaC = rC->Eta();
695         Double_t phiC = rC->Phi();
696         // massless jet, unit vector
697         Double_t pTC = rC->ChargedBgEnergy();
698         if(pTC<=0)pTC = 0.001; // for almost empty events
699         Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
700         Double_t pZC = pTC/TMath::Tan(thetaC);
701         Double_t pXC = pTC * TMath::Cos(phiC);
702         Double_t pYC = pTC * TMath::Sin(phiC);
703         Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
704         rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
705         rC->SetBgEnergy(0,0);
706         rC->SetEffArea(jetArea,0);
707       }
708     }
709     if(fTCARandomConesOutRan){
710       AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(0);
711       // same with random
712       if(rC){
713         Double_t etaC = rC->Eta();
714         Double_t phiC = rC->Phi();
715         // massless jet, unit vector
716         Double_t pTC = rC->ChargedBgEnergy();
717         if(pTC<=0)pTC = 0.001;// for almost empty events
718         Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
719         Double_t pZC = pTC/TMath::Tan(thetaC);
720         Double_t pXC = pTC * TMath::Cos(phiC);
721         Double_t pYC = pTC * TMath::Sin(phiC);
722         Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
723         rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
724         rC->SetBgEnergy(0,0);
725         rC->SetEffArea(jetArea,0);
726       }
727     }//find the random jets
728   }//inclusive Jets > 0
729
730  //Calculate delta pT
731  AliAODJet *randCone = (AliAODJet*)fTCARandomConesOut->At(0);
732  if(randCone){
733    //background is the backbround density per area and area=pi*0.4^2 -> backgroundCone is the background energy under the cone
734    Float_t backgroundCone = background * randCone->EffectiveAreaCharged();
735    //calculates difference between expected and measured energy density
736    Float_t ptSub = randCone->Pt() - backgroundCone;
737    fh1DeltapT->Fill(ptSub);// delta pT
738    fh1Rho->Fill(background);// background rho
739    fh1PtRandCone->Fill(randCone->Pt());// pT of random cone
740    fh1Area->Fill(randCone->EffectiveAreaCharged()); // area of random cone; should always be pi*0.4^2 = 0.5
741  }else{
742    if(fDebug)Printf("%s:%d No random Cone found",(char*)__FILE__,__LINE__);
743  }
744  
745
746  if (fDebug > 2){
747    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
748    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
749  }
750  PostData(1, fHistList);
751 }
752
753 void AliAnalysisTaskJetHBOM::Terminate(Option_t */*option*/)
754 {
755   //
756   // Terminate analysis
757   //
758     if (fDebug > 1) printf("AnalysisJetHBOM: Terminate() \n");
759
760     if(fMomResH1Fit) delete fMomResH1Fit;
761     if(fMomResH2Fit) delete fMomResH2Fit;
762     if(fMomResH3Fit) delete fMomResH3Fit;
763     
764 }
765
766
767 Int_t  AliAnalysisTaskJetHBOM::GetListOfTracks(TList *list,Int_t type){
768
769   //
770   // get list of tracks/particles for different types
771   //
772
773   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
774
775   Int_t iCount = 0;
776    if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
777     if(type!=kTrackAODextraonly) {
778       AliAODEvent *aod = 0;
779       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
780       else aod = AODEvent();
781       if(!aod){
782         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
783         return iCount;
784       }
785
786       for(int it = 0;it < aod->GetNumberOfTracks();++it){
787         AliAODTrack *tr = aod->GetTrack(it);
788         Bool_t bGood = false;
789         if(fFilterType == 0)bGood = true;
790         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
791         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
792         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
793           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());     
794           continue;
795         }
796         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
797           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());     
798           continue;
799         }
800         if(tr->Pt()<fTrackPtCut){
801           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());      
802           continue;
803         }
804         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());        
805         list->Add(tr);
806         iCount++;
807       }
808     }
809     if(type==kTrackAODextra || type==kTrackAODextraonly) {
810       AliAODEvent *aod = 0;
811       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
812       else aod = AODEvent();
813       
814       if(!aod){
815         return iCount;
816       }
817       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
818       if(!aodExtraTracks)return iCount;
819       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
820         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
821         if (!track) continue;
822
823         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
824         if(!trackAOD)continue;
825         Bool_t bGood = false;
826         if(fFilterType == 0)bGood = true;
827         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
828         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
829         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
830         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
831         if(trackAOD->Pt()<fTrackPtCut) continue;
832         list->Add(trackAOD);
833         iCount++;
834       }
835     }
836   }
837   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
838     AliMCEvent* mcEvent = MCEvent();
839     if(!mcEvent)return iCount;
840     // we want to have alivpartilces so use get track
841     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
842       if(!mcEvent->IsPhysicalPrimary(it))continue;
843       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
844       if(type == kTrackKineAll){
845         if(part->Pt()<fTrackPtCut)continue;
846         list->Add(part);
847         iCount++;
848       }
849       else if(type == kTrackKineCharged){
850         if(part->Particle()->GetPDG()->Charge()==0)continue;
851         if(part->Pt()<fTrackPtCut)continue;
852         list->Add(part);
853         iCount++;
854       }
855     }
856   }
857   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
858     AliAODEvent *aod = 0;
859     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
860     else aod = AODEvent();
861     if(!aod)return iCount;
862     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
863     if(!tca)return iCount;
864     for(int it = 0;it < tca->GetEntriesFast();++it){
865       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
866       if(!part->IsPhysicalPrimary())continue;
867       if(type == kTrackAODMCAll){
868         if(part->Pt()<fTrackPtCut)continue;
869         list->Add(part);
870         iCount++;
871       }
872       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
873         if(part->Charge()==0)continue;
874         if(part->Pt()<fTrackPtCut)continue;
875         if(kTrackAODMCCharged){
876           list->Add(part);
877         }
878         else {
879           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
880           list->Add(part);
881         }
882         iCount++;
883       }
884     }
885   }// AODMCparticle
886   list->Sort();
887   return iCount;
888 }
889
890 void AliAnalysisTaskJetHBOM::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
891
892   //
893   // set mom res profiles
894   //
895
896   fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
897   fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
898   fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
899 }
900
901 void AliAnalysisTaskJetHBOM:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
902   //
903   // set tracking efficiency histos
904   //
905
906   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
907   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
908   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
909 }
910
911 Double_t AliAnalysisTaskJetHBOM::GetMomentumSmearing(Int_t cat, Double_t pt) {
912
913   //
914   // Get smearing on generated momentum
915   //
916
917   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
918
919   TProfile *fMomRes = 0x0;
920   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
921   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
922   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
923
924   if(!fMomRes) {
925     return 0.;
926   }
927
928
929   Double_t smear = 0.;
930
931   if(pt>20.) {
932     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
933     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
934     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
935   }
936   else {
937
938     Int_t bin = fMomRes->FindBin(pt);
939
940     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
941
942   }
943
944   if(fMomRes) delete fMomRes;
945   
946   return smear;
947 }
948
949 void AliAnalysisTaskJetHBOM::FitMomentumResolution() {
950   //
951   // Fit linear function on momentum resolution at high pT
952   //
953
954   if(!fMomResH1Fit && fMomResH1) {
955     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
956     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
957     fMomResH1Fit ->SetRange(5.,100.);
958   }
959
960   if(!fMomResH2Fit && fMomResH2) {
961     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
962     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
963     fMomResH2Fit ->SetRange(5.,100.);
964   }
965
966   if(!fMomResH3Fit && fMomResH3) {
967     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
968     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
969     fMomResH3Fit ->SetRange(5.,100.);
970   }
971
972 }
973