9 #include "AliAnalysisTask.h"
10 #include "AliAnalysisManager.h"
13 #include "AliMCEvent.h"
14 #include "AliESDEvent.h"
15 #include "AliESDInputHandler.h"
16 #include "AliESDtrackCuts.h"
17 #include "AliCentrality.h"
19 #include "AliESDtrack.h"
20 #include "AliESDpid.h"
22 #include "AliEbyEFluctuationAnalysisTask.h"
24 // Event by event charge fluctuation analysis
25 // Authors: Satyajit Jena and Panos Cristakoglou
29 ClassImp(AliEbyEFluctuationAnalysisTask)
31 //________________________________________________________________________
32 AliEbyEFluctuationAnalysisTask::AliEbyEFluctuationAnalysisTask(const char *name)
33 : AliAnalysisTaskSE(name), fESD(0), fOutputList(0),
34 fHistEventStats(0), fHistCentrality(0),
35 fHistNMultMC(0), fHistNPlusNMinusMC(0),
37 fAnalysisType("ESD"), fAnalysisMode("TPC"),
38 fCentralityEstimator("V0M"), fCentralityBins20(kFALSE),
39 fVxMax(3.), fVyMax(3.), fVzMax(10.),
43 fOutputList = new TList();
44 fOutputList->SetOwner();
46 for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
47 fHistNMult[iBin-1] = NULL;
48 fHistNPlusNMinus[iBin-1] = NULL;
51 // Define input and output slots here
52 // Input slot #0 works with a TChain
53 // DefineInput(0, TChain::Class());
54 // Output slot #0 id reserved by the base class for AOD
55 // Output slot #1 writes into a TH1 container
56 DefineOutput(2, TList::Class());
59 //________________________________________________________________________
60 void AliEbyEFluctuationAnalysisTask::UserCreateOutputObjects() {
64 printf("in UserCreateOutputObjects\n");
67 TString gCutName[4] = {"Total","Offline trigger",
69 fHistEventStats = new TH1F("fHistEventStats",
70 "Event statistics;;N_{events}",
72 for(Int_t i = 1; i <= 4; i++)
73 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
74 fOutputList->Add(fHistEventStats);
77 if(fAnalysisType == "ESD") {
78 fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
79 nCentralityBins,0.5,nCentralityBins+0.5);
80 fOutputList->Add(fHistCentrality);
83 for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
84 histName = "fHistNMult"; histName += "Centrality"; histName += iBin;
85 fHistNMult[iBin-1] = new TH1F(histName.Data(),
88 fOutputList->Add(fHistNMult[iBin-1]);
90 for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
91 histName = "fHistNPlusNMinus"; histName += "Centrality"; histName += iBin;
92 fHistNPlusNMinus[iBin-1] = new TH2F(histName.Data(),
94 2000,0.5,2000.5,2000,0.5,2000.5);
95 fOutputList->Add(fHistNPlusNMinus[iBin-1]);
98 else if(fAnalysisType == "MC") {
99 TString histName = "fHistNMultMC";
100 fHistNMultMC = new TH1F(histName.Data(),
103 fOutputList->Add(fHistNMultMC);
105 histName = "fHistNPlusNMinusMC";
106 fHistNPlusNMinusMC = new TH2F(histName.Data(),
108 3000,0.5,3000.5,3000,0.5,3000.5);
109 fOutputList->Add(fHistNPlusNMinusMC);
112 PostData(2, fOutputList);
115 //________________________________________________________________________
116 void AliEbyEFluctuationAnalysisTask::UserExec(Option_t *) {
118 // Called for each event
120 Int_t nPlus = 0, nMinus = 0;
124 if(fAnalysisType == "ESD") {
125 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
127 printf("ERROR: fESD not available\n");
131 fHistEventStats->Fill(1); //all events
133 //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
135 fHistEventStats->Fill(2); //triggered + centrality
136 Printf("Event accepted");
139 AliCentrality *centrality = fESD->GetCentrality();
140 Int_t nCentrality = 0;
141 if(fCentralityBins20)
142 nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());
144 nCentrality = centrality->GetCentralityClass10(fCentralityEstimator.Data());
146 if((nCentrality < 0)||(nCentrality > 19)) return;
148 if(fAnalysisMode == "TPC") {
149 const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
151 if(vertex->GetNContributors() > 0) {
152 if(vertex->GetZRes() != 0) {
153 fHistEventStats->Fill(3); //events with a proper vertex
154 if(TMath::Abs(vertex->GetXv()) < fVxMax) {
155 if(TMath::Abs(vertex->GetYv()) < fVyMax) {
156 if(TMath::Abs(vertex->GetZv()) < fVzMax) {
157 fHistEventStats->Fill(4); //analyzed events
159 //Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",
160 //centrality->GetCentralityPercentile(fCentralityEstimator.Data()),
161 //nCentrality,fESD->GetNumberOfTracks());
163 Int_t nAcceptedTracks = 0;
164 TObjArray *gTrackArray = 0;
166 gTrackArray = fESDtrackCuts->GetAcceptedTracks(fESD,kTRUE);
168 nAcceptedTracks = gTrackArray->GetEntries();
170 AliESDtrack* track = 0;
171 // Track loop to fill a pT spectrum
172 for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++) {
173 track = dynamic_cast<AliESDtrack *>(gTrackArray->At(iTracks));
175 printf("ERROR: Could not receive track %d\n", iTracks);
179 Short_t gCharge = track->Charge();
180 if(gCharge > 0) nPlus += 1;
181 if(gCharge < 0) nMinus += 1;
185 }//TObjArray valid object
186 //if((nCentrality >= 1)&&(nCentrality <= 20)) {
188 fHistCentrality->Fill(nCentrality);
189 fHistNPlusNMinus[nCentrality-1]->Fill(nPlus,nMinus);
190 fHistNMult[nCentrality-1]->Fill(nPlus+nMinus);
196 }//number of contributors
199 //}//physics selection
200 }//ESD analysis level
203 if(fAnalysisType == "MC") {
204 AliMCEvent* mcEvent = MCEvent();
206 Printf("ERROR: Could not retrieve MC event");
209 AliStack* stack = mcEvent->Stack();
211 Printf("ERROR: Could not retrieve MC stack");
215 fHistEventStats->Fill(1);
216 fHistEventStats->Fill(2);
217 fHistEventStats->Fill(3);
218 fHistEventStats->Fill(4); //analyzed events
219 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
220 AliVParticle* track = mcEvent->GetTrack(iTracks);
222 Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
226 if(!stack->IsPhysicalPrimary(iTracks)) continue;
227 if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;
228 if(TMath::Abs(track->Eta()) > 0.8) continue;
230 Short_t gCharge = track->Charge();
231 if(gCharge > 0) nPlus += 1;
232 if(gCharge < 0) nMinus += 1;
234 fHistNPlusNMinusMC->Fill(nPlus,nMinus);
235 fHistNMultMC->Fill(nPlus+nMinus);
240 PostData(2, fOutputList);
243 //________________________________________________________________________
244 void AliEbyEFluctuationAnalysisTask::Terminate(Option_t *) {
245 // Draw result to the screen
246 // Called once at the end of the query
248 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
250 printf("ERROR: Output list not available\n");
255 //___________________________________________________________
258 AliEbyEFluctuationAnalysisTask::HasTPCPID(AliESDtrack *track) const
264 /* check PID signal */
265 if (track->GetTPCsignal() <= 0. ||
266 track->GetTPCsignalN() == 0 ||
267 !track->GetInnerParam()) return kFALSE;
271 //___________________________________________________________
274 AliEbyEFluctuationAnalysisTask::HasTOFPID(AliESDtrack *track) const
280 /* check TOF matched track */
281 if (!(track->GetStatus() & AliESDtrack::kTOFout)||
282 !(track->GetStatus() & AliESDtrack::kTIME)) return kFALSE;
283 /* check integrated length */
284 if (track->GetIntegratedLength() < 350.) return kFALSE;
288 //___________________________________________________________
291 AliEbyEFluctuationAnalysisTask::MakeTPCPID(AliESDtrack *track, Double_t *nSigma) const
295 * returns measured dEdx if PID available, otherwise -1.
296 * fills nSigma array with TPC nsigmas for e, mu, pi, K, p
300 if (!HasTPCPID(track)) return -1.;
303 Double_t ptpc = track->GetInnerParam() ? track->GetInnerParam()->P() : 0.;
304 Double_t dEdx = track->GetTPCsignal();
305 Double_t dEdxN = track->GetTPCsignalN();
307 /* loop over particles */
308 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
309 Double_t bethe = fESDpid->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
310 Double_t diff = dEdx - bethe;
311 Double_t sigma = fESDpid->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
312 nSigma[ipart] = diff / sigma;
318 //___________________________________________________________
321 AliEbyEFluctuationAnalysisTask::MakeTOFPID(AliESDtrack *track, Double_t *nSigma) const
325 * returns measured beta if PID available, otherwise -1.
326 * fills nSigma array with TOF nsigmas for e, mu, pi, K, p
330 if (!HasTOFPID(track)) return -1.;
333 Double_t p = track->P();
334 Double_t time = track->GetTOFsignal() - fESDpid->GetTOFResponse().GetStartTime(p);
335 Double_t length = track->GetIntegratedLength();
337 track->GetIntegratedTimes(timei);
339 /* loop over particles */
340 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
341 Double_t timez = time - timei[ipart];
342 Double_t sigma = fESDpid->GetTOFResponse().GetExpectedSigma(p, timei[ipart], AliPID::ParticleMass(ipart));
343 nSigma[ipart] = timez / sigma;
346 return length / (time * 2.99792457999999984e-02);
349 //___________________________________________________________
352 AliEbyEFluctuationAnalysisTask::MakePID(AliESDtrack *track, Bool_t *pidFlag) const
357 * fills pidFlag array with PID flags for e, mu, pi, K, p
361 (better put them as static variables so they can be changed from outside) */
362 Double_t fgTPCPIDmomcut[AliPID::kSPECIES] = {0., 0., 0.5, 0.5, 0.7};
363 Double_t fgTPCPIDsigmacut[AliPID::kSPECIES] = {0., 0., 2., 2., 2.};
364 Double_t fgTPCPIDcompcut[AliPID::kSPECIES] = {0., 0., 3., 3., 3.};
365 Double_t fgTOFPIDmomcut[AliPID::kSPECIES] = {0., 0., 1.5, 1.5, 2.0};
366 Double_t fgTOFPIDsigmacut[AliPID::kSPECIES] = {0., 0., 2., 2., 2.};
368 /* make pid and check if available */
369 Double_t p = track->P();
370 Double_t pt = track->Pt();
371 Double_t ptpc = track->GetInnerParam() ? track->GetInnerParam()->P() : 0.;
372 Double_t nsigmaTPC[AliPID::kSPECIES];
373 Double_t nsigmaTOF[AliPID::kSPECIES];
374 Double_t dEdx = MakeTPCPID(track, nsigmaTPC);
375 Double_t beta = MakeTOFPID(track, nsigmaTPC);
376 Bool_t hasTPCPID = dEdx > 0.;
377 Bool_t hasTOFPID = beta > 0.;
381 fHistoTPCdEdx->Fill(ptpc, dEdx);
383 fHistoTOFbeta->Fill(p, beta);
385 /* loop over species */
386 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
389 pidFlag[ipart] = kFALSE;
393 fHistoNSigmaTPC[ipart]->Fill(pt, nsigmaTPC[ipart]);
395 fHistoNSigmaTOF[ipart]->Fill(pt, nsigmaTOF[ipart]);
397 /* combined PID cuts */
398 if (hasTPCPID && hasTOFPID) {
399 if (pt < fgTOFPIDmomcut[ipart] &&
400 TMath::Abs(nsigmaTOF[ipart]) < fgTOFPIDsigmacut[ipart] &&
401 TMath::Abs(nsigmaTPC[ipart]) < fgTPCPIDcompcut[ipart])
402 pidFlag[ipart] = kTRUE;
405 /* TPC-only PID cuts */
406 else if (hasTPCPID && !hasTOFPID) {
407 if (pt < fgTPCPIDmomcut[ipart] &&
408 TMath::Abs(nsigmaTPC[ipart]) < fgTPCPIDsigmacut[ipart])
409 pidFlag[ipart] = kTRUE;
411 /* TOF-only PID cuts */
412 else if (!hasTPCPID && hasTOFPID) {
413 if (pt < fgTOFPIDmomcut[ipart] &&
414 TMath::Abs(nsigmaTOF[ipart]) < fgTOFPIDsigmacut[ipart])
415 pidFlag[ipart] = kTRUE;
418 } /* end of loop over species */
422 //___________________________________________________________
425 AliEbyEFluctuationAnalysisTask::InitPID(AliESDEvent *event)
431 /* create ESDpid object if not there yet */
433 /* instance object */
434 Bool_t mcFlag = kFALSE; /*** WARNING: check whether is MC ***/
435 fESDpid = new AliESDpid(mcFlag);
437 fESDpid->SetOADBPath("$ALICE_ROOT/OADB");
441 Int_t recoPass = 2; /*** WARNING: need to set the recoPass somehow better ***/
442 fESDpid->InitialiseEvent(event, recoPass); /* warning: this apparently sets TOF time
443 * resolution to some hardcoded value,
444 * therefore we have to set correct
445 * resolution value after this call */
447 /* set TOF resolution */
448 Double_t tofReso = 85.; /* ps */ /*** WARNING: need to set tofReso somehow better ***/
449 fESDpid->GetTOFResponse().SetTimeResolution(tofReso);