1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: Satyajit Jena. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //=========================================================================//
18 // AliEbyE Analysis for Particle Ratio Fluctuation //
19 // Deepika Rathee | Satyajit Jena //
20 // drathee@cern.ch | sjena@cern.ch //
22 // V0.1 2013/03/25 Using THnSparse
23 // V0.2 2013/04/03 Cleanup
24 // V1.0 2013/04/10 Cleanup Bins for Less Memory
25 // V1.1 2013/04/15 Bins Added
26 // Todo: pp and pA, Mix Events
27 //=========================================================================//
37 #include "AliAnalysisTask.h"
38 #include "AliAnalysisManager.h"
39 #include "AliVEvent.h"
41 #include "AliESDEvent.h"
42 #include "AliAODEvent.h"
43 #include "AliAODMCParticle.h"
44 #include "AliAODMCHeader.h"
45 #include "AliPIDResponse.h"
46 #include "AliAODHeader.h"
47 #include "AliAODpidUtil.h"
48 #include "AliHelperPID.h"
51 #include "AliEbyEParticleRatioFluctuationTask.h"
53 ClassImp(AliEbyEParticleRatioFluctuationTask)
55 //-----------------------------------------------------------------------
56 AliEbyEParticleRatioFluctuationTask::AliEbyEParticleRatioFluctuationTask( const char *name ) : AliAnalysisTaskSE( name ),
59 fAnalysisData("PbPb"),
60 fCentralityEstimator("V0M"),
76 fHistoCorrelation(NULL) {
77 for(Int_t i = 0; i < 14; i++ ) fHistQA[i] = NULL;
78 DefineOutput(1, TList::Class()); //! Connect Outpput....
81 AliEbyEParticleRatioFluctuationTask::~AliEbyEParticleRatioFluctuationTask() {
83 if (fThnList) delete fThnList;
84 if (fHelperPID) delete fHelperPID;
87 //---------------------------------------------------------------------------------
88 void AliEbyEParticleRatioFluctuationTask::UserCreateOutputObjects() {
89 fThnList = new TList();
90 fThnList->SetOwner(kTRUE);
92 fEventCounter = new TH1D("fEventCounter","EventCounter", 300, 0.5,300.5);
93 if (isQA) fThnList->Add(fEventCounter);
95 fHistQA[0] = new TH2F("fHistQAvx", "Histo Vx Selected;Centrality;Vx", 100,0,100, 5000, -5., 5.);
96 fHistQA[1] = new TH2F("fHistQAvy", "Histo Vy Selected;Centrality;Vy", 100,0,100, 5000, -5., 5.);
97 fHistQA[2] = new TH2F("fHistQAvz", "Histo Vz Selected;Centrality;Vz", 100,0,100, 5000, -25., 25.);
98 fHistQA[3] = new TH2F("fHistQAvxA", "Histo Vx;Centrality;Vx", 100,0,100, 5000, -5., 5.);
99 fHistQA[4] = new TH2F("fHistQAvyA", "Histo Vy;Centrality;Vy", 100,0,100, 5000, -5., 5.);
100 fHistQA[5] = new TH2F("fHistQAvzA", "Histo Vz;Centrality;Vz", 100,0,100, 5000, -25., 25.);
102 fHistQA[6] = new TH2F("fHistQADcaXyA", "Histo DCAxy;Centrality;DCAxy",100,0,100, 600, -15., 15.);
103 fHistQA[7] = new TH2F("fHistQADcaZA", "Histo DCAz;Centrality;DCAz ",100,0,100, 600, -15., 15.);
104 fHistQA[8] = new TH2F("fHistQAPtA","p_{T} distribution;Centrality;p_{T}",100,0,100,1000,0,10);
105 fHistQA[9] = new TH2F("fHistQAEtaA","#eta distribution;Centrality;#eta",100,0,100,240,-1.2,1.2);
107 fHistQA[10] = new TH2F("fHistQADcaXy", "Histo DCAxy after Selected;Centrality;DCAxy", 100,0,100,600, -15., 15.);
108 fHistQA[11] = new TH2F("fHistQADcaZ", "Histo DCAz Selected;Centrality;DCAz", 100,0,100,600, -15., 15.);
109 fHistQA[12] = new TH2F("fHistQAPt","p_{T} distribution Selected;Centrality;p_{T}",100,0,100,1000,0,10);
110 fHistQA[13] = new TH2F("fHistQAEta","#eta distribution Selected;Centrality;#eta",100,0,100, 240,-1.2,1.2);
112 if (isQA) for(Int_t i = 0; i < 14; i++) fThnList->Add(fHistQA[i]);
114 Int_t fgSparseDataBins[kNSparseData] = {100, 5000, 5000, 2500, 2500, 3000, 1500, 1500, 1000, 500, 500, 500, 250, 250};
115 Double_t fgSparseDataMin[kNSparseData] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
116 Double_t fgSparseDataMax[kNSparseData] = {100.,5000.,5000.,2500.,2500.,3000.,1500.,1500.,1000.,500.,500.,500.,250.,250.};
118 const Char_t *fgkSparseDataTitle[] = {"centrality","RefMult","N_{ch}", "N_{+}","N_{-}","N_{#pi}", "N_{#pi^{+}}","N_{#pi^{-}}","N_{K}","N_{K^{+}}", "N_{K^{-}}","N_{p}","N_{p}","N_{#bar{p}}"};
120 fHistoCorrelation = new THnSparseI("fThnCorr", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax);
121 for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++)
122 fHistoCorrelation->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]);
124 if(!isQA) fThnList->Add(fHistoCorrelation);
127 fThnList->Add(fHelperPID);
129 PostData(1, fThnList);
132 //----------------------------------------------------------------------------------
133 void AliEbyEParticleRatioFluctuationTask::UserExec( Option_t * ){
135 if (isQA) fEventCounter->Fill(1);
137 AliAODEvent* event = dynamic_cast<AliAODEvent*>(InputEvent());
139 Printf("ERROR 01: AOD not found ");
144 Float_t gRefMul = -1;
146 AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(event->GetHeader());
147 if(!aodHeader) AliFatal("Not a standard AOD");
148 gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()); // V0M`
149 gRefMul = aodHeader->GetRefMultiplicity();
150 if (gCent < 0 || gCent > 100) return;
151 if (isQA) fEventCounter->Fill(2);
153 if (!AcceptEvent(event,gCent)) return;
159 for (Int_t i = 0; i < 2; i++) {
161 for (Int_t ii = 0; ii < 3; ii++) {
167 if(fAnalysisType == "AOD") {
169 fEventCounter->Fill(5);
170 fEventCounter->Fill(50+gCent);
173 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
174 AliAODTrack* track = dynamic_cast<AliAODTrack *>(event->GetTrack(itrk));
175 if (!track) continue;
177 if (!AcceptTrack(track, gCent)) continue;
179 Int_t a = fHelperPID->GetParticleSpecies((AliVTrack*) track,kTRUE);
181 if(a < 0 || a > 2) continue;
182 icharge = track->Charge() > 0 ? 0 : 1;
187 else if(fAnalysisType == "MCAOD") {
188 TClonesArray *arrayMC= 0;
189 arrayMC = dynamic_cast<TClonesArray*> (event->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
191 Printf("Error: MC particles branch not found!\n");
194 AliAODMCHeader *mcHdr=0;
195 mcHdr=(AliAODMCHeader*)event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
197 Printf("MC header branch not found!\n");
202 fEventCounter->Fill(5);
203 fEventCounter->Fill(50+gCent);
206 Int_t nMC = arrayMC->GetEntries();
207 for (Int_t iMC = 0; iMC < nMC; iMC++) {
208 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
209 if(!AcceptMCTrack(partMC, gCent)) continue;
210 Int_t a = fHelperPID->GetMCParticleSpecie((AliVEvent*) event,(AliVTrack*)partMC,1);
211 if(a < 0 || a > 2) continue;
212 icharge = partMC->Charge() > 0 ? 0 : 1;
218 printf(" No Event Type is Defined ");
222 if( (gCharge[0] + gCharge[1]) != 0 ) {
224 fEventCounter->Fill(6);
225 fEventCounter->Fill(160 + gCent);
228 Double_t vsparse[kNSparseData];
230 vsparse[1] = gRefMul;
231 vsparse[2] = gCharge[0] + gCharge[1];
232 vsparse[3] = gCharge[0];
233 vsparse[4] = gCharge[1];
234 vsparse[5] = gPid[0][0] + gPid[0][1];
235 vsparse[6] = gPid[0][0];
236 vsparse[7] = gPid[0][1];
237 vsparse[8] = gPid[1][0] + gPid[1][0];
238 vsparse[9] = gPid[1][0];
239 vsparse[10] = gPid[1][1];
240 vsparse[11] = gPid[2][0] + gPid[2][1];
241 vsparse[12] = gPid[2][0];
242 vsparse[13] = gPid[2][1];
243 fHistoCorrelation->Fill(vsparse);
247 if(fDebug && isQA) Printf(" %6d %6d %6d %6d %6d %6d %6d %6d %6d %6d %6d %6d %6d",
248 (Int_t)fEventCounter->GetBinContent(1),
249 (Int_t)fEventCounter->GetBinContent(2),
250 (Int_t)fEventCounter->GetBinContent(3),
251 (Int_t)gCent, (Int_t)gRefMul,
252 gCharge[0], gCharge[1],
253 gPid[0][0], gPid[0][1], gPid[1][0],
254 gPid[1][1], gPid[2][0], gPid[2][1]);
256 PostData(1, fThnList);
259 void AliEbyEParticleRatioFluctuationTask::Terminate( Option_t * ){
260 Info("AliEbyEParticleRatioFluctuationTask"," Task Successfully finished");
263 //___________________________________________________________
264 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptEvent(AliAODEvent *event, Int_t cent) const {
266 const AliAODVertex *vertex = event->GetPrimaryVertex();
269 vertex->GetCovarianceMatrix(fCov);
270 if(vertex->GetNContributors() > 0) {
274 fEventCounter->Fill(3);
275 fHistQA[3]->Fill(cent,vertex->GetX());
276 fHistQA[4]->Fill(cent,vertex->GetY());
277 fHistQA[5]->Fill(cent,vertex->GetZ());
280 if(TMath::Abs(vertex->GetX()) < fVxMax) {
281 if(TMath::Abs(vertex->GetY()) < fVyMax) {
282 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
285 fEventCounter->Fill(4);
286 fHistQA[0]->Fill(cent,vertex->GetX());
287 fHistQA[1]->Fill(cent,vertex->GetY());
288 fHistQA[2]->Fill(cent,vertex->GetZ());
297 AliCentrality *centrality = event->GetCentrality();
298 if (centrality->GetQuality() != 0) ver = kFALSE;
303 //___________________________________________________________
304 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptTrack(AliAODTrack *track, Int_t cent) const {
305 if(!track) return kFALSE;
306 if (track->Charge() == 0 ) return kFALSE;
309 fHistQA[6]->Fill(cent,track->DCA());
310 fHistQA[7]->Fill(cent,track->ZAtDCA());
311 fHistQA[8]->Fill(cent,track->Pt());
312 fHistQA[9]->Fill(cent,track->Eta());
315 if (!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
317 if (track->Eta() < fEtaLowerLimit ||
318 track->Eta() > fEtaHigherLimit) return kFALSE;
319 if (track->Pt() < fPtLowerLimit ||
320 track->Pt() > fPtHigherLimit) return kFALSE;
321 if ( track->DCA() > fDCAxy ) return kFALSE;
322 if ( track->ZAtDCA() > fDCAz ) return kFALSE;
325 fHistQA[10]->Fill(cent,track->DCA());
326 fHistQA[11]->Fill(cent,track->ZAtDCA());
327 fHistQA[12]->Fill(cent,track->Pt());
328 fHistQA[13]->Fill(cent,track->Eta());
335 //___________________________________________________________
336 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptMCTrack(AliAODMCParticle *track, Int_t cent) const {
337 if(!track) return kFALSE;
338 if(!track->IsPhysicalPrimary()) return kFALSE;
339 if (track->Charge() == 0 ) return kFALSE;
341 fHistQA[8]->Fill(cent,track->Pt());
342 fHistQA[9]->Fill(cent,track->Eta());
345 if (track->Eta() < fEtaLowerLimit ||
346 track->Eta() > fEtaHigherLimit) return kFALSE;
347 if (track->Pt() < fPtLowerLimit ||
348 track->Pt() > fPtHigherLimit) return kFALSE;
351 fHistQA[12]->Fill(cent,track->Pt());
352 fHistQA[13]->Fill(cent,track->Eta());