Proper usage of the bayesian approach (Alis Rodriguez Manso)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskBF.cxx
CommitLineData
3feee083 1#include "TChain.h"\r
2#include "TList.h"\r
3#include "TCanvas.h"\r
4#include "TLorentzVector.h"\r
5#include "TGraphErrors.h"\r
6#include "TH1F.h"\r
7#include "TH2F.h"\r
6e28a437 8#include "TArrayF.h"\r
6fbd0608 9#include "TF1.h"\r
10#include "TRandom.h"\r
3feee083 11\r
12#include "AliAnalysisTaskSE.h"\r
13#include "AliAnalysisManager.h"\r
14\r
15#include "AliESDVertex.h"\r
16#include "AliESDEvent.h"\r
17#include "AliESDInputHandler.h"\r
18#include "AliAODEvent.h"\r
19#include "AliAODTrack.h"\r
20#include "AliAODInputHandler.h"\r
6e28a437 21#include "AliGenEventHeader.h"\r
22#include "AliGenHijingEventHeader.h"\r
3feee083 23#include "AliMCEventHandler.h"\r
24#include "AliMCEvent.h"\r
25#include "AliStack.h"\r
26#include "AliESDtrackCuts.h"\r
27\r
3c3bbd53 28#include "TH2D.h" \r
29#include "AliPID.h" \r
30#include "AliPIDResponse.h" \r
31#include "AliPIDCombined.h" \r
32\r
3feee083 33#include "AliAnalysisTaskBF.h"\r
34#include "AliBalance.h"\r
35\r
36\r
37// Analysis task for the BF code\r
6e28a437 38// Authors: Panos.Christakoglou@nikhef.nl\r
3feee083 39\r
40ClassImp(AliAnalysisTaskBF)\r
41\r
42//________________________________________________________________________\r
43AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name) \r
44: AliAnalysisTaskSE(name), \r
45 fBalance(0),\r
10a747d2 46 fRunShuffling(kFALSE),\r
1900bf4e 47 fShuffledBalance(0),\r
3feee083 48 fList(0),\r
49 fListBF(0),\r
6432ac6a 50 fListBFS(0),\r
3c3bbd53 51 fHistListPIDQA(0),\r
3feee083 52 fHistEventStats(0),\r
a0ce3c86 53 fHistCentStats(0),\r
a839b0a3 54 fHistTriggerStats(0),\r
3feee083 55 fHistTrackStats(0),\r
56 fHistVx(0),\r
57 fHistVy(0),\r
58 fHistVz(0),\r
59 fHistClus(0),\r
60 fHistDCA(0),\r
61 fHistChi2(0),\r
62 fHistPt(0),\r
63 fHistEta(0),\r
64 fHistPhi(0),\r
4a1e5d10 65 fHistPhiBefore(0),\r
66 fHistPhiAfter(0),\r
3feee083 67 fHistV0M(0),\r
a0ce3c86 68 fHistRefTracks(0),\r
3c3bbd53 69 fHistdEdxVsPTPCbeforePID(NULL),\r
70 fHistBetavsPTOFbeforePID(NULL), \r
71 fHistProbTPCvsPtbeforePID(NULL), \r
72 fHistProbTOFvsPtbeforePID(NULL), \r
da31bdf4 73 fHistProbTPCTOFvsPtbeforePID(NULL),\r
3c3bbd53 74 fHistNSigmaTPCvsPtbeforePID(NULL), \r
75 fHistNSigmaTOFvsPtbeforePID(NULL), \r
76 fHistdEdxVsPTPCafterPID(NULL),\r
77 fHistBetavsPTOFafterPID(NULL), \r
78 fHistProbTPCvsPtafterPID(NULL), \r
79 fHistProbTOFvsPtafterPID(NULL), \r
da31bdf4 80 fHistProbTPCTOFvsPtafterPID(NULL),\r
3c3bbd53 81 fHistNSigmaTPCvsPtafterPID(NULL), \r
82 fHistNSigmaTOFvsPtafterPID(NULL), \r
83 fPIDResponse(0x0),\r
84 fPIDCombined(0x0),\r
85 fParticleOfInterest(kPion),\r
86 fPidDetectorConfig(kTPCTOF),\r
87 fUsePID(kFALSE),\r
88 fUsePIDnSigma(kTRUE),\r
89 fUsePIDPropabilities(kFALSE), \r
90 fPIDNSigma(3.),\r
91 fMinAcceptedPIDProbability(0.8),\r
3feee083 92 fESDtrackCuts(0),\r
32a4f715 93 fCentralityEstimator("V0M"),\r
c735aca6 94 fUseCentrality(kFALSE),\r
3feee083 95 fCentralityPercentileMin(0.), \r
96 fCentralityPercentileMax(5.),\r
6e28a437 97 fImpactParameterMin(0.),\r
98 fImpactParameterMax(20.),\r
99 fUseMultiplicity(kFALSE),\r
100 fNumberOfAcceptedTracksMin(0),\r
101 fNumberOfAcceptedTracksMax(10000),\r
a5ab9ca3 102 fHistNumberOfAcceptedTracks(0),\r
3feee083 103 fUseOfflineTrigger(kFALSE),\r
104 fVxMax(0.3),\r
105 fVyMax(0.3),\r
106 fVzMax(10.),\r
6432ac6a 107 nAODtrackCutBit(128),\r
3feee083 108 fPtMin(0.3),\r
109 fPtMax(1.5),\r
110 fEtaMin(-0.8),\r
1900bf4e 111 fEtaMax(-0.8),\r
a839b0a3 112 fDCAxyCut(-1),\r
113 fDCAzCut(-1),\r
114 fTPCchi2Cut(-1),\r
6fbd0608 115 fNClustersTPCCut(-1),\r
a1390237 116 fAcceptanceParameterization(0),\r
985da7aa 117 fDifferentialV2(0),\r
118 fUseFlowAfterBurner(kFALSE),\r
d9f855f3 119 fExcludeResonancesInMC(kFALSE),\r
120 fUseMCPdgCode(kFALSE),\r
121 fPDGCodeToBeAnalyzed(-1) {\r
3feee083 122 // Constructor\r
3feee083 123 // Define input and output slots here\r
124 // Input slot #0 works with a TChain\r
125 DefineInput(0, TChain::Class());\r
126 // Output slot #0 writes into a TH1 container\r
6432ac6a 127 DefineOutput(1, TList::Class());\r
128 DefineOutput(2, TList::Class());\r
3feee083 129 DefineOutput(3, TList::Class());\r
3c3bbd53 130 DefineOutput(4, TList::Class());\r
3feee083 131}\r
132\r
6432ac6a 133//________________________________________________________________________\r
134AliAnalysisTaskBF::~AliAnalysisTaskBF() {\r
135\r
136 // delete fBalance; \r
137 // delete fShuffledBalance; \r
138 // delete fList;\r
139 // delete fListBF; \r
140 // delete fListBFS;\r
141\r
142 // delete fHistEventStats; \r
143 // delete fHistTrackStats; \r
144 // delete fHistVx; \r
145 // delete fHistVy; \r
146 // delete fHistVz; \r
147\r
148 // delete fHistClus;\r
149 // delete fHistDCA;\r
150 // delete fHistChi2;\r
151 // delete fHistPt;\r
152 // delete fHistEta;\r
153 // delete fHistPhi;\r
154 // delete fHistV0M;\r
155}\r
3feee083 156\r
157//________________________________________________________________________\r
158void AliAnalysisTaskBF::UserCreateOutputObjects() {\r
159 // Create histograms\r
160 // Called once\r
161 if(!fBalance) {\r
162 fBalance = new AliBalance();\r
163 fBalance->SetAnalysisLevel("ESD");\r
6432ac6a 164 //fBalance->SetNumberOfBins(-1,16);\r
165 fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
3feee083 166 }\r
10a747d2 167 if(fRunShuffling) {\r
168 if(!fShuffledBalance) {\r
169 fShuffledBalance = new AliBalance();\r
170 fShuffledBalance->SetAnalysisLevel("ESD");\r
171 //fShuffledBalance->SetNumberOfBins(-1,16);\r
172 fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
173 }\r
1900bf4e 174 }\r
3feee083 175\r
176 //QA list\r
177 fList = new TList();\r
178 fList->SetName("listQA");\r
179 fList->SetOwner();\r
180\r
181 //Balance Function list\r
182 fListBF = new TList();\r
183 fListBF->SetName("listBF");\r
184 fListBF->SetOwner();\r
185\r
10a747d2 186 if(fRunShuffling) {\r
187 fListBFS = new TList();\r
188 fListBFS->SetName("listBFShuffled");\r
189 fListBFS->SetOwner();\r
190 }\r
6432ac6a 191\r
3c3bbd53 192 //PID QA list\r
193 if(fUsePID) {\r
194 fHistListPIDQA = new TList();\r
195 fHistListPIDQA->SetName("listQAPID");\r
196 fHistListPIDQA->SetOwner();\r
197 }\r
198\r
3feee083 199 //Event stats.\r
200 TString gCutName[4] = {"Total","Offline trigger",\r
201 "Vertex","Analyzed"};\r
202 fHistEventStats = new TH1F("fHistEventStats",\r
203 "Event statistics;;N_{events}",\r
204 4,0.5,4.5);\r
205 for(Int_t i = 1; i <= 4; i++)\r
206 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
207 fList->Add(fHistEventStats);\r
208\r
a0ce3c86 209 TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
210 fHistCentStats = new TH2F("fHistCentStats",\r
211 "Centrality statistics;;Cent percentile",\r
212 9,-0.5,8.5,220,-5,105);\r
213 for(Int_t i = 1; i <= 9; i++)\r
214 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
215 fList->Add(fHistCentStats);\r
216\r
a839b0a3 217 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
218 fList->Add(fHistTriggerStats);\r
219\r
a0ce3c86 220 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
3feee083 221 fList->Add(fHistTrackStats);\r
222\r
a5ab9ca3 223 fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);\r
224 fList->Add(fHistNumberOfAcceptedTracks);\r
225\r
3feee083 226 // Vertex distributions\r
227 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
228 fList->Add(fHistVx);\r
229 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
230 fList->Add(fHistVy);\r
231 fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
232 fList->Add(fHistVz);\r
233\r
234 // QA histograms\r
235 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
236 fList->Add(fHistClus);\r
237 fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
238 fList->Add(fHistChi2);\r
239 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
240 fList->Add(fHistDCA);\r
241 fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
242 fList->Add(fHistPt);\r
243 fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
244 fList->Add(fHistEta);\r
245 fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
246 fList->Add(fHistPhi);\r
4a1e5d10 247 fHistPhiBefore = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());\r
248 fList->Add(fHistPhiBefore);\r
249 fHistPhiAfter = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());\r
250 fList->Add(fHistPhiAfter);\r
3feee083 251 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
252 fList->Add(fHistV0M);\r
a0ce3c86 253 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
254 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
255 for(Int_t i = 1; i <= 6; i++)\r
256 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
257 fList->Add(fHistRefTracks);\r
3feee083 258\r
6432ac6a 259 // Balance function histograms\r
6432ac6a 260 // Initialize histograms if not done yet\r
10a747d2 261 if(!fBalance->GetHistNp(0)){\r
6432ac6a 262 AliWarning("Histograms not yet initialized! --> Will be done now");\r
263 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
264 fBalance->InitHistograms();\r
10a747d2 265 }\r
266\r
267 if(fRunShuffling) {\r
268 if(!fShuffledBalance->GetHistNp(0)) {\r
269 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
270 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
271 fShuffledBalance->InitHistograms();\r
272 }\r
3feee083 273 }\r
6432ac6a 274\r
275 for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
276 fListBF->Add(fBalance->GetHistNp(a));\r
277 fListBF->Add(fBalance->GetHistNn(a));\r
278 fListBF->Add(fBalance->GetHistNpn(a));\r
279 fListBF->Add(fBalance->GetHistNnn(a));\r
280 fListBF->Add(fBalance->GetHistNpp(a));\r
281 fListBF->Add(fBalance->GetHistNnp(a));\r
282\r
10a747d2 283 if(fRunShuffling) {\r
284 fListBFS->Add(fShuffledBalance->GetHistNp(a));\r
285 fListBFS->Add(fShuffledBalance->GetHistNn(a));\r
286 fListBFS->Add(fShuffledBalance->GetHistNpn(a));\r
287 fListBFS->Add(fShuffledBalance->GetHistNnn(a));\r
288 fListBFS->Add(fShuffledBalance->GetHistNpp(a));\r
289 fListBFS->Add(fShuffledBalance->GetHistNnp(a));\r
290 } \r
291 }\r
6432ac6a 292\r
293 if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
3feee083 294\r
3c3bbd53 295 //====================PID========================//\r
296 if(fUsePID) {\r
7a32a667 297 fPIDCombined = new AliPIDCombined();\r
298 fPIDCombined->SetDefaultTPCPriors();\r
299\r
300 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
3c3bbd53 301 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
302 \r
7a32a667 303 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
3c3bbd53 304 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
305 \r
da31bdf4 306 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
3c3bbd53 307 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
308 \r
da31bdf4 309 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
3c3bbd53 310 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
da31bdf4 311\r
312 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
313 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
3c3bbd53 314 \r
7a32a667 315 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
3c3bbd53 316 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
317 \r
7a32a667 318 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
3c3bbd53 319 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
320 \r
7a32a667 321 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
3c3bbd53 322 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
323 \r
7a32a667 324 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
3c3bbd53 325 fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
326 \r
7a32a667 327 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
3c3bbd53 328 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
329 \r
7a32a667 330 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
3c3bbd53 331 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition \r
332 \r
da31bdf4 333 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
334 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
335\r
7a32a667 336 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
3c3bbd53 337 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition \r
338 \r
7a32a667 339 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
3c3bbd53 340 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
341 }\r
342 //====================PID========================//\r
343\r
3feee083 344 // Post output data.\r
6432ac6a 345 PostData(1, fList);\r
346 PostData(2, fListBF);\r
10a747d2 347 if(fRunShuffling) PostData(3, fListBFS);\r
3c3bbd53 348 if(fUsePID) PostData(4, fHistListPIDQA); //PID\r
a839b0a3 349}\r
3feee083 350\r
351//________________________________________________________________________\r
352void AliAnalysisTaskBF::UserExec(Option_t *) {\r
353 // Main loop\r
354 // Called for each event\r
355 TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
356\r
a839b0a3 357 AliESDtrack *track_TPC = NULL;\r
358\r
6e28a437 359 Int_t gNumberOfAcceptedTracks = 0;\r
74b04dd4 360 Float_t fCentrality = 0.;\r
361\r
362 // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
363 vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled\r
364 vector<Double_t> *chargeVector[9]; // original charge\r
365 for(Int_t i = 0; i < 9; i++){\r
366 chargeVectorShuffle[i] = new vector<Double_t>;\r
367 chargeVector[i] = new vector<Double_t>;\r
368 }\r
369\r
370 Double_t v_charge;\r
371 Double_t v_y;\r
372 Double_t v_eta;\r
373 Double_t v_phi;\r
374 Double_t v_p[3];\r
375 Double_t v_pt;\r
376 Double_t v_E;\r
377\r
7a32a667 378 if(fUsePID) {\r
379 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
380 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
381 }\r
382 \r
3feee083 383 //ESD analysis\r
384 if(gAnalysisLevel == "ESD") {\r
385 AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
386 if (!gESD) {\r
387 Printf("ERROR: gESD not available");\r
388 return;\r
389 }\r
390\r
a839b0a3 391 // store offline trigger bits\r
392 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
393\r
6432ac6a 394 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
3feee083 395 fHistEventStats->Fill(1); //all events\r
396 Bool_t isSelected = kTRUE;\r
397 if(fUseOfflineTrigger)\r
398 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
399 if(isSelected) {\r
400 fHistEventStats->Fill(2); //triggered events\r
401\r
c735aca6 402 if(fUseCentrality) {\r
403 //Centrality stuff\r
404 AliCentrality *centrality = gESD->GetCentrality();\r
74b04dd4 405 \r
406 fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
c735aca6 407 \r
408 // take only events inside centrality class\r
409 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
410 fCentralityPercentileMax,\r
411 fCentralityEstimator.Data()))\r
412 return;\r
413 \r
3feee083 414 // centrality QA (V0M)\r
415 fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
c735aca6 416 }\r
3feee083 417 \r
c735aca6 418 const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
419 if(vertex) {\r
420 if(vertex->GetNContributors() > 0) {\r
421 if(vertex->GetZRes() != 0) {\r
422 fHistEventStats->Fill(3); //events with a proper vertex\r
423 if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
424 if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
425 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
426 fHistEventStats->Fill(4); //analayzed events\r
427 fHistVx->Fill(vertex->GetXv());\r
428 fHistVy->Fill(vertex->GetYv());\r
429 fHistVz->Fill(vertex->GetZv());\r
430 \r
431 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
432 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
433 AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
434 if (!track) {\r
435 Printf("ERROR: Could not receive track %d", iTracks);\r
436 continue;\r
437 } \r
3feee083 438 \r
c735aca6 439 // take only TPC only tracks\r
440 track_TPC = new AliESDtrack();\r
441 if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
442 \r
443 //ESD track cuts\r
444 if(fESDtrackCuts) \r
445 if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
446 \r
447 // fill QA histograms\r
448 Float_t b[2];\r
449 Float_t bCov[3];\r
450 track_TPC->GetImpactParameters(b,bCov);\r
451 if (bCov[0]<=0 || bCov[2]<=0) {\r
452 AliDebug(1, "Estimated b resolution lower or equal zero!");\r
453 bCov[0]=0; bCov[2]=0;\r
454 }\r
455 \r
456 Int_t nClustersTPC = -1;\r
457 nClustersTPC = track_TPC->GetTPCNclsIter1(); // TPC standalone\r
458 //nClustersTPC = track->GetTPCclusters(0); // global track\r
459 Float_t chi2PerClusterTPC = -1;\r
460 if (nClustersTPC!=0) {\r
461 chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
462 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
463 }\r
74b04dd4 464\r
3c3bbd53 465 //===========================PID===============================// \r
466 if(fUsePID) {\r
467 Double_t prob[AliPID::kSPECIES]={0.};\r
da31bdf4 468 Double_t probTPC[AliPID::kSPECIES]={0.};\r
469 Double_t probTOF[AliPID::kSPECIES]={0.};\r
470 Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
471\r
3c3bbd53 472 Double_t nSigma = 0.;\r
da31bdf4 473 UInt_t detUsedTPC = 0;\r
474 UInt_t detUsedTOF = 0;\r
475 UInt_t detUsedTPCTOF = 0;\r
476\r
3c3bbd53 477 //Decide what detector configuration we want to use\r
478 switch(fPidDetectorConfig) {\r
05d631ac 479 case kTPCpid:\r
3c3bbd53 480 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
7a32a667 481 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
da31bdf4 482 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
483 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
484 prob[iSpecies] = probTPC[iSpecies];\r
3c3bbd53 485 break;\r
05d631ac 486 case kTOFpid:\r
3c3bbd53 487 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
7a32a667 488 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
da31bdf4 489 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
490 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
491 prob[iSpecies] = probTOF[iSpecies];\r
3c3bbd53 492 break;\r
05d631ac 493 case kTPCTOF:\r
3c3bbd53 494 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
da31bdf4 495 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
496 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
497 prob[iSpecies] = probTPCTOF[iSpecies];\r
3c3bbd53 498 break;\r
499 default:\r
500 break;\r
501 }//end switch: define detector mask\r
502 \r
3c3bbd53 503 //Filling the PID QA\r
504 Double_t tofTime = -999., length = 999., tof = -999.;\r
505 Double_t c = TMath::C()*1.E-9;// m/ns\r
506 Double_t beta = -999.;\r
507 Double_t nSigmaTOFForParticleOfInterest = -999.;\r
508 if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
509 (track->IsOn(AliESDtrack::kTIME)) ) { \r
510 tofTime = track->GetTOFsignal();//in ps\r
511 length = track->GetIntegratedLength();\r
512 tof = tofTime*1E-3; // ns \r
513 \r
514 if (tof <= 0) {\r
515 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
516 continue;\r
517 }\r
518 if (length <= 0){\r
519 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
520 continue;\r
521 }\r
522 \r
523 length = length*0.01; // in meters\r
524 tof = tof*c;\r
525 beta = length/tof;\r
526 \r
527 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
528 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
da31bdf4 529 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
3c3bbd53 530 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
531 }//TOF signal \r
532 \r
533 \r
534 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
535 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
da31bdf4 536 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
3c3bbd53 537 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
da31bdf4 538 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
3c3bbd53 539 //end of QA-before pid\r
540 \r
da31bdf4 541 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
05d631ac 542 //Make the decision based on the n-sigma\r
543 if(fUsePIDnSigma) {\r
544 if(nSigma > fPIDNSigma) continue;}\r
545 \r
546 //Make the decision based on the bayesian\r
547 else if(fUsePIDPropabilities) {\r
548 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
7a32a667 549 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r
05d631ac 550 }\r
7a32a667 551 \r
552 //Fill QA after the PID\r
553 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
da31bdf4 554 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
7a32a667 555 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
556 \r
557 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
da31bdf4 558 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
559 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
7a32a667 560 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
3c3bbd53 561 }\r
562 \r
3c3bbd53 563 PostData(4, fHistListPIDQA);\r
564 }\r
565 //===========================PID===============================//\r
74b04dd4 566 v_charge = track_TPC->Charge();\r
567 v_y = track_TPC->Y();\r
568 v_eta = track_TPC->Eta();\r
569 v_phi = track_TPC->Phi() * TMath::RadToDeg();\r
570 v_E = track_TPC->E();\r
571 v_pt = track_TPC->Pt();\r
572 track_TPC->PxPyPz(v_p);\r
c735aca6 573 fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
574 fHistDCA->Fill(b[1],b[0]);\r
575 fHistChi2->Fill(chi2PerClusterTPC);\r
74b04dd4 576 fHistPt->Fill(v_pt);\r
577 fHistEta->Fill(v_eta);\r
578 fHistPhi->Fill(v_phi);\r
579\r
c735aca6 580 // fill charge vector\r
74b04dd4 581 chargeVector[0]->push_back(v_charge);\r
582 chargeVector[1]->push_back(v_y);\r
583 chargeVector[2]->push_back(v_eta);\r
584 chargeVector[3]->push_back(v_phi);\r
585 chargeVector[4]->push_back(v_p[0]);\r
586 chargeVector[5]->push_back(v_p[1]);\r
587 chargeVector[6]->push_back(v_p[2]);\r
588 chargeVector[7]->push_back(v_pt);\r
589 chargeVector[8]->push_back(v_E);\r
590\r
591 if(fRunShuffling) {\r
592 chargeVectorShuffle[0]->push_back(v_charge);\r
593 chargeVectorShuffle[1]->push_back(v_y);\r
594 chargeVectorShuffle[2]->push_back(v_eta);\r
595 chargeVectorShuffle[3]->push_back(v_phi);\r
596 chargeVectorShuffle[4]->push_back(v_p[0]);\r
597 chargeVectorShuffle[5]->push_back(v_p[1]);\r
598 chargeVectorShuffle[6]->push_back(v_p[2]);\r
599 chargeVectorShuffle[7]->push_back(v_pt);\r
600 chargeVectorShuffle[8]->push_back(v_E);\r
c735aca6 601 }\r
602 \r
603 delete track_TPC;\r
604 \r
605 } //track loop\r
606 }//Vz cut\r
607 }//Vy cut\r
608 }//Vx cut\r
609 }//proper vertex resolution\r
610 }//proper number of contributors\r
611 }//vertex object valid\r
3feee083 612 }//triggered event \r
613 }//ESD analysis\r
614 \r
3feee083 615 //AOD analysis (vertex and track cuts also here!!!!)\r
616 else if(gAnalysisLevel == "AOD") {\r
617 AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE\r
618 if(!gAOD) {\r
619 Printf("ERROR: gAOD not available");\r
620 return;\r
621 }\r
622\r
a839b0a3 623 AliAODHeader *aodHeader = gAOD->GetHeader();\r
624\r
625 // store offline trigger bits\r
626 fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());\r
6e28a437 627 \r
6432ac6a 628 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
3feee083 629 fHistEventStats->Fill(1); //all events\r
630 Bool_t isSelected = kTRUE;\r
631 if(fUseOfflineTrigger)\r
632 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
633 if(isSelected) {\r
634 fHistEventStats->Fill(2); //triggered events\r
3feee083 635 \r
636 //Centrality stuff (centrality in AOD header)\r
c735aca6 637 if(fUseCentrality) {\r
74b04dd4 638 fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
639 \r
640 // QA for centrality estimators\r
641 fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
642 fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));\r
643 fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));\r
644 fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));\r
645 fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));\r
646 fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));\r
647 fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
648 fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
649 fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
c735aca6 650 \r
74b04dd4 651 // take only events inside centrality class\r
652 if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
653 return;\r
c735aca6 654 \r
74b04dd4 655 // centrality QA (V0M)\r
656 fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());\r
c735aca6 657 \r
74b04dd4 658 // centrality QA (reference tracks)\r
659 fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());\r
660 fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());\r
661 fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());\r
662 fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());\r
663 fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));\r
664 fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));\r
665 fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));\r
666 fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));\r
667 fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));\r
c735aca6 668 }\r
a0ce3c86 669\r
c735aca6 670 const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
671 \r
672 if(vertex) {\r
74b04dd4 673 Double32_t fCov[6];\r
674 vertex->GetCovarianceMatrix(fCov);\r
c735aca6 675 \r
74b04dd4 676 if(vertex->GetNContributors() > 0) {\r
677 if(fCov[5] != 0) {\r
678 fHistEventStats->Fill(3); //events with a proper vertex\r
679 if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
680 if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
681 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
682 fHistEventStats->Fill(4); //analyzed events\r
683 fHistVx->Fill(vertex->GetX());\r
684 fHistVy->Fill(vertex->GetY());\r
685 fHistVz->Fill(vertex->GetZ());\r
c735aca6 686 \r
74b04dd4 687 //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());\r
688 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {\r
689 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));\r
690 if (!aodTrack) {\r
691 Printf("ERROR: Could not receive track %d", iTracks);\r
692 continue;\r
693 }\r
c735aca6 694 \r
74b04dd4 695 // AOD track cuts\r
c735aca6 696 \r
74b04dd4 697 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
698 // take only TPC only tracks \r
699 fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
700 if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
c735aca6 701 \r
74b04dd4 702 v_charge = aodTrack->Charge();\r
703 v_y = aodTrack->Y();\r
704 v_eta = aodTrack->Eta();\r
705 v_phi = aodTrack->Phi() * TMath::RadToDeg();\r
706 v_E = aodTrack->E();\r
707 v_pt = aodTrack->Pt();\r
708 aodTrack->PxPyPz(v_p);\r
c735aca6 709 \r
74b04dd4 710 Float_t DCAxy = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
711 Float_t DCAz = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
c735aca6 712 \r
713 \r
74b04dd4 714 // Kinematics cuts from ESD track cuts\r
715 if( v_pt < fPtMin || v_pt > fPtMax) continue;\r
716 if( v_eta < fEtaMin || v_eta > fEtaMax) continue;\r
c735aca6 717 \r
74b04dd4 718 // Extra DCA cuts (for systematic studies [!= -1])\r
719 if( fDCAxyCut != -1 && fDCAxyCut != -1){\r
720 if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
721 continue; // 2D cut\r
722 }\r
723 }\r
c735aca6 724 \r
74b04dd4 725 // Extra TPC cuts (for systematic studies [!= -1])\r
726 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
727 continue;\r
728 }\r
729 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
730 continue;\r
731 }\r
c735aca6 732 \r
74b04dd4 733 // fill QA histograms\r
734 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
735 fHistDCA->Fill(DCAz,DCAxy);\r
736 fHistChi2->Fill(aodTrack->Chi2perNDF());\r
737 fHistPt->Fill(v_pt);\r
738 fHistEta->Fill(v_eta);\r
739 fHistPhi->Fill(v_phi);\r
740\r
741 // fill charge vector\r
742 chargeVector[0]->push_back(v_charge);\r
743 chargeVector[1]->push_back(v_y);\r
744 chargeVector[2]->push_back(v_eta);\r
745 chargeVector[3]->push_back(v_phi);\r
746 chargeVector[4]->push_back(v_p[0]);\r
747 chargeVector[5]->push_back(v_p[1]);\r
748 chargeVector[6]->push_back(v_p[2]);\r
749 chargeVector[7]->push_back(v_pt);\r
750 chargeVector[8]->push_back(v_E);\r
751\r
752 if(fRunShuffling) {\r
753 chargeVectorShuffle[0]->push_back(v_charge);\r
754 chargeVectorShuffle[1]->push_back(v_y);\r
755 chargeVectorShuffle[2]->push_back(v_eta);\r
756 chargeVectorShuffle[3]->push_back(v_phi);\r
757 chargeVectorShuffle[4]->push_back(v_p[0]);\r
758 chargeVectorShuffle[5]->push_back(v_p[1]);\r
759 chargeVectorShuffle[6]->push_back(v_p[2]);\r
760 chargeVectorShuffle[7]->push_back(v_pt);\r
761 chargeVectorShuffle[8]->push_back(v_E);\r
762 }\r
763 \r
764 gNumberOfAcceptedTracks += 1;\r
c735aca6 765 \r
74b04dd4 766 } //track loop\r
767 }//Vz cut\r
768 }//Vy cut\r
769 }//Vx cut\r
770 }//proper vertex resolution\r
c735aca6 771 }//proper number of contributors\r
772 }//vertex object valid\r
3feee083 773 }//triggered event \r
774 }//AOD analysis\r
775\r
6e28a437 776 //MC-ESD analysis\r
777 if(gAnalysisLevel == "MCESD") {\r
3feee083 778 AliMCEvent* mcEvent = MCEvent(); \r
779 if (!mcEvent) {\r
780 Printf("ERROR: mcEvent not available");\r
781 return;\r
782 }\r
1900bf4e 783\r
6e28a437 784 AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
785 if (!gESD) {\r
786 Printf("ERROR: gESD not available");\r
787 return;\r
788 }\r
789\r
790 // store offline trigger bits\r
791 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
1900bf4e 792\r
6e28a437 793 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
794 fHistEventStats->Fill(1); //all events\r
795 Bool_t isSelected = kTRUE;\r
796 if(fUseOfflineTrigger)\r
797 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
798 if(isSelected) {\r
799 fHistEventStats->Fill(2); //triggered events\r
1900bf4e 800\r
6e28a437 801 if(fUseCentrality) {\r
802 //Centrality stuff\r
803 AliCentrality *centrality = gESD->GetCentrality();\r
74b04dd4 804\r
805 fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
6e28a437 806 \r
807 // take only events inside centrality class\r
808 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
809 fCentralityPercentileMax,\r
810 fCentralityEstimator.Data()))\r
811 return;\r
812 \r
813 // centrality QA (V0M)\r
814 fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
10a747d2 815 }\r
6e28a437 816 \r
817 const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
818 if(vertex) {\r
819 if(vertex->GetNContributors() > 0) {\r
820 if(vertex->GetZRes() != 0) {\r
821 fHistEventStats->Fill(3); //events with a proper vertex\r
822 if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
823 if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
824 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
825 fHistEventStats->Fill(4); //analayzed events\r
826 fHistVx->Fill(vertex->GetXv());\r
827 fHistVy->Fill(vertex->GetYv());\r
828 fHistVz->Fill(vertex->GetZv());\r
829 \r
830 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
831 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
832 AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
833 if (!track) {\r
834 Printf("ERROR: Could not receive track %d", iTracks);\r
835 continue;\r
836 } \r
837 \r
838 Int_t label = TMath::Abs(track->GetLabel());\r
839 if(label > mcEvent->GetNumberOfTracks()) continue;\r
840 if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
841 \r
842 AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
843 if(!mcTrack) continue;\r
844\r
845 // take only TPC only tracks\r
846 track_TPC = new AliESDtrack();\r
847 if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
848 \r
849 //ESD track cuts\r
850 if(fESDtrackCuts) \r
851 if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
852 \r
853 // fill QA histograms\r
854 Float_t b[2];\r
855 Float_t bCov[3];\r
856 track_TPC->GetImpactParameters(b,bCov);\r
857 if (bCov[0]<=0 || bCov[2]<=0) {\r
858 AliDebug(1, "Estimated b resolution lower or equal zero!");\r
859 bCov[0]=0; bCov[2]=0;\r
860 }\r
861 \r
862 Int_t nClustersTPC = -1;\r
863 nClustersTPC = track_TPC->GetTPCNclsIter1(); // TPC standalone\r
864 //nClustersTPC = track->GetTPCclusters(0); // global track\r
865 Float_t chi2PerClusterTPC = -1;\r
866 if (nClustersTPC!=0) {\r
867 chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
868 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
869 }\r
870 \r
74b04dd4 871 v_charge = track_TPC->Charge();\r
872 v_y = track_TPC->Y();\r
873 v_eta = track_TPC->Eta();\r
874 v_phi = track_TPC->Phi() * TMath::RadToDeg();\r
875 v_E = track_TPC->E();\r
876 v_pt = track_TPC->Pt();\r
877 track_TPC->PxPyPz(v_p);\r
878\r
6e28a437 879 fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
880 fHistDCA->Fill(b[1],b[0]);\r
881 fHistChi2->Fill(chi2PerClusterTPC);\r
74b04dd4 882 fHistPt->Fill(v_pt);\r
883 fHistEta->Fill(v_eta);\r
884 fHistPhi->Fill(v_phi);\r
885\r
6e28a437 886 // fill charge vector\r
74b04dd4 887 chargeVector[0]->push_back(v_charge);\r
888 chargeVector[1]->push_back(v_y);\r
889 chargeVector[2]->push_back(v_eta);\r
890 chargeVector[3]->push_back(v_phi);\r
891 chargeVector[4]->push_back(v_p[0]);\r
892 chargeVector[5]->push_back(v_p[1]);\r
893 chargeVector[6]->push_back(v_p[2]);\r
894 chargeVector[7]->push_back(v_pt);\r
895 chargeVector[8]->push_back(v_E);\r
896\r
897 if(fRunShuffling) {\r
898 chargeVectorShuffle[0]->push_back(v_charge);\r
899 chargeVectorShuffle[1]->push_back(v_y);\r
900 chargeVectorShuffle[2]->push_back(v_eta);\r
901 chargeVectorShuffle[3]->push_back(v_phi);\r
902 chargeVectorShuffle[4]->push_back(v_p[0]);\r
903 chargeVectorShuffle[5]->push_back(v_p[1]);\r
904 chargeVectorShuffle[6]->push_back(v_p[2]);\r
905 chargeVectorShuffle[7]->push_back(v_pt);\r
906 chargeVectorShuffle[8]->push_back(v_E);\r
6e28a437 907 }\r
908 \r
909 delete track_TPC;\r
63d15f4b 910 gNumberOfAcceptedTracks += 1;\r
6e28a437 911 \r
912 } //track loop\r
913 }//Vz cut\r
914 }//Vy cut\r
915 }//Vx cut\r
916 }//proper vertex resolution\r
917 }//proper number of contributors\r
918 }//vertex object valid\r
919 }//triggered event \r
920 }//MC-ESD analysis\r
a839b0a3 921\r
6e28a437 922 //MC analysis\r
923 else if(gAnalysisLevel == "MC") {\r
924 AliMCEvent* mcEvent = MCEvent(); \r
925 if (!mcEvent) {\r
926 Printf("ERROR: mcEvent not available");\r
927 return;\r
928 }\r
6a3f819e 929 fHistEventStats->Fill(1); //total events\r
930 fHistEventStats->Fill(2); //offline trigger\r
6e28a437 931\r
932 Double_t gReactionPlane = 0., gImpactParameter = 0.;\r
933 if(fUseCentrality) {\r
934 //Get the MC header\r
935 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
936 if (headerH) {\r
937 //Printf("=====================================================");\r
938 //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());\r
939 //Printf("Impact parameter: %lf",headerH->ImpactParameter());\r
940 //Printf("=====================================================");\r
941 gReactionPlane = headerH->ReactionPlaneAngle();\r
942 gImpactParameter = headerH->ImpactParameter();\r
4a1e5d10 943 fCentrality = gImpactParameter;\r
6e28a437 944 }\r
7ef516b2 945 fCentrality = gImpactParameter;\r
946\r
74b04dd4 947 // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW)\r
6e28a437 948 if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))\r
949 return;\r
950 }\r
951 \r
952 AliGenEventHeader *header = mcEvent->GenEventHeader();\r
953 if(!header) return;\r
954 \r
955 TArrayF gVertexArray;\r
956 header->PrimaryVertex(gVertexArray);\r
957 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
958 //gVertexArray.At(0),\r
959 //gVertexArray.At(1),\r
960 //gVertexArray.At(2));\r
961 fHistEventStats->Fill(3); //events with a proper vertex\r
962 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
963 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
964 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
965 fHistEventStats->Fill(4); //analayzed events\r
966 fHistVx->Fill(gVertexArray.At(0));\r
967 fHistVy->Fill(gVertexArray.At(1));\r
968 fHistVz->Fill(gVertexArray.At(2));\r
969 \r
970 Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());\r
971 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {\r
972 AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));\r
973 if (!track) {\r
974 Printf("ERROR: Could not receive particle %d", iTracks);\r
975 continue;\r
976 }\r
977 \r
d9f855f3 978 //exclude non stable particles\r
979 if(!mcEvent->IsPhysicalPrimary(iTracks)) continue;\r
980\r
74b04dd4 981 v_eta = track->Eta();\r
982 v_pt = track->Pt();\r
983 \r
984 if( v_pt < fPtMin || v_pt > fPtMax) \r
6e28a437 985 continue;\r
74b04dd4 986 if( v_eta < fEtaMin || v_eta > fEtaMax) \r
6e28a437 987 continue;\r
988 \r
d9f855f3 989 //analyze one set of particles\r
990 if(fUseMCPdgCode) {\r
991 TParticle *particle = track->Particle();\r
992 if(!particle) continue;\r
993 \r
994 Int_t gPdgCode = particle->GetPdgCode();\r
995 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
996 continue;\r
997 }\r
74b04dd4 998 \r
a1390237 999 //Use the acceptance parameterization\r
6fbd0608 1000 if(fAcceptanceParameterization) {\r
1001 Double_t gRandomNumber = gRandom->Rndm();\r
1002 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1003 continue;\r
1004 }\r
74b04dd4 1005 \r
a1390237 1006 //Exclude resonances\r
1007 if(fExcludeResonancesInMC) {\r
1008 TParticle *particle = track->Particle();\r
1009 if(!particle) continue;\r
1010 \r
1011 Bool_t kExcludeParticle = kFALSE;\r
1012 Int_t gMotherIndex = particle->GetFirstMother();\r
1013 if(gMotherIndex != -1) {\r
1014 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1015 if(motherTrack) {\r
1016 TParticle *motherParticle = motherTrack->Particle();\r
1017 if(motherParticle) {\r
1018 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
25cd739b 1019 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1020 if(pdgCodeOfMother == 113) {\r
a1390237 1021 kExcludeParticle = kTRUE;\r
1022 }\r
1023 }\r
1024 }\r
1025 }\r
1026 \r
1027 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1028 if(kExcludeParticle) continue;\r
1029 }\r
74b04dd4 1030\r
1031 v_charge = track->Charge();\r
1032 v_y = track->Y();\r
4a1e5d10 1033 v_phi = track->Phi();\r
74b04dd4 1034 v_E = track->E();\r
1035 track->PxPyPz(v_p);\r
4a1e5d10 1036 //Printf("phi (before): %lf",v_phi);\r
74b04dd4 1037\r
7ef516b2 1038 fHistPt->Fill(v_pt);\r
1039 fHistEta->Fill(v_eta);\r
1040 fHistPhi->Fill(v_phi);\r
1041\r
985da7aa 1042 //Flow after burner\r
1043 if(fUseFlowAfterBurner) {\r
1044 Double_t precisionPhi = 0.001;\r
1045 Int_t maxNumberOfIterations = 100;\r
1046\r
1047 Double_t phi0 = v_phi;\r
1048 Double_t gV2 = fDifferentialV2->Eval(v_pt);\r
1049\r
1050 for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
1051 Double_t phiprev = v_phi;\r
1052 Double_t fl = v_phi - phi0 + gV2*TMath::Sin(2.*(v_phi - gReactionPlane));\r
1053 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(v_phi - gReactionPlane)); \r
1054 v_phi -= fl/fp;\r
1055 if (TMath::AreEqualAbs(phiprev,v_phi,precisionPhi)) break;\r
1056 }\r
4a1e5d10 1057 //Printf("phi (after): %lf\n",v_phi);\r
3c3bbd53 1058 Double_t v_DeltaphiBefore = phi0 - gReactionPlane;\r
1059 if(v_DeltaphiBefore < 0) v_DeltaphiBefore += 2*TMath::Pi();\r
1060 fHistPhiBefore->Fill(v_DeltaphiBefore);\r
1061\r
1062 Double_t v_DeltaphiAfter = v_phi - gReactionPlane;\r
1063 if(v_DeltaphiAfter < 0) v_DeltaphiAfter += 2*TMath::Pi();\r
1064 fHistPhiAfter->Fill(v_DeltaphiAfter);\r
985da7aa 1065 }\r
4a1e5d10 1066 \r
1067 v_phi *= TMath::RadToDeg();\r
985da7aa 1068\r
6e28a437 1069 // fill charge vector\r
74b04dd4 1070 chargeVector[0]->push_back(v_charge);\r
1071 chargeVector[1]->push_back(v_y);\r
1072 chargeVector[2]->push_back(v_eta);\r
1073 chargeVector[3]->push_back(v_phi);\r
1074 chargeVector[4]->push_back(v_p[0]);\r
1075 chargeVector[5]->push_back(v_p[1]);\r
1076 chargeVector[6]->push_back(v_p[2]);\r
1077 chargeVector[7]->push_back(v_pt);\r
1078 chargeVector[8]->push_back(v_E);\r
1079 \r
6e28a437 1080 if(fRunShuffling) {\r
74b04dd4 1081 chargeVectorShuffle[0]->push_back(v_charge);\r
1082 chargeVectorShuffle[1]->push_back(v_y);\r
1083 chargeVectorShuffle[2]->push_back(v_eta);\r
1084 chargeVectorShuffle[3]->push_back(v_phi);\r
1085 chargeVectorShuffle[4]->push_back(v_p[0]);\r
1086 chargeVectorShuffle[5]->push_back(v_p[1]);\r
1087 chargeVectorShuffle[6]->push_back(v_p[2]);\r
1088 chargeVectorShuffle[7]->push_back(v_pt);\r
1089 chargeVectorShuffle[8]->push_back(v_E);\r
6e28a437 1090 }\r
63d15f4b 1091 gNumberOfAcceptedTracks += 1;\r
1092 \r
6e28a437 1093 } //track loop\r
1094 }//Vz cut\r
1095 }//Vy cut\r
1096 }//Vx cut\r
3feee083 1097 }//MC analysis\r
1098 \r
d9f855f3 1099 //multiplicity cut (used in pp)\r
a5ab9ca3 1100 if(fUseMultiplicity) {\r
d9f855f3 1101 if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
6e28a437 1102 return;\r
a5ab9ca3 1103 }\r
d9f855f3 1104 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks);\r
6e28a437 1105 \r
10a747d2 1106 // calculate balance function\r
63d15f4b 1107 if(fUseMultiplicity) \r
1108 fBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVector);\r
1109 else \r
1110 fBalance->CalculateBalance(fCentrality,chargeVector);\r
1900bf4e 1111\r
10a747d2 1112 if(fRunShuffling) {\r
1113 // shuffle charges\r
74b04dd4 1114 random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
63d15f4b 1115 if(fUseMultiplicity) \r
1116 fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVectorShuffle);\r
1117 else \r
1118 fShuffledBalance->CalculateBalance(fCentrality,chargeVectorShuffle);\r
10a747d2 1119 }\r
3feee083 1120} \r
1121\r
1122//________________________________________________________________________\r
1123void AliAnalysisTaskBF::FinishTaskOutput(){\r
1124 //Printf("END BF");\r
1125\r
3feee083 1126 if (!fBalance) {\r
1127 Printf("ERROR: fBalance not available");\r
1128 return;\r
1900bf4e 1129 } \r
10a747d2 1130 if(fRunShuffling) {\r
1131 if (!fShuffledBalance) {\r
1132 Printf("ERROR: fShuffledBalance not available");\r
1133 return;\r
1134 }\r
3feee083 1135 }\r
3feee083 1136\r
3feee083 1137}\r
1138\r
1139//________________________________________________________________________\r
1140void AliAnalysisTaskBF::Terminate(Option_t *) {\r
1141 // Draw result to the screen\r
1142 // Called once at the end of the query\r
1143\r
1144 // not implemented ...\r
6432ac6a 1145\r
3feee083 1146}\r