+++ /dev/null
-// **************************************
-// Task used for the correction of detector effects for background fluctuations in jet spectra by the HBOM method
-// *******************************************
-
-
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-
-#include <TROOT.h>
-#include <TRandom3.h>
-#include <TSystem.h>
-#include <TInterpreter.h>
-#include <TChain.h>
-#include <TRefArray.h>
-#include <TFile.h>
-#include <TKey.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TH3F.h>
-#include <TProfile.h>
-#include <TF1.h>
-#include <TList.h>
-#include <TLorentzVector.h>
-#include <TClonesArray.h>
-#include "TDatabasePDG.h"
-
-#include "AliAnalysisTaskJetHBOM.h"
-#include "AliAnalysisManager.h"
-#include "AliESDEvent.h"
-#include "AliAODEvent.h"
-#include "AliAODHandler.h"
-#include "AliAODExtension.h"
-#include "AliAODTrack.h"
-#include "AliAODJet.h"
-#include "AliAODMCParticle.h"
-#include "AliMCEventHandler.h"
-#include "AliMCEvent.h"
-#include "AliStack.h"
-#include "AliGenPythiaEventHeader.h"
-#include "AliGenCocktailEventHeader.h"
-#include "AliInputEventHandler.h"
-#include "AliAODJetEventBackground.h"
-
-
-using std::cout;
-using std::endl;
-using std::vector;
-
-ClassImp(AliAnalysisTaskJetHBOM)
-
-AliAnalysisTaskJetHBOM::~AliAnalysisTaskJetHBOM(){
- //
- // Destructor
- //
-
- delete fRef;
- fRef = 0;
- delete fRandom;
- fRandom = 0;
-
- if(fTCARandomConesOut)fTCARandomConesOut->Delete();
- delete fTCARandomConesOut;
- fTCARandomConesOut = 0;
-
-}
-
-AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM():
- AliAnalysisTaskSE(),
- fAOD(0x0),
- fAODExtension(0x0),
- fRef(0),
- fUseAODTrackInput(kFALSE),
- fUseAODMCInput(kFALSE),
- fEventSelection(kFALSE),
- fFilterMask(0),
- fFilterMaskBestPt(0),
- fFilterType(0),
- fJetTypes(kJet),
- fTrackTypeRec(kTrackUndef),
- fTrackTypeGen(kTrackUndef),
- fNSkipLeadingRan(0),
- fNSkipLeadingCone(0),
- fNRandomCones(0),
- randCone_pos(0),
- randCone_Eta(0),
- randCone_Phi(0),
- fNHBOM(0),
- fTrackEtaWindow(0.9),
- fTrackPtCut(0.),
- fJetOutputMinPt(0.150),
- fMaxTrackPtInJet(100.),
- fVtxZCut(8),
- fVtxR2Cut(1),
- fCentCutUp(0),
- fCentCutLo(0),
- fNonStdBranch(""),
- fBackgroundBranch(""),
- fNonStdFile(""),
- fMomResH1(0x0),
- fMomResH2(0x0),
- fMomResH3(0x0),
- fMomResH1Fit(0x0),
- fMomResH2Fit(0x0),
- fMomResH3Fit(0x0),
- fhEffH1(0x0),
- fhEffH2(0x0),
- fhEffH3(0x0),
- fUseTrMomentumSmearing(kFALSE),
- fUseDiceEfficiency(kFALSE),
- fRparam(1.0),
- fAlgorithm(fastjet::kt_algorithm),
- fStrategy(fastjet::Best),
- fRecombScheme(fastjet::BIpt_scheme),
- fAreaType(fastjet::active_area),
- fGhostArea(0.01),
- fActiveAreaRepeats(1),
- fGhostEtamax(1.5),
- background(0),
- fTCARandomConesOut(0x0),
- fRandom(0),
- fh1Xsec(0x0),
- fh1Trials(0x0),
- fh1PtHard(0x0),
- fh1PtHardNoW(0x0),
- fh1PtHardTrials(0x0),
- fh1Nch(0x0),
- fh1CentralityPhySel(0x0),
- fh1Centrality(0x0),
- fh1DeltapT(0x0),
- fh1Rho(0x0),
- fh1RhoSigma(0x0),
- fh1PtRandCone(0x0),
- fh1efficiencyPt(0x0),
- fh2efficiencyPhi(0x0),
- fh1ZPhySel(0x0),
- fh1Z(0x0),
- fHistList(0x0)
-{
-
-}
-
-AliAnalysisTaskJetHBOM::AliAnalysisTaskJetHBOM(const char* name):
- AliAnalysisTaskSE(name),
- fAOD(0x0),
- fAODExtension(0x0),
- fRef(0),
- fUseAODTrackInput(kFALSE),
- fUseAODMCInput(kFALSE),
- fEventSelection(kFALSE),
- fFilterMask(0),
- fFilterMaskBestPt(0),
- fFilterType(0),
- fJetTypes(kJet),
- fTrackTypeRec(kTrackUndef),
- fTrackTypeGen(kTrackUndef),
- fNSkipLeadingRan(0),
- fNSkipLeadingCone(0),
- fNRandomCones(0),
- randCone_pos(0),
- randCone_Eta(0),
- randCone_Phi(0),
- fNHBOM(0),
- fTrackEtaWindow(0.9),
- fTrackPtCut(0.),
- fJetOutputMinPt(0.150),
- fMaxTrackPtInJet(100.),
- fVtxZCut(8),
- fVtxR2Cut(1),
- fCentCutUp(0),
- fCentCutLo(0),
- fNonStdBranch(""),
- fBackgroundBranch(""),
- fNonStdFile(""),
- fMomResH1(0x0),
- fMomResH2(0x0),
- fMomResH3(0x0),
- fMomResH1Fit(0x0),
- fMomResH2Fit(0x0),
- fMomResH3Fit(0x0),
- fhEffH1(0x0),
- fhEffH2(0x0),
- fhEffH3(0x0),
- fUseTrMomentumSmearing(kFALSE),
- fUseDiceEfficiency(kFALSE),
- fRparam(1.0),
- fAlgorithm(fastjet::kt_algorithm),
- fStrategy(fastjet::Best),
- fRecombScheme(fastjet::BIpt_scheme),
- fAreaType(fastjet::active_area),
- fGhostArea(0.01),
- fActiveAreaRepeats(1),
- fGhostEtamax(1.5),
- background(0),
- fTCARandomConesOut(0x0),
- fRandom(0),
- fh1Xsec(0x0),
- fh1Trials(0x0),
- fh1PtHard(0x0),
- fh1PtHardNoW(0x0),
- fh1PtHardTrials(0x0),
- fh1Nch(0x0),
- fh1CentralityPhySel(0x0),
- fh1Centrality(0x0),
- fh1DeltapT(0x0),
- fh1Rho(0x0),
- fh1RhoSigma(0x0),
- fh1PtRandCone(0x0),
- fh1efficiencyPt(0x0),
- fh2efficiencyPhi(0x0),
- fh1ZPhySel(0x0),
- fh1Z(0x0),
- fHistList(0x0)
-{
- DefineOutput(1, TList::Class());
-}
-
-
-
-Bool_t AliAnalysisTaskJetHBOM::Notify()
-{
- //
- // Implemented Notify() to read the cross sections
- // and number of trials from pyxsec.root
- //
- return kTRUE;
-}
-
-void AliAnalysisTaskJetHBOM::UserCreateOutputObjects()
-{
-
- //
- // Create the output container
- //
-
- fRandom = new TRandom3(0);
-
-
- // Connect the AOD
-
-
- if (fDebug > 1) printf("AnalysisTaskJetHBOM::UserCreateOutputObjects() \n");
-
-
-
- if(fNonStdBranch.Length()!=0)
- {
- // only create the output branch if we have a name
- // Create a new branch for jets...
- // -> cleared in the UserExec....
- // here we can also have the case that the brnaches are written to a separate file
-
- // create the branch for the random cones with the same R
- TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
- if(fUseDiceEfficiency || fUseTrMomentumSmearing)
- cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone);
-
- //create array for the random cones; Until now only one cone per event is used
- if(!AODEvent()->FindListObject(cName.Data())){
- fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
- fTCARandomConesOut->SetName(cName.Data());
- AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
- }
-
- if(fNonStdFile.Length()!=0){
- //
- // case that we have an AOD extension we need to fetch the jets from the extended output
- // we identify the extension aod event by looking for the branchname
- AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
- // case that we have an AOD extension we need can fetch the background maybe from the extended output
- fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
- }
- }
-
- // FitMomentumResolution();
-
-
- if(!fHistList)fHistList = new TList();
- fHistList->SetOwner();
- PostData(1, fHistList); // post data in any case once
-
- Bool_t oldStatus = TH1::AddDirectoryStatus();
- TH1::AddDirectory(kFALSE);
-
- //
- // Histogram
-
- const Int_t nBinPt = 100;
- Double_t binLimitsPt[nBinPt+1];
- for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
- if(iPt == 0){
- binLimitsPt[iPt] = 0.0;
- }
- else {// 1.0
- binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
- }
- }
-
- const Int_t nBinPhi = 90;
- Double_t binLimitsPhi[nBinPhi+1];
- for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
- if(iPhi==0){
- binLimitsPhi[iPhi] = -1.*TMath::Pi();
- }
- else{
- binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
- }
- }
-
-
-
- const Int_t nBinEta = 40;
- Double_t binLimitsEta[nBinEta+1];
- for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
- if(iEta==0){
- binLimitsEta[iEta] = -2.0;
- }
- else{
- binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
- }
- }
-
- const int nChMax = 5000;
-
- fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
- fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
-
- fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
- fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
-
-
- fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
- fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
- fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
-
- fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
-
- fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
- fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
-
- fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
- fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
-
- fh1DeltapT = new TH1F("fh1DeltapT","DeltapT",100,-50,50);
- fh1Rho = new TH1F("fh1Rho","Rho",100,0,200);
- fh1RhoSigma = new TH1F("fh1RhoSigma","SigmaRho",40,0,40);
- fh1PtRandCone = new TH1F("fh1PtRandCone","pt",100,0,200);
-
- const Int_t saveLevel = 3; // large save level more histos
- if(saveLevel>0){
- fHistList->Add(fh1Xsec);
- fHistList->Add(fh1Trials);
-
- fHistList->Add(fh1Nch);
- fHistList->Add(fh1Centrality);
- fHistList->Add(fh1CentralityPhySel);
- fHistList->Add(fh1Z);
- fHistList->Add(fh1ZPhySel);
- fHistList->Add(fh1DeltapT);
- fHistList->Add(fh1Rho);
- fHistList->Add(fh1RhoSigma);
- fHistList->Add(fh1PtRandCone);
- }
-
- // =========== Switch on Sumw2 for all histos ===========
- for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
- TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
- if (h1){
- h1->Sumw2();
- continue;
- }
- THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
- if(hn)hn->Sumw2();
- }
- TH1::AddDirectory(oldStatus);
-}
-
-void AliAnalysisTaskJetHBOM::Init()
-{
- //
- // Initialization
- //
-
- if (fDebug > 1) printf("AnalysisTaskJetHBOM::Init() \n");
-
- FitMomentumResolution();
-
-}
-
-void AliAnalysisTaskJetHBOM::UserExec(Option_t */*option*/)
-{
-
- // handle and reset the output jet branch
-
- if(fTCARandomConesOut)fTCARandomConesOut->Delete();
-
- //
- // Execute analysis for current event
- //
- AliESDEvent *fESD = 0;
- if(fUseAODTrackInput){
- fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
- if(!fAOD){
- Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
- return;
- }
- // fetch the header
- }
- else{
- // assume that the AOD is in the general output...
- fAOD = AODEvent();
- if(!fAOD){
- Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
- return;
- }
- if(fDebug>0){
- fESD = dynamic_cast<AliESDEvent*> (InputEvent());
- }
- }
-
- //Check if information is provided detector level effects
- if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
- if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE;
-
- Bool_t selectEvent = false;
- Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
-
- Float_t cent = 0;
- Float_t zVtx = 0;
-
- if(fAOD){
- const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
- TString vtxTitle(vtxAOD->GetTitle());
- zVtx = vtxAOD->GetZ();
-
- cent = ((AliVAODHeader*)fAOD->GetHeader())->GetCentrality();
- if(physicsSelection){
- fh1CentralityPhySel->Fill(cent);
- fh1ZPhySel->Fill(zVtx);
- }
- // zVertex and centrality selection
- if(fEventSelection){
- if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
- Float_t yvtx = vtxAOD->GetY();
- Float_t xvtx = vtxAOD->GetX();
- Float_t r2 = yvtx*yvtx+xvtx*xvtx;
- if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){//usual fVtxZCut=10 and fVtxR2Cut=1 // apply vertex cut later on
- if(physicsSelection){
- selectEvent = true;
- }
- }
- }
- //centrality selection
- if(fCentCutUp>0){
- if(cent<fCentCutLo||cent>fCentCutUp){
- selectEvent = false;
- }
- }
- }else{
- selectEvent = true;
- }
- }
-
-
- if(!selectEvent){
- PostData(1, fHistList);
- return;
- }
- fh1Centrality->Fill(cent);
- fh1Z->Fill(zVtx);
- fh1Trials->Fill("#sum{ntrials}",1);
-
-
- if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
-
- // ==== General variables needed
-
-
-
- // we simply fetch the tracks/mc particles as a list of AliVParticles
-
- //reconstructed particles
- TList recParticles;
- Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
- Float_t nCh = recParticles.GetEntries();
- fh1Nch->Fill(nCh);
- if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
- //nT = GetListOfTracks(&genParticles,fTrackTypeGen);
- //if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
-
-
- //apply efficiency fNHBOM times
- if(fNHBOM>0){
- for(int particle=0;particle<recParticles.GetEntries();particle++){
- // hier Effizienzen laden und überprüfen ob das Teilchen nachgewiesen wird.
- AliVParticle *vp = (AliVParticle*)recParticles.At(particle);
-
- Double_t pT = vp->Pt();
- Double_t phi = vp->Phi();
-
- //load efficiency
- Double_t efficiencyPt = fh1efficiencyPt->GetBinContent(fh1efficiencyPt->FindBin(pT));
- Double_t efficiencyPhi = fh2efficiencyPhi->GetBinContent(fh2efficiencyPhi->FindBin(phi,pT));
- Double_t eff = efficiencyPt*efficiencyPhi; //over all efficiency
-
- // if ran<eff -> particle is detected eff^fNHBOM = efficiency to detect particle fNHBOM times
- Double_t ran = fRandom->Rndm();
- if(ran>TMath::Power(eff,fNHBOM)){
- recParticles.Remove(vp);
- }
- }
- }
-
-
- // find the jets....
-
- vector<fastjet::PseudoJet> inputParticlesRec;
- vector<fastjet::PseudoJet> inputParticlesRecRan;
-
- //randomize particles
- AliAODJet vTmpRan(1,0,0,1);
- for(int i = 0; i < recParticles.GetEntries(); i++){
- AliVParticle *vp = (AliVParticle*)recParticles.At(i);
-
- // Carefull energy is not well determined in real data, should not matter for p_T scheme?
- // we take total momentum here
-
- //Add particles to fastjet in case we are not running toy model
- fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
- jInp.set_user_index(i);
- inputParticlesRec.push_back(jInp);
-
- // the randomized input changes eta and phi, but keeps the p_T
- Double_t pT = vp->Pt();
- Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
- Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
-
- Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
- Double_t pZ = pT/TMath::Tan(theta);
-
- Double_t pX = pT * TMath::Cos(phi);
- Double_t pY = pT * TMath::Sin(phi);
- Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
- fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
-
- jInpRan.set_user_index(i);
- inputParticlesRecRan.push_back(jInpRan);
- vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
-
- }// recparticles
-
- if(inputParticlesRec.size()==0){
- if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
- PostData(1, fHistList);
- return;
- }
-
- // run fast jet
- // employ setters for these...
-
- fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
- fastjet::AreaType areaType = fastjet::active_area;
- fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
- fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
- fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
-
- //range where to compute background
- Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
- phiMin = 0;
- phiMax = 2*TMath::Pi();
- rapMax = fGhostEtamax - fRparam;
- rapMin = - fGhostEtamax + fRparam;
- fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
-
-
- const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
- const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
-
-
- // loop over all jets and fill information, first one is the leading jet
-
- if(inclusiveJets.size()>0){
-
- //background estimates:all bckg jets wo the 2 hardest
- vector<fastjet::PseudoJet> jets2=sortedJets;
- 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
- Double_t bkg1=0;
- Double_t sigma1=0.;
- Double_t meanarea1=0.;
- clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
- fh1RhoSigma->Fill(sigma1);// fluctuation of the background
- background = bkg1;//sets background variable of the task to the correct value
-
-
- // generate random cones
- if(fTCARandomConesOut){
- // create a random jet within the acceptance
- Double_t etaMax = fTrackEtaWindow - fRparam;//0.9 - 0.4
- Int_t nCone = 0;
- Double_t pTC = 1; // small number
- Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
- Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
- // use fixed position for random Cones
- if(randCone_pos){
- etaC = randCone_Eta;
- phiC = randCone_Phi;
- }
- // massless jet
- Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
- Double_t pZC = pTC/TMath::Tan(thetaC);
- Double_t pXC = pTC * TMath::Cos(phiC);
- Double_t pYC = pTC * TMath::Sin(phiC);
- Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
- AliAODJet tmpRecC (pXC,pYC,pZC, pC);
-
- tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
- if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
-
- // loop over the reconstructed particles and add up the pT in the random cones
- // maybe better to loop over randomized particles not in the real jets...
- // but this by definition brings dow average energy in the whole event
- AliAODJet vTmpRanR(1,0,0,1);
- for(int i = 0; i < recParticles.GetEntries(); i++){
- AliVParticle *vp = (AliVParticle*)recParticles.At(i);
- //add up energy in cone
- if(fTCARandomConesOut){
- AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(0);
- if(jC&&jC->DeltaR(vp)<fRparam){
- if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
- jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
- }
- }// add up energy in cone
-
- }// loop over recparticles
- } //fTCARandomConesOut
- Float_t jetArea = fRparam*fRparam*TMath::Pi();
- if(fTCARandomConesOut){
- // rescale the momentum vectors for the random cones
- AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(0);
- if(rC){
- Double_t etaC = rC->Eta();
- Double_t phiC = rC->Phi();
- // massless jet, unit vector
- Double_t pTC = rC->ChargedBgEnergy();
- if(pTC<=0)pTC = 0.001; // for almost empty events
- Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
- Double_t pZC = pTC/TMath::Tan(thetaC);
- Double_t pXC = pTC * TMath::Cos(phiC);
- Double_t pYC = pTC * TMath::Sin(phiC);
- Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
- rC->SetPxPyPzE(pXC,pYC,pZC, pC);
- rC->SetBgEnergy(0,0);
- rC->SetEffArea(jetArea,0);
- }
- }
- }//inclusive Jets > 0
-
- //Calculate delta pT
- AliAODJet *randCone = (AliAODJet*)fTCARandomConesOut->At(0);
- if(randCone){
- //background is the backbround density per area and area=pi*0.4^2 -> backgroundCone is the background energy under the cone
- Float_t backgroundCone = background * randCone->EffectiveAreaCharged();
- //calculates difference between expected and measured energy density
- Float_t ptSub = randCone->Pt() - backgroundCone;
- fh1DeltapT->Fill(ptSub);// delta pT
- fh1Rho->Fill(background);// background rho
- fh1PtRandCone->Fill(randCone->Pt());// pT of random cone
- }else{
- if(fDebug)Printf("%s:%d No random Cone found",(char*)__FILE__,__LINE__);
- }
-
-
- if (fDebug > 2){
- if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
- }
- PostData(1, fHistList);
-}
-
-void AliAnalysisTaskJetHBOM::Terminate(Option_t */*option*/)
-{
- //
- // Terminate analysis
- //
- if (fDebug > 1) printf("AnalysisJetHBOM: Terminate() \n");
-
- if(fMomResH1Fit) delete fMomResH1Fit;
- if(fMomResH2Fit) delete fMomResH2Fit;
- if(fMomResH3Fit) delete fMomResH3Fit;
-
-}
-
-
-Int_t AliAnalysisTaskJetHBOM::GetListOfTracks(TList *list,Int_t type){
-
- //
- // get list of tracks/particles for different types
- //
-
- if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
-
- Int_t iCount = 0;
- if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
- if(type!=kTrackAODextraonly) {
- AliAODEvent *aod = 0;
- if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
- else aod = AODEvent();
- if(!aod){
- if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
- return iCount;
- }
-
- for(int it = 0;it < aod->GetNumberOfTracks();++it){
- AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
- if(!tr) AliFatal("Not a standard AOD");
- Bool_t bGood = false;
- if(fFilterType == 0)bGood = true;
- else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
- else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
- if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
- if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
- continue;
- }
- if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
- if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
- continue;
- }
- if(tr->Pt()<fTrackPtCut){
- if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
- continue;
- }
- if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
- list->Add(tr);
- iCount++;
- }
- }
- if(type==kTrackAODextra || type==kTrackAODextraonly) {
- AliAODEvent *aod = 0;
- if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
- else aod = AODEvent();
-
- if(!aod){
- return iCount;
- }
- TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
- if(!aodExtraTracks)return iCount;
- for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
- AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
- if (!track) continue;
-
- AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
- if(!trackAOD)continue;
- Bool_t bGood = false;
- if(fFilterType == 0)bGood = true;
- else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
- else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
- if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
- if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
- if(trackAOD->Pt()<fTrackPtCut) continue;
- list->Add(trackAOD);
- iCount++;
- }
- }
- }
- else if (type == kTrackKineAll||type == kTrackKineCharged){
- AliMCEvent* mcEvent = MCEvent();
- if(!mcEvent)return iCount;
- // we want to have alivpartilces so use get track
- for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
- if(!mcEvent->IsPhysicalPrimary(it))continue;
- AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
- if(type == kTrackKineAll){
- if(part->Pt()<fTrackPtCut)continue;
- list->Add(part);
- iCount++;
- }
- else if(type == kTrackKineCharged){
- if(part->Particle()->GetPDG()->Charge()==0)continue;
- if(part->Pt()<fTrackPtCut)continue;
- list->Add(part);
- iCount++;
- }
- }
- }
- else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
- AliAODEvent *aod = 0;
- if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
- else aod = AODEvent();
- if(!aod)return iCount;
- TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
- if(!tca)return iCount;
- for(int it = 0;it < tca->GetEntriesFast();++it){
- AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
- if(!part->IsPhysicalPrimary())continue;
- if(type == kTrackAODMCAll){
- if(part->Pt()<fTrackPtCut)continue;
- list->Add(part);
- iCount++;
- }
- else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
- if(part->Charge()==0)continue;
- if(part->Pt()<fTrackPtCut)continue;
- if(kTrackAODMCCharged){
- list->Add(part);
- }
- else {
- if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
- list->Add(part);
- }
- iCount++;
- }
- }
- }// AODMCparticle
- list->Sort();
- return iCount;
-}
-
-void AliAnalysisTaskJetHBOM::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
-
- //
- // set mom res profiles
- //
-
- fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
- fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
- fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
-}
-
-void AliAnalysisTaskJetHBOM:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
- //
- // set tracking efficiency histos
- //
-
- fhEffH1 = (TH1*)h1->Clone("fhEffH1");
- fhEffH2 = (TH1*)h2->Clone("fhEffH2");
- fhEffH3 = (TH1*)h3->Clone("fhEffH3");
-}
-
-Double_t AliAnalysisTaskJetHBOM::GetMomentumSmearing(Int_t cat, Double_t pt) {
-
- //
- // Get smearing on generated momentum
- //
-
- //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
-
- TProfile *fMomRes = 0x0;
- if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
- if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
- if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
-
- if(!fMomRes) {
- return 0.;
- }
-
-
- Double_t smear = 0.;
-
- if(pt>20.) {
- if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
- if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
- if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
- }
- else {
-
- Int_t bin = fMomRes->FindBin(pt);
-
- smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
-
- }
-
- if(fMomRes) delete fMomRes;
-
- return smear;
-}
-
-void AliAnalysisTaskJetHBOM::FitMomentumResolution() {
- //
- // Fit linear function on momentum resolution at high pT
- //
-
- if(!fMomResH1Fit && fMomResH1) {
- fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
- fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
- fMomResH1Fit ->SetRange(5.,100.);
- }
-
- if(!fMomResH2Fit && fMomResH2) {
- fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
- fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
- fMomResH2Fit ->SetRange(5.,100.);
- }
-
- if(!fMomResH3Fit && fMomResH3) {
- fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
- fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
- fMomResH3Fit ->SetRange(5.,100.);
- }
-
-}
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff