//-*- Mode: C++ -*- #include "TMath.h" #include "TAxis.h" #include "TSystem.h" #include "TProfile.h" #include "TH2F.h" #include "TH3F.h" #include "TFile.h" #include "TPRegexp.h" #include "AliStack.h" #include "AliMCEvent.h" #include "AliMCParticle.h" #include "AliESDtrackCuts.h" #include "AliESDInputHandler.h" #include "AliESDpid.h" #include "AliCentrality.h" #include "AliTracker.h" #include "AliAODInputHandler.h" #include "AliAODEvent.h" #include "AliAODTrack.h" #include "AliAODMCParticle.h" #include "AliAnalysisNetParticleQA.h" using namespace std; /** * Class for for NetParticle QA * -- Create input for QA * Authors: Jochen Thaeder * Michael Weber */ ClassImp(AliAnalysisNetParticleQA) /* * --------------------------------------------------------------------------------- * Constructor / Destructor * --------------------------------------------------------------------------------- */ //________________________________________________________________________ AliAnalysisNetParticleQA::AliAnalysisNetParticleQA() : fHelper(NULL), fPdgCode(2212), fESD(NULL), fESDTrackCuts(NULL), fPIDResponse(NULL), fAOD(NULL), fArrayMC(NULL), fAODtrackCutBit(1024), fIsMC(kFALSE), fMCEvent(NULL), fStack(NULL), fCentralityBin(-1.), fNTracks(0), fHnQA(NULL) { // Constructor AliLog::SetClassDebugLevel("AliAnalysisNetParticleQA",10); } //________________________________________________________________________ AliAnalysisNetParticleQA::~AliAnalysisNetParticleQA() { // Destructor return; } /* * --------------------------------------------------------------------------------- * Public Methods * --------------------------------------------------------------------------------- */ //________________________________________________________________________ Int_t AliAnalysisNetParticleQA::Initialize(AliAnalysisNetParticleHelper* helper, AliESDtrackCuts* cuts, Bool_t isMC, Int_t trackCutBit) { // -- Initialize // -- ESD track cuts // ------------------- fESDTrackCuts = cuts; // -- Get Helper // --------------- fHelper = helper; // -- Is MC // ---------- fIsMC = isMC; // -- AOD track filter bit // ------------------------- fAODtrackCutBit = trackCutBit; // -- Get particle species / pdgCode // ------------------------- if (fIsMC) fPdgCode = AliPID::ParticleCode(fHelper->GetParticleSpecies()); // -- Create THnSparse Histograms // -------------------------------- CreateHistograms(); return 0; } //________________________________________________________________________ void AliAnalysisNetParticleQA::CreateHistograms() { // -- Add histograms to outlist // ------------------------------------------------------------------ // -- Create THnSparseF - QA // ------------------------------------------------------------------ Int_t binHnQA[15] = {AliAnalysisNetParticleHelper::fgkfHistNBinsCent, AliAnalysisNetParticleHelper::fgkfHistNBinsEta, // cent | eta AliAnalysisNetParticleHelper::fgkfHistNBinsRap, AliAnalysisNetParticleHelper::fgkfHistNBinsPhi, // y | phi AliAnalysisNetParticleHelper::fgkfHistNBinsPt, AliAnalysisNetParticleHelper::fgkfHistNBinsPt, // pt | pInner AliAnalysisNetParticleHelper::fgkfHistNBinsSign, 500, 50, 50, // sign | TPCsignal | nSigmaTPC | nSigmaTOF 50, 50, 50, 50, // DCAr | DCAz | nSigmaDCAr | nSigmaDCAz 3}; // MCisProbe Double_t minHnQA[15] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[0], AliAnalysisNetParticleHelper::fgkfHistRangeEta[0], AliAnalysisNetParticleHelper::fgkfHistRangeRap[0], AliAnalysisNetParticleHelper::fgkfHistRangePhi[0], AliAnalysisNetParticleHelper::fgkfHistRangePt[0], AliAnalysisNetParticleHelper::fgkfHistRangePt[0], AliAnalysisNetParticleHelper::fgkfHistRangeSign[0], 30, -5.0, -5.0, -10., -10., -10., -10., -0.5}; Double_t maxHnQA[15] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[1], AliAnalysisNetParticleHelper::fgkfHistRangeEta[1], AliAnalysisNetParticleHelper::fgkfHistRangeRap[1], AliAnalysisNetParticleHelper::fgkfHistRangePhi[1], AliAnalysisNetParticleHelper::fgkfHistRangePt[1], AliAnalysisNetParticleHelper::fgkfHistRangePt[1], AliAnalysisNetParticleHelper::fgkfHistRangeSign[1], 500, 5.0, 5.0, 10., 10., 10., 10., 2.5}; fHnQA = new THnSparseF("fHnQA", "cent:eta:y:phi:pt:pInner:sign:TPCsignal:nSigmaTPC:nSigmaTOF:DCAr:DCAz:nSigmaDCAr:nSigmaDCAz:MCisProbe", 15, binHnQA, minHnQA, maxHnQA); fHnQA->Sumw2(); fHnQA->GetAxis(0)->SetTitle("centrality"); // 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins fHnQA->GetAxis(1)->SetTitle("#eta"); // eta [-0.9, 0.9] fHnQA->GetAxis(2)->SetTitle("#it{y}"); // rapidity [-0.5, 0.5] fHnQA->GetAxis(3)->SetTitle("#varphi"); // phi [ 0. , 2Pi] fHnQA->GetAxis(4)->SetTitle("#it{p}_{T} (GeV/#it{c})"); // pt [ 0.1, 3.0] fHnQA->GetAxis(5)->SetTitle("#it{p}_{Inner} (GeV/#it{c})"); // pInner [ 0.1, 3.0] fHnQA->GetAxis(6)->SetTitle("sign"); // -1 | 0 | +1 fHnQA->GetAxis(7)->SetTitle("TPC signal"); // TPCsignal [30, 500] fHnQA->GetAxis(8)->SetTitle("n #sigma TPC"); // nSigma TPC [-5, 5] fHnQA->GetAxis(9)->SetTitle("n #sigma TOF"); // nSigma TOF [-5, 5] fHnQA->GetAxis(10)->SetTitle("DCAr"); // DCAr [-10, 10] fHnQA->GetAxis(11)->SetTitle("DCAz"); // DCAz [-10, 10] fHnQA->GetAxis(12)->SetTitle("n #sigma #sqrt(Cdd)/DCAr"); // nSigma DCAr [-10, 10] fHnQA->GetAxis(13)->SetTitle("n #sigma #sqrt(Czz)/DCAz"); // nSigma DCAz [-10, 10] fHnQA->GetAxis(14)->SetTitle("MCisProbe"); // 0 | 1 (isProbeParticle) | 2 (isProbeParticle wrong sign) fHelper->BinLogAxis(fHnQA, 4); fHelper->BinLogAxis(fHnQA, 5); return; } //________________________________________________________________________ Int_t AliAnalysisNetParticleQA::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) { // -- Setup Event ResetEvent(); // -- Get ESD objects if(esdHandler){ fPIDResponse = esdHandler->GetPIDResponse(); // NEEDED ? fESD = esdHandler->GetEvent(); fNTracks = fESD->GetNumberOfTracks(); } // -- Get AOD objects else if(aodHandler){ fPIDResponse = aodHandler->GetPIDResponse(); fAOD = aodHandler->GetEvent(); fNTracks = fAOD->GetNumberOfTracks(); if (fIsMC) { fArrayMC = dynamic_cast(fAOD->FindListObject(AliAODMCParticle::StdBranchName())); if (!fArrayMC) AliFatal("No array of MC particles found !!!"); // MW no AliFatal use return values } } // -- Get MC objects if (fIsMC && mcEvent) { fMCEvent = mcEvent; if (fMCEvent) fStack = fMCEvent->Stack(); } // -- Get CentralityBin fCentralityBin = fHelper->GetCentralityBin(); return 0; } //________________________________________________________________________ void AliAnalysisNetParticleQA::ResetEvent() { // -- Reset event // -- Reset ESD Event fESD = NULL; // -- Reset AOD Event fAOD = NULL; // -- Reset MC Event if (fIsMC) fMCEvent = NULL; } //________________________________________________________________________ Int_t AliAnalysisNetParticleQA::Process() { // -- Process NetParticle QA // -- Fill ESD/AOC tracks ProcessTracks(); return 0; } //________________________________________________________________________ Int_t AliAnalysisNetParticleQA::ProcessTracks() { // -- Process ESD/AOD tracks and fill QA histogram // -- Get ranges for AOD particles Float_t etaRange[2]; fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]); Float_t ptRange[2]; fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]); // -- Track Loop for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) { AliVTrack *track = (fESD) ? static_cast(fESD->GetTrack(idxTrack)) : static_cast(fAOD->GetTrack(idxTrack)); // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- // -- Check track cuts // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- // -- Check if track is accepted for basic parameters if (!fHelper->IsTrackAcceptedBasicCharged(track)) continue; // -- Check if accepted - ESD if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast(track))) continue; // -- Check if accepted - AOD if (fAOD){ AliAODTrack * trackAOD = dynamic_cast(track); if (!trackAOD) { AliError("Pointer to dynamic_cast(track) = ZERO"); continue; } if (!trackAOD->TestFilterBit(fAODtrackCutBit)) continue; // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs) if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1])) continue; } // -- Check if accepted in rapidity window -- for identified particles // for charged particles -- eta check is done in the trackcuts Double_t yP; if (fHelper->GetUsePID() && !fHelper->IsTrackAcceptedRapidity(track, yP)) continue; // -- Check if accepted with thighter DCA cuts // -- returns kTRUE for AODs for now : MW if (!fHelper->IsTrackAcceptedDCA(track)) continue; // -- Check if accepted by PID from TPC or TPC+TOF Double_t pid[2]; if (!fHelper->IsTrackAcceptedPID(track, pid)) continue; // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- // -- Fill Track // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- // -- Check if probe particle Int_t isProbeParticle = 0; if (fIsMC) { Int_t label = TMath::Abs(track->GetLabel()); AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast(fArrayMC->At(label)); if (particle) { if (TMath::Abs(particle->PdgCode()) == fPdgCode) { ++isProbeParticle; if (particle->PdgCode() != (track->Charge()*fPdgCode)) ++isProbeParticle; } } } // -- Get dca r/z Float_t dca[] = {0.,0.}; // dca_xy, dca_z, Float_t cov[] = {0.,0.,0.}; // sigma_xy, sigma_xy_z, sigma_z if (fESD) (static_cast(track))->GetImpactParameters(dca,cov); Float_t dcaRoverCdd = ( TMath::Sqrt(cov[0]) != 0. ) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99; Float_t dcaZoverCzz = ( TMath::Sqrt(cov[2]) != 0. ) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99; if (!fHelper->GetUsePID()) yP = track->Eta(); Double_t aTrack[15] = { Double_t(fCentralityBin), // 0 centrality track->Eta(), // 1 eta yP, // 2 rapidity track->Phi(), // 3 phi track->Pt(), // 4 pt track->GetTPCmomentum(), // 5 pInner track->Charge(), // 6 sign track->GetTPCsignal(), // 7 TPC dE/dx pid[0], // 8 n Sigma TPC pid[1], // 9 n Sigma TOF dca[0], // 10 dca r dca[1], // 11 dca z dcaRoverCdd, // 12 sqrt(cov[dd]) dcaZoverCzz, // 13 sqrt(cov[zz]) isProbeParticle // 14 isProbe }; fHnQA->Fill(aTrack); } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) { return 0; }