From 7099c89a17b457da231bd7113b101bd3b878f5bb Mon Sep 17 00:00:00 2001 From: esicking Date: Tue, 11 Oct 2011 13:03:18 +0000 Subject: [PATCH] Complete the set of correction steps, available for ESD and AODs. Histos which are not needed anymore were removed. --- PWG4/JetTasks/AliAnalysisTaskMinijet.cxx | 1540 ++++++++++++++-------- PWG4/JetTasks/AliAnalysisTaskMinijet.h | 175 ++- 2 files changed, 1108 insertions(+), 607 deletions(-) diff --git a/PWG4/JetTasks/AliAnalysisTaskMinijet.cxx b/PWG4/JetTasks/AliAnalysisTaskMinijet.cxx index 9d483e73652..4641de6d1bf 100644 --- a/PWG4/JetTasks/AliAnalysisTaskMinijet.cxx +++ b/PWG4/JetTasks/AliAnalysisTaskMinijet.cxx @@ -4,6 +4,8 @@ #include #include #include +#include +#include #include #include #include "TRandom.h" @@ -19,6 +21,7 @@ #include "AliAODEvent.h" #include "AliAODTrack.h" #include "AliAODMCParticle.h" +#include "AliAODMCHeader.h" #include "AliStack.h" #include "AliMCEvent.h" @@ -26,6 +29,12 @@ #include "AliGenEventHeader.h" #include "AliAnalysisManager.h" +#include "AliInputEventHandler.h" + +#include +#include +#include +using namespace std; #include "AliAnalysisTaskMinijet.h" @@ -44,11 +53,13 @@ ClassImp(AliAnalysisTaskMinijet) : AliAnalysisTaskSE(name), fUseMC(kFALSE), fMcOnly(kFALSE), + fBSign(0), + fAnalysePrimOnly(kFALSE),// not used + fPtMin(0.2), + fPtMax(50.0), fCuts(0), - fRadiusCut(0.7), fTriggerPtCut(0.7), fAssociatePtCut(0.4), - fLeadingOrRandom(1), fMode(0), fVertexZCut(10.), fEtaCut(0.9), @@ -58,19 +69,30 @@ ClassImp(AliAnalysisTaskMinijet) fESDEvent(0), fAODEvent(0), fNMcPrimAccept(0), + fNRecAccept(0), fVzEvent(0), + fMeanPtRec(0), + fLeadingPtRec(0), fHists(0), + fStep(0), fHistPt(0), fHistPtMC(0), fNmcNch(0), fPNmcNch(0), fChargedPi0(0) { + //Constructor - for(Int_t i = 0;i< 4;i++){ + for(Int_t i = 0;i< 6;i++){ + fMapSingleTrig[i] = 0; + fMapPair[i] = 0; + fMapEvent[i] = 0; + fMapAll[i] = 0; + fVertexZ[i] = 0; - + + fNcharge[i] = 0; fPt[i] = 0; fEta[i] = 0; fPhi[i] = 0; @@ -91,22 +113,17 @@ ClassImp(AliAnalysisTaskMinijet) fDPhiDEtaEventAxis[i] = 0; fDPhiDEtaEventAxisSeeds[i]= 0; fTriggerNch[i] = 0; + fTriggerNchSeeds[i] = 0; fTriggerTracklet[i] = 0; fNch07Nch[i] = 0; - fPNch07Nch[i] = 0; + fPNch07Nch[i] = 0; + fNch07Tracklet[i] = 0; fNchTracklet[i] = 0; - fPNch07Tracklet[i] = 0; - fDPhiEventAxis[i] = 0; - for(Int_t j=0;j<150;j++){ - fDPhiEventAxisNchBin[i][j] = 0; - fDPhiEventAxisNchBinTrig[i][j] = 0; - - fDPhiEventAxisTrackletBin[i][j] = 0; - fDPhiEventAxisTrackletBinTrig[i][j] = 0; - } + fPNch07Tracklet[i] = 0; + fDPhiEventAxis[i] = 0; } DefineOutput(1, TList::Class()); } @@ -126,29 +143,120 @@ void AliAnalysisTaskMinijet::UserCreateOutputObjects() // Called once if(fDebug) Printf("In User Create Output Objects."); + fStep = new TH1F("fStep", "fStep", 10, -0.5, 9.5); fHistPt = new TH1F("fHistPt", "P_{T} distribution REC", 150, 0.1, 3.1); fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); fHistPt->SetMarkerStyle(kFullCircle); - if (fUseMC) { fHistPtMC = new TH1F("fHistPtMC", "P_{T} distribution MC", 150, 0.1, 3.1); fHistPtMC->GetXaxis()->SetTitle("P_{T} (GeV/c)"); fHistPtMC->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); fHistPtMC->SetMarkerStyle(kFullCircle); - + fNmcNch = new TH2F("fNmcNch", "fNmcNch", 100,-0.5,99.5,100,-0.5,99.5); fPNmcNch = new TProfile("pNmcNch", "pNmcNch", 100,-0.5,99.5); } fChargedPi0 = new TH2F("fChargedPi0", "fChargedPi0", 200, -0.5, 199.5, 200, -0.5, 199.5); + + + //---------------------- + //bins for pt in THnSpare + Double_t ptMin = 0.0, ptMax = 100.; + Int_t nPtBins = 39; + Double_t binsPt[] = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, + 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, + 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0}; + + // Int_t nPtBins = 100; + // Double_t ptMin = 1.e-2, ptMax = 100.; + // Double_t *binsPt = 0; + // binsPt = CreateLogAxis(nPtBins,ptMin,ptMax); + + + Double_t ptMin2 = 0.0, ptMax2 = 100.; + Int_t nPtBins2 = 9; + Double_t binsPt2[] = {0.1, 0.4, 0.7, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 100.0}; + + // Int_t nPtBins2 = 10; + // Double_t ptMin2 = 0.4, ptMax2 = 100.; + // Double_t *binsPt2 = 0; + // binsPt2 = CreateLogAxis(nPtBins2,ptMin2,ptMax2); + + //3 dim matrix + Int_t binsEffHisto[3] = {Int_t(fEtaCut*20), nPtBins, 150 }; + Double_t minEffHisto[3] = {-fEtaCut, ptMin, -0.5 }; + Double_t maxEffHisto[3] = {fEtaCut, ptMax, 149.5 }; + + //5 dim matrix + Int_t binsEffHisto5[6] = { nPtBins2, nPtBins2, Int_t(fEtaCut*10), 180, 150 , 2 }; + Double_t minEffHisto5[6] = { ptMin2, ptMin2, -2*fEtaCut, -0.5*TMath::Pi(), -0.5 , -0.5 }; + Double_t maxEffHisto5[6] = { ptMax2, ptMax2, 2*fEtaCut, 1.5*TMath::Pi(), 149.5 , 1.5 }; + - TString labels[4]={"ESDrec", "ESDmc", "AODrec", "AODmc"}; + //4 dim matrix + Int_t binsEvent[4] = { 150, 20, 50, nPtBins }; + Double_t minEvent[4] = { -0.5, -10, 0, ptMin }; + Double_t maxEvent[4] = { 149.5, 10, 10, ptMax }; + + //3 dim matrix + Int_t binsAll[3] = {Int_t(fEtaCut*20), nPtBins2, 150 }; + Double_t minAll[3] = {-fEtaCut, ptMin2, -0.5 }; + Double_t maxAll[3] = {fEtaCut, ptMax2, 149.5 }; + + //-------------------- + TString dataType[2] ={"ESD", "AOD"}; + TString labels[6]={Form("%sAllAllMcNmc",dataType[fMode].Data()), + Form("%sTrigAllMcNmc",dataType[fMode].Data()), + Form("%sTrigAllMcNrec",dataType[fMode].Data()), + Form("%sTrigVtxMcNrec",dataType[fMode].Data()), + Form("%sTrigVtxRecMcPropNrec",dataType[fMode].Data()), + Form("%sTrigVtxRecNrec",dataType[fMode].Data())}; + + + for(Int_t i=0;i<6;i++){ + + fMapSingleTrig[i] = new THnSparseD(Form("fMapSingleTrig%s", labels[i].Data()),"eta:pt:Nrec", + 3,binsEffHisto,minEffHisto,maxEffHisto); + fMapSingleTrig[i]->SetBinEdges(1,binsPt); + fMapSingleTrig[i]->GetAxis(0)->SetTitle("#eta"); + fMapSingleTrig[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); + fMapSingleTrig[i]->GetAxis(2)->SetTitle("N_{rec}"); + fMapSingleTrig[i]->Sumw2(); + + fMapPair[i] = new THnSparseD(Form("fMapPair%s", labels[i].Data()),"pt_trig:pt_assoc:DeltaEta:DeltaPhi:Nrec:likesign", + 6,binsEffHisto5,minEffHisto5,maxEffHisto5); + fMapPair[i]->SetBinEdges(0,binsPt2); + fMapPair[i]->SetBinEdges(1,binsPt2); + fMapPair[i]->GetAxis(0)->SetTitle("p_{T, trig} (GeV/c)"); + fMapPair[i]->GetAxis(1)->SetTitle("p_{T, assoc} (GeV/c)"); + fMapPair[i]->GetAxis(2)->SetTitle("#Delta #eta"); + fMapPair[i]->GetAxis(3)->SetTitle("#Delta #phi"); + fMapPair[i]->GetAxis(4)->SetTitle("N_{rec}"); + fMapPair[i]->GetAxis(5)->SetTitle("Like-sign or Unlike-sign"); + fMapPair[i]->Sumw2(); + + + fMapEvent[i] = new THnSparseD(Form("fMapEvent%s", labels[i].Data()),"Nrec:vertexZ:meanPt:leadingPt", + 4,binsEvent,minEvent,maxEvent); + fMapEvent[i]->GetAxis(0)->SetTitle("N_{rec}"); + fMapEvent[i]->GetAxis(1)->SetTitle("z_{vertex} (cm)"); + fMapEvent[i]->GetAxis(2)->SetTitle("meanPt"); + fMapEvent[i]->SetBinEdges(3,binsPt); + fMapEvent[i]->GetAxis(3)->SetTitle("leadingPt"); + fMapEvent[i]->Sumw2(); + + fMapAll[i] = new THnSparseD(Form("fMapAll%s", labels[i].Data()),"eta:pt:Nrec", + 3,binsAll,minAll,maxAll); + fMapAll[i]->SetBinEdges(1,binsPt2); + fMapAll[i]->GetAxis(0)->SetTitle("#eta"); + fMapAll[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); + fMapAll[i]->GetAxis(2)->SetTitle("N_{rec}"); + fMapAll[i]->Sumw2(); - // for(Int_t i=2*fMode+fMcOnly;i<1+2*fMode+fUseMC;i++){ - for(Int_t i=0;i<4;i++){ fVertexZ[i] = new TH1F(Form("fVertexZ%s",labels[i].Data()), Form("fVertexZ%s",labels[i].Data()) , @@ -156,6 +264,9 @@ void AliAnalysisTaskMinijet::UserCreateOutputObjects() fPt[i] = new TH1F(Form("fPt%s",labels[i].Data()), Form("fPt%s",labels[i].Data()) , 200, 0., 50); + fNcharge[i] = new TH1F(Form("fNcharge%s",labels[i].Data()), + Form("fNcharge%s",labels[i].Data()) , + 250, -0.5, 249.5); fEta[i] = new TH1F (Form("fEta%s",labels[i].Data()), Form("fEta%s",labels[i].Data()) , 100, -2., 2); @@ -168,36 +279,34 @@ void AliAnalysisTaskMinijet::UserCreateOutputObjects() fDcaZ[i] = new TH1F(Form("fDcaZ%s",labels[i].Data()), Form("fDcaZ%s",labels[i].Data()) , 200, -10,10); - fPtSeed[i] = new TH1F(Form("fPSeedt%s",labels[i].Data()), - Form("fPtSeed%s",labels[i].Data()) , - 500, 0., 50); + Form("fPtSeed%s",labels[i].Data()) , + 500, 0., 50); fEtaSeed[i] = new TH1F (Form("fEtaSeed%s",labels[i].Data()), - Form("fEtaSeed%s",labels[i].Data()) , - 100, -2., 2); + Form("fEtaSeed%s",labels[i].Data()) , + 100, -2., 2); fPhiSeed[i] = new TH1F(Form("fPhiSeed%s",labels[i].Data()), - Form("fPhiSeed%s",labels[i].Data()) , - 360, 0.,2*TMath::Pi()); + Form("fPhiSeed%s",labels[i].Data()) , + 360, 0.,2*TMath::Pi()); fPtOthers[i] = new TH1F(Form("fPOtherst%s",labels[i].Data()), - Form("fPtOthers%s",labels[i].Data()) , - 500, 0., 50); + Form("fPtOthers%s",labels[i].Data()) , + 500, 0., 50); fEtaOthers[i] = new TH1F (Form("fEtaOthers%s",labels[i].Data()), - Form("fEtaOthers%s",labels[i].Data()) , - 100, -2., 2); + Form("fEtaOthers%s",labels[i].Data()) , + 100, -2., 2); fPhiOthers[i] = new TH1F(Form("fPhiOthers%s",labels[i].Data()), - Form("fPhiOthers%s",labels[i].Data()) , - 360, 0.,2*TMath::Pi()); + Form("fPhiOthers%s",labels[i].Data()) , + 360, 0.,2*TMath::Pi()); fPtEtaOthers[i] = new TH2F(Form("fPtEtaOthers%s",labels[i].Data()), Form("fPtEtaOthers%s",labels[i].Data()) , 500, 0., 50, 100, -1., 1); - fPhiEta[i] = new TH2F(Form("fPhiEta%s",labels[i].Data()), Form("fPhiEta%s",labels[i].Data()) , 180, 0., 2*TMath::Pi(), 100, -1.,1.); fDPhiDEtaEventAxis[i] = new TH2F(Form("fDPhiDEtaEventAxis%s",labels[i].Data()), Form("fDPhiDEtaEventAxis%s",labels[i].Data()) , - 180, -0.5* TMath::Pi(), 1.5*TMath::Pi(), 200, -4.,4.); + 180, -0.5* TMath::Pi(), 1.5*TMath::Pi(), 9, -1.8,1.8); fTriggerNch[i] = new TH1F(Form("fTriggerNch%s",labels[i].Data()), Form("fTriggerNch%s",labels[i].Data()) , 250, -0.5, 249.5); @@ -211,8 +320,8 @@ void AliAnalysisTaskMinijet::UserCreateOutputObjects() Form("fNch07Nch%s",labels[i].Data()) , 250, -2.5, 247.5,250, -2.5, 247.5); fPNch07Nch[i] = new TProfile(Form("pNch07Nch%s",labels[i].Data()), - Form("pNch07Nch%s",labels[i].Data()) , - 250, -2.5, 247.5); + Form("pNch07Nch%s",labels[i].Data()) , + 250, -2.5, 247.5); fNch07Tracklet[i] = new TH2F(Form("fNch07Tracklet%s",labels[i].Data()), Form("fNch07Tracklet%s",labels[i].Data()) , 250, -2.5, 247.5,250, -2.5, 247.5); @@ -220,55 +329,37 @@ void AliAnalysisTaskMinijet::UserCreateOutputObjects() Form("fNchTracklet%s",labels[i].Data()) , 250, -2.5, 247.5,250, -2.5, 247.5); fPNch07Tracklet[i] = new TProfile(Form("pNch07Tracklet%s",labels[i].Data()), - Form("pNch07Tracklet%s",labels[i].Data()) , - 250, -2.5, 247.5); + Form("pNch07Tracklet%s",labels[i].Data()) , + 250, -2.5, 247.5); fDPhiEventAxis[i] = new TH1F(Form("fDPhiEventAxis%s", labels[i].Data()), Form("fDPhiEventAxis%s", labels[i].Data()) , 180, -0.5*TMath::Pi(), 1.5*TMath::Pi()); - for(Int_t j=0;j<150;j++){ - fDPhiEventAxisNchBin[i][j] = new TH1F(Form("fDPhiEventAxisNchBin%s%02d", - labels[i].Data(),j), - Form("fDPhiEventAxisNchBin%s%02d", - labels[i].Data(),j) , - 180, -0.5*TMath::Pi(), 1.5*TMath::Pi()); - fDPhiEventAxisNchBinTrig[i][j] = new TH1F(Form("fDPhiEventAxisNchBinTrig%s%02d", - labels[i].Data(),j), - Form("fDPhiEventAxisNchBinTrig%s%02d", - labels[i].Data(),j) , - 180, -0.5*TMath::Pi(), 1.5*TMath::Pi()); - - fDPhiEventAxisTrackletBin[i][j] = new TH1F(Form("fDPhiEventAxisTrackletBin%s%02d", - labels[i].Data(),j), - Form("fDPhiEventAxisTrackletBin%s%02d", - labels[i].Data(),j) , - 180, -0.5*TMath::Pi(), 1.5*TMath::Pi()); - - fDPhiEventAxisTrackletBinTrig[i][j] = new TH1F(Form("fDPhiEventAxisTrackletBinTrig%s%02d", - labels[i].Data(),j), - Form("fDPhiEventAxisTrackletBinTrig%s%02d", - labels[i].Data(),j) , - 180, -0.5*TMath::Pi(), 1.5*TMath::Pi()); - } } fHists = new TList(); fHists->SetOwner(); + fHists->Add(fStep); fHists->Add(fHistPt); + if(fUseMC){ fHists->Add(fHistPtMC); fHists->Add(fNmcNch); fHists->Add(fPNmcNch); } fHists->Add(fChargedPi0); - - //for(Int_t i=2*fMode+fMcOnly;i<1+2*fMode+fUseMC;i++){ - for(Int_t i=0;i<4;i++){ + + for(Int_t i=0;i<6;i++){ + fHists->Add(fMapSingleTrig[i]); + fHists->Add(fMapPair[i]); + fHists->Add(fMapEvent[i]); + fHists->Add(fMapAll[i]); fHists->Add(fVertexZ[i]); fHists->Add(fPt[i]); + fHists->Add(fNcharge[i]); fHists->Add(fEta[i]); fHists->Add(fPhi[i]); fHists->Add(fDcaXY[i]); @@ -291,14 +382,8 @@ void AliAnalysisTaskMinijet::UserCreateOutputObjects() fHists->Add(fNchTracklet[i]); fHists->Add(fPNch07Tracklet[i]); fHists->Add(fDPhiEventAxis[i]); - for(Int_t j=0;j<150;j++){ - fHists->Add(fDPhiEventAxisNchBin[i][j]); - fHists->Add(fDPhiEventAxisNchBinTrig[i][j]); - - fHists->Add(fDPhiEventAxisTrackletBin[i][j]); - fHists->Add(fDPhiEventAxisTrackletBinTrig[i][j]); - } } + PostData(1, fHists); } @@ -306,137 +391,239 @@ void AliAnalysisTaskMinijet::UserCreateOutputObjects() //________________________________________________________________________ void AliAnalysisTaskMinijet::UserExec(Option_t *) { - // Main loop - // Called for each event - - //Printf("Event starts"); + // Main loop, called for each event + // Kinematics-only, ESD and AOD can be processed. + // Data is read (LoopESD, LoopAOD...) and then analysed (Analyse). + // - in case of MC with full detector simulation, all correction steps(0-5) can be processed + // - for Data, only step 5 is performed + // - for kinematics-only, only step 0 is processed + // step 5 = Triggered events, reconstructed accepted vertex, reconstructed tracks, reconstructed multiplicity, + // step 4 = Triggered events, reconstructed accepted vertex, reconstructed tracks with MC properties, reconstructed multiplicity + // step 3 = Triggered events, reconstructed accepted vertex, mc primary particles, reconstructed multiplicity, + // step 2 = Triggered events, all mc primary particles, reconstructed multiplicity + // step 1 = Triggered events, all mc primary particles, true multiplicity + // step 1 = Triggered events, all mc primary particles, true multiplicity + // step 0 = All events, all mc primary particles, true multiplicity + + if(fDebug) Printf("UserExec: Event starts"); AliVEvent *event = InputEvent(); if (!event) { Error("UserExec", "Could not retrieve event"); return; } + fBSign= event->GetMagneticField(); //get events, either ESD or AOD fESDEvent = dynamic_cast (InputEvent()); fAODEvent = dynamic_cast (InputEvent()); + vector pt; + vector eta; + vector phi; + vector charge; + vector px; + vector py; + vector pz; + vector theta; + - //arrays for track properties - Float_t *pt = 0; - Float_t *eta = 0; - Float_t *phi = 0; - Short_t *charge=0; //number of accepted tracks and tracklets Int_t ntracks = 0; //return value for reading functions for ESD and AOD - Int_t *nTracksTracklets = 0; // [0]=nAccepted,1=ntracklets,2=nall(also neutral in case of mc, - //for real nall=ncharged) - - - //read data and analyse. Implemented for ESD, ESDmc, AOD, AODmc - //------------------------------------------------------------- - if (fESDEvent && fMode==0) {//ESDs - //ESD MC reading and analysis - //------ - if (fUseMC){ - ntracks = LoopESDMC(&pt, &eta, &phi, &charge, &nTracksTracklets); //read mc particles - if(ntracks>0){ - Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 1);//analyse + //Int_t ntracksRemove = 0; //return value for reading functions for ESD and AOD + vector nTracksTracklets; // [0]=nAccepted,1=ntracklets,2=nall(also neutral in case of mc, + //for real nall=ncharged) + + + if(!fAODEvent && !fESDEvent)return; + + //=================== AOD =============== + + if(fAODEvent){//AOD loop + + //reset global values + fNMcPrimAccept=0;// number of accepted primaries + fNRecAccept=0; // number of accepted tracks + + if(CheckEvent(true)){//step 5 = TrigVtxRecNrec, step 4 = TrigVtxRecMcPropNrec ,step 3 = TrigVtxMcNrec + + if(!fMcOnly){ + //step 5 = TrigVtxRecNrec + ntracks = LoopAOD(pt, eta, phi, charge, nTracksTracklets, 5);//read tracks + if(pt.size()){ //(internally ntracks=fNRecAccept) + Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 5);//analyse + } + + if(fUseMC){ + // step 4 = TrigVtxRecMcPropNrec + ntracks = LoopAODRecMcProp(pt, eta, phi, charge, nTracksTracklets, 4);//read tracks + if(pt.size()){//(internally ntracks=fNRecAccept) + Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 4);//analyse + } + } } - if(pt || eta || phi || charge || nTracksTracklets) - CleanArrays(pt, eta, phi, charge, nTracksTracklets);// clean up array memory - } - //------ - //ESD reading and analysis - //------ - if(!fMcOnly){ - ntracks = LoopESD(&pt, &eta, &phi,&charge, &nTracksTracklets); //read tracks - if(ntracks>0){ - if(fDebug){ - Printf("ntracks=%d", nTracksTracklets[0]); - Printf("ntracklets=%d", nTracksTracklets[1]); + if(fUseMC){ + // step 3 = TrigVtxMcNrec + ntracks = LoopAODMC(pt, eta, phi, charge, nTracksTracklets, 3);//read tracks + if(pt.size()){//(internally ntracks=fNRecAccept) + Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);//analyse } - Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 0); //analyse tracks } - if(pt || eta || phi || charge || nTracksTracklets) - CleanArrays(pt, eta, phi, charge, nTracksTracklets); // clean up array memory - } - //------ - } - - if (fAODEvent && fMode==1) {//AODs - - //AOD MCreading and analysis - //------ - if (fUseMC){ - ntracks = LoopAODMC(&pt, &eta, &phi, &charge, &nTracksTracklets);//read tracks - if(ntracks>0){ - Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);//analyse + + + }//check event (true) + + //reset values + fNMcPrimAccept=0;// number of accepted primaries + fNRecAccept=0; // number of accepted tracks + + if(CheckEvent(false)){// all events, with and without reconstucted vertex + ntracks = LoopAOD (pt, eta, phi, charge, nTracksTracklets, 2);//need to compute Nrec once more for all events + ntracks = LoopAODMC(pt, eta, phi, charge, nTracksTracklets, 1);//read tracks + if(pt.size()){ + Analyse(pt, eta, phi, charge, fNRecAccept, nTracksTracklets[1],nTracksTracklets[2], 2); // step 2 = TrigAllMcNrec + + Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1); // step 1 = TrigAllMcNmc } - if(pt || eta || phi || charge || nTracksTracklets) - CleanArrays(pt, eta, phi, charge, nTracksTracklets);// clean up array memory + + //step 0 (include not triggered events) is not possible with AODs generated before October 2011 + //step 0 can be implemented for new AODs + } - //------ + }//AOD loop + + + //=================== ESD =============== + + + if(fESDEvent){//ESD loop + + //reset values + fNMcPrimAccept=0;// number of accepted primaries + fNRecAccept=0; // number of accepted tracks + + // instead of task->SelectCollisionCandidate(mask) in AddTask macro + Bool_t isSelected = (((AliInputEventHandler*) + (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())) + ->IsEventSelected() & AliVEvent::kMB); + - //AOD reading and analysis - //------ - if(!fMcOnly){ - ntracks = LoopAOD(&pt, &eta, &phi, &charge, &nTracksTracklets);//read tracks - if(ntracks>0){ - if(fDebug){ - Printf("ntracks=%d", nTracksTracklets[0]); - Printf("ntracklets=%d", nTracksTracklets[1]); + if(fDebug){ + Printf("IsSelected = %d", isSelected); + Printf("CheckEvent(true)= %d", CheckEvent(true)); + Printf("CheckEvent(false)= %d", CheckEvent(false)); + } + + //check trigger + if(isSelected){ // has offline trigger + + if(CheckEvent(true)){//step 5 = TrigVtxRecNrec, step 4 = TrigVtxRecMcPropNrec ,step 3 = TrigVtxMcNrec + + if(!fMcOnly){ + //step 5 = TrigVtxRecNrec + ntracks = LoopESD(pt, eta, phi, charge,nTracksTracklets, 5);//read tracks + if(pt.size()){ //(internally ntracks=fNRecAccept) + Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 5);//analyse + } + + if(fUseMC){ + // step 4 = TrigVtxRecMcPropNrec + ntracks = LoopESDRecMcProp(pt, eta, phi, charge, nTracksTracklets, 4);//read tracks + if(pt.size()){//(internally ntracks=fNRecAccept) + Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 4);//analyse + } + } } - Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 2);//analyse + + if(fUseMC){ + // step 3 = TrigVtxMcNrec + ntracks = LoopESDMC(pt, eta, phi, charge, nTracksTracklets, 3);//read tracks + if(pt.size()){//(internally ntracks=fNRecAccept) + Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);//analyse + } + } + + + }//check event (true) + + //reset values + fNMcPrimAccept=0;// number of accepted primaries + fNRecAccept=0; // number of accepted tracks + + if(CheckEvent(false)){// all events, with and without reconstucted vertex + ntracks = LoopESD (pt, eta, phi, charge, nTracksTracklets, 2);//need to compute Nrec once more for all events + ntracks = LoopESDMC(pt, eta, phi, charge, nTracksTracklets, 1);//read tracks + if(pt.size()){ + Analyse(pt, eta, phi, charge, fNRecAccept, nTracksTracklets[1],nTracksTracklets[2], 2); // step 2 = TrigAllMcNrec + + Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1); // step 1 = TrigAllMcNmc + + Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0); //first part of step 0 // step 0 = AllAllMcNmc + } + + + } + }// triggered event + + else { // not selected by physics selection task = not triggered + if(CheckEvent(false)){ + ntracks = LoopESDMC(pt, eta, phi, charge, nTracksTracklets, 0);//read tracks + if(pt.size())Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0); //second part of step 0 // step 0 = AllAllMcNmc } - if(pt || eta || phi || charge || nTracksTracklets) - CleanArrays(pt, eta, phi, charge, nTracksTracklets);// clean up array memory } - //------ - } - //------------------------------------------------------------- - - - // Post output data. - // PostData(1, fHists); + + + }//ESD loop + + + + } //________________________________________________________________________ -Int_t AliAnalysisTaskMinijet::LoopESD(Float_t **ptArray, Float_t ** etaArray, - Float_t ** phiArray, Short_t ** chargeArray, - Int_t **nTracksTracklets ) +Int_t AliAnalysisTaskMinijet::LoopESD( vector &ptArray, vector &etaArray, + vector &phiArray, vector &chargeArray, + vector &nTracksTracklets, const Int_t step) { // gives back the number of esd tracks and pointer to arrays with track // properties (pt, eta, phi) // Uses TPC tracks with SPD vertex from now on + ptArray.clear(); + etaArray.clear(); + phiArray.clear(); + chargeArray.clear(); + nTracksTracklets.clear(); + + const AliESDVertex* vtxSPD = fESDEvent->GetPrimaryVertexSPD(); // uses track or SPD vertexer + fVertexZ[step]->Fill(vtxSPD->GetZ()); + // Retreive the number of all tracks for this event. Int_t ntracks = fESDEvent->GetNumberOfTracks(); - if(fDebug)Printf("all ESD tracks: %d", ntracks); - - const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD(); + if(fDebug>1) Printf("all ESD tracks: %d", ntracks); //first loop to check how many tracks are accepted //------------------ Int_t nAcceptedTracks=0; for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { - AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks); + AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks); if (!esdTrack) { - Error("UserExec", "Could not receive track %d", iTracks); + Error("LoopESD", "Could not receive track %d", iTracks); continue; } - //if(!fCuts->AcceptTrack(esdTrack)) continue; - - + // use TPC only tracks with non default SPD vertex AliESDtrack *track = AliESDtrackCuts:: GetTPCOnlyTrack(dynamic_cast(fESDEvent),esdTrack->GetID()); if(!track) continue; - if(!fCuts->AcceptTrack(track)) continue;// apply TPC track cuts before constrain to SPD vertex + if(!fCuts->AcceptTrack(track)) { + delete track; + continue;// apply TPC track cuts before constrain to SPD vertex + } if(track->Pt()>0.){ // only constrain tracks above threshold AliExternalTrackParam exParam; @@ -451,59 +638,102 @@ Int_t AliAnalysisTaskMinijet::LoopESD(Float_t **ptArray, Float_t ** etaArray, track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(), exParam.GetCovariance()); } - else continue;// only if tracks have pt<=0 - + else{ + delete track; + continue;// only if tracks have pt<=0 + } - if (TMath::Abs(track->Eta())Pt()>0.2 && track->Pt()<200.) + if (TMath::Abs(track->Eta())Pt()>fPtMin && track->Pt()Pt()); + etaArray.push_back(track->Eta()); + phiArray.push_back(track->Phi()); + chargeArray.push_back(track->Charge()); ++nAcceptedTracks; + fHistPt->Fill(track->Pt()); + } // TPC only track memory needs to be cleaned up if(track) delete track; } + + //need to be checked + if(nAcceptedTracks==0) return 0; + + //tracklet loop + Int_t ntrackletsAccept=0; + AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity()); + Int_t ntracklets = mult->GetNumberOfTracklets(); + for(Int_t i=0;i< ntracklets;i++){ + if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))GetZ(); // needed for correction map + if(step==5 || step ==2) fNRecAccept=nAcceptedTracks; + return fNRecAccept; - //check if event is pile up or no tracks are accepted, return to user exec - if(fESDEvent->IsPileupFromSPD(3,0,8)) return 0; - if(nAcceptedTracks==0) return 0; - //accept event only, if vertex is good and is within fVertexZcut region - const AliESDVertex* vertexESD = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer +} - if(!vertexESD) return 0; - if(vertexESD->GetNContributors()<=0)return 0; - Float_t fVz= vertexESD->GetZ(); - if(TMath::Abs(fVz)>fVertexZCut) return 0; - fVertexZ[0]->Fill(fVz); +//________________________________________________________________________ +Int_t AliAnalysisTaskMinijet::LoopESDRecMcProp( vector &ptArray, vector &etaArray, + vector &phiArray, vector &chargeArray, + vector &nTracksTracklets, const Int_t step) +{ + // gives back the number of esd tracks and pointer to arrays with track + // properties (pt, eta, phi) of mc particles if available + // Uses TPC tracks with SPD vertex from now on - //variables for DCA QA check per track - Float_t fXY = 0.; - Float_t fZ = 0.; + ptArray.clear(); + etaArray.clear(); + phiArray.clear(); + chargeArray.clear(); + nTracksTracklets.clear(); - // Track loop - Int_t iAcceptedTrack=0; + + AliMCEvent *mcEvent = (AliMCEvent*) MCEvent(); + if (!mcEvent) { + Error("LoopESDRecMcProp", "Could not retrieve MC event"); + return 0; + } + AliStack* stack = MCEvent()->Stack(); + if(!stack) return 0; + + + // Retreive the number of all tracks for this event. + Int_t ntracks = fESDEvent->GetNumberOfTracks(); + if(fDebug>1)Printf("all ESD tracks: %d", ntracks); + + const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD(); + fVertexZ[step]->Fill(vtxSPD->GetZ()); + + //track loop + Int_t nAcceptedTracks=0; for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { + + AliVParticle *vtrack = fESDEvent->GetTrack(iTracks); AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks); if (!esdTrack) { - Error("UserExec", "Could not receive track %d", iTracks); + Error("LoopESDRecMcProp", "Could not receive track %d", iTracks); continue; } - //if(!fCuts->AcceptTrack(esdTrack)) continue; - - // create a tpc only track with SPD vertex + + // use TPC only tracks with non default SPD vertex AliESDtrack *track = AliESDtrackCuts:: GetTPCOnlyTrack(dynamic_cast(fESDEvent),esdTrack->GetID()); if(!track) continue; - if(!fCuts->AcceptTrack(track)) continue; // apply TPC track cuts before constrain to SPD vertex - + if(!fCuts->AcceptTrack(track)) { + delete track; + continue;// apply TPC track cuts before constrain to SPD vertex + } if(track->Pt()>0.){ // only constrain tracks above threshold AliExternalTrackParam exParam; @@ -518,74 +748,89 @@ Int_t AliAnalysisTaskMinijet::LoopESD(Float_t **ptArray, Float_t ** etaArray, track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(), exParam.GetCovariance()); } - else continue; - - - if (TMath::Abs(track->Eta())Pt()>0.2 && track->Pt()<200.){ - (*ptArray)[iAcceptedTrack] = track->Pt(); - (*etaArray)[iAcceptedTrack] = track->Eta(); - (*phiArray)[iAcceptedTrack] = track->Phi(); - (*chargeArray)[iAcceptedTrack++] = track->Charge(); - fHistPt->Fill(track->Pt()); - - fXY = 0.; - fZ = 0.; - track->GetImpactParameters(fXY,fZ); - fDcaXY[0]->Fill(fXY); - fDcaZ[0]->Fill(fZ); + else{ + delete track; + continue;// only if tracks have pt<=0 + } + //count tracks, if available, use mc particle properties + if(vtrack->GetLabel()<0){ + if (TMath::Abs(track->Eta())Pt()>fPtMin && track->Pt()Pt()); + etaArray.push_back(track->Eta()); + phiArray.push_back(track->Phi()); + chargeArray.push_back(track->Charge()); + ++nAcceptedTracks; + } } - + else{ + TParticle *partOfTrack = stack->Particle(vtrack->GetLabel()); + if (TMath::Abs(partOfTrack->Eta())Pt()>fPtMin && partOfTrack->Pt()Pt()); + etaArray.push_back(partOfTrack->Eta()); + phiArray.push_back(partOfTrack->Phi()); + chargeArray.push_back(vtrack->Charge()); + ++nAcceptedTracks; + } + } + // TPC only track memory needs to be cleaned up if(track) delete track; } - + if(nAcceptedTracks==0) return 0; //tracklet loop Int_t ntrackletsAccept=0; AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity()); Int_t ntracklets = mult->GetNumberOfTracklets(); for(Int_t i=0;i< ntracklets;i++){ - if(mult->GetDeltaPhi(i)<0.05){ + if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))GenEventHeader(); + TArrayF mcV; + header->PrimaryVertex(mcV); + fVzEvent= mcV[2]; + + return fNRecAccept; // give back reconstructed value - if(fUseMC){ - //Printf("Number of MC particles from ESDMC = %d",fNMcPrimAccept); - //Printf("Number of tracks from ESD = %d",nAcceptedTracks); - fNmcNch->Fill(fNMcPrimAccept,nAcceptedTracks); - fPNmcNch->Fill(fNMcPrimAccept,nAcceptedTracks); - return fNMcPrimAccept; // also possible to use reconstructed Nch -> return nAcceptedTracks; - } - else{ - fVzEvent=fVz; // needed for correction maps - return nAcceptedTracks; - } } + + + //________________________________________________________________________ -Int_t AliAnalysisTaskMinijet::LoopESDMC(Float_t **ptArray, Float_t ** etaArray, - Float_t ** phiArray, Short_t ** chargeArray, - Int_t ** nTracksTracklets) +Int_t AliAnalysisTaskMinijet::LoopESDMC(vector &ptArray, vector &etaArray, + vector &phiArray, vector &chargeArray, + vector &nTracksTracklets, const Int_t step) { // gives back the number of charged prim MC particle and pointer to arrays // with particle properties (pt, eta, phi) - // Printf("Event starts: Loop ESD MC"); + ptArray.clear(); + etaArray.clear(); + phiArray.clear(); + chargeArray.clear(); + nTracksTracklets.clear(); + + fNMcPrimAccept=0; AliMCEvent *mcEvent = (AliMCEvent*) MCEvent(); if (!mcEvent) { - Error("UserExec", "Could not retrieve MC event"); + Error("LoopESDMC", "Could not retrieve MC event"); return 0; } @@ -593,8 +838,14 @@ Int_t AliAnalysisTaskMinijet::LoopESDMC(Float_t **ptArray, Float_t ** etaArray, if(!stack) return 0; Int_t ntracks = mcEvent->GetNumberOfTracks(); - if(fDebug)Printf("MC particles: %d", ntracks); + if(fDebug>1)Printf("MC particles: %d", ntracks); + //vertex + AliGenEventHeader* header = MCEvent()->GenEventHeader(); + TArrayF mcV; + header->PrimaryVertex(mcV); + Float_t vzMC = mcV[2]; + fVertexZ[step]->Fill(vzMC); //---------------------------------- //first track loop to check how many chared primary tracks are available @@ -605,78 +856,52 @@ Int_t AliAnalysisTaskMinijet::LoopESDMC(Float_t **ptArray, Float_t ** etaArray, for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { AliMCParticle *track = dynamic_cast(mcEvent->GetTrack(iTracks)); if (!track) { - Error("UserExec", "Could not receive track %d", iTracks); + Error("LoopESDMC", "Could not receive track %d", iTracks); continue; } if(//count also charged particles in case of fSelectParticles==2 (only neutral) !SelectParticlePlusCharged( - track->Charge(), - track->PdgCode(), - stack->IsPhysicalPrimary(track->Label()) - ) + track->Charge(), + track->PdgCode(), + stack->IsPhysicalPrimary(track->Label()) + ) ) continue; - - // Printf("fSelectParticles= %d\n", fSelectParticles); //count number of pseudo tracklets - if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.035) ++nPseudoTracklets; + if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.0) ++nPseudoTracklets; //0.035 //same cuts as on ESDtracks - if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue; + if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()Pt()>fPtMax) continue; //count all primaries ++nAllPrimaries; //count charged primaries if (track->Charge() != 0) ++nChargedPrimaries; - // Printf("PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) ); + if(fDebug>2) Printf("PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) ); } //---------------------------------- - // Printf("All in acceptance=%d",nAllPrimaries); - // Printf("Charged in acceptance =%d",nChargedPrimaries); + if(fDebug>2){ + Printf("All in acceptance=%d",nAllPrimaries); + Printf("Charged in acceptance =%d",nChargedPrimaries); + } fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries); - if(fSelectParticles!=2){ - *ptArray = new Float_t[nAllPrimaries]; - *etaArray = new Float_t[nAllPrimaries]; - *phiArray = new Float_t[nAllPrimaries]; - *chargeArray = new Short_t[nAllPrimaries]; - } - else{ - *ptArray = new Float_t[nAllPrimaries-nChargedPrimaries]; - *etaArray = new Float_t[nAllPrimaries-nChargedPrimaries]; - *phiArray = new Float_t[nAllPrimaries-nChargedPrimaries]; - *chargeArray = new Short_t[nAllPrimaries-nChargedPrimaries]; - } - - *nTracksTracklets = new Int_t[3]; - if(nAllPrimaries==0) return 0; if(nChargedPrimaries==0) return 0; - - //vertex cut - AliGenEventHeader* header = MCEvent()->GenEventHeader(); - TArrayF mcV; - header->PrimaryVertex(mcV); - if(TMath::Abs(mcV[0])<1e-8 && TMath::Abs(mcV[0])<1e-8 && TMath::Abs(mcV[0])<1e-8) return 0; - Float_t vzMC = mcV[2]; - if(TMath::Abs(vzMC)>fVertexZCut) return 0; - fVertexZ[1]->Fill(vzMC); - - //track loop - Int_t iChargedPiK=0; + //Int_t iChargedPiK=0; for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { AliMCParticle *track = dynamic_cast(mcEvent->GetTrack(iTracks)); if (!track) { - Error("UserExec", "Could not receive track %d", iTracks); + Error("LoopESDMC", "Could not receive track %d", iTracks); continue; } @@ -690,219 +915,310 @@ Int_t AliAnalysisTaskMinijet::LoopESDMC(Float_t **ptArray, Float_t ** etaArray, //same cuts as on ESDtracks - if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 - || track->Pt()>200) continue; + if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()Pt()>fPtMax) continue; - // Printf("After: PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) ); + if(fDebug>2) Printf("After: PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) ); fHistPtMC->Fill(track->Pt()); //fills arrays with track properties - (*ptArray)[iChargedPiK] = track->Pt(); - (*etaArray)[iChargedPiK] = track->Eta(); - (*phiArray)[iChargedPiK] = track->Phi(); - (*chargeArray)[iChargedPiK++] = track->Charge(); - + ptArray.push_back(track->Pt()); + etaArray.push_back(track->Eta()); + phiArray.push_back(track->Phi()); + chargeArray.push_back(track->Charge()); + + } //track loop - (*nTracksTracklets)[0] = nChargedPrimaries; - (*nTracksTracklets)[1] = nPseudoTracklets; + nTracksTracklets.push_back(nChargedPrimaries); + nTracksTracklets.push_back(nPseudoTracklets); if(fSelectParticles!=2){ - (*nTracksTracklets)[2] = nAllPrimaries; + nTracksTracklets.push_back(nAllPrimaries); } else{ - (*nTracksTracklets)[2] = nAllPrimaries-nChargedPrimaries; // only neutral + nTracksTracklets.push_back(nAllPrimaries-nChargedPrimaries); // only neutral } - // Printf("Number of MC particles = %d",nChargedPrimaries); + fNMcPrimAccept=nChargedPrimaries; - return nChargedPrimaries; + + if(step==1){ + fNmcNch->Fill(fNMcPrimAccept,fNRecAccept); + fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept); + } + + fVzEvent= mcV[2]; + return fNRecAccept; } //________________________________________________________________________ -Int_t AliAnalysisTaskMinijet::LoopAOD(Float_t **ptArray, Float_t ** etaArray, - Float_t ** phiArray, Short_t ** chargeArray, - Int_t ** nTracksTracklets) +Int_t AliAnalysisTaskMinijet::LoopAOD( vector &ptArray, vector &etaArray, + vector &phiArray, vector &chargeArray, + vector &nTracksTracklets, const Int_t step) { // gives back the number of AOD tracks and pointer to arrays with track // properties (pt, eta, phi) + ptArray.clear(); + etaArray.clear(); + phiArray.clear(); + chargeArray.clear(); + nTracksTracklets.clear(); + + TClonesArray *mcArray=0x0; + if(fAnalysePrimOnly){ + mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()); + } + + + AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex() + Double_t vzAOD=vertex->GetZ(); + fVertexZ[step]->Fill(vzAOD); // Retreive the number of tracks for this event. Int_t ntracks = fAODEvent->GetNumberOfTracks(); - if(fDebug) Printf("AOD tracks: %d", ntracks); + if(fDebug>1) Printf("AOD tracks: %d", ntracks); Int_t nAcceptedTracks=0; for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks); if (!track) { - Error("UserExec", "Could not receive track %d", iTracks); + Error("LoopAOD", "Could not receive track %d", iTracks); continue; } - // filter bit needs to be ajusted to TPC only tracks with SPD vertex as soon as AOD are available (so far ESD track cuts ITS TPC 2010) - if(track->TestFilterBit(16) && TMath::Abs(track->Eta())Pt()>0.2 && track->Pt()<200.) nAcceptedTracks++; + + AliVParticle *vtrack = fAODEvent->GetTrack(iTracks); + + //use only tracks from primaries + if(fAnalysePrimOnly){ + if(vtrack->GetLabel()<0)continue; + if(!(static_cast(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue; + } + + if(track->TestFilterBit(128) && TMath::Abs(track->Eta())Pt()>fPtMin && track->Pt()DCA()); + //save track properties in vector + ptArray.push_back(track->Pt()); + etaArray.push_back(track->Eta()); + phiArray.push_back(track->Phi()); + chargeArray.push_back(track->Charge()); + fHistPt->Fill(track->Pt()); + + } + } + //need to check this option for MC + if(nAcceptedTracks==0) return 0; + + + //tracklet loop + Int_t ntrackletsAccept=0; + AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets(); + for(Int_t i=0;iGetNumberOfTracklets();++i){ + if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))GetPrimaryVertex(); +} + +//________________________________________________________________________ +Int_t AliAnalysisTaskMinijet::LoopAODRecMcProp( vector &ptArray, vector &etaArray, + vector &phiArray, vector &chargeArray, + vector &nTracksTracklets, const Int_t step) +{ + // gives back the number of AOD tracks and pointer to arrays with track + // properties (pt, eta, phi) + + + ptArray.clear(); + etaArray.clear(); + phiArray.clear(); + chargeArray.clear(); + nTracksTracklets.clear(); - // TODO: need to check how to implement IsPileupFromSPD(3,0.8) - // function of esd event - // first solution: exclude this call in esd loop for comparison (QA) - if(!vertex) return 0; - Double_t vzAOD=vertex->GetZ(); - if(vertex->GetNContributors()<=0) return 0; - if(TMath::Abs(vzAOD)<1e-9) return 0; + // Retreive the number of tracks for this event. + Int_t ntracks = fAODEvent->GetNumberOfTracks(); + if(fDebug>1) Printf("AOD tracks: %d", ntracks); - if(TMath::Abs(vzAOD)>fVertexZCut) return 0; - fVertexZ[2]->Fill(vzAOD); + + //get array of mc particles + TClonesArray *mcArray = (TClonesArray*)fAODEvent-> + FindListObject(AliAODMCParticle::StdBranchName()); + if(!mcArray){ + Printf("No MC particle branch found"); + return kFALSE; + } + + AliAODVertex* vtx= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex() + Double_t vzAOD=vtx->GetZ(); + fVertexZ[step]->Fill(vzAOD); - // Track loop to fill a pT spectrum - Int_t iAcceptedTracks=0; + Int_t nAcceptedTracks=0; for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks); + + AliVParticle *vtrack = fAODEvent->GetTrack(iTracks); + if (!track) { - Error("UserExec", "Could not receive track %d", iTracks); + Error("LoopAODRecMcProp", "Could not receive track %d", iTracks); continue; } - if(!track->TestFilterBit(16) || TMath::Abs(track->Eta())>fEtaCut - || track->Pt()<0.2 || track->Pt()>200.) continue; - fHistPt->Fill(track->Pt()); + + //use only tracks from primaries + if(fAnalysePrimOnly){ + if(vtrack->GetLabel()<0)continue; + if(!(static_cast(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue; + } - //fills arrays with track properties - (*ptArray)[iAcceptedTracks] = track->Pt(); - (*etaArray)[iAcceptedTracks] = track->Eta(); - (*phiArray)[iAcceptedTracks] = track->Phi(); - (*chargeArray)[iAcceptedTracks++] = track->Charge(); + if(track->TestFilterBit(128) && TMath::Abs(track->Eta())Pt()>fPtMin && track->Pt()GetLabel()<0){ //fake tracks + // Printf("Fake track"); + // continue; + ptArray.push_back(track->Pt()); + etaArray.push_back(track->Eta()); + phiArray.push_back(track->Phi()); + chargeArray.push_back(track->Charge()); + } + else{//mc properties + AliAODMCParticle *partOfTrack = (AliAODMCParticle*)mcArray->At(vtrack->GetLabel()); + + ptArray.push_back(partOfTrack->Pt()); + etaArray.push_back(partOfTrack->Eta()); + phiArray.push_back(partOfTrack->Phi()); + chargeArray.push_back(vtrack->Charge());//partOfTrack? + } - } //track loop + } + } + //need to check this option for MC + if(nAcceptedTracks==0) return 0; //tracklet loop Int_t ntrackletsAccept=0; AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets(); for(Int_t i=0;iGetNumberOfTracklets();++i){ - if(TMath::Abs(mult->GetDeltaPhi(i))<0.05){ + if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))Fill(fNMcPrimAccept,nAcceptedTracks); - return fNMcPrimAccept; - } - else - return nAcceptedTracks; -} + //check vertex (mc) + AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent-> + FindListObject(AliAODMCHeader::StdBranchName()); + Float_t vzMC = aodMCheader->GetVtxZ(); + + fVzEvent= vzMC; + return fNRecAccept;//this is the rec value from step 5 + +} + + //________________________________________________________________________ -Int_t AliAnalysisTaskMinijet::LoopAODMC(Float_t **ptArray, Float_t ** etaArray, - Float_t ** phiArray, Short_t ** chargeArray, - Int_t ** nTracksTracklets) +Int_t AliAnalysisTaskMinijet::LoopAODMC( vector &ptArray, vector &etaArray, + vector &phiArray, vector &chargeArray, + vector &nTracksTracklets, const Int_t step) { // gives back the number of AOD MC particles and pointer to arrays with particle // properties (pt, eta, phi) + ptArray.clear(); + etaArray.clear(); + phiArray.clear(); + chargeArray.clear(); + nTracksTracklets.clear(); + + //check vertex + AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent-> + FindListObject(AliAODMCHeader::StdBranchName()); + Float_t vzMC = aodMCheader->GetVtxZ(); + fVertexZ[step]->Fill(vzMC); + + //retreive MC particles from event TClonesArray *mcArray = (TClonesArray*)fAODEvent-> FindListObject(AliAODMCParticle::StdBranchName()); if(!mcArray){ - Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__); + Printf("No MC particle branch found"); return kFALSE; } Int_t ntracks = mcArray->GetEntriesFast(); - if(fDebug)Printf("MC particles: %d", ntracks); + if(fDebug>1)Printf("MC particles: %d", ntracks); // Track loop: chek how many particles will be accepted - Float_t vzMC=0.; + //Float_t vzMC=0.; Int_t nChargedPrim=0; Int_t nAllPrim=0; Int_t nPseudoTracklets=0; for (Int_t it = 0; it < ntracks; it++) { AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it); if (!track) { - Error("UserExec", "Could not receive track %d", it); + Error("LoopAODMC", "Could not receive particle %d", it); continue; } if(!SelectParticlePlusCharged( - track->Charge(), - track->PdgCode(), - track->IsPhysicalPrimary() - ) + track->Charge(), + track->PdgCode(), + track->IsPhysicalPrimary() + ) ) continue; - if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.035)++nPseudoTracklets; - if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue; - // if(TMath::Abs(track->Eta())<1e-9 && TMath::Abs(track->Phi())<1e-9)continue; - - //same cuts as in ESD filter - //if(!track->TestBit(16))continue; //cuts set in ESD filter during AOD generation - + if(TMath::Abs(track->Eta())Pt()>0.0)++nPseudoTracklets; //0.035 + if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()Pt()>fPtMax) continue; + nAllPrim++; if(track->Charge()!=0) nChargedPrim++; - // Printf("eta=%f,phi=%f,pt=%f",track->Eta(),track->Phi(),track->Pt()); - // Printf("prim=%d,pdg=%d,charge=%d",track->IsPhysicalPrimary(),track->PdgCode(),track->Charge()); - - - if(nAllPrim==1) vzMC= track->Zv(); // check only one time. (only one vertex per event allowed) - } - - //generate array with size of number of accepted tracks - fChargedPi0->Fill(nAllPrim,nChargedPrim); - - if(fSelectParticles!=2){ - *ptArray = new Float_t[nAllPrim]; - *etaArray = new Float_t[nAllPrim]; - *phiArray = new Float_t[nAllPrim]; - *chargeArray = new Short_t[nAllPrim]; - } - else{ - *ptArray = new Float_t[nAllPrim-nChargedPrim]; - *etaArray = new Float_t[nAllPrim-nChargedPrim]; - *phiArray = new Float_t[nAllPrim-nChargedPrim]; - *chargeArray = new Short_t[nAllPrim-nChargedPrim]; } - - *nTracksTracklets = new Int_t[3]; - //Printf("nAllPrim=%d", nAllPrim); - //Printf("nChargedPrim=%d", nChargedPrim); - if(nAllPrim==0) return 0; if(nChargedPrim==0) return 0; - - if(TMath::Abs(vzMC)>fVertexZCut) return 0; - fVertexZ[3]->Fill(vzMC); - + //generate array with size of number of accepted tracks + fChargedPi0->Fill(nAllPrim,nChargedPrim); + // Track loop: fill arrays for accepted tracks - Int_t iChargedPrim=0; + // Int_t iChargedPrim=0; for (Int_t it = 0; it < ntracks; it++) { AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it); if (!track) { - Error("UserExec", "Could not receive track %d", it); + Error("LoopAODMC", "Could not receive particle %d", it); continue; } @@ -914,89 +1230,134 @@ Int_t AliAnalysisTaskMinijet::LoopAODMC(Float_t **ptArray, Float_t ** etaArray, ) continue; - - if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue; - //if(TMath::Abs(track->Eta())<1e-8 && TMath::Abs(track->Phi())<1e-8)continue; - - //Printf("eta=%f,phi=%f,pt=%f",track->Eta(),track->Phi(),track->Pt()); - //Printf("prim=%d,pdg=%d,charge=%d",track->IsPhysicalPrimary(),track->PdgCode(),track->Charge()); - - //if(track->TestBit(16))continue; - + if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()Pt()>fPtMax) continue; + fHistPtMC->Fill(track->Pt()); - (*ptArray)[iChargedPrim] = track->Pt(); - (*etaArray)[iChargedPrim] = track->Eta(); - (*phiArray)[iChargedPrim] = track->Phi(); - (*chargeArray)[iChargedPrim++] = track->Charge(); - + ptArray.push_back(track->Pt()); + etaArray.push_back(track->Eta()); + phiArray.push_back(track->Phi()); + chargeArray.push_back(track->Charge()); } - (*nTracksTracklets)[0] = nChargedPrim; - (*nTracksTracklets)[1] = nPseudoTracklets; + nTracksTracklets.push_back(nChargedPrim); + nTracksTracklets.push_back(nPseudoTracklets); if(fSelectParticles!=2){ - (*nTracksTracklets)[2] = nAllPrim; + nTracksTracklets.push_back(nAllPrim); } else{ - (*nTracksTracklets)[2] = nAllPrim-nChargedPrim; // only neutral + nTracksTracklets.push_back(nAllPrim-nChargedPrim); // only neutral } + + - fNMcPrimAccept=nChargedPrim; - return nChargedPrim; - // Printf("ntracks=%d", nChargedPrim); + fVzEvent= vzMC; + fNMcPrimAccept = nChargedPrim; + if(step==1){ // step 1 = Trig All Mc Nmc + fNmcNch->Fill( fNMcPrimAccept,fNRecAccept); + fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept); + } + return fNRecAccept; // rec value from step 5 or step 2 + } //________________________________________________________________________ -void AliAnalysisTaskMinijet::Analyse(const Float_t *pt, const Float_t *eta, const Float_t *phi, - const Short_t *charge, Int_t ntracksCharged, - Int_t ntracklets, const Int_t nAll, Int_t mode) +void AliAnalysisTaskMinijet::Analyse(const vector &pt, + const vector &eta, + const vector &phi, + const vector &charge, + const Int_t ntracksCharged, + const Int_t ntracklets, + const Int_t nAll, + const Int_t step) { // analyse track properties (comming from either ESDs or AODs) in order to compute // mini jet activity (chared tracks) as function of charged multiplicity - // ntracks and ntracklets are already the number of accepted tracks and tracklets + fStep->Fill(step); if(fDebug){ - Printf("In Analysis\n"); - Printf("nAll=%d",nAll); - Printf("nCharged=%d",ntracksCharged); + Printf("Analysis Step=%d", step); + if(fDebug>2){ + Printf("nAll=%d",nAll); + Printf("nCharged=%d",ntracksCharged); + for (Int_t i = 0; i < nAll; i++) { + Printf("pt[%d]=%f",i,pt[i]); + } + } + } + + //calculation of mean pt for all tracks in case of step==0 + if(step==5 || step ==2){ // rec step + Double_t meanPt=0.; + Double_t leadingPt=0.; + for (Int_t i = 0; i < nAll; i++) { + if(pt[i]<0.01)continue; + meanPt+=pt[i]; + if(leadingPtFill(propEvent); + + fNcharge[step]->Fill(ntracksCharged); - Float_t ptEventAxis=0; // pt leading - Float_t etaEventAxis=0; // eta leading - Float_t phiEventAxis=0; // phi leading + Float_t ptEventAxis=0; // pt event axis + Float_t etaEventAxis=0; // eta event axis + Float_t phiEventAxis=0; // phi event axis + Short_t chargeEventAxis=0; // charge event axis Float_t ptOthers = 0; // pt others // for all other tracks around event axis -> see loop Float_t etaOthers = 0; // eta others Float_t phiOthers = 0; // phi others + Short_t chargeOthers = 0; // charge others - Int_t *pindexInnerEta = new Int_t[nAll+1]; // Fix for coverity defect complaining in TMath:Sort in line 1022 - Float_t *ptInnerEta = new Float_t[nAll+1]; // ToDo: Proper fix needs to be investigated + Int_t *pindexInnerEta = new Int_t[nAll]; + Float_t *ptInnerEta = new Float_t[nAll]; + + for (Int_t i = 0; i < nAll; i++) { + + if(pt[i]<0.01)continue; + + //fill single particle correction for first step of pair correction + Double_t propAll[3] = {eta[i],pt[i],ntracksCharged }; + fMapAll[step]->Fill(propAll); + + //filling of simple check plots - fPt[mode] -> Fill( pt[i]); - fEta[mode] -> Fill(eta[i]); - fPhi[mode] -> Fill(phi[i]); - fPhiEta[mode]-> Fill(phi[i], eta[i]); + if(pt[i]<0.7) continue; + fPt[step] -> Fill( pt[i]); + fEta[step] -> Fill(eta[i]); + fPhi[step] -> Fill(phi[i]); + fPhiEta[step]-> Fill(phi[i], eta[i]); pindexInnerEta[i]=0; //set all values to zero //fill new array for eta check ptInnerEta[i]= pt[i]; - } + + // define event axis: leading or random track (with pt>fTriggerPtCut) // --------------------------------------- Int_t highPtTracks=0; Int_t highPtTracksInnerEta=0; + // Double_t highPtTracksInnerEtaCorr=0; Int_t mult09=0; //count high pt tracks and high pt tracks in acceptance for seeds for(Int_t i=0;iFill(ntracksCharged, highPtTracksInnerEta); - fPNch07Nch[mode]->Fill(ntracksCharged, highPtTracksInnerEta); + // plot of multiplicity distributions + fNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta); + fPNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta); + if(ntracklets){ - fNch07Tracklet[mode]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds - fNchTracklet[mode]->Fill(ntracklets, ntracksCharged); - fPNch07Tracklet[mode]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds + fNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds + fNchTracklet[step]->Fill(ntracklets, ntracksCharged); + fPNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds } //analysis can only be performed with event axis, defined by high pt track - + if(highPtTracks>0 && highPtTracksInnerEta>0){ - //Printf("%s:%d",(char*)__FILE__,__LINE__); - //check setter of event axis - //default option: random=1, - //second option:leading=0 - - // Int_t axis=-1; - // if(fLeadingOrRandom==0)axis=0; - // else if (fLeadingOrRandom==1)axis= Int_t((highPtTracksInnerEta)*gRandom->Rndm()); - // else Printf("Wrong settings for event axis."); - - // if(fDebug){ - // Printf("Axis tracks has pT=%f, test=%f", ptInnerEta[pindexInnerEta[axis]], pt[pindexInnerEta[axis]]); - // Printf("Axis tracks has eta=%f", eta[pindexInnerEta[axis]]); - // } - - //--------------------------------------- - //Printf("Number of seeds = %d",highPtTracksInnerEta); - - + // build pairs in two track loops // loop over all possible trigger particles (defined by pt_trig and eta_acceptance) - for(Int_t axis=0;axis1){ // require at least two tracks (leading and prob. accosicates) - + for(Int_t axis=0;(axis=fTriggerPtCut); axis++){ + //EventAxisRandom track properties ptEventAxis = pt [pindexInnerEta[axis]]; etaEventAxis = eta[pindexInnerEta[axis]]; phiEventAxis = phi[pindexInnerEta[axis]]; - fPtSeed[mode] -> Fill( ptEventAxis); - fEtaSeed[mode] -> Fill(etaEventAxis); - fPhiSeed[mode] -> Fill(phiEventAxis); + chargeEventAxis = charge[pindexInnerEta[axis]]; + fPtSeed[step] -> Fill( ptEventAxis); + fEtaSeed[step] -> Fill(etaEventAxis); + fPhiSeed[step] -> Fill(phiEventAxis); - //track loop for event propoerties around event axis with pt>triggerPtCut - //loop only over already accepted tracks except event axis - // if(ptEventAxis>fTriggerPtCut && ptEventAxis < 0.8){ - if(ptEventAxis>fTriggerPtCut){ - for (Int_t iTrack = 0; iTrack < nAll; iTrack++) { - - if(pindexInnerEta[axis]==iTrack)continue; // no double counting + Double_t prop[3] = {etaEventAxis,ptEventAxis,ntracksCharged }; + fMapSingleTrig[step]->Fill(prop); - if(fSelectParticlesAssoc==1){ - if(charge[iTrack]==0)continue; - } - if(fSelectParticlesAssoc==2){ - if(charge[iTrack]!=0)continue; - } - // Printf("Charge=%d", charge[iTrack]); + //associated tracks + for (Int_t iTrack = axis+1; iTrack < nAll; iTrack++) { + if(pt[pindexInnerEta[iTrack]]0.4 && ptOthers<0.5){ // only tracks which fullfill associate pt cut - if(ptOthers>fAssociatePtCut){ // only tracks which fullfill associate pt cut - - //plot only properties of tracks with pt>ptassoc - fPtOthers[mode] -> Fill( ptOthers); - fEtaOthers[mode] -> Fill(etaOthers); - fPhiOthers[mode] -> Fill(phiOthers); - fPtEtaOthers[mode] -> Fill(ptOthers, etaOthers); - - Float_t dPhi = phiOthers-phiEventAxis; - if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi(); - else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi(); - Float_t dEta=etaOthers-etaEventAxis; - fDPhiDEtaEventAxis[mode]->Fill(dPhi, dEta); + ptOthers = pt [pindexInnerEta[iTrack]]; + etaOthers = eta[pindexInnerEta[iTrack]]; + phiOthers = phi[pindexInnerEta[iTrack]]; + chargeOthers = charge[pindexInnerEta[iTrack]]; + + + //plot only properties of tracks with pt>ptassoc + fPtOthers[step] -> Fill( ptOthers); + fEtaOthers[step] -> Fill(etaOthers); + fPhiOthers[step] -> Fill(phiOthers); + fPtEtaOthers[step] -> Fill(ptOthers, etaOthers); - fDPhiEventAxis[mode]->Fill(dPhi); - if(ntracksCharged<150){ // in case of MC data, ntracksCharged represents Nmcprim for ESD and MC loop - fDPhiEventAxisNchBin[mode][ntracksCharged]->Fill(dPhi); - if(ptOthers>fTriggerPtCut && TMath::Abs(etaOthers)Fill(dPhi); - } - - if(ntracklets<150){ - fDPhiEventAxisTrackletBin[mode][ntracklets]->Fill(dPhi); - if(ptOthers>fTriggerPtCut && TMath::Abs(etaOthers)Fill(dPhi); - } - - }//tracks fulfill assoc track cut - - }// end of inner track loop - - - // fill histogram with number of tracks (pt>fAssociatePtCut) around event axis - // how often is there a trigger particle at a certain Nch bin + // if(fDebug>2)Printf("%f, %f", pt[pindexInnerEta[axis]], pt[pindexInnerEta[iTrack]]); - - }//if track pt is at least trigger pt - - }//loop over all highPtInnerEta tracks + Float_t dPhi = phiOthers-phiEventAxis; + if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi(); + else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi(); + Float_t dEta=etaOthers-etaEventAxis; - } //if there are more than 1 track + fDPhiDEtaEventAxis[step]->Fill(dPhi, dEta); + fDPhiEventAxis[step]->Fill(dPhi); + + //check outliers + if(ptEventAxis< 0.4 || ptEventAxis > 100) Printf("particles out of range pt"); + if(ntracksCharged<0 || ntracksCharged>150) Printf("particles out of range ncharge"); + if(TMath::Abs(dEta)>2*fEtaCut) { + Printf("particles out of range dEta"); + Printf("eta1=%f, eta2=%f", etaOthers, etaEventAxis); + Printf("step=%d",step); + } + if(dPhi<-0.5*TMath::Pi() || dPhi>1.5*TMath::Pi()){ + Printf("particles out of range dPhi"); + Printf("phi1=%f, phi2=%f", phiOthers, phiEventAxis); + } + + Bool_t isLikeSign = CheckLikeSign(chargeEventAxis, chargeOthers); + + Double_t prop6[6] = {ptEventAxis,ptOthers,dEta,dPhi,ntracksCharged, isLikeSign }; + fMapPair[step]->Fill(prop6); + + }// end of inner track loop + + } //end of outer track loop + + fTriggerNch[step]->Fill(ntracksCharged,highPtTracksInnerEta); + fTriggerNchSeeds[step]->Fill(ntracksCharged,highPtTracksInnerEta); + fTriggerTracklet[step]->Fill(ntracklets); - fTriggerNch[mode]->Fill(ntracksCharged,highPtTracksInnerEta); - fTriggerNchSeeds[mode]->Fill(ntracksCharged,highPtTracksInnerEta); - fTriggerTracklet[mode]->Fill(ntracklets); }//if there is at least one high pt track @@ -1152,21 +1490,22 @@ void AliAnalysisTaskMinijet::Analyse(const Float_t *pt, const Float_t *eta, cons delete[] ptInnerEta; ptInnerEta=0; } - - PostData(1, fHists); } + + //________________________________________________________________________ void AliAnalysisTaskMinijet::Terminate(Option_t*) { //terminate function is called at the end //can be used to draw histograms etc. + } //________________________________________________________________________ -Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, const Bool_t prim) +const Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(Short_t charge, Int_t pdg, Bool_t prim) { //selection of mc particle //fSelectParticles=0: use charged primaries and pi0 and k0 @@ -1200,7 +1539,7 @@ Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, c } //________________________________________________________________________ -Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim) +const Bool_t AliAnalysisTaskMinijet::SelectParticle(Short_t charge, Int_t pdg, Bool_t prim) { //selection of mc particle //fSelectParticles=0: use charged primaries and pi0 and k0 @@ -1216,8 +1555,7 @@ Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t if(!prim) return false; } - } - + } else if (fSelectParticles==1){ if (charge==0 || !prim){ return false; @@ -1233,34 +1571,166 @@ Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t } //________________________________________________________________________ -void AliAnalysisTaskMinijet::CleanArrays(const Float_t* pt, - const Float_t* eta, - const Float_t* phi, - const Short_t * charge, - const Int_t* nTracksTracklets) +const Bool_t AliAnalysisTaskMinijet::CheckEvent(const Bool_t recVertex) { - //clean up of memory used for arrays of track properties + // This function tests the quality of an event (ESD/AOD) (rec and/or mc part) + // recVertex=false: check if Mc events and stack is there, Nmc>0 + // recVertex=false: " + check if there is a good, reconstructed SPD vertex + // defined by |z|0, no PileUpFromSPD(3,0,8) + + + if(fMode==0){//esd + + //mc + if(fUseMC){ + + //mc event + AliMCEvent *mcEvente = (AliMCEvent*) MCEvent(); + if (!mcEvente) { + Error("CheckEvent", "Could not retrieve MC event"); + return false; + } + + //stack + AliStack* stackg = MCEvent()->Stack(); + if(!stackg) return false; + Int_t ntracksg = mcEvente->GetNumberOfTracks(); + if(ntracksg<0) return false; + + //vertex + //AliGenEventHeader* headerg = MCEvent()->GenEventHeader(); + //TArrayF mcVg; + //headerg->PrimaryVertex(mcVg); + // if(TMath::Abs(mcVg[0])<1e-8 && TMath::Abs(mcVg[0])<1e-8 && + // TMath::Abs(mcVg[0])<1e-8) return false; + // Float_t vzMCg = mcVg[2]; + // if(TMath::Abs(vzMCg)>fVertexZCut) return false; + } + + //rec + if(recVertex==true){ + if(fESDEvent->IsPileupFromSPD(3,0,8))return false; + + //rec vertex + // const AliESDVertex* vertexESDg = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer + // if(!vertexESDg) return false; + // if(vertexESDg->GetNContributors()<=0)return false; + // Float_t fVzg= vertexESDg->GetZ(); + // if(TMath::Abs(fVzg)>fVertexZCut) return false; + + //rec spd vertex + const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD(); + if(!vtxSPD) return false; + if(vtxSPD->GetNContributors()<=0)return false; + Float_t fVzSPD= vtxSPD->GetZ(); + if(TMath::Abs(fVzSPD)>fVertexZCut) return false; + } + return true; - if(pt){ - delete[] pt; - pt=0; - } - if(eta){ - delete[] eta; - eta=0; } - if(phi){ - delete[] phi; - phi=0; + + + else if(fMode==1){ //aod + + if(fUseMC){ + //mc + // AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent-> + // FindListObject(AliAODMCHeader::StdBranchName()); + // Float_t vzMC = aodMCheader->GetVtxZ(); + // if(TMath::Abs(vzMC)>fVertexZCut) return false; + + //retreive MC particles from event + TClonesArray *mcArray = (TClonesArray*)fAODEvent-> + FindListObject(AliAODMCParticle::StdBranchName()); + if(!mcArray){ + Printf("No MC particle branch found"); + return false; + } + } + + //rec + if(recVertex==true){ + if(fAODEvent->IsPileupFromSPD(3,0.8))return false; + + // AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex(); + // if(!vertex) return false; + // if(vertex->GetNContributors()<=0) return false; + // Double_t vzAOD=vertex->GetZ(); + // if(TMath::Abs(vzAOD)<1e-9) return false; + // if(TMath::Abs(vzAOD)>fVertexZCut) return false; + + AliAODVertex* vertexSPD= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD(); + if(!vertexSPD) return false; + if(vertexSPD->GetNContributors()<=0) return false; + Double_t vzSPD=vertexSPD->GetZ(); + //if(TMath::Abs(vzSPD)<1e-9) return false; + if(TMath::Abs(vzSPD)>fVertexZCut) return false; + } + return true; + } - if(charge){ - delete[] charge; - charge=0; + + else { + Printf("Wrong mode!"); + return false; } - if(nTracksTracklets){ - delete[] nTracksTracklets; - nTracksTracklets=0; + +} + +//_____________________________________________________________________________ +const Double_t * AliAnalysisTaskMinijet::CreateLogAxis(const Int_t nbins, + const Double_t xmin, + const Double_t xmax) +{ + // returns pointer to an array with can be used to build a logarithmic axis + // it is user responsibility to delete the array + + Double_t logxmin = TMath::Log10(xmin); + Double_t logxmax = TMath::Log10(xmax); + Double_t binwidth = (logxmax-logxmin)/nbins; + + Double_t *xbins = new Double_t[nbins+1]; + + xbins[0] = xmin; + for (Int_t i=1;i<=nbins;i++) { + xbins[i] = xmin + TMath::Power(10,logxmin+i*binwidth); } + return xbins; } +//_____________________________________________________________________________ +const Bool_t AliAnalysisTaskMinijet::CheckLikeSign(const Short_t chargeEventAxis, + const Short_t chargeOthers) +{ + // compute if charge of two particles/tracks has same sign or different sign + + if(fDebug>2) Printf("Charge1=%d, Charge2=%d",chargeEventAxis,chargeOthers); + + if(chargeEventAxis<0){ + if(chargeOthers<0){ + return true; + } + else if(chargeOthers>0){ + return false; + } + } + + else if(chargeEventAxis>0){ + if(chargeOthers>0){ + return true; + } + else if(chargeOthers<0){ + return false; + } + } + + else{ + Printf("Error: Charge not lower nor higher as zero"); + return false; + } + + Printf("Error: Check values of Charge"); + return false; +} + diff --git a/PWG4/JetTasks/AliAnalysisTaskMinijet.h b/PWG4/JetTasks/AliAnalysisTaskMinijet.h index d24d0d1e694..d92ed37064e 100644 --- a/PWG4/JetTasks/AliAnalysisTaskMinijet.h +++ b/PWG4/JetTasks/AliAnalysisTaskMinijet.h @@ -9,58 +9,82 @@ class TList; class TH1F; class TH2F; +class TH3F; class TProfile; - +class THnSparse; class AliESDtrackCuts; - #include "AliAnalysisTaskSE.h" +#include class AliAnalysisTaskMinijet : public AliAnalysisTaskSE { public: AliAnalysisTaskMinijet(const char *name=""); virtual ~AliAnalysisTaskMinijet(); - virtual void UserCreateOutputObjects(); - virtual void UserExec(Option_t* option); - virtual void Terminate(Option_t *); + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t* option); + virtual void Terminate(Option_t *); + virtual void SetCuts(AliESDtrackCuts* cuts){fCuts = cuts;} - Int_t LoopESD (Float_t **pt, Float_t **eta, Float_t **phi, Short_t **charge, Int_t **nTracksTracklets); - Int_t LoopESDMC(Float_t **pt, Float_t **eta, Float_t **phi, Short_t **charge, Int_t **nTracksTracklets); - Int_t LoopAOD (Float_t **pt, Float_t **eta, Float_t **phi, Short_t **charge, Int_t **nTracksTracklets); - Int_t LoopAODMC(Float_t **pt, Float_t **eta, Float_t **phi, Short_t **charge, Int_t **nTracksTracklets); - void Analyse (const Float_t* pt, const Float_t* eta, const Float_t* phi, const Short_t *charge, Int_t ntacks, Int_t ntacklets=0, const Int_t nAll=0, Int_t mode=0); - void CleanArrays(const Float_t *pt, const Float_t *eta, const Float_t *phi, const Short_t *charge, const Int_t *nTracksTracklets=0); - Bool_t SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, const Bool_t prim); - Bool_t SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim); - - - void UseMC(Bool_t useMC=kTRUE, Bool_t mcOnly=kFALSE) {fUseMC = useMC; fMcOnly=mcOnly;} + void SetUseMC(Bool_t useMC=kTRUE, Bool_t mcOnly=kFALSE) {fUseMC = useMC; fMcOnly=mcOnly;} + void SetAnalyseOnlyPrimaries(Bool_t analysePrimOnly) {fAnalysePrimOnly = analysePrimOnly;} // not used anymore + void SetPtRangeMultiplicity(Float_t ptMin, Float_t ptMax) {fPtMin = ptMin; fPtMax = ptMax; } + void SetTriggerPtCut(Float_t triggerPtCut) {fTriggerPtCut = triggerPtCut;} + void SetAssociatePtCut(Float_t associatePtCut) {fAssociatePtCut = associatePtCut;} + void SetModeEsdAod(Int_t mode) {fMode = mode;} + void SetMaxVertexZ(Float_t vertexZCut) {fVertexZCut = vertexZCut;} + void SetMaxEta(Float_t etaCut) {fEtaCut = etaCut;} + void SetMaxEtaSeed(Float_t etaCutSeed) {fEtaCutSeed = etaCutSeed;} + void SetSelectParticles(Int_t selectParticles) {fSelectParticles = selectParticles;} + void SetSelectParticlesAssoc(Int_t selectParticlesAssoc) {fSelectParticlesAssoc = selectParticlesAssoc;} - virtual void SetCuts(AliESDtrackCuts* cuts) {fCuts = cuts;} - void SetRadiusCut(Float_t radiusCut) {fRadiusCut = radiusCut;} - void SetTriggerPtCut(Float_t triggerPtCut) {fTriggerPtCut = triggerPtCut;} - void SetAssociatePtCut(Float_t associatePtCut) {fAssociatePtCut = associatePtCut;} - void SetEventAxis(Int_t leadingOrRandom) {fLeadingOrRandom = leadingOrRandom;} - void SetMode(Int_t mode) {fMode = mode;} - void SetMaxVertexZ(Float_t vertexZCut) {fVertexZCut = vertexZCut;} - void SetMaxEta(Float_t etaCut) {fEtaCut = etaCut;} - void SetMaxEtaSeed(Float_t etaCutSeed) {fEtaCutSeed = etaCutSeed;} - - void SelectParticles(Int_t selectParticles) {fSelectParticles = selectParticles;} - void SelectParticlesAssoc(Int_t selectParticlesAssoc) {fSelectParticlesAssoc = selectParticlesAssoc;} + private: + Int_t LoopESD (std::vector &pt, std::vector &eta, + std::vector &phi, std::vector &charge, + std::vector &nTracksTracklets, const Int_t step); + Int_t LoopESDRecMcProp(std::vector &pt, std::vector &eta, + std::vector &phi, std::vector &charge, + std::vector &nTracksTracklets, const Int_t step); + Int_t LoopESDMC (std::vector &pt, std::vector &eta, + std::vector &phi, std::vector &charge, + std::vector &nTracksTracklets, const Int_t step); + + Int_t LoopAOD (std::vector &pt, std::vector &eta, + std::vector &phi, std::vector &charge, + std::vector &nTracksTracklets, const Int_t step); + Int_t LoopAODRecMcProp(std::vector &pt, std::vector &eta, + std::vector &phi, std::vector &charge, + std::vector &nTracksTracklets, const Int_t step); + Int_t LoopAODMC (std::vector &pt, std::vector &eta, + std::vector &phi, std::vector &charge, + std::vector &nTracksTracklets, const Int_t step); + + void Analyse (const std::vector &pt, + const std::vector &eta, + const std::vector &phi, + const std::vector &charge, + const Int_t ntacks, const Int_t ntacklets=0, + const Int_t nAll=0, const Int_t step=0); + + const Bool_t SelectParticlePlusCharged(Short_t charge, Int_t pdg, Bool_t prim); + const Bool_t SelectParticle(Short_t charge, Int_t pdg, Bool_t prim); + const Bool_t CheckEvent(const Bool_t recVertex); + const Double_t* CreateLogAxis(const Int_t nbins, const Double_t xmin, const Double_t xmax); + const Bool_t CheckLikeSign(const Short_t chargeEventAxis, const Short_t chargeOthers); - private: Bool_t fUseMC; // flag for Monte Carlo usages Bool_t fMcOnly; // flag defines, if only MC data is used in analysis or also reconstructed data + Double_t fBSign; // magnetic field + Bool_t fAnalysePrimOnly; // flag for analysis of primaries only (also in reconstructed data) + Float_t fPtMin; // set lower limit for pt acceptance for mutliplicity defintion + Float_t fPtMax; // set upper limit for pt acceptance for mutliplicity defintion AliESDtrackCuts* fCuts; // List of cuts for ESDs - Float_t fRadiusCut; // radius cut Float_t fTriggerPtCut; // cut on particle pt used as event axis Float_t fAssociatePtCut; // cut on particle pt used for correlations - Int_t fLeadingOrRandom; // event axis:leading track or random track Int_t fMode; // ESD(=0) of AOD(=1) reading Float_t fVertexZCut; // vertex cut Float_t fEtaCut; // eta acceptance cut @@ -71,51 +95,58 @@ class AliAnalysisTaskMinijet : public AliAnalysisTaskSE { AliESDEvent *fESDEvent; //! esd event AliAODEvent *fAODEvent; //! aod event Int_t fNMcPrimAccept; // global variable for mc multiplucity + Int_t fNRecAccept; // global variable for rec multiplucity Float_t fVzEvent; // global variable for rec vertex position + Double_t fMeanPtRec; // global variable for rec mean pt + Double_t fLeadingPtRec; // global variable for rec mean pt - TList *fHists; // output list - TH1F *fHistPt; // Pt spectrum ESD - TH1F *fHistPtMC; // Pt spectrum MC - TH2F *fNmcNch; // N mc - N ch rec + TList *fHists; // output list + TH1F *fStep; // how many events have passed which correction step + TH1F *fHistPt; // Pt spectrum ESD + TH1F *fHistPtMC; // Pt spectrum MC + + TH2F *fNmcNch; // N mc - N ch rec TProfile *fPNmcNch; // N mc - N ch rec - TH2F *fChargedPi0; // charged versus charged+Pi0 - TH1F * fVertexZ[4]; // z of vertex - - TH1F * fPt[4]; // pt - TH1F * fEta[4]; // et - TH1F * fPhi[4]; // phi - TH1F * fDcaXY[4]; // dca xy direction - TH1F * fDcaZ[4]; // dca z direction - - TH1F * fPtSeed[4]; // pt of seed (event axis) - TH1F * fEtaSeed[4]; // eta of seed - TH1F * fPhiSeed[4]; // phi of seed - - TH1F * fPtOthers[4]; // pt of all other particels used in dEtadPhi - TH1F * fEtaOthers[4]; // eta of all other particels used in dEtadPhi - TH1F * fPhiOthers[4]; // phi of all other particels used in dEtadPhi - TH2F * fPtEtaOthers[4]; // pt-eta of all other particels used in dEtadPhi - - - TH2F * fPhiEta[4]; // eta - phi - TH2F * fDPhiDEtaEventAxis[4]; // correlation dEta-dPhi towards event axis - TH2F * fDPhiDEtaEventAxisSeeds[4]; // correlation dEta-dPhi towards event axis of trigger particles - TH1F * fTriggerNch[4]; // number of triggers with accepted-track number - TH2F * fTriggerNchSeeds[4]; // number of triggers with accepted-track number - TH1F * fTriggerTracklet[4]; // number of triggers with accepted-tracklet number - TH2F * fNch07Nch[4]; // nCharged with pT>fTriggerPtCut vs nCharged - TProfile * fPNch07Nch[4]; // nCharged with pT>fTriggerPtCut vs nCharged - TH2F * fNch07Tracklet[4]; // nCharged with pT>fTriggerPtCut vs nTracklet - TH2F * fNchTracklet[4]; // nCharged vs nTracklet - TProfile * fPNch07Tracklet[4]; // nCharged with pT>fTriggerPtCut vs nTracklet - - TH1F * fDPhiEventAxis[4]; // delta phi of associate tracks to event axis - TH1F * fDPhiEventAxisNchBin[4][150];// delta phi of associate tracks to event axis per Nch bin - TH1F * fDPhiEventAxisNchBinTrig[4][150];// "" for all possoble trigger particles - - TH1F * fDPhiEventAxisTrackletBin[4][150]; // delta phi of associate tracks to event axis per Nch bin - TH1F * fDPhiEventAxisTrackletBinTrig[4][150]; // "" for all possible trigger particles + TH2F *fChargedPi0; // charged versus charged+Pi0 + THnSparse *fMapSingleTrig[6]; //! multi-dim histo for trigger track properties + THnSparse *fMapPair[6]; //! multi-dim histo for pair properties + THnSparse *fMapEvent[6]; //! multi-dim histo for event properties + THnSparse *fMapAll[6]; //! multi-dim histo for properties of all analysed tracks + + TH1F * fVertexZ[6]; // z of vertex + TH1F * fNcharge[6]; // pt + TH1F * fPt[6]; // pt + TH1F * fEta[6]; // eta + TH1F * fPhi[6]; // phi + TH1F * fDcaXY[6]; // dca xy direction + TH1F * fDcaZ[6]; // dca z direction + + TH1F * fPtSeed[6]; // pt of seed (event axis) + TH1F * fEtaSeed[6]; // eta of seed + TH1F * fPhiSeed[6]; // phi of seed + + TH1F * fPtOthers[6]; // pt of all other particels used in dEtadPhi + TH1F * fEtaOthers[6]; // eta of all other particels used in dEtadPhi + TH1F * fPhiOthers[6]; // phi of all other particels used in dEtadPhi + TH2F * fPtEtaOthers[6]; // pt-eta of all other particels used in dEtadPhi + + + TH2F * fPhiEta[6]; // eta - phi + TH2F * fDPhiDEtaEventAxis[6]; // correlation dEta-dPhi towards event axis + TH2F * fDPhiDEtaEventAxisSeeds[6]; // correlation dEta-dPhi towards event axis of trigger particles + TH1F * fTriggerNch[6]; // number of triggers with accepted-track number + TH2F * fTriggerNchSeeds[6]; // number of triggers with accepted-track number + TH1F * fTriggerTracklet[6]; // number of triggers with accepted-tracklet number + TH2F * fNch07Nch[6]; // nCharged with pT>fTriggerPtCut vs nCharged + TProfile * fPNch07Nch[6]; // nCharged with pT>fTriggerPtCut vs nCharged + + TH2F * fNch07Tracklet[6]; // nCharged with pT>fTriggerPtCut vs nTracklet + TH2F * fNchTracklet[6]; // nCharged vs nTracklet + TProfile * fPNch07Tracklet[6]; // nCharged with pT>fTriggerPtCut vs nTracklet + + TH1F * fDPhiEventAxis[6]; // delta phi of associate tracks to event axis + AliAnalysisTaskMinijet(const AliAnalysisTaskMinijet&); // not implemented AliAnalysisTaskMinijet& operator=(const AliAnalysisTaskMinijet&); // not implemented -- 2.43.0