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