end-of-line normalization
[u/mrichter/AliRoot.git] / PWGCF / EBYE / Fluctuations / AliEbyEFluctuationAnalysisTaskTrain.cxx
CommitLineData
a65a7e70 1#include "TChain.h"
2#include "TTree.h"
3#include "TH1F.h"
4#include "TH2F.h"
5#include "TCanvas.h"
6#include "TParticle.h"
7#include "TObjArray.h"
8
9#include "AliAnalysisTask.h"
10#include "AliAnalysisManager.h"
11
12#include "AliStack.h"
13#include "AliMCEvent.h"
14#include "AliESDEvent.h"
15#include "AliESDInputHandler.h"
16//#include "AliESDtrackCuts.h"
17#include "AliCentrality.h"
18
19#include "AliEbyEFluctuationAnalysisTaskTrain.h"
20
21// Event by event charge fluctuation analysis
22// Authors: Satyajit Jena and Panos Cristakoglou
23//
24
25ClassImp(AliEbyEFluctuationAnalysisTaskTrain)
26
27//________________________________________________________________________
28AliEbyEFluctuationAnalysisTaskTrain::AliEbyEFluctuationAnalysisTaskTrain(const char *name)
29 : AliAnalysisTaskSE(name), fESD(0), fOutputList(0),
30 fHistEventStats(0), fHistCentrality(0),
31 fHistNMult(0), fHistNPlusNMinus(0),
32 fHistNMultMC(0), fHistNPlusNMinusMC(0),
33 fAnalysisType("ESD"), fAnalysisMode("TPC"),
34 fCentralityEstimator("V0M"),
35 fCentralityPercentileMin(0.), fCentralityPercentileMax(5.),
36 fVxMax(3.), fVyMax(3.), fVzMax(10.),
37 fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0),
38 fMaxDCAxy(3.0), fMaxDCAz(3.0) {
39 // Constructor
40
41 // Define input and output slots here
42 // Input slot #0 works with a TChain
43 DefineInput(0, TChain::Class());
44 // Output slot #0 id reserved by the base class for AOD
45 // Output slot #1 writes into a TH1 container
46 DefineOutput(1, TList::Class());
47}
48
49//________________________________________________________________________
50void AliEbyEFluctuationAnalysisTaskTrain::UserCreateOutputObjects() {
51 // Create histograms
52 // Called once
53
54 fOutputList = new TList();
55 fOutputList->SetOwner();
56
57 //Event stats.
58 TString gCutName[4] = {"Total","Offline trigger",
59 "Vertex","Analyzed"};
60 fHistEventStats = new TH1F("fHistEventStats",
61 "Event statistics;;N_{events}",
62 4,0.5,4.5);
63 for(Int_t i = 1; i <= 4; i++)
64 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
65 fOutputList->Add(fHistEventStats);
66
67 //ESD analysis
68 if(fAnalysisType.CompareTo("ESD") == 0 ) {
69 fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
70 20,0.5,20.5);
71 fOutputList->Add(fHistCentrality);
72
73 TString histName;
74 histName = "fHistNMult";
75 fHistNMult = new TH1F(histName.Data(),
76 ";N_{mult.}",
77 800,0,4000);
78 fOutputList->Add(fHistNMult);
79
80 histName = "fHistNPlusNMinus";
81 fHistNPlusNMinus = new TH2F(histName.Data(),
82 ";N_{+};N_{-}",
83 2000,0.5,2000.5,2000,0.5,2000.5);
84 fOutputList->Add(fHistNPlusNMinus);
85 }//ESD analysis
86 else if(fAnalysisType.CompareTo("MC") == 0) {
87 TString histName = "fHistNMultMC";
88 fHistNMultMC = new TH1F(histName.Data(),
89 ";N_{mult.}",
90 800,0,8000);
91 fOutputList->Add(fHistNMultMC);
92
93 histName = "fHistNPlusNMinusMC";
94 fHistNPlusNMinusMC = new TH2F(histName.Data(),
95 ";N_{+};N_{-}",
96 3000,0.5,3000.5,3000,0.5,3000.5);
97 fOutputList->Add(fHistNPlusNMinusMC);
98 }//MC analysis
99
100 PostData(1, fOutputList);
101}
102
103//________________________________________________________________________
104void AliEbyEFluctuationAnalysisTaskTrain::UserExec(Option_t *) {
105 // Main loop
106 // Called for each event
107
108 Int_t nPlus = 0, nMinus = 0;
109
110 // Post output data.
111 //ESD analysis
112 if(fAnalysisType.CompareTo("ESD") == 0) {
113 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
114 if (!fESD) {
115 printf("ERROR: fESD not available\n");
116 return;
117 }
118
119 fHistEventStats->Fill(1); //all events
120
121 //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
122 //if(isSelected) {
123 fHistEventStats->Fill(2); //triggered + centrality
124
125
126 //Centrality stuff
127 AliCentrality *centrality = fESD->GetCentrality();
128 Int_t nCentrality = 0;
129
130 if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
131 fCentralityPercentileMax,
132 fCentralityEstimator.Data())) {
133 if(fAnalysisMode.CompareTo("TPC") == 0 ) {
134 const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
135 if(vertex) {
136 if(vertex->GetNContributors() > 0) {
137 if(vertex->GetZRes() != 0) {
138 fHistEventStats->Fill(3); //events with a proper vertex
139 if(TMath::Abs(vertex->GetXv()) < fVxMax) {
140 if(TMath::Abs(vertex->GetYv()) < fVyMax) {
141 if(TMath::Abs(vertex->GetZv()) < fVzMax) {
142 fHistEventStats->Fill(4); //analyzed events
143
144 Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",
145 centrality->GetCentralityPercentile(fCentralityEstimator.Data()),
146 nCentrality,fESD->GetNumberOfTracks());
147
148
149 Int_t nAcceptedTracks = 0;
150
151 nAcceptedTracks = fESD->GetNumberOfTracks();
152 AliESDtrack* track = 0;
153
154 Float_t dcaXY = 0.0, dcaZ = 0.0;
155
156 for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++)
157 {
158 track = dynamic_cast<AliESDtrack *>(fESD->GetTrack(iTracks));
159 if (!track) {
160 printf("ERROR: Could not receive track %d\n", iTracks);
161 continue;
162 }
163
164 AliESDtrack *tpcOnlyTrack = new AliESDtrack();
165
166 if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {
167 delete tpcOnlyTrack;
168 continue;
169 }
170
171 tpcOnlyTrack->GetImpactParameters(dcaXY,dcaZ);
172
173 Int_t nClustersITS = track->GetITSclusters(0x0);
174 Int_t nClustersTPC = track->GetTPCclusters(0x0);
175
176 Float_t chi2PerClusterITS = -1;
177 if (nClustersITS!=0)
178 chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
179 Float_t chi2PerClusterTPC = -1;
180 if (nClustersTPC!=0)
181 chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
182
183 //
184
185 if(nClustersTPC < fMinNumberOfTPCClusters) continue;
186 if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) continue;
187 if(TMath::Abs(dcaXY) > fMaxDCAxy) continue;
188 if(TMath::Abs(dcaZ) > fMaxDCAz) continue;
189
190 Short_t gCharge = tpcOnlyTrack->Charge();
191 if(gCharge > 0) nPlus += 1;
192 if(gCharge < 0) nMinus += 1;
193 delete tpcOnlyTrack;
194 } // track loop
195
196 if( nPlus == 0 && nMinus == 0) return;
197
198 // cout<<"=========NPlus: "<<nPlus<<"=========NMinus: "<<nMinus<<endl;
199 fHistCentrality->Fill(nCentrality);
200 fHistNPlusNMinus->Fill(nPlus,nMinus);
201 fHistNMult->Fill(nPlus+nMinus);
202
203
204 }//Vz cut
205 }//Vy cut
206 }//Vx cut
207 }//Vz resolution
208 }//number of contributors
209 }//valid vertex
210 }//TPC analysis mode
211 }//centrality
212 //}//physics selection
213 }//ESD analysis level
214
215 //MC analysis
216 if(fAnalysisType.CompareTo("MC") == 0 ) {
217 AliMCEvent* mcEvent = MCEvent();
218 if (!mcEvent) {
219 Printf("ERROR: Could not retrieve MC event");
220 return;
221 }
222 AliStack* stack = mcEvent->Stack();
223 if (!stack) {
224 Printf("ERROR: Could not retrieve MC stack");
225 return;
226 }
227
228 fHistEventStats->Fill(1);
229 fHistEventStats->Fill(2);
230 fHistEventStats->Fill(3);
231 fHistEventStats->Fill(4); //analyzed events
232 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
233 AliVParticle* track = mcEvent->GetTrack(iTracks);
234 if (!track) {
235 Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
236 continue;
237 }
238
239 if(!stack->IsPhysicalPrimary(iTracks)) continue;
240 if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;
241 if(TMath::Abs(track->Eta()) > 0.8) continue;
242
243 Short_t gCharge = track->Charge();
244 if(gCharge > 0) nPlus += 1;
245 if(gCharge < 0) nMinus += 1;
246 }//particle loop
247 fHistNPlusNMinusMC->Fill(nPlus,nMinus);
248 fHistNMultMC->Fill(nPlus+nMinus);
249
250 }//MC analysis level
251
252
253}
254
255//________________________________________________________________________
256void AliEbyEFluctuationAnalysisTaskTrain::Terminate(Option_t *) {
257 // Draw result to the screen
258 // Called once at the end of the query
259
260 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
261 if (!fOutputList) {
262 printf("ERROR: Output list not available\n");
263 return;
264 }
265}