1 // **************************************
2 // Task used for the correction of detector effects for background fluctuations in jet spectra by the HBOM method
3 // *******************************************
6 /**************************************************************************
7 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9 * Author: The ALICE Off-line Project. *
10 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
25 #include <TInterpreter.h>
27 #include <TRefArray.h>
36 #include <TLorentzVector.h>
37 #include <TClonesArray.h>
38 #include "TDatabasePDG.h"
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"
52 #include "AliGenPythiaEventHeader.h"
53 #include "AliGenCocktailEventHeader.h"
54 #include "AliInputEventHandler.h"
55 #include "AliAODJetEventBackground.h"
62 ClassImp(AliAnalysisTaskJetHBOM)
64 AliAnalysisTaskJetHBOM::~AliAnalysisTaskJetHBOM(){
74 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
75 delete fTCARandomConesOut;
76 fTCARandomConesOut = 0;
80 AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM():
85 fUseAODTrackInput(kFALSE),
86 fUseAODMCInput(kFALSE),
87 fEventSelection(kFALSE),
92 fTrackTypeRec(kTrackUndef),
93 fTrackTypeGen(kTrackUndef),
101 fTrackEtaWindow(0.9),
103 fJetOutputMinPt(0.150),
104 fMaxTrackPtInJet(100.),
110 fBackgroundBranch(""),
121 fUseTrMomentumSmearing(kFALSE),
122 fUseDiceEfficiency(kFALSE),
124 fAlgorithm(fastjet::kt_algorithm),
125 fStrategy(fastjet::Best),
126 fRecombScheme(fastjet::BIpt_scheme),
127 fAreaType(fastjet::active_area),
129 fActiveAreaRepeats(1),
132 fTCARandomConesOut(0x0),
138 fh1PtHardTrials(0x0),
140 fh1CentralityPhySel(0x0),
146 fh1efficiencyPt(0x0),
147 fh2efficiencyPhi(0x0),
155 AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM(const char* name):
156 AliAnalysisTaskSE(name),
160 fUseAODTrackInput(kFALSE),
161 fUseAODMCInput(kFALSE),
162 fEventSelection(kFALSE),
164 fFilterMaskBestPt(0),
167 fTrackTypeRec(kTrackUndef),
168 fTrackTypeGen(kTrackUndef),
170 fNSkipLeadingCone(0),
176 fTrackEtaWindow(0.9),
178 fJetOutputMinPt(0.150),
179 fMaxTrackPtInJet(100.),
185 fBackgroundBranch(""),
196 fUseTrMomentumSmearing(kFALSE),
197 fUseDiceEfficiency(kFALSE),
199 fAlgorithm(fastjet::kt_algorithm),
200 fStrategy(fastjet::Best),
201 fRecombScheme(fastjet::BIpt_scheme),
202 fAreaType(fastjet::active_area),
204 fActiveAreaRepeats(1),
207 fTCARandomConesOut(0x0),
213 fh1PtHardTrials(0x0),
215 fh1CentralityPhySel(0x0),
221 fh1efficiencyPt(0x0),
222 fh2efficiencyPhi(0x0),
227 DefineOutput(1, TList::Class());
232 Bool_t AliAnalysisTaskJetHBOM::Notify()
235 // Implemented Notify() to read the cross sections
236 // and number of trials from pyxsec.root
241 void AliAnalysisTaskJetHBOM::UserCreateOutputObjects()
245 // Create the output container
248 fRandom = new TRandom3(0);
254 if (fDebug > 1) printf("AnalysisTaskJetHBOM::UserCreateOutputObjects() \n");
258 if(fNonStdBranch.Length()!=0)
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
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);
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());
277 if(fNonStdFile.Length()!=0){
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);
287 // FitMomentumResolution();
290 if(!fHistList)fHistList = new TList();
291 fHistList->SetOwner();
292 PostData(1, fHistList); // post data in any case once
294 Bool_t oldStatus = TH1::AddDirectoryStatus();
295 TH1::AddDirectory(kFALSE);
300 const Int_t nBinPt = 100;
301 Double_t binLimitsPt[nBinPt+1];
302 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
304 binLimitsPt[iPt] = 0.0;
307 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
311 const Int_t nBinPhi = 90;
312 Double_t binLimitsPhi[nBinPhi+1];
313 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
315 binLimitsPhi[iPhi] = -1.*TMath::Pi();
318 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
324 const Int_t nBinEta = 40;
325 Double_t binLimitsEta[nBinEta+1];
326 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
328 binLimitsEta[iEta] = -2.0;
331 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
335 const int nChMax = 5000;
337 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
338 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
340 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
341 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
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);
348 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
350 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
351 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
353 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
354 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
356 fh1DeltapT = new TH1F("fh1DeltapT","DeltapT",100,-50,50);
357 fh1Rho = new TH1F("fh1Rho","Rho",100,0,200);
358 fh1RhoSigma = new TH1F("fh1RhoSigma","SigmaRho",40,0,40);
359 fh1PtRandCone = new TH1F("fh1PtRandCone","pt",100,0,200);
361 const Int_t saveLevel = 3; // large save level more histos
363 fHistList->Add(fh1Xsec);
364 fHistList->Add(fh1Trials);
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);
373 fHistList->Add(fh1RhoSigma);
374 fHistList->Add(fh1PtRandCone);
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));
384 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
387 TH1::AddDirectory(oldStatus);
390 void AliAnalysisTaskJetHBOM::Init()
396 if (fDebug > 1) printf("AnalysisTaskJetHBOM::Init() \n");
398 FitMomentumResolution();
402 void AliAnalysisTaskJetHBOM::UserExec(Option_t */*option*/)
405 // handle and reset the output jet branch
407 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
410 // Execute analysis for current event
412 AliESDEvent *fESD = 0;
413 if(fUseAODTrackInput){
414 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
416 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
422 // assume that the AOD is in the general output...
425 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
429 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
433 //Check if information is provided detector level effects
434 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
435 if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE;
437 Bool_t selectEvent = false;
438 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
444 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
445 TString vtxTitle(vtxAOD->GetTitle());
446 zVtx = vtxAOD->GetZ();
448 cent = fAOD->GetHeader()->GetCentrality();
449 if(physicsSelection){
450 fh1CentralityPhySel->Fill(cent);
451 fh1ZPhySel->Fill(zVtx);
453 // zVertex and centrality selection
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){
465 //centrality selection
467 if(cent<fCentCutLo||cent>fCentCutUp){
478 PostData(1, fHistList);
481 fh1Centrality->Fill(cent);
483 fh1Trials->Fill("#sum{ntrials}",1);
486 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
488 // ==== General variables needed
492 // we simply fetch the tracks/mc particles as a list of AliVParticles
494 //reconstructed particles
496 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
497 Float_t nCh = recParticles.GetEntries();
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());
504 //apply efficiency fNHBOM times
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);
510 Double_t pT = vp->Pt();
511 Double_t phi = vp->Phi();
514 Double_t efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
515 Double_t efficiencyPhi = fh2efficiencyPhi->GetBinContent(fh2efficiencyPhi->FindBin(phi,pT));
516 Double_t eff = efficiencyPt*efficiencyPhi; //over all efficiency
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);
529 vector<fastjet::PseudoJet> inputParticlesRec;
530 vector<fastjet::PseudoJet> inputParticlesRecRan;
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);
537 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
538 // we take total momentum here
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);
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();
550 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
551 Double_t pZ = pT/TMath::Tan(theta);
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);
558 jInpRan.set_user_index(i);
559 inputParticlesRecRan.push_back(jInpRan);
560 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
564 if(inputParticlesRec.size()==0){
565 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
566 PostData(1, fHistList);
571 // employ setters for these...
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);
579 //range where to compute background
580 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
582 phiMax = 2*TMath::Pi();
583 rapMax = fGhostEtamax - fRparam;
584 rapMin = - fGhostEtamax + fRparam;
585 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
588 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
589 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
592 // loop over all jets and fill information, first one is the leading jet
594 if(inclusiveJets.size()>0){
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
601 Double_t meanarea1=0.;
602 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
603 fh1RhoSigma->Fill(sigma1);// fluctuation of the background
604 background = bkg1;//sets background variable of the task to the correct value
607 // generate random cones
608 if(fTCARandomConesOut){
609 // create a random jet within the acceptance
610 Double_t etaMax = fTrackEtaWindow - fRparam;//0.9 - 0.4
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
615 // use fixed position for random Cones
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);
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);
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);
644 }// add up energy in cone
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);
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);
668 }//inclusive Jets > 0
671 AliAODJet *randCone = (AliAODJet*)fTCARandomConesOut->At(0);
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
681 if(fDebug)Printf("%s:%d No random Cone found",(char*)__FILE__,__LINE__);
686 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
688 PostData(1, fHistList);
691 void AliAnalysisTaskJetHBOM::Terminate(Option_t */*option*/)
694 // Terminate analysis
696 if (fDebug > 1) printf("AnalysisJetHBOM: Terminate() \n");
698 if(fMomResH1Fit) delete fMomResH1Fit;
699 if(fMomResH2Fit) delete fMomResH2Fit;
700 if(fMomResH3Fit) delete fMomResH3Fit;
705 Int_t AliAnalysisTaskJetHBOM::GetListOfTracks(TList *list,Int_t type){
708 // get list of tracks/particles for different types
711 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
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();
720 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
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());
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());
738 if(tr->Pt()<fTrackPtCut){
739 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
742 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
747 if(type==kTrackAODextra || type==kTrackAODextraonly) {
748 AliAODEvent *aod = 0;
749 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
750 else aod = AODEvent();
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;
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;
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;
787 else if(type == kTrackKineCharged){
788 if(part->Particle()->GetPDG()->Charge()==0)continue;
789 if(part->Pt()<fTrackPtCut)continue;
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;
810 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
811 if(part->Charge()==0)continue;
812 if(part->Pt()<fTrackPtCut)continue;
813 if(kTrackAODMCCharged){
817 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
828 void AliAnalysisTaskJetHBOM::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
831 // set mom res profiles
834 fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
835 fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
836 fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
839 void AliAnalysisTaskJetHBOM:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
841 // set tracking efficiency histos
844 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
845 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
846 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
849 Double_t AliAnalysisTaskJetHBOM::GetMomentumSmearing(Int_t cat, Double_t pt) {
852 // Get smearing on generated momentum
855 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
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");
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);
876 Int_t bin = fMomRes->FindBin(pt);
878 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
882 if(fMomRes) delete fMomRes;
887 void AliAnalysisTaskJetHBOM::FitMomentumResolution() {
889 // Fit linear function on momentum resolution at high pT
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.);
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.);
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.);