#include "TChain.h" #include "TTree.h" #include "TH1F.h" #include "TH2F.h" #include "TH3D.h" #include "TCanvas.h" #include "TParticle.h" #include "TObjArray.h" #include "AliAnalysisTask.h" #include "AliAnalysisManager.h" #include "AliStack.h" #include "AliMCEvent.h" #include "AliESDEvent.h" #include "AliESDInputHandler.h" #include "AliESDtrackCuts.h" #include "AliCentrality.h" #include "AliLog.h" #include "AliAnalysisTaskEfficiencyBF.h" // --------------------------------------------------------------------- // // Task for calculating the efficiency of the Balance Function // for single particles and pairs // // Authors: Panos Christakoglou, Michael Weber // // --------------------------------------------------------------------- ClassImp(AliAnalysisTaskEfficiencyBF) //________________________________________________________________________ AliAnalysisTaskEfficiencyBF::AliAnalysisTaskEfficiencyBF(const char *name) : AliAnalysisTaskSE(name), fESD(0), fQAList(0), fOutputList(0), fHistEventStats(0), fHistCentrality(0), fHistNMult(0), fHistGeneratedEtaPtPhiPlus(0), fHistFindableEtaPtPhiPlus(0), fHistReconstructedEtaPtPhiPlus(0), fHistSurvivedEtaPtPhiPlus(0), fHistGeneratedEtaPtPhiMinus(0), fHistFindableEtaPtPhiMinus(0), fHistReconstructedEtaPtPhiMinus(0), fHistSurvivedEtaPtPhiMinus(0), fHistGeneratedEtaPtPlusControl(0), fHistFindableEtaPtPlusControl(0), fHistReconstructedEtaPtPlusControl(0), fHistSurvivedEtaPtPlusControl(0), fHistGeneratedEtaPtMinusControl(0), fHistFindableEtaPtMinusControl(0), fHistReconstructedEtaPtMinusControl(0), fHistSurvivedEtaPtMinusControl(0), fHistGeneratedEtaPtPlusPlus(0), fHistFindableEtaPtPlusPlus(0), fHistReconstructedEtaPtPlusPlus(0), fHistSurvivedEtaPtPlusPlus(0), fHistGeneratedEtaPtMinusMinus(0), fHistFindableEtaPtMinusMinus(0), fHistReconstructedEtaPtMinusMinus(0), fHistSurvivedEtaPtMinusMinus(0), fHistGeneratedEtaPtPlusMinus(0), fHistFindableEtaPtPlusMinus(0), fHistReconstructedEtaPtPlusMinus(0), fHistSurvivedEtaPtPlusMinus(0), fHistGeneratedPhiEtaPlusPlus(0), fHistFindablePhiEtaPlusPlus(0), fHistReconstructedPhiEtaPlusPlus(0), fHistSurvivedPhiEtaPlusPlus(0), fHistGeneratedPhiEtaMinusMinus(0), fHistFindablePhiEtaMinusMinus(0), fHistReconstructedPhiEtaMinusMinus(0), fHistSurvivedPhiEtaMinusMinus(0), fHistGeneratedPhiEtaPlusMinus(0), fHistFindablePhiEtaPlusMinus(0), fHistReconstructedPhiEtaPlusMinus(0), fHistSurvivedPhiEtaPlusMinus(0), fESDtrackCuts(0), fAnalysisMode(0), fCentralityEstimator("V0M"), fCentralityPercentileMin(0.0), fCentralityPercentileMax(5.0), fVxMax(3.0), fVyMax(3.0), fVzMax(10.), fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), fMaxDCAxy(3.0), fMaxDCAz(3.0), fMinPt(0.3), fMaxPt(1.5), fMinEta(-0.8), fMaxEta(0.8),fEtaRangeMin(0.0), fEtaRangeMax(1.6), fPtRangeMin(0.3), fPtRangeMax(1.5), fPhiRangeMin(0.0),fPhiRangeMax(360.),fdPhiRangeMax(180.), fEtaBin(100),fdEtaBin(64),fPtBin(49),fPhiBin(100),fdPhiBin(90) { // Define input and output slots here // Input slot #0 works with a TChain DefineInput(0, TChain::Class()); // Output slot #0 id reserved by the base class for AOD // Output slot #1 writes into a TH1 container DefineOutput(1, TList::Class()); DefineOutput(2, TList::Class()); } //________________________________________________________________________ void AliAnalysisTaskEfficiencyBF::UserCreateOutputObjects() { // Create histograms // Called once // global switch disabling the reference // (to avoid "Replacing existing TH1" if several wagons are created in train) Bool_t oldStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); fQAList = new TList(); fQAList->SetName("QAList"); fQAList->SetOwner(); fOutputList = new TList(); fOutputList->SetName("OutputList"); fOutputList->SetOwner(); //Event stats. TString gCutName[4] = {"Total","Offline trigger", "Vertex","Analyzed"}; fHistEventStats = new TH1F("fHistEventStats", "Event statistics;;N_{events}", 4,0.5,4.5); for(Int_t i = 1; i <= 4; i++) fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data()); fQAList->Add(fHistEventStats); //ESD analysis fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events", 20,0.5,20.5); fQAList->Add(fHistCentrality); //multiplicity (good MC tracks) TString histName; histName = "fHistNMult"; fHistNMult = new TH1F(histName.Data(), ";N_{mult.}", 200,0,20000); fQAList->Add(fHistNMult); //eta vs pt for MC positives fHistGeneratedEtaPtPhiPlus = new TH3D("fHistGeneratedEtaPtPhiPlus", "Generated positive primaries;#eta;p_{T} (GeV/c);#phi", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax); fOutputList->Add(fHistGeneratedEtaPtPhiPlus); fHistFindableEtaPtPhiPlus = new TH3D("fHistFindableEtaPtPhiPlus", "Findable positive primaries;#eta;p_{T} (GeV/c);#phi", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax); fOutputList->Add(fHistFindableEtaPtPhiPlus); fHistReconstructedEtaPtPhiPlus = new TH3D("fHistReconstructedEtaPtPhiPlus", "Reconstructed positive primaries;#eta;p_{T} (GeV/c);#phi", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax); fOutputList->Add(fHistReconstructedEtaPtPhiPlus); fHistSurvivedEtaPtPhiPlus = new TH3D("fHistSurvivedEtaPtPhiPlus", "Survived positive primaries;#eta;p_{T} (GeV/c);#phi", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax); fOutputList->Add(fHistSurvivedEtaPtPhiPlus); //eta vs pt for MC negatives fHistGeneratedEtaPtPhiMinus = new TH3D("fHistGeneratedEtaPtPhiMinus", "Generated positive primaries;#eta;p_{T} (GeV/c);#phi", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax); fOutputList->Add(fHistGeneratedEtaPtPhiMinus); fHistFindableEtaPtPhiMinus = new TH3D("fHistFindableEtaPtPhiMinus", "Findable positive primaries;#eta;p_{T} (GeV/c);#phi", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax); fOutputList->Add(fHistFindableEtaPtPhiMinus); fHistReconstructedEtaPtPhiMinus = new TH3D("fHistReconstructedEtaPtPhiMinus", "Reconstructed positive primaries;#eta;p_{T} (GeV/c);#phi", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax); fOutputList->Add(fHistReconstructedEtaPtPhiMinus); fHistSurvivedEtaPtPhiMinus = new TH3D("fHistSurvivedEtaPtPhiMinus", "Survived positive primaries;#eta;p_{T} (GeV/c);#phi", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax); fOutputList->Add(fHistSurvivedEtaPtPhiMinus); //eta vs pt for MC positives (control) fHistGeneratedEtaPtPlusControl = new TH2F("fHistGeneratedEtaPtPlusControl", "Generated positive primaries;#eta;p_{T} (GeV/c)", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistGeneratedEtaPtPlusControl); fHistFindableEtaPtPlusControl = new TH2F("fHistFindableEtaPtPlusControl", "Findable positive primaries;#eta;p_{T} (GeV/c)", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistFindableEtaPtPlusControl); fHistReconstructedEtaPtPlusControl = new TH2F("fHistReconstructedEtaPtPlusControl", "Reconstructed positive primaries;#eta;p_{T} (GeV/c)", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistReconstructedEtaPtPlusControl); fHistSurvivedEtaPtPlusControl = new TH2F("fHistSurvivedEtaPtPlusControl", "Survived positive primaries;#eta;p_{T} (GeV/c)", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistSurvivedEtaPtPlusControl); //eta vs pt for MC negatives (control) fHistGeneratedEtaPtMinusControl = new TH2F("fHistGeneratedEtaPtMinusControl", "Generated positive primaries;#eta;p_{T} (GeV/c)", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistGeneratedEtaPtMinusControl); fHistFindableEtaPtMinusControl = new TH2F("fHistFindableEtaPtMinusControl", "Findable positive primaries;#eta;p_{T} (GeV/c)", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistFindableEtaPtMinusControl); fHistReconstructedEtaPtMinusControl = new TH2F("fHistReconstructedEtaPtMinusControl", "Reconstructed positive primaries;#eta;p_{T} (GeV/c)", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistReconstructedEtaPtMinusControl); fHistSurvivedEtaPtMinusControl = new TH2F("fHistSurvivedEtaPtMinusControl", "Survived positive primaries;#eta;p_{T} (GeV/c)", fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistSurvivedEtaPtMinusControl); //eta vs pt for MC ++ fHistGeneratedEtaPtPlusPlus = new TH2F("fHistGeneratedEtaPtPlusPlus", "Generated ++ primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistGeneratedEtaPtPlusPlus); fHistFindableEtaPtPlusPlus = new TH2F("fHistFindableEtaPtPlusPlus", "Findable ++ primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistFindableEtaPtPlusPlus); fHistReconstructedEtaPtPlusPlus = new TH2F("fHistReconstructedEtaPtPlusPlus", "Reconstructed ++ primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistReconstructedEtaPtPlusPlus); fHistSurvivedEtaPtPlusPlus = new TH2F("fHistSurvivedEtaPtPlusPlus", "Survived ++ primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistSurvivedEtaPtPlusPlus); //eta vs pt for MC -- fHistGeneratedEtaPtMinusMinus = new TH2F("fHistGeneratedEtaPtMinusMinus", "Generated -- primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistGeneratedEtaPtMinusMinus); fHistFindableEtaPtMinusMinus = new TH2F("fHistFindableEtaPtMinusMinus", "Findable -- primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistFindableEtaPtMinusMinus); fHistReconstructedEtaPtMinusMinus = new TH2F("fHistReconstructedEtaPtMinusMinus", "Reconstructed -- primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistReconstructedEtaPtMinusMinus); fHistSurvivedEtaPtMinusMinus = new TH2F("fHistSurvivedEtaPtMinusMinus", "Survived -- primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistSurvivedEtaPtMinusMinus); //eta vs pt for MC +- fHistGeneratedEtaPtPlusMinus = new TH2F("fHistGeneratedEtaPtPlusMinus", "Generated +- primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistGeneratedEtaPtPlusMinus); fHistFindableEtaPtPlusMinus = new TH2F("fHistFindableEtaPtPlusMinus", "Findable +- primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistFindableEtaPtPlusMinus); fHistReconstructedEtaPtPlusMinus = new TH2F("fHistReconstructedEtaPtPlusMinus", "Reconstructed +- primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistReconstructedEtaPtPlusMinus); fHistSurvivedEtaPtPlusMinus = new TH2F("fHistSurvivedEtaPtPlusMinus", "Survived +- primaries;#Delta#eta;p_{T} (GeV/c)", fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax); fOutputList->Add(fHistSurvivedEtaPtPlusMinus); //=============================// //phi vs eta for MC ++ fHistGeneratedPhiEtaPlusPlus = new TH2F("fHistGeneratedPhiEtaPlusPlus", "Generated ++ primaries;#Delta#phi", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistGeneratedPhiEtaPlusPlus); fHistFindablePhiEtaPlusPlus = new TH2F("fHistFindablePhiEtaPlusPlus", "Findable ++ primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistFindablePhiEtaPlusPlus); fHistReconstructedPhiEtaPlusPlus = new TH2F("fHistReconstructedPhiEtaPlusPlus", "Reconstructed ++ primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistReconstructedPhiEtaPlusPlus); fHistSurvivedPhiEtaPlusPlus = new TH2F("fHistSurvivedPhiEtaPlusPlus", "Survived ++ primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistSurvivedPhiEtaPlusPlus); //phi vs eta for MC -- fHistGeneratedPhiEtaMinusMinus = new TH2F("fHistGeneratedPhiEtaMinusMinus", "Generated -- primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistGeneratedPhiEtaMinusMinus); fHistFindablePhiEtaMinusMinus = new TH2F("fHistFindablePhiEtaMinusMinus", "Findable -- primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistFindablePhiEtaMinusMinus); fHistReconstructedPhiEtaMinusMinus = new TH2F("fHistReconstructedPhiEtaMinusMinus", "Reconstructed -- primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistReconstructedPhiEtaMinusMinus); fHistSurvivedPhiEtaMinusMinus = new TH2F("fHistSurvivedPhiEtaMinusMinus", "Survived -- primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistSurvivedPhiEtaMinusMinus); //phi vs eta for MC +- fHistGeneratedPhiEtaPlusMinus = new TH2F("fHistGeneratedPhiEtaPlusMinus", "Generated +- primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistGeneratedPhiEtaPlusMinus); fHistFindablePhiEtaPlusMinus = new TH2F("fHistFindablePhiEtaPlusMinus", "Findable +- primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistFindablePhiEtaPlusMinus); fHistReconstructedPhiEtaPlusMinus = new TH2F("fHistReconstructedPhiEtaPlusMinus", "Reconstructed +- primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistReconstructedPhiEtaPlusMinus); fHistSurvivedPhiEtaPlusMinus = new TH2F("fHistSurvivedPhiEtaPlusMinus", "Survived +- primaries;#Delta#phi;#Delta#eta", fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax); fOutputList->Add(fHistSurvivedPhiEtaPlusMinus); //=============================// fQAList->Print(); fOutputList->Print(); PostData(1, fQAList); PostData(2, fOutputList); TH1::AddDirectory(oldStatus); } //________________________________________________________________________ void AliAnalysisTaskEfficiencyBF::UserExec(Option_t *) { // Main loop // Called for each event // Post output data. //ESD analysis fESD = dynamic_cast(InputEvent()); if (!fESD) { printf("ERROR: fESD not available\n"); return; } AliMCEvent* mcEvent = MCEvent(); if (!mcEvent) { AliError("ERROR: Could not retrieve MC event"); return; } AliStack* stack = mcEvent->Stack(); if (!stack) { AliError("ERROR: Could not retrieve MC stack"); return; } // arrays for 2 particle histograms Int_t nMCLabelCounter = 0; const Int_t maxMCLabelCounter = 20000; Double_t eta[maxMCLabelCounter]; Double_t pt[maxMCLabelCounter]; Double_t phi[maxMCLabelCounter]; Int_t level[maxMCLabelCounter]; Int_t charge[maxMCLabelCounter]; //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fESD->GetNumberOfTracks())); fHistEventStats->Fill(1); //all events //Centrality stuff AliCentrality *centrality = fESD->GetCentrality(); Int_t nCentrality = 0; nCentrality = (Int_t)(centrality->GetCentralityPercentile(fCentralityEstimator.Data())/10.); //Printf("Centrality: %lf",centrality->GetCentralityPercentile(fCentralityEstimator.Data())); if(centrality->IsEventInCentralityClass(fCentralityPercentileMin, fCentralityPercentileMax, fCentralityEstimator.Data())) { fHistEventStats->Fill(2); //triggered + centrality fHistCentrality->Fill(nCentrality+1); //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax); if(fAnalysisMode.CompareTo("TPC") == 0 ) { const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC(); if(vertex) { if(vertex->GetNContributors() > 0) { if(vertex->GetZRes() != 0) { fHistEventStats->Fill(3); //events with a proper vertex if(TMath::Abs(vertex->GetX()) < fVxMax) { if(TMath::Abs(vertex->GetY()) < fVyMax) { if(TMath::Abs(vertex->GetZ()) < fVzMax) { fHistEventStats->Fill(4); //analyzed events Int_t nMCParticles = mcEvent->GetNumberOfTracks(); TArrayI labelMCArray(nMCParticles); for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) { AliMCParticle *mcTrack = (AliMCParticle*) mcEvent->GetTrack(iTracks); if (!mcTrack) { AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks)); continue; } //exclude particles generated out of the acceptance Double_t vz = mcTrack->Zv(); if (TMath::Abs(vz) > 50.) continue; //acceptance if(TMath::Abs(mcTrack->Eta()) > fMaxEta) continue; if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt)) continue; if((mcTrack->Phi() > fPhiRangeMax)||(mcTrack->Phi() < fPhiRangeMin)) continue; TParticle* particle = mcTrack->Particle(); if(!particle) continue; if(!stack->IsPhysicalPrimary(iTracks)) continue; if(iTracks <= stack->GetNprimary()) { Short_t gMCCharge = mcTrack->Charge(); Double_t phiRad = particle->Phi(); Double_t phiDeg = phiRad*TMath::RadToDeg(); if(gMCCharge > 0) fHistGeneratedEtaPtPhiPlus->Fill(particle->Eta(), particle->Pt(), phiDeg); else if(gMCCharge < 0) fHistGeneratedEtaPtPhiMinus->Fill(particle->Eta(), particle->Pt(), phiDeg); // findable tracks --> DOES NOT WORK???? // Loop over Track References Bool_t labelTPC = kTRUE; AliTrackReference* trackRef = 0; for (Int_t iTrackRef = 0; iTrackRef < mcTrack->GetNumberOfTrackReferences(); iTrackRef++) { trackRef = mcTrack->GetTrackReference(iTrackRef); if(trackRef) { Int_t detectorId = trackRef->DetectorId(); if (detectorId == AliTrackReference::kTPC) { labelTPC = kTRUE; break; } } }//loop over track references if(labelTPC) { labelMCArray.AddAt(iTracks,nMCLabelCounter); if(nMCLabelCounter >= maxMCLabelCounter){ AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter)); break; } //fill the arrays for 2 particle analysis eta[nMCLabelCounter] = particle->Eta(); pt[nMCLabelCounter] = particle->Pt(); phi[nMCLabelCounter] = particle->Phi()*TMath::RadToDeg(); charge[nMCLabelCounter] = gMCCharge; // findable = generated in this case! level[nMCLabelCounter] = 1; nMCLabelCounter += 1; if(gMCCharge > 0) fHistFindableEtaPtPhiPlus->Fill(particle->Eta(), particle->Pt(), phiDeg); else if(gMCCharge < 0) fHistFindableEtaPtPhiMinus->Fill(particle->Eta(), particle->Pt(), phiDeg); } }//primaries }//loop over MC particles fHistNMult->Fill(nMCLabelCounter); // not used so far //Float_t dcaXY = 0.0, dcaZ = 0.0; //ESD track loop Int_t nGoodTracks = fESD->GetNumberOfTracks(); TArrayI labelArray(nGoodTracks); Int_t labelCounter = 0; for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { AliESDtrack* track = fESD->GetTrack(iTracks); //AliESDtrack* track = fESDtrackCuts->GetTPCOnlyTrack(fESD,iTracks); if(!track) continue; AliESDtrack *tpcOnlyTrack = new AliESDtrack(); if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) { delete tpcOnlyTrack; continue; } Int_t label = TMath::Abs(track->GetTPCLabel()); if(IsLabelUsed(labelArray,label)) continue; labelArray.AddAt(label,labelCounter); labelCounter += 1; Int_t mcGoods = nMCLabelCounter; for (Int_t k = 0; k < mcGoods; k++) { Int_t mcLabel = labelMCArray.At(k); if (mcLabel != TMath::Abs(label)) continue; if(mcLabel != label) continue; if(label > stack->GetNtrack()) continue; TParticle *particle = stack->Particle(label); if(!particle) continue; //acceptance if(TMath::Abs(particle->Eta()) > fMaxEta) continue; if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt)) continue; if((particle->Phi() > fPhiRangeMax)||(particle->Phi() < fPhiRangeMin)) continue; if(!stack->IsPhysicalPrimary(label)) continue; if(label <= stack->GetNprimary()) { // reconstructed level[k] = 2; Short_t gCharge = track->Charge(); Double_t phiRad = particle->Phi(); Double_t phiDeg = phiRad*TMath::RadToDeg(); if(gCharge > 0) fHistReconstructedEtaPtPhiPlus->Fill(particle->Eta(), particle->Pt(), phiDeg); else if(gCharge < 0) fHistReconstructedEtaPtPhiMinus->Fill(particle->Eta(), particle->Pt(), phiDeg); // track cuts + analysis kinematic cuts if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta()) < fMaxEta && track->Pt() > fMinPt && track->Pt() < fMaxPt ){ // survived level[k] = 3; if(gCharge > 0) fHistSurvivedEtaPtPhiPlus->Fill(particle->Eta(), particle->Pt(), phiDeg); else if(gCharge < 0) fHistSurvivedEtaPtPhiMinus->Fill(particle->Eta(), particle->Pt(), phiDeg); }//track cuts }//primary particles }//findable track loop }//ESD track loop labelMCArray.Reset(); labelArray.Reset(); }//Vz cut }//Vy cut }//Vx cut }//Vz resolution }//number of contributors }//valid vertex }//TPC analysis mode }//centrality // Here comes the 2 particle analysis // loop over all good MC particles for (Int_t i = 0; i < nMCLabelCounter ; i++) { // control 1D histograms (charge might be different?) if(charge[i] > 0){ if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]); if(level[i] > 1) fHistReconstructedEtaPtPlusControl->Fill(eta[i],pt[i]); if(level[i] > 2) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]); } else if(charge[i] < 0){ if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]); if(level[i] > 1) fHistReconstructedEtaPtMinusControl->Fill(eta[i],pt[i]); if(level[i] > 2) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]); } for (Int_t j = i+1; j < nMCLabelCounter ; j++) { if(charge[i] > 0 && charge[j] > 0 ){ if(level[i] > 0 && level[j] > 0) { fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]); if (TMath::Abs(phi[i]-phi[j]) < 180) fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); } if(level[i] > 1 && level[j] > 1) { fHistReconstructedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]); if (TMath::Abs(phi[i]-phi[j]) < 180) fHistReconstructedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); } if(level[i] > 2 && level[j] > 2) { fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]); if (TMath::Abs(phi[i]-phi[j]) < 180) fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); } } else if(charge[i] < 0 && charge[j] < 0 ){ if(level[i] > 0 && level[j] > 0) { fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]); if (TMath::Abs(phi[i]-phi[j]) < 180) fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); } if(level[i] > 1 && level[j] > 1) { fHistReconstructedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]); fHistReconstructedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); } if(level[i] > 2 && level[j] > 2) { fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]); if (TMath::Abs(phi[i]-phi[j]) < 180) fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); } } else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){ if(level[i] > 0 && level[j] > 0) { fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]); if (TMath::Abs(phi[i]-phi[j]) < 180) fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); } if(level[i] > 1 && level[j] > 1) { fHistReconstructedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]); fHistReconstructedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); } if(level[i] > 2 && level[j] > 2) { fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]); if (TMath::Abs(phi[i]-phi[j]) < 180) fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j])); } } } } //Checking Entries for generated and survived /*TH1F* GeneratedEtaPt = (TH1F*)fHistGeneratedEtaPtPlusMinus->ProjectionX("GeneratedEtaPt",5,15); Int_t xGeneratedPt = fHistGeneratedEtaPtPlusMinus->GetNbinsX(); for (Int_t h=1; h < xGeneratedPt+1; h++){ Double_t binEntriesGenerated = GeneratedEtaPt->GetBinContent(h); Printf("binEntriesGenerated: %lf - xGeneratedPt: %d",binEntriesGenerated,h); } TH1F* GeneratedPhiEta = (TH1F*)fHistGeneratedPhiEtaPlusMinus->ProjectionY("GeneratedPhiEta",5,15); Int_t yGeneratedPhi = fHistGeneratedPhiEtaPlusMinus->GetNbinsY(); for (Int_t h=1; h < yGeneratedPhi+1; h++){ Double_t binEntriesGenerated = GeneratedPhiEta->GetBinContent(h); Printf("binEntriesGenerated: %lf - yGeneratedPhi: %d",binEntriesGenerated,h); }*/ /*TH1F* SurvivedEtaPt = (TH1F*)fHistSurvivedEtaPtPlusMinus->ProjectionX("SurvivedEtaPt",5,15); Int_t xSurvivedPt = fHistSurvivedEtaPtPlusMinus->GetNbinsX(); for (Int_t h=1; h < xSurvivedPt+1; h++){ Double_t binEntriesSurvived = SurvivedEtaPt->GetBinContent(h); Printf("binEntriesSurvived: %lf - xSurvivedPt: %d",binEntriesSurvived,h); } TH1F* SurvivedPhiEta = (TH1F*)fHistSurvivedPhiEtaPlusMinus->ProjectionY("SurvivedPhiEta",5,15); Int_t ySurvivedPhi = fHistSurvivedPhiEtaPlusMinus->GetNbinsY(); for (Int_t h=1; h < ySurvivedPhi+1; h++){ Double_t binEntriesSurvived = SurvivedPhiEta->GetBinContent(h); Printf("binEntriesSurvived: %lf - ySurvivedPhi: %d",binEntriesSurvived,h); }*/ } //________________________________________________________________________ void AliAnalysisTaskEfficiencyBF::Terminate(Option_t *) { // Draw result to the screen // Called once at the end of the query fOutputList = dynamic_cast (GetOutputData(1)); if (!fOutputList) { printf("ERROR: Output list not available\n"); return; } } //____________________________________________________________________// Bool_t AliAnalysisTaskEfficiencyBF::IsLabelUsed(TArrayI labelArray, Int_t label) { //Checks if the label is used already Bool_t status = kFALSE; for(Int_t i = 0; i < labelArray.GetSize(); i++) { if(labelArray.At(i) == label) status = kTRUE; } return status; }