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