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