coverity fixes
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetHBOM.cxx
CommitLineData
bcf7a476 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
c64cb1f6 58using std::cout;
59using std::endl;
60using std::vector;
61
bcf7a476 62ClassImp(AliAnalysisTaskJetHBOM)
63
64AliAnalysisTaskJetHBOM::~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
80AliAnalysisTaskJetHBOM::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),
cfcde656 97 randCone_pos(0),
98 randCone_Eta(0),
99 randCone_Phi(0),
bcf7a476 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
157AliAnalysisTaskJetHBOM::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),
cfcde656 174 randCone_pos(0),
175 randCone_Eta(0),
176 randCone_Phi(0),
bcf7a476 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
236Bool_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
245void 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
401void AliAnalysisTaskJetHBOM::Init()
402{
403 //
404 // Initialization
405 //
406
407 if (fDebug > 1) printf("AnalysisTaskJetHBOM::Init() \n");
408
409 FitMomentumResolution();
410
411}
412
413void 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
cfcde656 627 if(fTCARandomConesOut){
bcf7a476 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
cfcde656 635 // use fixed position for random Cones
636 if(randCone_pos){
637 etaC = randCone_Eta;
638 phiC = randCone_Phi;
639 }
bcf7a476 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
753void 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
767Int_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
890void 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
901void 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
911Double_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
949void 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