4 #include "TLorentzVector.h"
5 #include "TGraphErrors.h"
16 #include "AliAnalysisTaskSE.h"
17 #include "AliAnalysisManager.h"
19 #include "AliESDVertex.h"
20 #include "AliESDEvent.h"
21 #include "AliESDInputHandler.h"
22 #include "AliAODEvent.h"
23 #include "AliAODTrack.h"
24 #include "AliAODInputHandler.h"
25 #include "AliAODMCParticle.h"
26 #include "AliCollisionGeometry.h"
27 #include "AliGenEventHeader.h"
28 #include "AliMCEventHandler.h"
29 #include "AliMCEvent.h"
31 #include "AliESDtrackCuts.h"
32 #include "AliEventplane.h"
35 #include "AliAnalysisUtils.h"
37 #include "AliEventPoolManager.h"
40 #include "AliPIDResponse.h"
41 #include "AliPIDCombined.h"
43 #include "AliAnalysisTaskBFPsi.h"
44 #include "AliBalancePsi.h"
45 #include "AliAnalysisTaskTriggeredBF.h"
48 // Analysis task for the BF vs Psi code
49 // Authors: Panos.Christakoglou@nikhef.nl
54 ClassImp(AliAnalysisTaskBFPsi)
56 //________________________________________________________________________
57 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
58 : AliAnalysisTaskSE(name),
60 fArrayMC(0), //+++++++++++++
62 fRunShuffling(kFALSE),
65 fRunMixingEventPlane(kFALSE),
76 fHistCentStatsUsed(0),
82 fHistTPCvsVZEROMultiplicity(0),
100 fHistdEdxVsPTPCbeforePID(NULL),
101 fHistBetavsPTOFbeforePID(NULL),
102 fHistProbTPCvsPtbeforePID(NULL),
103 fHistProbTOFvsPtbeforePID(NULL),
104 fHistProbTPCTOFvsPtbeforePID(NULL),
105 fHistNSigmaTPCvsPtbeforePID(NULL),
106 fHistNSigmaTOFvsPtbeforePID(NULL),
107 fHistBetaVsdEdXbeforePID(NULL), //+++++++
108 fHistNSigmaTPCTOFvsPtbeforePID(NULL), //++++++
109 fHistdEdxVsPTPCafterPID(NULL),
110 fHistBetavsPTOFafterPID(NULL),
111 fHistProbTPCvsPtafterPID(NULL),
112 fHistProbTOFvsPtafterPID(NULL),
113 fHistProbTPCTOFvsPtafterPID(NULL),
114 fHistNSigmaTPCvsPtafterPID(NULL),
115 fHistNSigmaTOFvsPtafterPID(NULL),
116 fHistBetaVsdEdXafterPID(NULL), //+++++++
117 fHistNSigmaTPCTOFvsPtafterPID(NULL), //+++++++
118 fHistdEdxVsPTPCbeforePIDelectron(NULL), //+++++++
119 fHistNSigmaTPCvsPtbeforePIDelectron(NULL), //+++++++
120 fHistdEdxVsPTPCafterPIDelectron(NULL), //+++++++
121 fHistNSigmaTPCvsPtafterPIDelectron(NULL), //+++++++
122 fCentralityArrayBinsForCorrections(kCENTRALITY),
123 fCentralityWeights(0x0),
126 fParticleOfInterest(kPion),
127 fPidDetectorConfig(kTPCTOF),
129 fUsePIDnSigma(kTRUE),
130 fUsePIDPropabilities(kFALSE),
132 fMinAcceptedPIDProbability(0.8),
133 fElectronRejection(kFALSE),
134 fElectronOnlyRejection(kFALSE),
135 fElectronRejectionNSigma(-1.),
136 fElectronRejectionMinPt(0.),
137 fElectronRejectionMaxPt(1000.),
139 fCentralityEstimator("V0M"),
140 fUseCentrality(kFALSE),
141 fCentralityPercentileMin(0.),
142 fCentralityPercentileMax(5.),
143 fImpactParameterMin(0.),
144 fImpactParameterMax(20.),
145 fMultiplicityEstimator("V0A"),
146 fUseMultiplicity(kFALSE),
147 fNumberOfAcceptedTracksMin(0),
148 fNumberOfAcceptedTracksMax(10000),
149 fHistNumberOfAcceptedTracks(0),
150 fHistMultiplicity(0),
151 fUseOfflineTrigger(kFALSE),
152 fCheckFirstEventInChunk(kFALSE),
153 fCheckPileUp(kFALSE),
154 fCheckPrimaryFlagAOD(kFALSE),
155 fUseMCforKinematics(kFALSE),
159 fnAODtrackCutBit(128),
169 fNClustersTPCCut(-1),
170 fAcceptanceParameterization(0),
172 fUseFlowAfterBurner(kFALSE),
173 fExcludeResonancesInMC(kFALSE),
174 fExcludeElectronsInMC(kFALSE),
175 fUseMCPdgCode(kFALSE),
176 fPDGCodeToBeAnalyzed(-1),
177 fEventClass("EventPlane"),
179 fHistVZEROAGainEqualizationMap(0),
180 fHistVZEROCGainEqualizationMap(0),
181 fHistVZEROChannelGainEqualizationMap(0) {
183 // Define input and output slots here
184 // Input slot #0 works with a TChain
186 //======================================================correction
187 for (Int_t i=0; i<kCENTRALITY; i++){
188 fHistCorrectionPlus[i] = NULL;
189 fHistCorrectionMinus[i] = NULL;
190 fCentralityArrayForCorrections[i] = -1.;
192 //=====================================================correction
194 DefineInput(0, TChain::Class());
195 // Output slot #0 writes into a TH1 container
196 DefineOutput(1, TList::Class());
197 DefineOutput(2, TList::Class());
198 DefineOutput(3, TList::Class());
199 DefineOutput(4, TList::Class());
200 DefineOutput(5, TList::Class());
203 //________________________________________________________________________
204 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
207 // delete fShuffledBalance;
212 // delete fHistEventStats;
213 // delete fHistTrackStats;
224 // delete fHistEtaPhiPos;
225 // delete fHistEtaPhiNeg;
229 //________________________________________________________________________
230 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
234 // global switch disabling the reference
235 // (to avoid "Replacing existing TH1" if several wagons are created in train)
236 Bool_t oldStatus = TH1::AddDirectoryStatus();
237 TH1::AddDirectory(kFALSE);
240 fBalance = new AliBalancePsi();
241 fBalance->SetAnalysisLevel("ESD");
242 fBalance->SetEventClass(fEventClass);
243 //fBalance->SetNumberOfBins(-1,16);
244 //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
247 if(!fShuffledBalance) {
248 fShuffledBalance = new AliBalancePsi();
249 fShuffledBalance->SetAnalysisLevel("ESD");
250 fShuffledBalance->SetEventClass(fEventClass);
251 //fShuffledBalance->SetNumberOfBins(-1,16);
252 //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
257 fMixedBalance = new AliBalancePsi();
258 fMixedBalance->SetAnalysisLevel("ESD");
259 fMixedBalance->SetEventClass(fEventClass);
265 fList->SetName("listQA");
268 //Balance Function list
269 fListBF = new TList();
270 fListBF->SetName("listBF");
274 fListBFS = new TList();
275 fListBFS->SetName("listBFShuffled");
276 fListBFS->SetOwner();
280 fListBFM = new TList();
281 fListBFM->SetName("listTriggeredBFMixed");
282 fListBFM->SetOwner();
286 if(fUsePID || fElectronRejection) {
287 fHistListPIDQA = new TList();
288 fHistListPIDQA->SetName("listQAPID");
289 fHistListPIDQA->SetOwner();
293 TString gCutName[7] = {"Total","Offline trigger",
294 "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};
295 fHistEventStats = new TH2F("fHistEventStats",
296 "Event statistics;;Centrality percentile;N_{events}",
297 7,0.5,7.5,220,-5,105);
298 for(Int_t i = 1; i <= 7; i++)
299 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
300 fList->Add(fHistEventStats);
302 TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
303 fHistCentStats = new TH2F("fHistCentStats",
304 "Centrality statistics;;Cent percentile",
305 13,-0.5,12.5,220,-5,105);
306 for(Int_t i = 1; i <= 13; i++){
307 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
308 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
310 fList->Add(fHistCentStats);
312 fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);
313 fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());
314 fList->Add(fHistCentStatsUsed);
316 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
317 fList->Add(fHistTriggerStats);
319 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
320 fList->Add(fHistTrackStats);
322 fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
323 fList->Add(fHistNumberOfAcceptedTracks);
325 fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);
326 fList->Add(fHistMultiplicity);
328 // Vertex distributions
329 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
331 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
333 fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
336 //TPC vs VZERO multiplicity
337 fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
338 if(fMultiplicityEstimator == "V0A")
339 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");
340 else if(fMultiplicityEstimator == "V0C")
341 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");
343 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");
344 fList->Add(fHistTPCvsVZEROMultiplicity);
346 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
347 fList->Add(fHistVZEROSignal);
350 fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
351 fList->Add(fHistEventPlane);
354 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
355 fList->Add(fHistClus);
356 fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);
357 fList->Add(fHistChi2);
358 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
359 fList->Add(fHistDCA);
360 fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);
362 fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);
363 fList->Add(fHistEta);
364 fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);
365 fList->Add(fHistRapidity);
366 fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);
367 fList->Add(fHistPhi);
368 fHistEtaPhiPos = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);
369 fList->Add(fHistEtaPhiPos);
370 fHistEtaPhiNeg = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);
371 fList->Add(fHistEtaPhiNeg);
372 fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
373 fList->Add(fHistPhiBefore);
374 fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
375 fList->Add(fHistPhiAfter);
376 fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
377 fList->Add(fHistPhiPos);
378 fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);
379 fList->Add(fHistPhiNeg);
380 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
381 fList->Add(fHistV0M);
382 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
383 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
384 for(Int_t i = 1; i <= 6; i++)
385 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
386 fList->Add(fHistRefTracks);
388 // Balance function histograms
389 // Initialize histograms if not done yet (including the custom binning)
390 if(!fBalance->GetHistNp()){
391 AliInfo("Histograms not yet initialized! --> Will be done now");
392 fBalance->SetCustomBinning(fCustomBinning);
393 fBalance->InitHistograms();
397 if(!fShuffledBalance->GetHistNp()) {
398 AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");
399 fShuffledBalance->SetCustomBinning(fCustomBinning);
400 fShuffledBalance->InitHistograms();
405 if(!fMixedBalance->GetHistNp()) {
406 AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");
407 fMixedBalance->SetCustomBinning(fCustomBinning);
408 fMixedBalance->InitHistograms();
412 // QA histograms for different cuts
413 fList->Add(fBalance->GetQAHistHBTbefore());
414 fList->Add(fBalance->GetQAHistHBTafter());
415 fList->Add(fBalance->GetQAHistConversionbefore());
416 fList->Add(fBalance->GetQAHistConversionafter());
417 fList->Add(fBalance->GetQAHistPsiMinusPhi());
418 fList->Add(fBalance->GetQAHistResonancesBefore());
419 fList->Add(fBalance->GetQAHistResonancesRho());
420 fList->Add(fBalance->GetQAHistResonancesK0());
421 fList->Add(fBalance->GetQAHistResonancesLambda());
422 fList->Add(fBalance->GetQAHistQbefore());
423 fList->Add(fBalance->GetQAHistQafter());
425 //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
426 fListBF->Add(fBalance->GetHistNp());
427 fListBF->Add(fBalance->GetHistNn());
428 fListBF->Add(fBalance->GetHistNpn());
429 fListBF->Add(fBalance->GetHistNnn());
430 fListBF->Add(fBalance->GetHistNpp());
431 fListBF->Add(fBalance->GetHistNnp());
434 fListBFS->Add(fShuffledBalance->GetHistNp());
435 fListBFS->Add(fShuffledBalance->GetHistNn());
436 fListBFS->Add(fShuffledBalance->GetHistNpn());
437 fListBFS->Add(fShuffledBalance->GetHistNnn());
438 fListBFS->Add(fShuffledBalance->GetHistNpp());
439 fListBFS->Add(fShuffledBalance->GetHistNnp());
443 fListBFM->Add(fMixedBalance->GetHistNp());
444 fListBFM->Add(fMixedBalance->GetHistNn());
445 fListBFM->Add(fMixedBalance->GetHistNpn());
446 fListBFM->Add(fMixedBalance->GetHistNnn());
447 fListBFM->Add(fMixedBalance->GetHistNpp());
448 fListBFM->Add(fMixedBalance->GetHistNnp());
455 Int_t trackDepth = fMixingTracks;
456 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
459 Double_t* centbins = NULL;
460 Int_t nCentralityBins;
461 if(fBalance->IsUseVertexBinning()){
462 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);
465 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);
469 Double_t* multbins = NULL;
470 Int_t nMultiplicityBins;
471 multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);
474 Double_t* vtxbins = NULL;
476 if(fBalance->IsUseVertexBinning()){
477 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);
480 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);
483 // Event plane angle (Psi) bins
484 Double_t* psibins = NULL;
486 psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);
489 // run the event mixing also in bins of event plane (statistics!)
490 if(fRunMixingEventPlane){
491 if(fEventClass=="Multiplicity"){
492 if(multbins && vtxbins && psibins){
493 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
497 if(centbins && vtxbins && psibins){
498 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
503 if(fEventClass=="Multiplicity"){
504 if(multbins && vtxbins){
505 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
509 if(centbins && vtxbins){
510 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
515 if(centbins) delete [] centbins;
516 if(multbins) delete [] multbins;
517 if(vtxbins) delete [] vtxbins;
519 // check pool manager
521 AliError("Event Mixing required, but Pool Manager not initialized...");
526 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
528 //====================PID========================//
530 fPIDCombined = new AliPIDCombined();
531 fPIDCombined->SetDefaultTPCPriors();
533 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
534 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
536 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
537 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID);
539 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
540 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID);
542 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
543 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);
545 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
546 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);
548 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
549 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);
551 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
552 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID);
554 fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2);
555 fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++
557 fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500);
558 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++
560 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
561 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
563 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
564 fHistListPIDQA->Add(fHistBetavsPTOFafterPID);
566 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
567 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);
569 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
570 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID);
572 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
573 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);
575 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
576 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
578 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
579 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);
581 fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2);
582 fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++
584 fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500);
585 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++
588 // for electron rejection only TPC nsigma histograms
589 if(fElectronRejection) {
591 fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000);
592 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);
594 fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500);
595 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);
597 fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000);
598 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);
600 fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500);
601 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron);
603 //====================PID========================//
607 PostData(2, fListBF);
608 if(fRunShuffling) PostData(3, fListBFS);
609 if(fRunMixing) PostData(4, fListBFM);
610 if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID
612 AliInfo("Finished setting up the Output");
614 TH1::AddDirectory(oldStatus);
619 //________________________________________________________________________
620 void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename,
621 Int_t nCentralityBins,
622 Double_t *centralityArrayForCorrections) {
623 //Open files that will be used for correction
624 fCentralityArrayBinsForCorrections = nCentralityBins;
625 for (Int_t i=0; i<nCentralityBins; i++)
626 fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];
628 // No file specified -> run without corrections
629 if(!filename.Contains(".root")) {
630 AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));
634 //Open the input file
635 TFile *f = TFile::Open(filename);
637 AliInfo(Form("File %s not found --> run without corrections",filename.Data()));
641 //TString listEffName = "";
642 for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {
643 //Printf("iCent %d:",iCent);
644 TString histoName = "fHistCorrectionPlus";
645 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
646 fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));
647 if(!fHistCorrectionPlus[iCent]) {
648 AliError(Form("fHist %s not found!!!",histoName.Data()));
652 histoName = "fHistCorrectionMinus";
653 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
654 fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data()));
655 if(!fHistCorrectionMinus[iCent]) {
656 AliError(Form("fHist %s not found!!!",histoName.Data()));
659 }//loop over centralities: ONLY the PbPb case is covered
662 //________________________________________________________________________
663 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
665 // Called for each event
667 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
668 Int_t gNumberOfAcceptedTracks = 0;
669 Double_t lMultiplicityVar = -999.; //-1
670 Double_t gReactionPlane = -1.;
673 // get the event (for generator level: MCEvent())
674 AliVEvent* eventMain = NULL;
675 if(gAnalysisLevel == "MC") {
676 eventMain = dynamic_cast<AliVEvent*>(MCEvent());
679 eventMain = dynamic_cast<AliVEvent*>(InputEvent());
680 // for HBT like cuts need magnetic field sign
681 bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
684 AliError("eventMain not available");
688 // PID Response task active?
689 if(fUsePID || fElectronRejection) {
690 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
691 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
694 // check event cuts and fill event histograms
695 if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){
698 // get the reaction plane
699 if(fEventClass != "Multiplicity" && gAnalysisLevel!="AODnano") {
700 gReactionPlane = GetEventPlane(eventMain);
701 fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
702 if(gReactionPlane < 0){
707 // get the accepted tracks in main event
708 TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
709 gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
711 //multiplicity cut (used in pp)
712 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
714 // store charges of all accepted tracks,shuffle and reassign(two extra loops)
715 TObjArray* tracksShuffled = NULL;
717 tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
723 // 1. First get an event pool corresponding in mult (cent) and
724 // zvertex to the current event. Once initialized, the pool
725 // should contain nMix (reduced) events. This routine does not
726 // pre-scan the chain. The first several events of every chain
727 // will be skipped until the needed pools are filled to the
728 // specified depth. If the pool categories are not too rare, this
729 // should not be a problem. If they are rare, you could lose`
732 // 2. Collect the whole pool's content of tracks into one TObjArray
733 // (bgTracks), which is effectively a single background super-event.
735 // 3. The reduced and bgTracks arrays must both be passed into
736 // FillCorrelations(). Also nMix should be passed in, so a weight
737 // of 1./nMix can be applied.
739 AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
742 AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
748 if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){
751 Int_t nMix = pool->GetCurrentNEvents();
752 //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
754 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
755 //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
756 //if (pool->IsReady())
757 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
759 // Fill mixed-event histos here
760 for (Int_t jMix=0; jMix<nMix; jMix++)
762 TObjArray* tracksMixed = pool->GetEvent(jMix);
763 fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
767 // Update the Event pool
768 pool->UpdatePool(tracksMain);
774 // calculate balance function
775 fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
777 // calculate shuffled balance function
778 if(fRunShuffling && tracksShuffled != NULL) {
779 fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
783 //________________________________________________________________________
784 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
785 // Checks the Event cuts
786 // Fills Event statistics histograms
788 Bool_t isSelectedMain = kTRUE;
789 Float_t gRefMultiplicity = -1.;
790 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
792 AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
794 fHistEventStats->Fill(1,gRefMultiplicity); //all events
796 // check first event in chunk (is not needed for new reconstructions)
797 if(fCheckFirstEventInChunk){
799 if(ut.IsFirstEventInChunk(event))
801 fHistEventStats->Fill(6,gRefMultiplicity);
803 // check for pile-up event
806 ut.SetUseMVPlpSelection(kTRUE);
807 ut.SetUseOutOfBunchPileUp(kTRUE);
808 if(ut.IsPileUpEvent(event))
810 fHistEventStats->Fill(7,gRefMultiplicity);
813 // Event trigger bits
814 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
815 if(fUseOfflineTrigger)
816 isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
819 fHistEventStats->Fill(2,gRefMultiplicity); //triggered events
822 if(gAnalysisLevel == "MC") {
824 AliError("mcEvent not available");
829 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
831 TArrayF gVertexArray;
832 header->PrimaryVertex(gVertexArray);
833 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
834 //gVertexArray.At(0),
835 //gVertexArray.At(1),
836 //gVertexArray.At(2));
837 fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex
838 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
839 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
840 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
841 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
843 // get the reference multiplicty or centrality
844 gRefMultiplicity = GetRefMultiOrCentrality(event);
846 fHistVx->Fill(gVertexArray.At(0));
847 fHistVy->Fill(gVertexArray.At(1));
848 fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);
850 // take only events inside centrality class
852 if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){
853 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
854 return gRefMultiplicity;
857 // take events only within the same multiplicity class
858 else if(fUseMultiplicity) {
859 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
860 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
861 return gRefMultiplicity;
863 }//multiplicity range
871 // Event Vertex AOD, ESD, ESDMC
873 const AliVVertex *vertex = event->GetPrimaryVertex();
877 vertex->GetCovarianceMatrix(fCov);
878 if(vertex->GetNContributors() > 0) {
880 fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex
881 if(TMath::Abs(vertex->GetX()) < fVxMax) {
882 if(TMath::Abs(vertex->GetY()) < fVyMax) {
883 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
884 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
886 // get the reference multiplicty or centrality
887 gRefMultiplicity = GetRefMultiOrCentrality(event);
889 fHistVx->Fill(vertex->GetX());
890 fHistVy->Fill(vertex->GetY());
891 fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);
893 // take only events inside centrality class
895 if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){
897 // centrality weighting (optional for 2011 if central and semicentral triggers are used)
898 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity)){
899 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
903 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
904 return gRefMultiplicity;
907 // take events only within the same multiplicity class
908 else if(fUseMultiplicity) {
910 //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,
911 //fNumberOfAcceptedTracksMax,gRefMultiplicity);
913 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
914 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
915 return gRefMultiplicity;
917 }//multiplicity range
921 }//proper vertex resolution
922 }//proper number of contributors
923 }//vertex object valid
927 // in all other cases return -1 (event not accepted)
932 //________________________________________________________________________
933 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
934 // Checks the Event cuts
935 // Fills Event statistics histograms
937 Float_t gCentrality = -1.;
938 Double_t gMultiplicity = -1.;
939 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
942 // calculate centrality always (not only in centrality mode)
943 if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++
944 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
946 gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
948 // QA for centrality estimators
949 fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
950 fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));
951 fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));
952 fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
953 fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
954 fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL"));
955 fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
956 fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
957 fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));
958 fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));
959 fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
960 fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
961 fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
963 // Centrality estimator USED ++++++++++++++++++++++++++++++
964 fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));
966 // centrality QA (V0M)
967 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
969 // centrality QA (reference tracks)
970 fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
971 fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
972 fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
973 fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
974 fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
975 fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
976 fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
977 fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
978 fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
983 // calculate centrality always (not only in centrality mode)
984 else if(gAnalysisLevel == "AODnano" ) { //centrality via JF workaround
986 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
988 gCentrality = (Float_t) gROOT->ProcessLine(Form("100.0 + 100.0 * ((AliNanoAODHeader*) %p)->GetCentrality(\"%s\")", header,fCentralityEstimator.Data())) / 100 - 1.0;
991 fHistCentStatsUsed->Fill(0.,gCentrality);
996 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
997 AliCentrality *centrality = event->GetCentrality();
998 gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
1000 // QA for centrality estimators
1001 fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));
1002 fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));
1003 fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));
1004 fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));
1005 fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));
1006 fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));
1007 fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));
1008 fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));
1009 fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));
1010 fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));
1011 fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));
1012 fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));
1013 fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));
1015 // Centrality estimator USED ++++++++++++++++++++++++++++++
1016 fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
1018 // centrality QA (V0M)
1019 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
1022 else if(gAnalysisLevel == "MC"){
1023 Double_t gImpactParameter = 0.;
1024 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1026 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
1028 gImpactParameter = headerH->ImpactParameter();
1029 gCentrality = gImpactParameter;
1038 // calculate reference multiplicity always (not only in multiplicity mode)
1039 if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){
1040 AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
1042 gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
1043 fHistMultiplicity->Fill(gMultiplicity);
1047 else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){
1048 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
1049 if ((fMultiplicityEstimator == "V0M")||
1050 (fMultiplicityEstimator == "V0A")||
1051 (fMultiplicityEstimator == "V0C") ||
1052 (fMultiplicityEstimator == "TPC")) {
1053 gMultiplicity = GetReferenceMultiplicityFromAOD(event);
1054 if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);
1058 gMultiplicity = header->GetRefMultiplicity();
1059 if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);
1061 fHistMultiplicity->Fill(gMultiplicity);
1063 else if(gAnalysisLevel == "MC") {
1064 AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);
1065 //Calculating the multiplicity as the number of charged primaries
1066 //within \pm 0.8 in eta and pT > 0.1 GeV/c
1067 for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
1068 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
1070 AliError(Form("Could not receive particle %d", iParticle));
1074 //exclude non stable particles
1075 if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
1078 if (fMultiplicityEstimator == "V0M") {
1079 if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7))
1081 else if (fMultiplicityEstimator == "V0A") {
1082 if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}
1083 else if (fMultiplicityEstimator == "V0C") {
1084 if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}
1085 else if (fMultiplicityEstimator == "TPC") {
1086 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
1087 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
1090 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
1091 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
1095 if(track->Charge() == 0) continue;
1098 }//loop over primaries
1099 fHistMultiplicity->Fill(gMultiplicity);
1106 // decide what should be returned only here
1107 Double_t lReturnVal = -100;
1108 if(fEventClass=="Multiplicity"){
1109 lReturnVal = gMultiplicity;
1110 }else if(fEventClass=="Centrality"){
1111 lReturnVal = gCentrality;
1116 //________________________________________________________________________
1117 Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){
1118 //Function that returns the reference multiplicity from AODs (data or reco MC)
1119 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
1120 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
1121 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
1123 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
1125 Printf("ERROR: AOD header not available");
1128 Int_t gRunNumber = header->GetRunNumber();
1130 // Loop over tracks in event
1131 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1132 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1134 AliError(Form("Could not receive track %d", iTracks));
1139 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1141 if(aodTrack->Charge() == 0) continue;
1142 // Kinematics cuts from ESD track cuts
1143 if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;
1144 if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;
1146 //=================PID (so far only for electron rejection)==========================//
1147 if(fElectronRejection) {
1148 // get the electron nsigma
1149 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1151 // check only for given momentum range
1152 if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){
1153 //look only at electron nsigma
1154 if(!fElectronOnlyRejection) {
1155 //Make the decision based on the n-sigma of electrons
1156 if(nSigma < fElectronRejectionNSigma) continue;
1159 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1160 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1161 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1163 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1164 if(nSigma < fElectronRejectionNSigma
1165 && nSigmaPions > fElectronRejectionNSigma
1166 && nSigmaKaons > fElectronRejectionNSigma
1167 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1170 }//electron rejection
1172 gRefMultiplicityTPC += 1.0;
1175 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
1176 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
1177 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
1180 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
1181 else if(iChannel >= 32)
1182 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
1185 //Equalization of gain
1186 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
1188 gRefMultiplicityVZEROA /= gFactorA;
1189 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
1191 gRefMultiplicityVZEROC /= gFactorC;
1192 if((gFactorA != 0)&&(gFactorC != 0))
1193 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
1196 Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);
1198 fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
1200 if(fMultiplicityEstimator == "TPC")
1201 gRefMultiplicity = gRefMultiplicityTPC;
1202 else if(fMultiplicityEstimator == "V0M")
1203 gRefMultiplicity = gRefMultiplicityVZERO;
1204 else if(fMultiplicityEstimator == "V0A")
1205 gRefMultiplicity = gRefMultiplicityVZEROA;
1206 else if(fMultiplicityEstimator == "V0C")
1207 gRefMultiplicity = gRefMultiplicityVZEROC;
1209 return gRefMultiplicity;
1212 //________________________________________________________________________
1213 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
1214 // Get the event plane
1216 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1218 Float_t gVZEROEventPlane = -10.;
1219 Float_t gReactionPlane = -10.;
1220 Double_t qxTot = 0.0, qyTot = 0.0;
1222 //MC: from reaction plane
1223 if(gAnalysisLevel == "MC"){
1225 AliError("mcEvent not available");
1229 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1231 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
1233 gReactionPlane = headerH->ReactionPlaneAngle();
1234 //gReactionPlane *= TMath::RadToDeg();
1239 // AOD,ESD,ESDMC: from VZERO Event Plane
1242 AliEventplane *ep = event->GetEventplane();
1244 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
1245 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
1246 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
1247 gReactionPlane = gVZEROEventPlane;
1251 return gReactionPlane;
1254 //________________________________________________________________________
1255 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta,
1259 Double_t gCentrality) {
1260 // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality)
1262 Double_t correction = 1.;
1263 Int_t gCentralityInt = -1;
1265 for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){
1266 if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){
1272 // centrality not in array --> no correction
1273 if(gCentralityInt < 0){
1278 //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);
1280 if(fHistCorrectionPlus[gCentralityInt]){
1282 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1283 //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);
1286 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1287 //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);
1293 }//centrality in array
1295 if (correction == 0.) {
1296 AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt));
1303 //________________________________________________________________________
1304 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){
1305 // Returns TObjArray with tracks after all track cuts (only for AOD!)
1306 // Fills QA histograms
1308 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1310 //output TObjArray holding all good tracks
1311 TObjArray* tracksAccepted = new TObjArray;
1312 tracksAccepted->SetOwner(kTRUE);
1321 if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
1322 // Loop over tracks in event
1324 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1325 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1327 AliError(Form("Could not receive track %d", iTracks));
1333 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
1334 // take only TPC only tracks
1335 //fHistTrackStats->Fill(aodTrack->GetFilterMap());
1336 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1337 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1340 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1343 // additional check on kPrimary flag
1344 if(fCheckPrimaryFlagAOD){
1345 if(aodTrack->GetType() != AliAODTrack::kPrimary)
1350 vCharge = aodTrack->Charge();
1351 vEta = aodTrack->Eta();
1353 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1354 vPt = aodTrack->Pt();
1356 //===========================PID (so far only for electron rejection)===============================//
1357 if(fElectronRejection) {
1359 // get the electron nsigma
1360 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1362 //Fill QA before the PID
1363 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1364 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1365 //end of QA-before pid
1367 // check only for given momentum range
1368 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1370 //look only at electron nsigma
1371 if(!fElectronOnlyRejection){
1373 //Make the decision based on the n-sigma of electrons
1374 if(nSigma < fElectronRejectionNSigma) continue;
1378 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1379 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1380 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1382 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1383 if(nSigma < fElectronRejectionNSigma
1384 && nSigmaPions > fElectronRejectionNSigma
1385 && nSigmaKaons > fElectronRejectionNSigma
1386 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1390 //Fill QA after the PID
1391 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1392 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1395 //===========================end of PID (so far only for electron rejection)===============================//
1397 //+++++++++++++++++++++++++++++//
1398 //===========================PID===============================//
1400 Double_t prob[AliPID::kSPECIES]={0.};
1401 Double_t probTPC[AliPID::kSPECIES]={0.};
1402 Double_t probTOF[AliPID::kSPECIES]={0.};
1403 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1405 Double_t nSigma = 0.;
1406 Double_t nSigmaTPC = 0.; //++++
1407 Double_t nSigmaTOF = 0.; //++++
1408 Double_t nSigmaTPCTOF = 0.; //++++
1409 UInt_t detUsedTPC = 0;
1410 UInt_t detUsedTOF = 0;
1411 UInt_t detUsedTPCTOF = 0;
1413 nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
1414 nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
1415 nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);
1417 //Decide what detector configuration we want to use
1418 switch(fPidDetectorConfig) {
1420 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1422 detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
1423 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1424 prob[iSpecies] = probTPC[iSpecies];
1427 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1429 detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
1430 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1431 prob[iSpecies] = probTOF[iSpecies];
1434 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1435 nSigma = nSigmaTPCTOF;
1436 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);
1437 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1438 prob[iSpecies] = probTPCTOF[iSpecies];
1442 }//end switch: define detector mask
1444 //Filling the PID QA
1445 Double_t tofTime = -999., length = 999., tof = -999.;
1446 Double_t c = TMath::C()*1.E-9;// m/ns
1447 Double_t beta = -999.;
1448 if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) &&
1449 (aodTrack->IsOn(AliAODTrack::kTIME)) ) {
1450 tofTime = aodTrack->GetTOFsignal();//in ps
1451 length = aodTrack->GetIntegratedLength();
1452 tof = tofTime*1E-3; // ns
1455 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
1459 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
1463 length = length*0.01; // in meters
1467 fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1468 fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1469 fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);
1472 fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1473 fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
1474 fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC);
1476 fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1479 fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
1480 fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF);
1481 Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF);
1483 //end of QA-before pid
1485 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
1486 //Make the decision based on the n-sigma
1488 if(nSigma > fPIDNSigma) continue;
1490 fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF);
1491 fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC);
1492 fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF);
1494 //Make the decision based on the bayesian
1495 else if(fUsePIDPropabilities) {
1496 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
1497 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
1499 fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1500 fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
1501 fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1505 //Fill QA after the PID
1506 fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1507 fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1508 fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
1511 //===========================PID===============================//
1512 //+++++++++++++++++++++++++++++//
1515 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
1516 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
1519 // Kinematics cuts from ESD track cuts
1520 if( vPt < fPtMin || vPt > fPtMax) continue;
1521 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1523 // Extra DCA cuts (for systematic studies [!= -1])
1524 if( fDCAxyCut != -1 && fDCAzCut != -1){
1525 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1530 // Extra TPC cuts (for systematic studies [!= -1])
1531 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1534 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1538 // fill QA histograms
1539 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1540 fHistDCA->Fill(dcaZ,dcaXY);
1541 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1542 fHistPt->Fill(vPt,gCentrality);
1543 fHistEta->Fill(vEta,gCentrality);
1544 fHistRapidity->Fill(vY,gCentrality);
1545 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1546 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1547 fHistPhi->Fill(vPhi,gCentrality);
1548 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1549 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1551 //=======================================correction
1552 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1553 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1555 // add the track to the TObjArray
1556 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1562 else if(gAnalysisLevel == "AODnano") { // not fully supported yet (PID missing)
1563 // Loop over tracks in event
1565 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1566 AliVTrack* aodTrack = dynamic_cast<AliVTrack *>(event->GetTrack(iTracks));
1568 AliError(Form("Could not receive track %d", iTracks));
1572 // AOD track cuts (not needed)
1573 //if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1575 vCharge = aodTrack->Charge();
1576 vEta = aodTrack->Eta();
1578 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1579 vPt = aodTrack->Pt();
1582 // Kinematics cuts from ESD track cuts
1583 if( vPt < fPtMin || vPt > fPtMax) continue;
1584 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1587 // fill QA histograms
1588 fHistPt->Fill(vPt,gCentrality);
1589 fHistEta->Fill(vEta,gCentrality);
1590 fHistRapidity->Fill(vY,gCentrality);
1591 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1592 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1593 fHistPhi->Fill(vPhi,gCentrality);
1594 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1595 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1597 //=======================================correction
1598 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1599 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1601 // add the track to the TObjArray
1602 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1604 }// AOD nano analysis
1607 //==============================================================================================================
1608 else if(gAnalysisLevel == "MCAOD") {
1610 AliMCEvent* mcEvent = MCEvent();
1612 AliError("ERROR: Could not retrieve MC event");
1616 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
1617 AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
1619 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
1623 if(!aodTrack->IsPhysicalPrimary()) continue;
1625 vCharge = aodTrack->Charge();
1626 vEta = aodTrack->Eta();
1628 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1629 vPt = aodTrack->Pt();
1631 // Kinematics cuts from ESD track cuts
1632 if( vPt < fPtMin || vPt > fPtMax) continue;
1633 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1635 // Remove neutral tracks
1636 if( vCharge == 0 ) continue;
1638 //Exclude resonances
1639 if(fExcludeResonancesInMC) {
1641 Bool_t kExcludeParticle = kFALSE;
1642 Int_t gMotherIndex = aodTrack->GetMother();
1643 if(gMotherIndex != -1) {
1644 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1646 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1647 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1648 //if(pdgCodeOfMother == 113) {
1649 if(pdgCodeOfMother == 113 // rho0
1650 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1651 // || pdgCodeOfMother == 221 // eta
1652 // || pdgCodeOfMother == 331 // eta'
1653 // || pdgCodeOfMother == 223 // omega
1654 // || pdgCodeOfMother == 333 // phi
1655 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1656 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1657 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1658 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1659 || pdgCodeOfMother == 111 // pi0 Dalitz
1660 || pdgCodeOfMother == 22 // photon
1662 kExcludeParticle = kTRUE;
1667 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1668 if(kExcludeParticle) continue;
1671 //Exclude electrons with PDG
1672 if(fExcludeElectronsInMC) {
1674 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
1678 // fill QA histograms
1679 fHistPt->Fill(vPt,gCentrality);
1680 fHistEta->Fill(vEta,gCentrality);
1681 fHistRapidity->Fill(vY,gCentrality);
1682 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1683 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1684 fHistPhi->Fill(vPhi,gCentrality);
1685 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1686 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1688 //=======================================correction
1689 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1690 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1692 // add the track to the TObjArray
1693 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1697 //==============================================================================================================
1699 //==============================================================================================================
1700 else if(gAnalysisLevel == "MCAODrec") {
1702 /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1704 printf("ERROR: fAOD not available\n");
1708 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
1710 AliError("No array of MC particles found !!!");
1713 AliMCEvent* mcEvent = MCEvent();
1715 AliError("ERROR: Could not retrieve MC event");
1716 return tracksAccepted;
1719 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1720 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1722 AliError(Form("Could not receive track %d", iTracks));
1726 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1727 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1729 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1731 vCharge = aodTrack->Charge();
1732 vEta = aodTrack->Eta();
1734 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1735 vPt = aodTrack->Pt();
1737 //===========================use MC information for Kinematics===============================//
1738 if(fUseMCforKinematics){
1740 Int_t label = TMath::Abs(aodTrack->GetLabel());
1741 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1744 vCharge = AODmcTrack->Charge();
1745 vEta = AODmcTrack->Eta();
1746 vY = AODmcTrack->Y();
1747 vPhi = AODmcTrack->Phi();// * TMath::RadToDeg();
1748 vPt = AODmcTrack->Pt();
1751 AliDebug(1, "no MC particle for this track");
1755 //===========================end of use MC information for Kinematics========================//
1758 //===========================PID (so far only for electron rejection)===============================//
1759 if(fElectronRejection) {
1761 // get the electron nsigma
1762 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1764 //Fill QA before the PID
1765 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1766 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1767 //end of QA-before pid
1769 // check only for given momentum range
1770 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1772 //look only at electron nsigma
1773 if(!fElectronOnlyRejection){
1775 //Make the decision based on the n-sigma of electrons
1776 if(nSigma < fElectronRejectionNSigma) continue;
1780 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1781 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1782 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1784 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1785 if(nSigma < fElectronRejectionNSigma
1786 && nSigmaPions > fElectronRejectionNSigma
1787 && nSigmaKaons > fElectronRejectionNSigma
1788 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1792 //Fill QA after the PID
1793 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1794 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1797 //===========================end of PID (so far only for electron rejection)===============================//
1799 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
1800 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
1802 // Kinematics cuts from ESD track cuts
1803 if( vPt < fPtMin || vPt > fPtMax) continue;
1804 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1806 // Extra DCA cuts (for systematic studies [!= -1])
1807 if( fDCAxyCut != -1 && fDCAzCut != -1){
1808 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1813 // Extra TPC cuts (for systematic studies [!= -1])
1814 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1817 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1821 //Exclude resonances
1822 if(fExcludeResonancesInMC) {
1824 Bool_t kExcludeParticle = kFALSE;
1826 Int_t label = TMath::Abs(aodTrack->GetLabel());
1827 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1830 //if (AODmcTrack->IsPhysicalPrimary()){
1832 Int_t gMotherIndex = AODmcTrack->GetMother();
1833 if(gMotherIndex != -1) {
1834 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1836 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1837 if(pdgCodeOfMother == 113 // rho0
1838 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1839 // || pdgCodeOfMother == 221 // eta
1840 // || pdgCodeOfMother == 331 // eta'
1841 // || pdgCodeOfMother == 223 // omega
1842 // || pdgCodeOfMother == 333 // phi
1843 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1844 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1845 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1846 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1847 || pdgCodeOfMother == 111 // pi0 Dalitz
1848 || pdgCodeOfMother == 22 // photon
1850 kExcludeParticle = kTRUE;
1855 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1856 if(kExcludeParticle) continue;
1859 //Exclude electrons with PDG
1860 if(fExcludeElectronsInMC) {
1862 Int_t label = TMath::Abs(aodTrack->GetLabel());
1863 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1866 if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
1870 // fill QA histograms
1871 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1872 fHistDCA->Fill(dcaZ,dcaXY);
1873 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1874 fHistPt->Fill(vPt,gCentrality);
1875 fHistEta->Fill(vEta,gCentrality);
1876 fHistRapidity->Fill(vY,gCentrality);
1877 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1878 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1879 fHistPhi->Fill(vPhi,gCentrality);
1880 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1881 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1883 //=======================================correction
1884 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1885 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1887 // add the track to the TObjArray
1888 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1891 //==============================================================================================================
1893 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
1895 AliESDtrack *trackTPC = NULL;
1896 AliMCParticle *mcTrack = NULL;
1898 AliMCEvent* mcEvent = NULL;
1900 // for MC ESDs use also MC information
1901 if(gAnalysisLevel == "MCESD"){
1902 mcEvent = MCEvent();
1904 AliError("mcEvent not available");
1905 return tracksAccepted;
1909 // Loop over tracks in event
1910 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1911 AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));
1913 AliError(Form("Could not receive track %d", iTracks));
1917 // for MC ESDs use also MC information --> MC track not used further???
1918 if(gAnalysisLevel == "MCESD"){
1919 Int_t label = TMath::Abs(track->GetLabel());
1920 if(label > mcEvent->GetNumberOfTracks()) continue;
1921 if(label > mcEvent->GetNumberOfPrimaries()) continue;
1923 mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
1924 if(!mcTrack) continue;
1927 // take only TPC only tracks
1928 trackTPC = new AliESDtrack();
1929 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
1933 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
1935 // fill QA histograms
1938 trackTPC->GetImpactParameters(b,bCov);
1939 if (bCov[0]<=0 || bCov[2]<=0) {
1940 AliDebug(1, "Estimated b resolution lower or equal zero!");
1941 bCov[0]=0; bCov[2]=0;
1944 Int_t nClustersTPC = -1;
1945 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
1946 //nClustersTPC = track->GetTPCclusters(0); // global track
1947 Float_t chi2PerClusterTPC = -1;
1948 if (nClustersTPC!=0) {
1949 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
1950 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
1953 //===========================PID===============================//
1955 Double_t prob[AliPID::kSPECIES]={0.};
1956 Double_t probTPC[AliPID::kSPECIES]={0.};
1957 Double_t probTOF[AliPID::kSPECIES]={0.};
1958 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1960 Double_t nSigma = 0.;
1961 UInt_t detUsedTPC = 0;
1962 UInt_t detUsedTOF = 0;
1963 UInt_t detUsedTPCTOF = 0;
1965 //Decide what detector configuration we want to use
1966 switch(fPidDetectorConfig) {
1968 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1969 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
1970 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
1971 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1972 prob[iSpecies] = probTPC[iSpecies];
1975 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1976 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
1977 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
1978 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1979 prob[iSpecies] = probTOF[iSpecies];
1982 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1983 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
1984 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1985 prob[iSpecies] = probTPCTOF[iSpecies];
1989 }//end switch: define detector mask
1991 //Filling the PID QA
1992 Double_t tofTime = -999., length = 999., tof = -999.;
1993 Double_t c = TMath::C()*1.E-9;// m/ns
1994 Double_t beta = -999.;
1995 Double_t nSigmaTOFForParticleOfInterest = -999.;
1996 if ( (track->IsOn(AliESDtrack::kTOFin)) &&
1997 (track->IsOn(AliESDtrack::kTIME)) ) {
1998 tofTime = track->GetTOFsignal();//in ps
1999 length = track->GetIntegratedLength();
2000 tof = tofTime*1E-3; // ns
2003 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
2007 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
2011 length = length*0.01; // in meters
2015 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
2016 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
2017 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
2018 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
2022 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
2023 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
2024 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
2025 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
2026 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
2027 //end of QA-before pid
2029 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
2030 //Make the decision based on the n-sigma
2032 if(nSigma > fPIDNSigma) continue;}
2034 //Make the decision based on the bayesian
2035 else if(fUsePIDPropabilities) {
2036 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
2037 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
2040 //Fill QA after the PID
2041 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
2042 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
2043 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
2045 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
2046 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
2047 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
2048 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
2051 //===========================PID===============================//
2052 vCharge = trackTPC->Charge();
2054 vEta = trackTPC->Eta();
2055 vPhi = trackTPC->Phi();// * TMath::RadToDeg();
2056 vPt = trackTPC->Pt();
2057 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
2058 fHistDCA->Fill(b[1],b[0]);
2059 fHistChi2->Fill(chi2PerClusterTPC,gCentrality);
2060 fHistPt->Fill(vPt,gCentrality);
2061 fHistEta->Fill(vEta,gCentrality);
2062 fHistPhi->Fill(vPhi,gCentrality);
2063 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2064 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2065 fHistRapidity->Fill(vY,gCentrality);
2066 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2067 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2069 //=======================================correction
2070 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
2071 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2073 // add the track to the TObjArray
2074 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
2080 else if(gAnalysisLevel == "MC"){
2082 AliError("mcEvent not available");
2086 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
2088 // Loop over tracks in event
2089 for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
2090 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
2092 AliError(Form("Could not receive particle %d", iTracks));
2096 //exclude non stable particles
2097 if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
2099 vCharge = track->Charge();
2100 vEta = track->Eta();
2104 if( vPt < fPtMin || vPt > fPtMax)
2107 if( vEta < fEtaMin || vEta > fEtaMax) continue;
2110 if( vY < fEtaMin || vY > fEtaMax) continue;
2113 // Remove neutral tracks
2114 if( vCharge == 0 ) continue;
2116 //analyze one set of particles
2118 TParticle *particle = track->Particle();
2119 if(!particle) continue;
2121 Int_t gPdgCode = particle->GetPdgCode();
2122 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
2126 //Use the acceptance parameterization
2127 if(fAcceptanceParameterization) {
2128 Double_t gRandomNumber = gRandom->Rndm();
2129 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
2133 //Exclude resonances
2134 if(fExcludeResonancesInMC) {
2135 TParticle *particle = track->Particle();
2136 if(!particle) continue;
2138 Bool_t kExcludeParticle = kFALSE;
2139 Int_t gMotherIndex = particle->GetFirstMother();
2140 if(gMotherIndex != -1) {
2141 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
2143 TParticle *motherParticle = motherTrack->Particle();
2144 if(motherParticle) {
2145 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
2146 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
2147 if(pdgCodeOfMother == 113 // rho0
2148 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
2149 // || pdgCodeOfMother == 221 // eta
2150 // || pdgCodeOfMother == 331 // eta'
2151 // || pdgCodeOfMother == 223 // omega
2152 // || pdgCodeOfMother == 333 // phi
2153 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
2154 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
2155 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
2156 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
2157 || pdgCodeOfMother == 111 // pi0 Dalitz
2159 kExcludeParticle = kTRUE;
2165 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
2166 if(kExcludeParticle) continue;
2169 //Exclude electrons with PDG
2170 if(fExcludeElectronsInMC) {
2172 TParticle *particle = track->Particle();
2175 if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
2179 vPhi = track->Phi();
2180 //Printf("phi (before): %lf",vPhi);
2182 fHistPt->Fill(vPt,gCentrality);
2183 fHistEta->Fill(vEta,gCentrality);
2184 fHistPhi->Fill(vPhi,gCentrality);
2185 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2186 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2187 //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2188 fHistRapidity->Fill(vY,gCentrality);
2189 //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2190 //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2191 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2192 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2195 if(fUseFlowAfterBurner) {
2196 Double_t precisionPhi = 0.001;
2197 Int_t maxNumberOfIterations = 100;
2199 Double_t phi0 = vPhi;
2200 Double_t gV2 = fDifferentialV2->Eval(vPt);
2202 for (Int_t j = 0; j < maxNumberOfIterations; j++) {
2203 Double_t phiprev = vPhi;
2204 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2205 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2207 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
2209 //Printf("phi (after): %lf\n",vPhi);
2210 Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
2211 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
2212 fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
2214 Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
2215 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
2216 fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
2220 //vPhi *= TMath::RadToDeg();
2222 //=======================================correction
2223 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
2224 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2226 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
2231 return tracksAccepted;
2234 //________________________________________________________________________
2235 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
2236 // Clones TObjArray and returns it with tracks after shuffling the charges
2238 TObjArray* tracksShuffled = new TObjArray;
2239 tracksShuffled->SetOwner(kTRUE);
2241 vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks
2243 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2245 AliVParticle* track = (AliVParticle*) tracks->At(i);
2246 chargeVector->push_back(track->Charge());
2249 random_shuffle(chargeVector->begin(), chargeVector->end());
2251 for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
2252 AliVParticle* track = (AliVParticle*) tracks->At(i);
2253 //==============================correction
2254 Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
2255 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2256 tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
2259 delete chargeVector;
2261 return tracksShuffled;
2264 //________________________________________________________________________
2265 void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,
2266 const char* lhcPeriod) {
2267 //Function to setup the VZERO gain equalization
2268 //============Get the equilization map============//
2269 TFile *calibrationFile = TFile::Open(filename);
2270 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
2271 Printf("No calibration file found!!!");
2275 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
2277 Printf("Calibration TList not found!!!");
2281 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
2282 if(!fHistVZEROAGainEqualizationMap) {
2283 Printf("VZERO-A calibration object not found!!!");
2286 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
2287 if(!fHistVZEROCGainEqualizationMap) {
2288 Printf("VZERO-C calibration object not found!!!");
2292 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
2293 if(!fHistVZEROChannelGainEqualizationMap) {
2294 Printf("VZERO channel calibration object not found!!!");
2299 //________________________________________________________________________
2300 Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run,
2303 if(!fHistVZEROAGainEqualizationMap) return 1.0;
2305 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
2306 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2307 if(gRunNumber == run)
2308 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
2314 //________________________________________________________________________
2315 Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run,
2318 if(!fHistVZEROAGainEqualizationMap) return 1.0;
2320 TString gVZEROSide = side;
2321 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
2322 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2323 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
2324 if(gRunNumber == run) {
2325 if(gVZEROSide == "A")
2326 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
2327 else if(gVZEROSide == "C")
2328 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
2335 //____________________________________________________________________
2336 Bool_t AliAnalysisTaskBFPsi::AcceptEventCentralityWeight(Double_t centrality)
2338 // copied from AliAnalysisTaskPhiCorrelations
2340 // rejects "randomly" events such that the centrality gets flat
2341 // uses fCentralityWeights histogram
2343 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
2345 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
2346 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
2348 Bool_t result = kFALSE;
2349 if (centralityDigits < weight)
2352 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
2357 //________________________________________________________________________
2358 void AliAnalysisTaskBFPsi::FinishTaskOutput(){
2362 AliError("fBalance not available");
2366 if (!fShuffledBalance) {
2367 AliError("fShuffledBalance not available");
2374 //________________________________________________________________________
2375 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
2376 // Draw result to the screen
2377 // Called once at the end of the query
2379 // not implemented ...