4 #include "TLorentzVector.h"
\r
5 #include "TGraphErrors.h"
\r
11 #include "TArrayF.h"
\r
13 #include "TRandom.h"
\r
15 #include "AliAnalysisTaskSE.h"
\r
16 #include "AliAnalysisManager.h"
\r
18 #include "AliESDVertex.h"
\r
19 #include "AliESDEvent.h"
\r
20 #include "AliESDInputHandler.h"
\r
21 #include "AliAODEvent.h"
\r
22 #include "AliAODTrack.h"
\r
23 #include "AliAODInputHandler.h"
\r
24 #include "AliCollisionGeometry.h"
\r
25 #include "AliGenEventHeader.h"
\r
26 #include "AliMCEventHandler.h"
\r
27 #include "AliMCEvent.h"
\r
28 #include "AliStack.h"
\r
29 #include "AliESDtrackCuts.h"
\r
30 #include "AliEventplane.h"
\r
31 #include "AliTHn.h"
\r
33 #include "AliEventPoolManager.h"
\r
35 #include "AliPID.h"
\r
36 #include "AliPIDResponse.h"
\r
37 #include "AliPIDCombined.h"
\r
39 #include "AliAnalysisTaskBFPsi.h"
\r
40 #include "AliBalancePsi.h"
\r
41 #include "AliAnalysisTaskTriggeredBF.h"
\r
44 // Analysis task for the BF vs Psi code
\r
45 // Authors: Panos.Christakoglou@nikhef.nl
\r
50 ClassImp(AliAnalysisTaskBFPsi)
\r
52 //________________________________________________________________________
\r
53 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
\r
54 : AliAnalysisTaskSE(name),
\r
56 fRunShuffling(kFALSE),
\r
57 fShuffledBalance(0),
\r
59 fRunMixingEventPlane(kFALSE),
\r
60 fMixingTracks(50000),
\r
70 fHistTriggerStats(0),
\r
91 fHistdEdxVsPTPCbeforePID(NULL),
\r
92 fHistBetavsPTOFbeforePID(NULL),
\r
93 fHistProbTPCvsPtbeforePID(NULL),
\r
94 fHistProbTOFvsPtbeforePID(NULL),
\r
95 fHistProbTPCTOFvsPtbeforePID(NULL),
\r
96 fHistNSigmaTPCvsPtbeforePID(NULL),
\r
97 fHistNSigmaTOFvsPtbeforePID(NULL),
\r
98 fHistdEdxVsPTPCafterPID(NULL),
\r
99 fHistBetavsPTOFafterPID(NULL),
\r
100 fHistProbTPCvsPtafterPID(NULL),
\r
101 fHistProbTOFvsPtafterPID(NULL),
\r
102 fHistProbTPCTOFvsPtafterPID(NULL),
\r
103 fHistNSigmaTPCvsPtafterPID(NULL),
\r
104 fHistNSigmaTOFvsPtafterPID(NULL),
\r
107 fParticleOfInterest(kPion),
\r
108 fPidDetectorConfig(kTPCTOF),
\r
110 fUsePIDnSigma(kTRUE),
\r
111 fUsePIDPropabilities(kFALSE),
\r
113 fMinAcceptedPIDProbability(0.8),
\r
115 fCentralityEstimator("V0M"),
\r
116 fUseCentrality(kFALSE),
\r
117 fCentralityPercentileMin(0.),
\r
118 fCentralityPercentileMax(5.),
\r
119 fImpactParameterMin(0.),
\r
120 fImpactParameterMax(20.),
\r
121 fUseMultiplicity(kFALSE),
\r
122 fNumberOfAcceptedTracksMin(0),
\r
123 fNumberOfAcceptedTracksMax(10000),
\r
124 fHistNumberOfAcceptedTracks(0),
\r
125 fUseOfflineTrigger(kFALSE),
\r
129 nAODtrackCutBit(128),
\r
132 fPtMinForCorrections(0.3),//=================================correction
\r
133 fPtMaxForCorrections(1.5),//=================================correction
\r
134 fPtBinForCorrections(36), //=================================correction
\r
137 fEtaMinForCorrections(-0.8),//=================================correction
\r
138 fEtaMaxForCorrections(0.8),//=================================correction
\r
139 fEtaBinForCorrections(16), //=================================correction
\r
142 fPhiMinForCorrections(0.),//=================================correction
\r
143 fPhiMaxForCorrections(360.),//=================================correction
\r
144 fPhiBinForCorrections(100), //=================================correction
\r
148 fNClustersTPCCut(-1),
\r
149 fAcceptanceParameterization(0),
\r
150 fDifferentialV2(0),
\r
151 fUseFlowAfterBurner(kFALSE),
\r
152 fExcludeResonancesInMC(kFALSE),
\r
153 fUseMCPdgCode(kFALSE),
\r
154 fPDGCodeToBeAnalyzed(-1),
\r
155 fEventClass("EventPlane")
\r
158 // Define input and output slots here
\r
159 // Input slot #0 works with a TChain
\r
161 //======================================================correction
\r
162 for (Int_t i=0; i<kCENTRALITY; i++){
\r
163 fHistMatrixCorrectionPlus[i] = NULL;
\r
164 fHistMatrixCorrectionMinus[i] = NULL;
\r
166 //=====================================================correction
\r
168 DefineInput(0, TChain::Class());
\r
169 // Output slot #0 writes into a TH1 container
\r
170 DefineOutput(1, TList::Class());
\r
171 DefineOutput(2, TList::Class());
\r
172 DefineOutput(3, TList::Class());
\r
173 DefineOutput(4, TList::Class());
\r
174 DefineOutput(5, TList::Class());
\r
177 //________________________________________________________________________
\r
178 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
\r
180 // delete fBalance;
\r
181 // delete fShuffledBalance;
\r
183 // delete fListBF;
\r
184 // delete fListBFS;
\r
186 // delete fHistEventStats;
\r
187 // delete fHistTrackStats;
\r
188 // delete fHistVx;
\r
189 // delete fHistVy;
\r
190 // delete fHistVz;
\r
192 // delete fHistClus;
\r
193 // delete fHistDCA;
\r
194 // delete fHistChi2;
\r
196 // delete fHistEta;
\r
197 // delete fHistPhi;
\r
198 // delete fHistEtaPhiPos;
\r
199 // delete fHistEtaPhiNeg;
\r
200 // delete fHistV0M;
\r
203 //________________________________________________________________________
\r
204 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
\r
205 // Create histograms
\r
208 // global switch disabling the reference
\r
209 // (to avoid "Replacing existing TH1" if several wagons are created in train)
\r
210 Bool_t oldStatus = TH1::AddDirectoryStatus();
\r
211 TH1::AddDirectory(kFALSE);
\r
214 fBalance = new AliBalancePsi();
\r
215 fBalance->SetAnalysisLevel("ESD");
\r
216 fBalance->SetEventClass(fEventClass);
\r
217 //fBalance->SetNumberOfBins(-1,16);
\r
218 //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
\r
220 if(fRunShuffling) {
\r
221 if(!fShuffledBalance) {
\r
222 fShuffledBalance = new AliBalancePsi();
\r
223 fShuffledBalance->SetAnalysisLevel("ESD");
\r
224 fShuffledBalance->SetEventClass(fEventClass);
\r
225 //fShuffledBalance->SetNumberOfBins(-1,16);
\r
226 //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
\r
230 if(!fMixedBalance) {
\r
231 fMixedBalance = new AliBalancePsi();
\r
232 fMixedBalance->SetAnalysisLevel("ESD");
\r
233 fMixedBalance->SetEventClass(fEventClass);
\r
238 fList = new TList();
\r
239 fList->SetName("listQA");
\r
242 //Balance Function list
\r
243 fListBF = new TList();
\r
244 fListBF->SetName("listBF");
\r
245 fListBF->SetOwner();
\r
247 if(fRunShuffling) {
\r
248 fListBFS = new TList();
\r
249 fListBFS->SetName("listBFShuffled");
\r
250 fListBFS->SetOwner();
\r
254 fListBFM = new TList();
\r
255 fListBFM->SetName("listTriggeredBFMixed");
\r
256 fListBFM->SetOwner();
\r
261 fHistListPIDQA = new TList();
\r
262 fHistListPIDQA->SetName("listQAPID");
\r
263 fHistListPIDQA->SetOwner();
\r
267 TString gCutName[5] = {"Total","Offline trigger",
\r
268 "Vertex","Analyzed","sel. Centrality"};
\r
269 fHistEventStats = new TH2F("fHistEventStats",
\r
270 "Event statistics;;Centrality percentile;N_{events}",
\r
271 5,0.5,5.5,220,-5,105);
\r
272 for(Int_t i = 1; i <= 5; i++)
\r
273 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
\r
274 fList->Add(fHistEventStats);
\r
276 TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
\r
277 fHistCentStats = new TH2F("fHistCentStats",
\r
278 "Centrality statistics;;Cent percentile",
\r
279 9,-0.5,8.5,220,-5,105);
\r
280 for(Int_t i = 1; i <= 9; i++)
\r
281 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
\r
282 fList->Add(fHistCentStats);
\r
284 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
\r
285 fList->Add(fHistTriggerStats);
\r
287 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
\r
288 fList->Add(fHistTrackStats);
\r
290 fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
\r
291 fList->Add(fHistNumberOfAcceptedTracks);
\r
293 // Vertex distributions
\r
294 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
\r
295 fList->Add(fHistVx);
\r
296 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
\r
297 fList->Add(fHistVy);
\r
298 fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
\r
299 fList->Add(fHistVz);
\r
302 fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
\r
303 fList->Add(fHistEventPlane);
\r
306 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
\r
307 fList->Add(fHistClus);
\r
308 fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);
\r
309 fList->Add(fHistChi2);
\r
310 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
\r
311 fList->Add(fHistDCA);
\r
312 fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);
\r
313 fList->Add(fHistPt);
\r
314 fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);
\r
315 fList->Add(fHistEta);
\r
316 fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);
\r
317 fList->Add(fHistRapidity);
\r
318 fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);
\r
319 fList->Add(fHistPhi);
\r
320 fHistEtaPhiPos = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);
\r
321 fList->Add(fHistEtaPhiPos);
\r
322 fHistEtaPhiNeg = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);
\r
323 fList->Add(fHistEtaPhiNeg);
\r
324 fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
\r
325 fList->Add(fHistPhiBefore);
\r
326 fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
\r
327 fList->Add(fHistPhiAfter);
\r
328 fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
\r
329 fList->Add(fHistPhiPos);
\r
330 fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);
\r
331 fList->Add(fHistPhiNeg);
\r
332 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
\r
333 fList->Add(fHistV0M);
\r
334 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
\r
335 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
\r
336 for(Int_t i = 1; i <= 6; i++)
\r
337 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
\r
338 fList->Add(fHistRefTracks);
\r
340 // QA histograms for different cuts
\r
341 fList->Add(fBalance->GetQAHistHBTbefore());
\r
342 fList->Add(fBalance->GetQAHistHBTafter());
\r
343 fList->Add(fBalance->GetQAHistConversionbefore());
\r
344 fList->Add(fBalance->GetQAHistConversionafter());
\r
345 fList->Add(fBalance->GetQAHistPsiMinusPhi());
\r
346 fList->Add(fBalance->GetQAHistResonancesBefore());
\r
347 fList->Add(fBalance->GetQAHistResonancesRho());
\r
348 fList->Add(fBalance->GetQAHistResonancesK0());
\r
349 fList->Add(fBalance->GetQAHistResonancesLambda());
\r
350 fList->Add(fBalance->GetQAHistQbefore());
\r
351 fList->Add(fBalance->GetQAHistQafter());
\r
353 // Balance function histograms
\r
354 // Initialize histograms if not done yet
\r
355 if(!fBalance->GetHistNp()){
\r
356 AliWarning("Histograms not yet initialized! --> Will be done now");
\r
357 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
\r
358 fBalance->InitHistograms();
\r
361 if(fRunShuffling) {
\r
362 if(!fShuffledBalance->GetHistNp()) {
\r
363 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
\r
364 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
\r
365 fShuffledBalance->InitHistograms();
\r
369 //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
\r
370 fListBF->Add(fBalance->GetHistNp());
\r
371 fListBF->Add(fBalance->GetHistNn());
\r
372 fListBF->Add(fBalance->GetHistNpn());
\r
373 fListBF->Add(fBalance->GetHistNnn());
\r
374 fListBF->Add(fBalance->GetHistNpp());
\r
375 fListBF->Add(fBalance->GetHistNnp());
\r
377 if(fRunShuffling) {
\r
378 fListBFS->Add(fShuffledBalance->GetHistNp());
\r
379 fListBFS->Add(fShuffledBalance->GetHistNn());
\r
380 fListBFS->Add(fShuffledBalance->GetHistNpn());
\r
381 fListBFS->Add(fShuffledBalance->GetHistNnn());
\r
382 fListBFS->Add(fShuffledBalance->GetHistNpp());
\r
383 fListBFS->Add(fShuffledBalance->GetHistNnp());
\r
387 fListBFM->Add(fMixedBalance->GetHistNp());
\r
388 fListBFM->Add(fMixedBalance->GetHistNn());
\r
389 fListBFM->Add(fMixedBalance->GetHistNpn());
\r
390 fListBFM->Add(fMixedBalance->GetHistNnn());
\r
391 fListBFM->Add(fMixedBalance->GetHistNpp());
\r
392 fListBFM->Add(fMixedBalance->GetHistNnp());
\r
399 Int_t trackDepth = fMixingTracks;
\r
400 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
\r
403 Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
\r
404 Double_t* centbins = centralityBins;
\r
405 Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1;
\r
407 // multiplicity bins
\r
408 Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
\r
409 Double_t* multbins = multiplicityBins;
\r
410 Int_t nMultiplicityBins = sizeof(multiplicityBins) / sizeof(Double_t) - 1;
\r
413 Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
\r
414 Double_t* vtxbins = vertexBins;
\r
415 Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1;
\r
417 // Event plane angle (Psi) bins
\r
418 Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
\r
419 Double_t* psibins = psiBins;
\r
420 Int_t nPsiBins = sizeof(psiBins) / sizeof(Double_t) - 1;
\r
422 // run the event mixing also in bins of event plane (statistics!)
\r
423 if(fRunMixingEventPlane){
\r
424 if(fEventClass=="Multiplicity"){
\r
425 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
\r
428 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
\r
432 if(fEventClass=="Multiplicity"){
\r
433 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
\r
436 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
\r
441 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
\r
443 //====================PID========================//
\r
445 fPIDCombined = new AliPIDCombined();
\r
446 fPIDCombined->SetDefaultTPCPriors();
\r
448 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
\r
449 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition
\r
451 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
\r
452 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition
\r
454 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
\r
455 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition
\r
457 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
\r
458 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition
\r
460 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
\r
461 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition
\r
463 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
\r
464 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition
\r
466 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
\r
467 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition
\r
469 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
\r
470 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition
\r
472 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
\r
473 fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition
\r
475 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
\r
476 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition
\r
478 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
\r
479 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition
\r
481 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
\r
482 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition
\r
484 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
\r
485 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition
\r
487 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
\r
488 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition
\r
490 //====================PID========================//
\r
492 // Post output data.
\r
493 PostData(1, fList);
\r
494 PostData(2, fListBF);
\r
495 if(fRunShuffling) PostData(3, fListBFS);
\r
496 if(fRunMixing) PostData(4, fListBFM);
\r
497 if(fUsePID) PostData(5, fHistListPIDQA); //PID
\r
499 TH1::AddDirectory(oldStatus);
\r
503 //________________________________________________________________________
\r
504 void AliAnalysisTaskBFPsi::SetInputCorrection(const char* filename,
\r
505 const char* gCollSystem) {
\r
506 //Open files that will be used for correction
\r
507 TString gCollidingSystem = gCollSystem;
\r
509 //Open the input file
\r
510 TFile *f = TFile::Open(filename);
\r
512 Printf("File not found!!!");
\r
516 //Get the TDirectoryFile and TList
\r
517 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionEfficiencyAnalysis"));
\r
519 Printf("TDirectoryFile not found!!!");
\r
523 TString listEffName = "";
\r
524 for (Int_t iCent = 0; iCent < kCENTRALITY; iCent++) {
\r
525 listEffName = "listEfficiencyBF_";
\r
526 listEffName += (TString)((Int_t)(centralityArrayForPbPb[iCent]));
\r
527 listEffName += "_";
\r
528 listEffName += (TString)((Int_t)(centralityArrayForPbPb[iCent+1]));
\r
529 TList *list = (TList*)dir->Get(listEffName.Data());
\r
531 Printf("TList Efficiency not found!!!");
\r
535 TString histoName = "fHistMatrixCorrectionPlus";
\r
536 fHistMatrixCorrectionPlus[iCent]= dynamic_cast<TH3D *>(list->FindObject(histoName.Data()));
\r
537 if(!fHistMatrixCorrectionPlus[iCent]) {
\r
538 Printf("fHist not found!!!");
\r
542 histoName = "fHistMatrixCorrectionMinus";
\r
543 fHistMatrixCorrectionMinus[iCent] = dynamic_cast<TH3D *>(list->FindObject(histoName.Data()));
\r
544 if(!fHistMatrixCorrectionMinus[iCent]) {
\r
545 Printf("fHist not found!!!");
\r
548 }//loop over centralities: ONLY the PbPb case is covered
\r
550 if(fHistMatrixCorrectionPlus[0]){
\r
551 fEtaMinForCorrections = fHistMatrixCorrectionPlus[0]->GetXaxis()->GetXmin();
\r
552 fEtaMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetXaxis()->GetXmax();
\r
553 fEtaBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsX();
\r
555 fPtMinForCorrections = fHistMatrixCorrectionPlus[0]->GetYaxis()->GetXmin();
\r
556 fPtMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetYaxis()->GetXmax();
\r
557 fPtBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsY();
\r
559 fPhiMinForCorrections = fHistMatrixCorrectionPlus[0]->GetZaxis()->GetXmin();
\r
560 fPhiMaxForCorrections = fHistMatrixCorrectionPlus[0]->GetZaxis()->GetXmax();
\r
561 fPhiBinForCorrections = fHistMatrixCorrectionPlus[0]->GetNbinsZ();
\r
565 //________________________________________________________________________
\r
566 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
\r
568 // Called for each event
\r
570 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
571 Int_t gNumberOfAcceptedTracks = 0;
\r
572 Double_t gCentrality = -1.;
\r
573 Double_t gReactionPlane = -1.;
\r
574 Float_t bSign = 0.;
\r
576 // get the event (for generator level: MCEvent())
\r
577 AliVEvent* eventMain = NULL;
\r
578 if(gAnalysisLevel == "MC") {
\r
579 eventMain = dynamic_cast<AliVEvent*>(MCEvent());
\r
582 eventMain = dynamic_cast<AliVEvent*>(InputEvent());
\r
584 // for HBT like cuts need magnetic field sign
\r
585 bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
\r
588 AliError("eventMain not available");
\r
592 // PID Response task active?
\r
594 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
\r
595 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
\r
598 // check event cuts and fill event histograms
\r
599 if((gCentrality = IsEventAccepted(eventMain)) < 0){
\r
603 //Compute Multiplicity or Centrality variable
\r
604 Double_t lMultiplicityVar = GetRefMultiOrCentrality( eventMain );
\r
606 // get the reaction plane
\r
607 gReactionPlane = GetEventPlane(eventMain);
\r
608 fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
\r
609 if(gReactionPlane < 0){
\r
613 // get the accepted tracks in main event
\r
614 TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
\r
615 gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
\r
617 //multiplicity cut (used in pp)
\r
618 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
\r
619 if(fUseMultiplicity) {
\r
620 if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))
\r
624 // store charges of all accepted tracks, shuffle and reassign (two extra loops!)
\r
625 TObjArray* tracksShuffled = NULL;
\r
627 tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
\r
633 // 1. First get an event pool corresponding in mult (cent) and
\r
634 // zvertex to the current event. Once initialized, the pool
\r
635 // should contain nMix (reduced) events. This routine does not
\r
636 // pre-scan the chain. The first several events of every chain
\r
637 // will be skipped until the needed pools are filled to the
\r
638 // specified depth. If the pool categories are not too rare, this
\r
639 // should not be a problem. If they are rare, you could lose`
\r
642 // 2. Collect the whole pool's content of tracks into one TObjArray
\r
643 // (bgTracks), which is effectively a single background super-event.
\r
645 // 3. The reduced and bgTracks arrays must both be passed into
\r
646 // FillCorrelations(). Also nMix should be passed in, so a weight
\r
647 // of 1./nMix can be applied.
\r
649 AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
\r
652 AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
\r
656 //pool->SetDebug(1);
\r
658 if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){
\r
661 Int_t nMix = pool->GetCurrentNEvents();
\r
662 //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
\r
664 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
\r
665 //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
\r
666 //if (pool->IsReady())
\r
667 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
\r
669 // Fill mixed-event histos here
\r
670 for (Int_t jMix=0; jMix<nMix; jMix++)
\r
672 TObjArray* tracksMixed = pool->GetEvent(jMix);
\r
673 fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
\r
677 // Update the Event pool
\r
678 pool->UpdatePool(tracksMain);
\r
679 //pool->PrintInfo();
\r
681 }//pool NULL check
\r
684 // calculate balance function
\r
685 fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
\r
687 // calculate shuffled balance function
\r
688 if(fRunShuffling && tracksShuffled != NULL) {
\r
689 fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
\r
693 //________________________________________________________________________
\r
694 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
\r
695 // Checks the Event cuts
\r
696 // Fills Event statistics histograms
\r
698 Bool_t isSelectedMain = kTRUE;
\r
699 Float_t gCentrality = -1.;
\r
700 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
702 fHistEventStats->Fill(1,gCentrality); //all events
\r
704 // Event trigger bits
\r
705 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
\r
706 if(fUseOfflineTrigger)
\r
707 isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
\r
709 if(isSelectedMain) {
\r
710 fHistEventStats->Fill(2,gCentrality); //triggered events
\r
712 //Centrality stuff
\r
713 if(fUseCentrality) {
\r
714 if(gAnalysisLevel == "AOD") { //centrality in AOD header
\r
715 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
\r
717 gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
\r
719 // QA for centrality estimators
\r
720 fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
\r
721 fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
\r
722 fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
\r
723 fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("TKL"));
\r
724 fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
\r
725 fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
\r
726 fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
\r
727 fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
\r
728 fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
\r
730 // centrality QA (V0M)
\r
731 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
\r
733 // centrality QA (reference tracks)
\r
734 fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
\r
735 fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
\r
736 fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
\r
737 fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
\r
738 fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
\r
739 fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
\r
740 fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
\r
741 fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
\r
742 fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
\r
746 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
\r
747 AliCentrality *centrality = event->GetCentrality();
\r
748 gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
\r
750 // QA for centrality estimators
\r
751 fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));
\r
752 fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("FMD"));
\r
753 fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("TRK"));
\r
754 fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("TKL"));
\r
755 fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("CL0"));
\r
756 fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("CL1"));
\r
757 fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("V0MvsFMD"));
\r
758 fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("TKLvsV0M"));
\r
759 fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZEMvsZDC"));
\r
761 // centrality QA (V0M)
\r
762 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
\r
764 else if(gAnalysisLevel == "MC"){
\r
765 Double_t gImpactParameter = 0.;
\r
766 if(dynamic_cast<AliMCEvent*>(event)){
\r
767 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());
\r
769 gImpactParameter = headerH->ImpactParameter();
\r
770 gCentrality = gImpactParameter;
\r
780 if(gAnalysisLevel == "MC"){
\r
782 AliError("mcEvent not available");
\r
786 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
788 AliGenEventHeader *header = gMCEvent->GenEventHeader();
\r
790 TArrayF gVertexArray;
\r
791 header->PrimaryVertex(gVertexArray);
\r
792 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
\r
793 //gVertexArray.At(0),
\r
794 //gVertexArray.At(1),
\r
795 //gVertexArray.At(2));
\r
796 fHistEventStats->Fill(3,gCentrality); //events with a proper vertex
\r
797 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
\r
798 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
\r
799 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
\r
800 fHistEventStats->Fill(4,gCentrality); //analayzed events
\r
801 fHistVx->Fill(gVertexArray.At(0));
\r
802 fHistVy->Fill(gVertexArray.At(1));
\r
803 fHistVz->Fill(gVertexArray.At(2),gCentrality);
\r
805 // take only events inside centrality class
\r
806 if((fImpactParameterMin < gCentrality) && (fImpactParameterMax > gCentrality)){
\r
807 fHistEventStats->Fill(5,gCentrality); //events with correct centrality
\r
808 return gCentrality;
\r
809 }//centrality class
\r
817 // Event Vertex AOD, ESD, ESDMC
\r
819 const AliVVertex *vertex = event->GetPrimaryVertex();
\r
822 Double32_t fCov[6];
\r
823 vertex->GetCovarianceMatrix(fCov);
\r
824 if(vertex->GetNContributors() > 0) {
\r
826 fHistEventStats->Fill(3,gCentrality); //events with a proper vertex
\r
827 if(TMath::Abs(vertex->GetX()) < fVxMax) {
\r
828 if(TMath::Abs(vertex->GetY()) < fVyMax) {
\r
829 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
\r
830 fHistEventStats->Fill(4,gCentrality); //analyzed events
\r
831 fHistVx->Fill(vertex->GetX());
\r
832 fHistVy->Fill(vertex->GetY());
\r
833 fHistVz->Fill(vertex->GetZ(),gCentrality);
\r
835 // take only events inside centrality class
\r
836 if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){
\r
837 fHistEventStats->Fill(5,gCentrality); //events with correct centrality
\r
838 return gCentrality;
\r
839 }//centrality class
\r
843 }//proper vertex resolution
\r
844 }//proper number of contributors
\r
845 }//vertex object valid
\r
846 }//triggered event
\r
849 // in all other cases return -1 (event not accepted)
\r
854 //________________________________________________________________________
\r
855 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
\r
856 // Checks the Event cuts
\r
857 // Fills Event statistics histograms
\r
859 Float_t gCentrality = -1.;
\r
860 Double_t fMultiplicity = -100.;
\r
861 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
862 if(fEventClass == "Centrality"){
\r
865 if(gAnalysisLevel == "AOD") { //centrality in AOD header
\r
866 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
\r
868 gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
\r
872 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
\r
873 AliCentrality *centrality = event->GetCentrality();
\r
874 gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
\r
876 else if(gAnalysisLevel == "MC"){
\r
877 Double_t gImpactParameter = 0.;
\r
878 AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
880 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
\r
882 gImpactParameter = headerH->ImpactParameter();
\r
883 gCentrality = gImpactParameter;
\r
890 }//End if "Centrality"
\r
891 if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){
\r
892 AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
\r
894 fMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
\r
895 }//AliESDevent cast
\r
897 if(fEventClass=="Multiplicity"&&gAnalysisLevel != "ESD"){
\r
898 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
\r
900 fMultiplicity = header->GetRefMultiplicity();
\r
903 Double_t lReturnVal = -100;
\r
904 if(fEventClass=="Multiplicity"){
\r
905 lReturnVal = fMultiplicity;
\r
906 }else if(fEventClass=="Centrality"){
\r
907 lReturnVal = gCentrality;
\r
912 //________________________________________________________________________
\r
913 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
\r
914 // Get the event plane
\r
916 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
918 Float_t gVZEROEventPlane = -10.;
\r
919 Float_t gReactionPlane = -10.;
\r
920 Double_t qxTot = 0.0, qyTot = 0.0;
\r
922 //MC: from reaction plane
\r
923 if(gAnalysisLevel == "MC"){
\r
925 AliError("mcEvent not available");
\r
929 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
931 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
\r
933 gReactionPlane = headerH->ReactionPlaneAngle();
\r
934 //gReactionPlane *= TMath::RadToDeg();
\r
939 // AOD,ESD,ESDMC: from VZERO Event Plane
\r
942 AliEventplane *ep = event->GetEventplane();
\r
944 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
\r
945 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
\r
946 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
\r
947 gReactionPlane = gVZEROEventPlane;
\r
951 return gReactionPlane;
\r
954 //________________________________________________________________________
\r
955 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta,
\r
959 Double_t gCentrality) {
\r
960 // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality)
\r
962 Double_t correction = 1.;
\r
963 Int_t binEta = 0, binPt = 0, binPhi = 0;
\r
965 //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);
\r
966 if(fEtaBinForCorrections != 0) {
\r
967 Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;
\r
968 if(fEtaMaxForCorrections != fEtaMinForCorrections)
\r
969 binEta = (Int_t)(vEta/widthEta)+1;
\r
972 if(fPtBinForCorrections != 0) {
\r
973 Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;
\r
974 if(fPtMaxForCorrections != fPtMinForCorrections)
\r
975 binPt = (Int_t)(vPt/widthPt) + 1;
\r
978 if(fPhiBinForCorrections != 0) {
\r
979 Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;
\r
980 if(fPhiMaxForCorrections != fPhiMinForCorrections)
\r
981 binPhi = (Int_t)(vPhi/widthPhi)+ 1;
\r
984 Int_t gCentralityInt = 1;
\r
985 for (Int_t i=0; i<kCENTRALITY; i++){
\r
986 if((centralityArrayForPbPb[i] <= gCentrality)&&(gCentrality <= centralityArrayForPbPb[i+1]))
\r
987 gCentralityInt = i;
\r
990 if(fHistMatrixCorrectionPlus[gCentralityInt]){
\r
992 correction = fHistMatrixCorrectionPlus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));
\r
995 correction = fHistMatrixCorrectionMinus[gCentralityInt]->GetBinContent(fHistMatrixCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));
\r
1001 if (correction == 0.) {
\r
1002 AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt));
\r
1006 return correction;
\r
1009 //________________________________________________________________________
\r
1010 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){
\r
1011 // Returns TObjArray with tracks after all track cuts (only for AOD!)
\r
1012 // Fills QA histograms
\r
1014 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
1016 //output TObjArray holding all good tracks
\r
1017 TObjArray* tracksAccepted = new TObjArray;
\r
1018 tracksAccepted->SetOwner(kTRUE);
\r
1027 if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
\r
1028 // Loop over tracks in event
\r
1029 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1030 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
\r
1032 AliError(Form("Could not receive track %d", iTracks));
\r
1038 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
\r
1039 // take only TPC only tracks
\r
1040 //fHistTrackStats->Fill(aodTrack->GetFilterMap());
\r
1041 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
\r
1042 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
\r
1044 if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;
\r
1046 vCharge = aodTrack->Charge();
\r
1047 vEta = aodTrack->Eta();
\r
1048 vY = aodTrack->Y();
\r
1049 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
\r
1050 vPt = aodTrack->Pt();
\r
1052 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1053 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1056 // Kinematics cuts from ESD track cuts
\r
1057 if( vPt < fPtMin || vPt > fPtMax) continue;
\r
1058 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
1060 // Extra DCA cuts (for systematic studies [!= -1])
\r
1061 if( fDCAxyCut != -1 && fDCAzCut != -1){
\r
1062 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
\r
1063 continue; // 2D cut
\r
1067 // Extra TPC cuts (for systematic studies [!= -1])
\r
1068 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
\r
1071 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
\r
1075 // fill QA histograms
\r
1076 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
\r
1077 fHistDCA->Fill(dcaZ,dcaXY);
\r
1078 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
\r
1079 fHistPt->Fill(vPt,gCentrality);
\r
1080 fHistEta->Fill(vEta,gCentrality);
\r
1081 fHistRapidity->Fill(vY,gCentrality);
\r
1082 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
1083 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
1084 fHistPhi->Fill(vPhi,gCentrality);
\r
1085 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
1086 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
1088 //=======================================correction
\r
1089 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
1091 // add the track to the TObjArray
\r
1092 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
1097 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
\r
1099 AliESDtrack *trackTPC = NULL;
\r
1100 AliMCParticle *mcTrack = NULL;
\r
1102 AliMCEvent* mcEvent = NULL;
\r
1104 // for MC ESDs use also MC information
\r
1105 if(gAnalysisLevel == "MCESD"){
\r
1106 mcEvent = MCEvent();
\r
1108 AliError("mcEvent not available");
\r
1109 return tracksAccepted;
\r
1113 // Loop over tracks in event
\r
1114 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1115 AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));
\r
1117 AliError(Form("Could not receive track %d", iTracks));
\r
1121 // for MC ESDs use also MC information --> MC track not used further???
\r
1122 if(gAnalysisLevel == "MCESD"){
\r
1123 Int_t label = TMath::Abs(track->GetLabel());
\r
1124 if(label > mcEvent->GetNumberOfTracks()) continue;
\r
1125 if(label > mcEvent->GetNumberOfPrimaries()) continue;
\r
1127 mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
\r
1128 if(!mcTrack) continue;
\r
1131 // take only TPC only tracks
\r
1132 trackTPC = new AliESDtrack();
\r
1133 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
\r
1136 if(fESDtrackCuts)
\r
1137 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
\r
1139 // fill QA histograms
\r
1142 trackTPC->GetImpactParameters(b,bCov);
\r
1143 if (bCov[0]<=0 || bCov[2]<=0) {
\r
1144 AliDebug(1, "Estimated b resolution lower or equal zero!");
\r
1145 bCov[0]=0; bCov[2]=0;
\r
1148 Int_t nClustersTPC = -1;
\r
1149 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
\r
1150 //nClustersTPC = track->GetTPCclusters(0); // global track
\r
1151 Float_t chi2PerClusterTPC = -1;
\r
1152 if (nClustersTPC!=0) {
\r
1153 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
\r
1154 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
\r
1157 //===========================PID===============================//
\r
1159 Double_t prob[AliPID::kSPECIES]={0.};
\r
1160 Double_t probTPC[AliPID::kSPECIES]={0.};
\r
1161 Double_t probTOF[AliPID::kSPECIES]={0.};
\r
1162 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
\r
1164 Double_t nSigma = 0.;
\r
1165 UInt_t detUsedTPC = 0;
\r
1166 UInt_t detUsedTOF = 0;
\r
1167 UInt_t detUsedTPCTOF = 0;
\r
1169 //Decide what detector configuration we want to use
\r
1170 switch(fPidDetectorConfig) {
\r
1172 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
\r
1173 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
\r
1174 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
\r
1175 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1176 prob[iSpecies] = probTPC[iSpecies];
\r
1179 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
\r
1180 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
\r
1181 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
\r
1182 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1183 prob[iSpecies] = probTOF[iSpecies];
\r
1186 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
\r
1187 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
\r
1188 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1189 prob[iSpecies] = probTPCTOF[iSpecies];
\r
1193 }//end switch: define detector mask
\r
1195 //Filling the PID QA
\r
1196 Double_t tofTime = -999., length = 999., tof = -999.;
\r
1197 Double_t c = TMath::C()*1.E-9;// m/ns
\r
1198 Double_t beta = -999.;
\r
1199 Double_t nSigmaTOFForParticleOfInterest = -999.;
\r
1200 if ( (track->IsOn(AliESDtrack::kTOFin)) &&
\r
1201 (track->IsOn(AliESDtrack::kTIME)) ) {
\r
1202 tofTime = track->GetTOFsignal();//in ps
\r
1203 length = track->GetIntegratedLength();
\r
1204 tof = tofTime*1E-3; // ns
\r
1207 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
\r
1211 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
\r
1215 length = length*0.01; // in meters
\r
1217 beta = length/tof;
\r
1219 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
\r
1220 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
\r
1221 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
\r
1222 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
\r
1226 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
\r
1227 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
\r
1228 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
\r
1229 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
\r
1230 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
\r
1231 //end of QA-before pid
\r
1233 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
\r
1234 //Make the decision based on the n-sigma
\r
1235 if(fUsePIDnSigma) {
\r
1236 if(nSigma > fPIDNSigma) continue;}
\r
1238 //Make the decision based on the bayesian
\r
1239 else if(fUsePIDPropabilities) {
\r
1240 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
\r
1241 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
\r
1244 //Fill QA after the PID
\r
1245 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
\r
1246 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
\r
1247 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
\r
1249 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
\r
1250 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
\r
1251 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
\r
1252 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
\r
1255 //===========================PID===============================//
\r
1256 vCharge = trackTPC->Charge();
\r
1257 vY = trackTPC->Y();
\r
1258 vEta = trackTPC->Eta();
\r
1259 vPhi = trackTPC->Phi();// * TMath::RadToDeg();
\r
1260 vPt = trackTPC->Pt();
\r
1261 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
\r
1262 fHistDCA->Fill(b[1],b[0]);
\r
1263 fHistChi2->Fill(chi2PerClusterTPC,gCentrality);
\r
1264 fHistPt->Fill(vPt,gCentrality);
\r
1265 fHistEta->Fill(vEta,gCentrality);
\r
1266 fHistPhi->Fill(vPhi,gCentrality);
\r
1267 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
1268 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
1269 fHistRapidity->Fill(vY,gCentrality);
\r
1270 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
1271 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
1273 //=======================================correction
\r
1274 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
1276 // add the track to the TObjArray
\r
1277 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
1283 else if(gAnalysisLevel == "MC"){
\r
1285 AliError("mcEvent not available");
\r
1289 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
1291 // Loop over tracks in event
\r
1292 for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
\r
1293 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
\r
1295 AliError(Form("Could not receive particle %d", iTracks));
\r
1299 //exclude non stable particles
\r
1300 if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
\r
1302 vEta = track->Eta();
\r
1303 vPt = track->Pt();
\r
1306 if( vPt < fPtMin || vPt > fPtMax)
\r
1309 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
1311 else if (fUsePID){
\r
1312 if( vY < fEtaMin || vY > fEtaMax) continue;
\r
1315 //analyze one set of particles
\r
1316 if(fUseMCPdgCode) {
\r
1317 TParticle *particle = track->Particle();
\r
1318 if(!particle) continue;
\r
1320 Int_t gPdgCode = particle->GetPdgCode();
\r
1321 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
\r
1325 //Use the acceptance parameterization
\r
1326 if(fAcceptanceParameterization) {
\r
1327 Double_t gRandomNumber = gRandom->Rndm();
\r
1328 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
\r
1332 //Exclude resonances
\r
1333 if(fExcludeResonancesInMC) {
\r
1334 TParticle *particle = track->Particle();
\r
1335 if(!particle) continue;
\r
1337 Bool_t kExcludeParticle = kFALSE;
\r
1338 Int_t gMotherIndex = particle->GetFirstMother();
\r
1339 if(gMotherIndex != -1) {
\r
1340 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
\r
1342 TParticle *motherParticle = motherTrack->Particle();
\r
1343 if(motherParticle) {
\r
1344 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
\r
1345 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
\r
1346 if(pdgCodeOfMother == 113) {
\r
1347 kExcludeParticle = kTRUE;
\r
1353 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
\r
1354 if(kExcludeParticle) continue;
\r
1357 vCharge = track->Charge();
\r
1358 vPhi = track->Phi();
\r
1359 //Printf("phi (before): %lf",vPhi);
\r
1361 fHistPt->Fill(vPt,gCentrality);
\r
1362 fHistEta->Fill(vEta,gCentrality);
\r
1363 fHistPhi->Fill(vPhi,gCentrality);
\r
1364 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
1365 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
1366 //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
\r
1367 fHistRapidity->Fill(vY,gCentrality);
\r
1368 //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
\r
1369 //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
\r
1370 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
1371 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
1373 //Flow after burner
\r
1374 if(fUseFlowAfterBurner) {
\r
1375 Double_t precisionPhi = 0.001;
\r
1376 Int_t maxNumberOfIterations = 100;
\r
1378 Double_t phi0 = vPhi;
\r
1379 Double_t gV2 = fDifferentialV2->Eval(vPt);
\r
1381 for (Int_t j = 0; j < maxNumberOfIterations; j++) {
\r
1382 Double_t phiprev = vPhi;
\r
1383 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
\r
1384 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
\r
1386 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
\r
1388 //Printf("phi (after): %lf\n",vPhi);
\r
1389 Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
\r
1390 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
\r
1391 fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
\r
1393 Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
\r
1394 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
\r
1395 fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
\r
1399 //vPhi *= TMath::RadToDeg();
\r
1401 //=======================================correction
\r
1402 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge,gCentrality);
\r
1404 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
1406 }//MC event object
\r
1409 return tracksAccepted;
\r
1412 //________________________________________________________________________
\r
1413 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
\r
1414 // Clones TObjArray and returns it with tracks after shuffling the charges
\r
1416 TObjArray* tracksShuffled = new TObjArray;
\r
1417 tracksShuffled->SetOwner(kTRUE);
\r
1419 vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks
\r
1421 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
\r
1423 AliVParticle* track = (AliVParticle*) tracks->At(i);
\r
1424 chargeVector->push_back(track->Charge());
\r
1427 random_shuffle(chargeVector->begin(), chargeVector->end());
\r
1429 for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
\r
1430 AliVParticle* track = (AliVParticle*) tracks->At(i);
\r
1431 //==============================correction
\r
1432 Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
\r
1433 tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
\r
1436 delete chargeVector;
\r
1438 return tracksShuffled;
\r
1442 //________________________________________________________________________
\r
1443 void AliAnalysisTaskBFPsi::FinishTaskOutput(){
\r
1444 //Printf("END BF");
\r
1447 AliError("fBalance not available");
\r
1450 if(fRunShuffling) {
\r
1451 if (!fShuffledBalance) {
\r
1452 AliError("fShuffledBalance not available");
\r
1459 //________________________________________________________________________
\r
1460 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
\r
1461 // Draw result to the screen
\r
1462 // Called once at the end of the query
\r
1464 // not implemented ...
\r