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 "AliAODMCParticle.h"
\r
25 #include "AliCollisionGeometry.h"
\r
26 #include "AliGenEventHeader.h"
\r
27 #include "AliMCEventHandler.h"
\r
28 #include "AliMCEvent.h"
\r
29 #include "AliStack.h"
\r
30 #include "AliESDtrackCuts.h"
\r
31 #include "AliEventplane.h"
\r
32 #include "AliTHn.h"
\r
34 #include "AliAnalysisUtils.h"
\r
36 #include "AliEventPoolManager.h"
\r
38 #include "AliPID.h"
\r
39 #include "AliPIDResponse.h"
\r
40 #include "AliPIDCombined.h"
\r
42 #include "AliAnalysisTaskBFPsi.h"
\r
43 #include "AliBalancePsi.h"
\r
44 #include "AliAnalysisTaskTriggeredBF.h"
\r
47 // Analysis task for the BF vs Psi code
\r
48 // Authors: Panos.Christakoglou@nikhef.nl
\r
53 ClassImp(AliAnalysisTaskBFPsi)
\r
55 //________________________________________________________________________
\r
56 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
\r
57 : AliAnalysisTaskSE(name),
\r
58 fDebugLevel(kFALSE),
\r
59 fArrayMC(0), //+++++++++++++
\r
61 fRunShuffling(kFALSE),
\r
62 fShuffledBalance(0),
\r
64 fRunMixingEventPlane(kFALSE),
\r
65 fMixingTracks(50000),
\r
75 fHistCentStatsUsed(0),
\r
76 fHistTriggerStats(0),
\r
81 fHistTPCvsVZEROMultiplicity(0),
\r
82 fHistVZEROSignal(0),
\r
99 fHistdEdxVsPTPCbeforePID(NULL),
\r
100 fHistBetavsPTOFbeforePID(NULL),
\r
101 fHistProbTPCvsPtbeforePID(NULL),
\r
102 fHistProbTOFvsPtbeforePID(NULL),
\r
103 fHistProbTPCTOFvsPtbeforePID(NULL),
\r
104 fHistNSigmaTPCvsPtbeforePID(NULL),
\r
105 fHistNSigmaTOFvsPtbeforePID(NULL),
\r
106 fHistBetaVsdEdXbeforePID(NULL), //+++++++
\r
107 fHistNSigmaTPCTOFvsPtbeforePID(NULL), //++++++
\r
108 fHistdEdxVsPTPCafterPID(NULL),
\r
109 fHistBetavsPTOFafterPID(NULL),
\r
110 fHistProbTPCvsPtafterPID(NULL),
\r
111 fHistProbTOFvsPtafterPID(NULL),
\r
112 fHistProbTPCTOFvsPtafterPID(NULL),
\r
113 fHistNSigmaTPCvsPtafterPID(NULL),
\r
114 fHistNSigmaTOFvsPtafterPID(NULL),
\r
115 fHistBetaVsdEdXafterPID(NULL), //+++++++
\r
116 fHistNSigmaTPCTOFvsPtafterPID(NULL), //+++++++
\r
117 fCentralityArrayBinsForCorrections(kCENTRALITY),
\r
120 fParticleOfInterest(kPion),
\r
121 fPidDetectorConfig(kTPCTOF),
\r
123 fUsePIDnSigma(kTRUE),
\r
124 fUsePIDPropabilities(kFALSE),
\r
126 fMinAcceptedPIDProbability(0.8),
\r
127 fElectronRejection(kFALSE),
\r
128 fElectronOnlyRejection(kFALSE),
\r
129 fElectronRejectionNSigma(-1.),
\r
130 fElectronRejectionMinPt(0.),
\r
131 fElectronRejectionMaxPt(1000.),
\r
133 fCentralityEstimator("V0M"),
\r
134 fUseCentrality(kFALSE),
\r
135 fCentralityPercentileMin(0.),
\r
136 fCentralityPercentileMax(5.),
\r
137 fImpactParameterMin(0.),
\r
138 fImpactParameterMax(20.),
\r
139 fMultiplicityEstimator("V0A"),
\r
140 fUseMultiplicity(kFALSE),
\r
141 fNumberOfAcceptedTracksMin(0),
\r
142 fNumberOfAcceptedTracksMax(10000),
\r
143 fHistNumberOfAcceptedTracks(0),
\r
144 fHistMultiplicity(0),
\r
145 fUseOfflineTrigger(kFALSE),
\r
146 fCheckFirstEventInChunk(kFALSE),
\r
147 fCheckPileUp(kFALSE),
\r
148 fUseMCforKinematics(kFALSE),
\r
152 fnAODtrackCutBit(128),
\r
155 fPtMinForCorrections(0.3),//=================================correction
\r
156 fPtMaxForCorrections(1.5),//=================================correction
\r
157 fPtBinForCorrections(36), //=================================correction
\r
160 fEtaMinForCorrections(-0.8),//=================================correction
\r
161 fEtaMaxForCorrections(0.8),//=================================correction
\r
162 fEtaBinForCorrections(16), //=================================correction
\r
165 fPhiMinForCorrections(0.),//=================================correction
\r
166 fPhiMaxForCorrections(360.),//=================================correction
\r
167 fPhiBinForCorrections(100), //=================================correction
\r
171 fNClustersTPCCut(-1),
\r
172 fAcceptanceParameterization(0),
\r
173 fDifferentialV2(0),
\r
174 fUseFlowAfterBurner(kFALSE),
\r
175 fExcludeResonancesInMC(kFALSE),
\r
176 fExcludeElectronsInMC(kFALSE),
\r
177 fUseMCPdgCode(kFALSE),
\r
178 fPDGCodeToBeAnalyzed(-1),
\r
179 fEventClass("EventPlane"),
\r
180 fCustomBinning(""),
\r
181 fHistVZEROAGainEqualizationMap(0),
\r
182 fHistVZEROCGainEqualizationMap(0),
\r
183 fHistVZEROChannelGainEqualizationMap(0) {
\r
185 // Define input and output slots here
\r
186 // Input slot #0 works with a TChain
\r
188 //======================================================correction
\r
189 for (Int_t i=0; i<kCENTRALITY; i++){
\r
190 fHistCorrectionPlus[i] = NULL;
\r
191 fHistCorrectionMinus[i] = NULL;
\r
192 fCentralityArrayForCorrections[i] = -1.;
\r
194 //=====================================================correction
\r
196 DefineInput(0, TChain::Class());
\r
197 // Output slot #0 writes into a TH1 container
\r
198 DefineOutput(1, TList::Class());
\r
199 DefineOutput(2, TList::Class());
\r
200 DefineOutput(3, TList::Class());
\r
201 DefineOutput(4, TList::Class());
\r
202 DefineOutput(5, TList::Class());
\r
205 //________________________________________________________________________
\r
206 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
\r
208 // delete fBalance;
\r
209 // delete fShuffledBalance;
\r
211 // delete fListBF;
\r
212 // delete fListBFS;
\r
214 // delete fHistEventStats;
\r
215 // delete fHistTrackStats;
\r
216 // delete fHistVx;
\r
217 // delete fHistVy;
\r
218 // delete fHistVz;
\r
220 // delete fHistClus;
\r
221 // delete fHistDCA;
\r
222 // delete fHistChi2;
\r
224 // delete fHistEta;
\r
225 // delete fHistPhi;
\r
226 // delete fHistEtaPhiPos;
\r
227 // delete fHistEtaPhiNeg;
\r
228 // delete fHistV0M;
\r
231 //________________________________________________________________________
\r
232 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
\r
233 // Create histograms
\r
236 // global switch disabling the reference
\r
237 // (to avoid "Replacing existing TH1" if several wagons are created in train)
\r
238 Bool_t oldStatus = TH1::AddDirectoryStatus();
\r
239 TH1::AddDirectory(kFALSE);
\r
242 fBalance = new AliBalancePsi();
\r
243 fBalance->SetAnalysisLevel("ESD");
\r
244 fBalance->SetEventClass(fEventClass);
\r
245 //fBalance->SetNumberOfBins(-1,16);
\r
246 //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
\r
248 if(fRunShuffling) {
\r
249 if(!fShuffledBalance) {
\r
250 fShuffledBalance = new AliBalancePsi();
\r
251 fShuffledBalance->SetAnalysisLevel("ESD");
\r
252 fShuffledBalance->SetEventClass(fEventClass);
\r
253 //fShuffledBalance->SetNumberOfBins(-1,16);
\r
254 //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
\r
258 if(!fMixedBalance) {
\r
259 fMixedBalance = new AliBalancePsi();
\r
260 fMixedBalance->SetAnalysisLevel("ESD");
\r
261 fMixedBalance->SetEventClass(fEventClass);
\r
266 fList = new TList();
\r
267 fList->SetName("listQA");
\r
270 //Balance Function list
\r
271 fListBF = new TList();
\r
272 fListBF->SetName("listBF");
\r
273 fListBF->SetOwner();
\r
275 if(fRunShuffling) {
\r
276 fListBFS = new TList();
\r
277 fListBFS->SetName("listBFShuffled");
\r
278 fListBFS->SetOwner();
\r
282 fListBFM = new TList();
\r
283 fListBFM->SetName("listTriggeredBFMixed");
\r
284 fListBFM->SetOwner();
\r
288 if(fUsePID || fElectronRejection) {
\r
289 fHistListPIDQA = new TList();
\r
290 fHistListPIDQA->SetName("listQAPID");
\r
291 fHistListPIDQA->SetOwner();
\r
295 TString gCutName[7] = {"Total","Offline trigger",
\r
296 "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};
\r
297 fHistEventStats = new TH2F("fHistEventStats",
\r
298 "Event statistics;;Centrality percentile;N_{events}",
\r
299 7,0.5,7.5,220,-5,105);
\r
300 for(Int_t i = 1; i <= 7; i++)
\r
301 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
\r
302 fList->Add(fHistEventStats);
\r
304 TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
\r
305 fHistCentStats = new TH2F("fHistCentStats",
\r
306 "Centrality statistics;;Cent percentile",
\r
307 13,-0.5,12.5,220,-5,105);
\r
308 for(Int_t i = 1; i <= 13; i++){
\r
309 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
\r
310 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
\r
312 fList->Add(fHistCentStats);
\r
314 fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);
\r
315 fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());
\r
316 fList->Add(fHistCentStatsUsed);
\r
318 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
\r
319 fList->Add(fHistTriggerStats);
\r
321 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
\r
322 fList->Add(fHistTrackStats);
\r
324 fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
\r
325 fList->Add(fHistNumberOfAcceptedTracks);
\r
327 fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);
\r
328 fList->Add(fHistMultiplicity);
\r
330 // Vertex distributions
\r
331 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
\r
332 fList->Add(fHistVx);
\r
333 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
\r
334 fList->Add(fHistVy);
\r
335 fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
\r
336 fList->Add(fHistVz);
\r
338 //TPC vs VZERO multiplicity
\r
339 fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",3001,-0.5,30000.5,4001,-0.5,4000.5);
\r
340 if(fMultiplicityEstimator == "V0A")
\r
341 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");
\r
342 else if(fMultiplicityEstimator == "V0C")
\r
343 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");
\r
345 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");
\r
346 fList->Add(fHistTPCvsVZEROMultiplicity);
\r
348 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
\r
349 fList->Add(fHistVZEROSignal);
\r
352 fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
\r
353 fList->Add(fHistEventPlane);
\r
356 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
\r
357 fList->Add(fHistClus);
\r
358 fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);
\r
359 fList->Add(fHistChi2);
\r
360 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
\r
361 fList->Add(fHistDCA);
\r
362 fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);
\r
363 fList->Add(fHistPt);
\r
364 fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);
\r
365 fList->Add(fHistEta);
\r
366 fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);
\r
367 fList->Add(fHistRapidity);
\r
368 fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);
\r
369 fList->Add(fHistPhi);
\r
370 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
371 fList->Add(fHistEtaPhiPos);
\r
372 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
373 fList->Add(fHistEtaPhiNeg);
\r
374 fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
\r
375 fList->Add(fHistPhiBefore);
\r
376 fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
\r
377 fList->Add(fHistPhiAfter);
\r
378 fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
\r
379 fList->Add(fHistPhiPos);
\r
380 fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);
\r
381 fList->Add(fHistPhiNeg);
\r
382 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
\r
383 fList->Add(fHistV0M);
\r
384 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
\r
385 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
\r
386 for(Int_t i = 1; i <= 6; i++)
\r
387 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
\r
388 fList->Add(fHistRefTracks);
\r
390 // Balance function histograms
\r
391 // Initialize histograms if not done yet (including the custom binning)
\r
392 if(!fBalance->GetHistNp()){
\r
393 AliInfo("Histograms not yet initialized! --> Will be done now");
\r
394 fBalance->SetCustomBinning(fCustomBinning);
\r
395 fBalance->InitHistograms();
\r
398 if(fRunShuffling) {
\r
399 if(!fShuffledBalance->GetHistNp()) {
\r
400 AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");
\r
401 fShuffledBalance->SetCustomBinning(fCustomBinning);
\r
402 fShuffledBalance->InitHistograms();
\r
407 if(!fMixedBalance->GetHistNp()) {
\r
408 AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");
\r
409 fMixedBalance->SetCustomBinning(fCustomBinning);
\r
410 fMixedBalance->InitHistograms();
\r
414 // QA histograms for different cuts
\r
415 fList->Add(fBalance->GetQAHistHBTbefore());
\r
416 fList->Add(fBalance->GetQAHistHBTafter());
\r
417 fList->Add(fBalance->GetQAHistConversionbefore());
\r
418 fList->Add(fBalance->GetQAHistConversionafter());
\r
419 fList->Add(fBalance->GetQAHistPsiMinusPhi());
\r
420 fList->Add(fBalance->GetQAHistResonancesBefore());
\r
421 fList->Add(fBalance->GetQAHistResonancesRho());
\r
422 fList->Add(fBalance->GetQAHistResonancesK0());
\r
423 fList->Add(fBalance->GetQAHistResonancesLambda());
\r
424 fList->Add(fBalance->GetQAHistQbefore());
\r
425 fList->Add(fBalance->GetQAHistQafter());
\r
427 //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
\r
428 fListBF->Add(fBalance->GetHistNp());
\r
429 fListBF->Add(fBalance->GetHistNn());
\r
430 fListBF->Add(fBalance->GetHistNpn());
\r
431 fListBF->Add(fBalance->GetHistNnn());
\r
432 fListBF->Add(fBalance->GetHistNpp());
\r
433 fListBF->Add(fBalance->GetHistNnp());
\r
435 if(fRunShuffling) {
\r
436 fListBFS->Add(fShuffledBalance->GetHistNp());
\r
437 fListBFS->Add(fShuffledBalance->GetHistNn());
\r
438 fListBFS->Add(fShuffledBalance->GetHistNpn());
\r
439 fListBFS->Add(fShuffledBalance->GetHistNnn());
\r
440 fListBFS->Add(fShuffledBalance->GetHistNpp());
\r
441 fListBFS->Add(fShuffledBalance->GetHistNnp());
\r
445 fListBFM->Add(fMixedBalance->GetHistNp());
\r
446 fListBFM->Add(fMixedBalance->GetHistNn());
\r
447 fListBFM->Add(fMixedBalance->GetHistNpn());
\r
448 fListBFM->Add(fMixedBalance->GetHistNnn());
\r
449 fListBFM->Add(fMixedBalance->GetHistNpp());
\r
450 fListBFM->Add(fMixedBalance->GetHistNnp());
\r
457 Int_t trackDepth = fMixingTracks;
\r
458 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
\r
461 Double_t* centbins;
\r
462 Int_t nCentralityBins;
\r
463 if(fBalance->IsUseVertexBinning()){
\r
464 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);
\r
467 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);
\r
470 // multiplicity bins
\r
471 Double_t* multbins;
\r
472 Int_t nMultiplicityBins;
\r
473 multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);
\r
476 Double_t* vtxbins;
\r
478 if(fBalance->IsUseVertexBinning()){
\r
479 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);
\r
482 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);
\r
485 // Event plane angle (Psi) bins
\r
488 psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);
\r
491 // run the event mixing also in bins of event plane (statistics!)
\r
492 if(fRunMixingEventPlane){
\r
493 if(fEventClass=="Multiplicity"){
\r
494 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
\r
497 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
\r
501 if(fEventClass=="Multiplicity"){
\r
502 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
\r
505 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
\r
510 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
\r
512 //====================PID========================//
\r
514 fPIDCombined = new AliPIDCombined();
\r
515 fPIDCombined->SetDefaultTPCPriors();
\r
517 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
\r
518 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
\r
520 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
\r
521 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID);
\r
523 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
\r
524 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID);
\r
526 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
\r
527 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);
\r
529 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
\r
530 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);
\r
532 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
\r
533 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);
\r
535 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
\r
536 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID);
\r
538 fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2);
\r
539 fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++
\r
541 fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500);
\r
542 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++
\r
544 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
\r
545 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
\r
547 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
\r
548 fHistListPIDQA->Add(fHistBetavsPTOFafterPID);
\r
550 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
\r
551 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);
\r
553 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
\r
554 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID);
\r
556 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
\r
557 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);
\r
559 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
\r
560 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
\r
562 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
\r
563 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);
\r
565 fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2);
\r
566 fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++
\r
568 fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500);
\r
569 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++
\r
572 // for electron rejection only TPC nsigma histograms
\r
573 if(!fUsePID && fElectronRejection) {
\r
575 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
\r
576 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
\r
578 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
\r
579 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition
\r
581 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
\r
582 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
\r
584 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
\r
585 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
\r
587 //====================PID========================//
\r
589 // Post output data.
\r
590 PostData(1, fList);
\r
591 PostData(2, fListBF);
\r
592 if(fRunShuffling) PostData(3, fListBFS);
\r
593 if(fRunMixing) PostData(4, fListBFM);
\r
594 if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID
\r
596 AliInfo("Finished setting up the Output");
\r
598 TH1::AddDirectory(oldStatus);
\r
602 //________________________________________________________________________
\r
603 void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename,
\r
604 Int_t nCentralityBins,
\r
605 Double_t *centralityArrayForCorrections) {
\r
606 //Open files that will be used for correction
\r
607 fCentralityArrayBinsForCorrections = nCentralityBins;
\r
608 for (Int_t i=0; i<nCentralityBins; i++)
\r
609 fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];
\r
611 // No file specified -> run without corrections
\r
612 if(!filename.Contains(".root")) {
\r
613 AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));
\r
617 //Open the input file
\r
618 TFile *f = TFile::Open(filename);
\r
620 AliInfo(Form("File %s not found --> run without corrections",filename.Data()));
\r
624 //TString listEffName = "";
\r
625 for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {
\r
626 //Printf("iCent %d:",iCent);
\r
627 TString histoName = "fHistCorrectionPlus";
\r
628 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
\r
629 fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));
\r
630 if(!fHistCorrectionPlus[iCent]) {
\r
631 AliError(Form("fHist %s not found!!!",histoName.Data()));
\r
635 histoName = "fHistCorrectionMinus";
\r
636 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
\r
637 fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data()));
\r
638 if(!fHistCorrectionMinus[iCent]) {
\r
639 AliError(Form("fHist %s not found!!!",histoName.Data()));
\r
642 }//loop over centralities: ONLY the PbPb case is covered
\r
644 if(fHistCorrectionPlus[0]){
\r
645 fEtaMinForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmin();
\r
646 fEtaMaxForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmax();
\r
647 fEtaBinForCorrections = fHistCorrectionPlus[0]->GetNbinsX();
\r
649 fPtMinForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmin();
\r
650 fPtMaxForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmax();
\r
651 fPtBinForCorrections = fHistCorrectionPlus[0]->GetNbinsY();
\r
653 fPhiMinForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmin();
\r
654 fPhiMaxForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmax();
\r
655 fPhiBinForCorrections = fHistCorrectionPlus[0]->GetNbinsZ();
\r
659 //________________________________________________________________________
\r
660 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
\r
662 // Called for each event
\r
664 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
665 Int_t gNumberOfAcceptedTracks = 0;
\r
666 Double_t lMultiplicityVar = -999.; //-1
\r
667 Double_t gReactionPlane = -1.;
\r
668 Float_t bSign = 0.;
\r
670 // get the event (for generator level: MCEvent())
\r
671 AliVEvent* eventMain = NULL;
\r
672 if(gAnalysisLevel == "MC") {
\r
673 eventMain = dynamic_cast<AliVEvent*>(MCEvent());
\r
676 eventMain = dynamic_cast<AliVEvent*>(InputEvent());
\r
677 // for HBT like cuts need magnetic field sign
\r
678 bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
\r
681 AliError("eventMain not available");
\r
685 // PID Response task active?
\r
686 if(fUsePID || fElectronRejection) {
\r
687 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
\r
688 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
\r
691 // check event cuts and fill event histograms
\r
692 if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){
\r
695 // get the reaction plane
\r
696 if(fEventClass != "Multiplicity") {
\r
697 gReactionPlane = GetEventPlane(eventMain);
\r
698 fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
\r
699 if(gReactionPlane < 0){
\r
704 // get the accepted tracks in main event
\r
705 TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
\r
706 gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
\r
708 //multiplicity cut (used in pp)
\r
709 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
\r
711 // store charges of all accepted tracks,shuffle and reassign(two extra loops)
\r
712 TObjArray* tracksShuffled = NULL;
\r
714 tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
\r
720 // 1. First get an event pool corresponding in mult (cent) and
\r
721 // zvertex to the current event. Once initialized, the pool
\r
722 // should contain nMix (reduced) events. This routine does not
\r
723 // pre-scan the chain. The first several events of every chain
\r
724 // will be skipped until the needed pools are filled to the
\r
725 // specified depth. If the pool categories are not too rare, this
\r
726 // should not be a problem. If they are rare, you could lose`
\r
729 // 2. Collect the whole pool's content of tracks into one TObjArray
\r
730 // (bgTracks), which is effectively a single background super-event.
\r
732 // 3. The reduced and bgTracks arrays must both be passed into
\r
733 // FillCorrelations(). Also nMix should be passed in, so a weight
\r
734 // of 1./nMix can be applied.
\r
736 AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
\r
739 AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
\r
743 //pool->SetDebug(1);
\r
745 if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){
\r
748 Int_t nMix = pool->GetCurrentNEvents();
\r
749 //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
\r
751 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
\r
752 //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
\r
753 //if (pool->IsReady())
\r
754 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
\r
756 // Fill mixed-event histos here
\r
757 for (Int_t jMix=0; jMix<nMix; jMix++)
\r
759 TObjArray* tracksMixed = pool->GetEvent(jMix);
\r
760 fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
\r
764 // Update the Event pool
\r
765 pool->UpdatePool(tracksMain);
\r
766 //pool->PrintInfo();
\r
768 }//pool NULL check
\r
771 // calculate balance function
\r
772 fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
\r
774 // calculate shuffled balance function
\r
775 if(fRunShuffling && tracksShuffled != NULL) {
\r
776 fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
\r
780 //________________________________________________________________________
\r
781 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
\r
782 // Checks the Event cuts
\r
783 // Fills Event statistics histograms
\r
785 Bool_t isSelectedMain = kTRUE;
\r
786 Float_t gRefMultiplicity = -1.;
\r
787 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
789 AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
\r
791 fHistEventStats->Fill(1,gRefMultiplicity); //all events
\r
793 // check first event in chunk (is not needed for new reconstructions)
\r
794 if(fCheckFirstEventInChunk){
\r
795 AliAnalysisUtils ut;
\r
796 if(ut.IsFirstEventInChunk(event))
\r
798 fHistEventStats->Fill(6,gRefMultiplicity);
\r
800 // check for pile-up event
\r
802 AliAnalysisUtils ut;
\r
803 ut.SetUseMVPlpSelection(kTRUE);
\r
804 ut.SetUseOutOfBunchPileUp(kTRUE);
\r
805 if(ut.IsPileUpEvent(event))
\r
807 fHistEventStats->Fill(7,gRefMultiplicity);
\r
810 // Event trigger bits
\r
811 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
\r
812 if(fUseOfflineTrigger)
\r
813 isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
\r
815 if(isSelectedMain) {
\r
816 fHistEventStats->Fill(2,gRefMultiplicity); //triggered events
\r
819 if(gAnalysisLevel == "MC") {
\r
821 AliError("mcEvent not available");
\r
826 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
\r
828 TArrayF gVertexArray;
\r
829 header->PrimaryVertex(gVertexArray);
\r
830 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
\r
831 //gVertexArray.At(0),
\r
832 //gVertexArray.At(1),
\r
833 //gVertexArray.At(2));
\r
834 fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex
\r
835 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
\r
836 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
\r
837 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
\r
838 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
\r
840 // get the reference multiplicty or centrality
\r
841 gRefMultiplicity = GetRefMultiOrCentrality(event);
\r
843 fHistVx->Fill(gVertexArray.At(0));
\r
844 fHistVy->Fill(gVertexArray.At(1));
\r
845 fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);
\r
847 // take only events inside centrality class
\r
848 if(fUseCentrality) {
\r
849 if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){
\r
850 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
\r
851 return gRefMultiplicity;
\r
852 }//centrality class
\r
854 // take events only within the same multiplicity class
\r
855 else if(fUseMultiplicity) {
\r
856 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
\r
857 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
\r
858 return gRefMultiplicity;
\r
860 }//multiplicity range
\r
868 // Event Vertex AOD, ESD, ESDMC
\r
870 const AliVVertex *vertex = event->GetPrimaryVertex();
\r
873 Double32_t fCov[6];
\r
874 vertex->GetCovarianceMatrix(fCov);
\r
875 if(vertex->GetNContributors() > 0) {
\r
877 fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex
\r
878 if(TMath::Abs(vertex->GetX()) < fVxMax) {
\r
879 if(TMath::Abs(vertex->GetY()) < fVyMax) {
\r
880 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
\r
881 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
\r
883 // get the reference multiplicty or centrality
\r
884 gRefMultiplicity = GetRefMultiOrCentrality(event);
\r
886 fHistVx->Fill(vertex->GetX());
\r
887 fHistVy->Fill(vertex->GetY());
\r
888 fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);
\r
890 // take only events inside centrality class
\r
891 if(fUseCentrality) {
\r
892 if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){
\r
893 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
\r
894 return gRefMultiplicity;
\r
895 }//centrality class
\r
897 // take events only within the same multiplicity class
\r
898 else if(fUseMultiplicity) {
\r
900 //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,
\r
901 //fNumberOfAcceptedTracksMax,gRefMultiplicity);
\r
903 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
\r
904 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
\r
905 return gRefMultiplicity;
\r
907 }//multiplicity range
\r
911 }//proper vertex resolution
\r
912 }//proper number of contributors
\r
913 }//vertex object valid
\r
914 }//triggered event
\r
917 // in all other cases return -1 (event not accepted)
\r
922 //________________________________________________________________________
\r
923 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
\r
924 // Checks the Event cuts
\r
925 // Fills Event statistics histograms
\r
927 Float_t gCentrality = -1.;
\r
928 Double_t gMultiplicity = -1.;
\r
929 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
932 // calculate centrality always (not only in centrality mode)
\r
933 if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++
\r
934 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
\r
936 gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
\r
938 // QA for centrality estimators
\r
939 fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
\r
940 fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));
\r
941 fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));
\r
942 fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
\r
943 fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
\r
944 fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL"));
\r
945 fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
\r
946 fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
\r
947 fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));
\r
948 fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));
\r
949 fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
\r
950 fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
\r
951 fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
\r
953 // Centrality estimator USED ++++++++++++++++++++++++++++++
\r
954 fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));
\r
956 // centrality QA (V0M)
\r
957 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
\r
959 // centrality QA (reference tracks)
\r
960 fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
\r
961 fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
\r
962 fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
\r
963 fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
\r
964 fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
\r
965 fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
\r
966 fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
\r
967 fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
\r
968 fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
\r
973 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
\r
974 AliCentrality *centrality = event->GetCentrality();
\r
975 gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
\r
977 // QA for centrality estimators
\r
978 fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));
\r
979 fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));
\r
980 fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));
\r
981 fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));
\r
982 fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));
\r
983 fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));
\r
984 fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));
\r
985 fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));
\r
986 fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));
\r
987 fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));
\r
988 fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));
\r
989 fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));
\r
990 fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));
\r
992 // Centrality estimator USED ++++++++++++++++++++++++++++++
\r
993 fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
\r
995 // centrality QA (V0M)
\r
996 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
\r
999 else if(gAnalysisLevel == "MC"){
\r
1000 Double_t gImpactParameter = 0.;
\r
1001 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
1003 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
\r
1005 gImpactParameter = headerH->ImpactParameter();
\r
1006 gCentrality = gImpactParameter;
\r
1012 gCentrality = -1.;
\r
1015 // calculate reference multiplicity always (not only in multiplicity mode)
\r
1016 if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){
\r
1017 AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
\r
1019 gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
\r
1020 fHistMultiplicity->Fill(gMultiplicity);
\r
1021 }//AliESDevent cast
\r
1024 else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){
\r
1025 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
\r
1026 if ((fMultiplicityEstimator == "V0M")||
\r
1027 (fMultiplicityEstimator == "V0A")||
\r
1028 (fMultiplicityEstimator == "V0C") ||
\r
1029 (fMultiplicityEstimator == "TPC")) {
\r
1030 gMultiplicity = GetReferenceMultiplicityFromAOD(event);
\r
1031 if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);
\r
1035 gMultiplicity = header->GetRefMultiplicity();
\r
1036 if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);
\r
1038 fHistMultiplicity->Fill(gMultiplicity);
\r
1040 else if(gAnalysisLevel == "MC") {
\r
1041 AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
1042 //Calculating the multiplicity as the number of charged primaries
\r
1043 //within \pm 0.8 in eta and pT > 0.1 GeV/c
\r
1044 for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
\r
1045 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
\r
1047 AliError(Form("Could not receive particle %d", iParticle));
\r
1051 //exclude non stable particles
\r
1052 if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
\r
1054 //++++++++++++++++
\r
1055 if (fMultiplicityEstimator == "V0M") {
\r
1056 if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7))
\r
1058 else if (fMultiplicityEstimator == "V0A") {
\r
1059 if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}
\r
1060 else if (fMultiplicityEstimator == "V0C") {
\r
1061 if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}
\r
1062 else if (fMultiplicityEstimator == "TPC") {
\r
1063 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
\r
1064 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
\r
1067 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
\r
1068 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
\r
1070 //++++++++++++++++
\r
1072 if(track->Charge() == 0) continue;
\r
1074 gMultiplicity += 1;
\r
1075 }//loop over primaries
\r
1076 fHistMultiplicity->Fill(gMultiplicity);
\r
1079 gMultiplicity = -1;
\r
1083 // decide what should be returned only here
\r
1084 Double_t lReturnVal = -100;
\r
1085 if(fEventClass=="Multiplicity"){
\r
1086 lReturnVal = gMultiplicity;
\r
1087 }else if(fEventClass=="Centrality"){
\r
1088 lReturnVal = gCentrality;
\r
1090 return lReturnVal;
\r
1093 //________________________________________________________________________
\r
1094 Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){
\r
1095 //Function that returns the reference multiplicity from AODs (data or reco MC)
\r
1096 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
\r
1097 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
\r
1098 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
\r
1100 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
\r
1102 Printf("ERROR: AOD header not available");
\r
1105 Int_t gRunNumber = header->GetRunNumber();
\r
1107 // Loop over tracks in event
\r
1108 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1109 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
\r
1111 AliError(Form("Could not receive track %d", iTracks));
\r
1116 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
\r
1118 if(aodTrack->Charge() == 0) continue;
\r
1119 // Kinematics cuts from ESD track cuts
\r
1120 if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;
\r
1121 if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;
\r
1123 //=================PID (so far only for electron rejection)==========================//
\r
1124 if(fElectronRejection) {
\r
1125 // get the electron nsigma
\r
1126 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
\r
1128 // check only for given momentum range
\r
1129 if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){
\r
1130 //look only at electron nsigma
\r
1131 if(!fElectronOnlyRejection) {
\r
1132 //Make the decision based on the n-sigma of electrons
\r
1133 if(nSigma < fElectronRejectionNSigma) continue;
\r
1136 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
\r
1137 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
\r
1138 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
\r
1140 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
\r
1141 if(nSigma < fElectronRejectionNSigma
\r
1142 && nSigmaPions > fElectronRejectionNSigma
\r
1143 && nSigmaKaons > fElectronRejectionNSigma
\r
1144 && nSigmaProtons > fElectronRejectionNSigma ) continue;
\r
1147 }//electron rejection
\r
1149 gRefMultiplicityTPC += 1.0;
\r
1152 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
\r
1153 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
\r
1154 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
\r
1156 if(iChannel < 32)
\r
1157 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
\r
1158 else if(iChannel >= 32)
\r
1159 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
\r
1162 //Equalization of gain
\r
1163 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
\r
1165 gRefMultiplicityVZEROA /= gFactorA;
\r
1166 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
\r
1168 gRefMultiplicityVZEROC /= gFactorC;
\r
1169 if((gFactorA != 0)&&(gFactorC != 0))
\r
1170 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
\r
1173 Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);
\r
1175 fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
\r
1177 if(fMultiplicityEstimator == "TPC")
\r
1178 gRefMultiplicity = gRefMultiplicityTPC;
\r
1179 else if(fMultiplicityEstimator == "V0M")
\r
1180 gRefMultiplicity = gRefMultiplicityVZERO;
\r
1181 else if(fMultiplicityEstimator == "V0A")
\r
1182 gRefMultiplicity = gRefMultiplicityVZEROA;
\r
1183 else if(fMultiplicityEstimator == "V0C")
\r
1184 gRefMultiplicity = gRefMultiplicityVZEROC;
\r
1186 return gRefMultiplicity;
\r
1189 //________________________________________________________________________
\r
1190 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
\r
1191 // Get the event plane
\r
1193 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
1195 Float_t gVZEROEventPlane = -10.;
\r
1196 Float_t gReactionPlane = -10.;
\r
1197 Double_t qxTot = 0.0, qyTot = 0.0;
\r
1199 //MC: from reaction plane
\r
1200 if(gAnalysisLevel == "MC"){
\r
1202 AliError("mcEvent not available");
\r
1206 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
1208 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
\r
1210 gReactionPlane = headerH->ReactionPlaneAngle();
\r
1211 //gReactionPlane *= TMath::RadToDeg();
\r
1216 // AOD,ESD,ESDMC: from VZERO Event Plane
\r
1219 AliEventplane *ep = event->GetEventplane();
\r
1221 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
\r
1222 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
\r
1223 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
\r
1224 gReactionPlane = gVZEROEventPlane;
\r
1228 return gReactionPlane;
\r
1231 //________________________________________________________________________
\r
1232 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta,
\r
1236 Double_t gCentrality) {
\r
1237 // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality)
\r
1239 Double_t correction = 1.;
\r
1240 Int_t binEta = 0, binPt = 0, binPhi = 0;
\r
1242 //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);
\r
1243 if(fEtaBinForCorrections != 0) {
\r
1244 Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;
\r
1245 if(fEtaMaxForCorrections != fEtaMinForCorrections)
\r
1246 binEta = (Int_t)((vEta-fEtaMinForCorrections)/widthEta)+1;
\r
1249 if(fPtBinForCorrections != 0) {
\r
1250 Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;
\r
1251 if(fPtMaxForCorrections != fPtMinForCorrections)
\r
1252 binPt = (Int_t)((vPt-fPtMinForCorrections)/widthPt) + 1;
\r
1255 if(fPhiBinForCorrections != 0) {
\r
1256 Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;
\r
1257 if(fPhiMaxForCorrections != fPhiMinForCorrections)
\r
1258 binPhi = (Int_t)((vPhi-fPhiMinForCorrections)/widthPhi)+ 1;
\r
1261 Int_t gCentralityInt = -1;
\r
1262 for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){
\r
1263 if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){
\r
1264 gCentralityInt = i;
\r
1269 // centrality not in array --> no correction
\r
1270 if(gCentralityInt < 0){
\r
1275 //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);
\r
1277 if(fHistCorrectionPlus[gCentralityInt]){
\r
1278 if (vCharge > 0) {
\r
1279 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));
\r
1280 //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);
\r
1282 if (vCharge < 0) {
\r
1283 correction = fHistCorrectionMinus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));
\r
1284 //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);
\r
1290 }//centrality in array
\r
1292 if (correction == 0.) {
\r
1293 AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt));
\r
1297 return correction;
\r
1300 //________________________________________________________________________
\r
1301 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){
\r
1302 // Returns TObjArray with tracks after all track cuts (only for AOD!)
\r
1303 // Fills QA histograms
\r
1305 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
1307 //output TObjArray holding all good tracks
\r
1308 TObjArray* tracksAccepted = new TObjArray;
\r
1309 tracksAccepted->SetOwner(kTRUE);
\r
1318 if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
\r
1319 // Loop over tracks in event
\r
1321 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1322 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
\r
1324 AliError(Form("Could not receive track %d", iTracks));
\r
1330 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
\r
1331 // take only TPC only tracks
\r
1332 //fHistTrackStats->Fill(aodTrack->GetFilterMap());
\r
1333 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
\r
1334 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
\r
1337 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
\r
1339 vCharge = aodTrack->Charge();
\r
1340 vEta = aodTrack->Eta();
\r
1341 vY = aodTrack->Y();
\r
1342 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
\r
1343 vPt = aodTrack->Pt();
\r
1345 //===========================PID (so far only for electron rejection)===============================//
\r
1346 if(fElectronRejection) {
\r
1348 // get the electron nsigma
\r
1349 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
\r
1351 //Fill QA before the PID
\r
1352 fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1353 fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
\r
1354 //end of QA-before pid
\r
1356 // check only for given momentum range
\r
1357 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
\r
1359 //look only at electron nsigma
\r
1360 if(!fElectronOnlyRejection){
\r
1362 //Make the decision based on the n-sigma of electrons
\r
1363 if(nSigma < fElectronRejectionNSigma) continue;
\r
1367 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
\r
1368 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
\r
1369 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
\r
1371 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
\r
1372 if(nSigma < fElectronRejectionNSigma
\r
1373 && nSigmaPions > fElectronRejectionNSigma
\r
1374 && nSigmaKaons > fElectronRejectionNSigma
\r
1375 && nSigmaProtons > fElectronRejectionNSigma ) continue;
\r
1379 //Fill QA after the PID
\r
1380 fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1381 fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
\r
1384 //===========================end of PID (so far only for electron rejection)===============================//
\r
1386 //+++++++++++++++++++++++++++++//
\r
1387 //===========================PID===============================//
\r
1389 Double_t prob[AliPID::kSPECIES]={0.};
\r
1390 Double_t probTPC[AliPID::kSPECIES]={0.};
\r
1391 Double_t probTOF[AliPID::kSPECIES]={0.};
\r
1392 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
\r
1394 Double_t nSigma = 0.;
\r
1395 Double_t nSigmaTPC = 0.; //++++
\r
1396 Double_t nSigmaTOF = 0.; //++++
\r
1397 UInt_t detUsedTPC = 0;
\r
1398 UInt_t detUsedTOF = 0;
\r
1399 UInt_t detUsedTPCTOF = 0;
\r
1401 //Decide what detector configuration we want to use
\r
1402 switch(fPidDetectorConfig) {
\r
1404 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
\r
1405 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
\r
1406 detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
\r
1407 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1408 prob[iSpecies] = probTPC[iSpecies];
\r
1411 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
\r
1412 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
\r
1413 detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
\r
1414 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1415 prob[iSpecies] = probTOF[iSpecies];
\r
1418 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
\r
1419 nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest)); //++++++
\r
1420 nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest)); //++++++
\r
1421 nSigma = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);//++++++
\r
1422 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);
\r
1423 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1424 prob[iSpecies] = probTPCTOF[iSpecies];
\r
1428 }//end switch: define detector mask
\r
1430 //Filling the PID QA
\r
1431 Double_t tofTime = -999., length = 999., tof = -999.;
\r
1432 Double_t c = TMath::C()*1.E-9;// m/ns
\r
1433 Double_t beta = -999.;
\r
1434 Double_t nSigmaTOFForParticleOfInterest = -999.;
\r
1435 if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) &&
\r
1436 (aodTrack->IsOn(AliAODTrack::kTIME)) ) {
\r
1437 tofTime = aodTrack->GetTOFsignal();//in ps
\r
1438 length = aodTrack->GetIntegratedLength();
\r
1439 tof = tofTime*1E-3; // ns
\r
1442 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
\r
1446 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
\r
1450 length = length*0.01; // in meters
\r
1452 beta = length/tof;
\r
1454 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest);
\r
1455 fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
\r
1456 fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
\r
1457 fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOFForParticleOfInterest);
\r
1460 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest);
\r
1461 fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1462 fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
\r
1463 fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCForParticleOfInterest);
\r
1464 fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
\r
1466 //combined TPC&TOF
\r
1467 fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
\r
1468 Double_t nSigmaTPCTOFForParticleOfInterest = -999.;//++++++++
\r
1469 nSigmaTPCTOFForParticleOfInterest = TMath::Sqrt(nSigmaTPCForParticleOfInterest*nSigmaTPCForParticleOfInterest + nSigmaTOFForParticleOfInterest*nSigmaTOFForParticleOfInterest);//++++++++
\r
1470 fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOFForParticleOfInterest); //++++++++
\r
1472 //Printf("nSigma %lf",nSigma); //++++++++++++
\r
1473 //Printf("nSigmaTPCTOF %lf",nSigmaTPCTOFForParticleOfInterest); //++++++++++++
\r
1474 //end of QA-before pid
\r
1476 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
\r
1477 //Make the decision based on the n-sigma
\r
1478 if(fUsePIDnSigma) {
\r
1479 if(nSigma > fPIDNSigma) continue;
\r
1481 fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigma);
\r
1482 fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigma);
\r
1483 fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigma); //++++++++
\r
1485 //Make the decision based on the bayesian
\r
1486 else if(fUsePIDPropabilities) {
\r
1487 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
\r
1488 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
\r
1490 fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
\r
1491 fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
\r
1492 fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
\r
1496 //Printf("nSigmaAFTER %lf", nSigma); //++++++++++++
\r
1497 //Printf("nSigmaTPCTOFAFTER %lf", nSigmaTPCTOFForParticleOfInterest); //++++++++++++
\r
1498 //Fill QA after the PID
\r
1499 fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
\r
1500 fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1501 fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
\r
1504 //===========================PID===============================//
\r
1505 //+++++++++++++++++++++++++++++//
\r
1508 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1509 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1512 // Kinematics cuts from ESD track cuts
\r
1513 if( vPt < fPtMin || vPt > fPtMax) continue;
\r
1514 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
1516 // Extra DCA cuts (for systematic studies [!= -1])
\r
1517 if( fDCAxyCut != -1 && fDCAzCut != -1){
\r
1518 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
\r
1519 continue; // 2D cut
\r
1523 // Extra TPC cuts (for systematic studies [!= -1])
\r
1524 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
\r
1527 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
\r
1531 // fill QA histograms
\r
1532 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
\r
1533 fHistDCA->Fill(dcaZ,dcaXY);
\r
1534 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
\r
1535 fHistPt->Fill(vPt,gCentrality);
\r
1536 fHistEta->Fill(vEta,gCentrality);
\r
1537 fHistRapidity->Fill(vY,gCentrality);
\r
1538 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
1539 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
1540 fHistPhi->Fill(vPhi,gCentrality);
\r
1541 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
1542 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
1544 //=======================================correction
\r
1545 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
1546 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
1548 // add the track to the TObjArray
\r
1549 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
1553 //==============================================================================================================
\r
1554 else if(gAnalysisLevel == "MCAOD") {
\r
1556 AliMCEvent* mcEvent = MCEvent();
\r
1558 AliError("ERROR: Could not retrieve MC event");
\r
1562 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
\r
1563 AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
\r
1565 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
\r
1569 if(!aodTrack->IsPhysicalPrimary()) continue;
\r
1571 vCharge = aodTrack->Charge();
\r
1572 vEta = aodTrack->Eta();
\r
1573 vY = aodTrack->Y();
\r
1574 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
\r
1575 vPt = aodTrack->Pt();
\r
1577 // Kinematics cuts from ESD track cuts
\r
1578 if( vPt < fPtMin || vPt > fPtMax) continue;
\r
1579 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
1581 // Remove neutral tracks
\r
1582 if( vCharge == 0 ) continue;
\r
1584 //Exclude resonances
\r
1585 if(fExcludeResonancesInMC) {
\r
1587 Bool_t kExcludeParticle = kFALSE;
\r
1588 Int_t gMotherIndex = aodTrack->GetMother();
\r
1589 if(gMotherIndex != -1) {
\r
1590 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
\r
1592 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
\r
1593 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
\r
1594 //if(pdgCodeOfMother == 113) {
\r
1595 if(pdgCodeOfMother == 113 // rho0
\r
1596 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
\r
1597 // || pdgCodeOfMother == 221 // eta
\r
1598 // || pdgCodeOfMother == 331 // eta'
\r
1599 // || pdgCodeOfMother == 223 // omega
\r
1600 // || pdgCodeOfMother == 333 // phi
\r
1601 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
\r
1602 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
\r
1603 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
\r
1604 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
\r
1605 || pdgCodeOfMother == 111 // pi0 Dalitz
\r
1606 || pdgCodeOfMother == 22 // photon
\r
1608 kExcludeParticle = kTRUE;
\r
1613 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
\r
1614 if(kExcludeParticle) continue;
\r
1617 //Exclude electrons with PDG
\r
1618 if(fExcludeElectronsInMC) {
\r
1620 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
\r
1624 // fill QA histograms
\r
1625 fHistPt->Fill(vPt,gCentrality);
\r
1626 fHistEta->Fill(vEta,gCentrality);
\r
1627 fHistRapidity->Fill(vY,gCentrality);
\r
1628 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
1629 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
1630 fHistPhi->Fill(vPhi,gCentrality);
\r
1631 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
1632 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
1634 //=======================================correction
\r
1635 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
1636 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
1638 // add the track to the TObjArray
\r
1639 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
1643 //==============================================================================================================
\r
1645 //==============================================================================================================
\r
1646 else if(gAnalysisLevel == "MCAODrec") {
\r
1648 /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
\r
1650 printf("ERROR: fAOD not available\n");
\r
1654 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
\r
1656 AliError("No array of MC particles found !!!");
\r
1659 AliMCEvent* mcEvent = MCEvent();
\r
1661 AliError("ERROR: Could not retrieve MC event");
\r
1662 return tracksAccepted;
\r
1665 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1666 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
\r
1668 AliError(Form("Could not receive track %d", iTracks));
\r
1672 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
\r
1673 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
\r
1675 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
\r
1677 vCharge = aodTrack->Charge();
\r
1678 vEta = aodTrack->Eta();
\r
1679 vY = aodTrack->Y();
\r
1680 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
\r
1681 vPt = aodTrack->Pt();
\r
1683 //===========================use MC information for Kinematics===============================//
\r
1684 if(fUseMCforKinematics){
\r
1686 Int_t label = TMath::Abs(aodTrack->GetLabel());
\r
1687 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
\r
1690 vCharge = AODmcTrack->Charge();
\r
1691 vEta = AODmcTrack->Eta();
\r
1692 vY = AODmcTrack->Y();
\r
1693 vPhi = AODmcTrack->Phi();// * TMath::RadToDeg();
\r
1694 vPt = AODmcTrack->Pt();
\r
1697 AliDebug(1, "no MC particle for this track");
\r
1701 //===========================end of use MC information for Kinematics========================//
\r
1704 //===========================PID (so far only for electron rejection)===============================//
\r
1705 if(fElectronRejection) {
\r
1707 // get the electron nsigma
\r
1708 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
\r
1710 //Fill QA before the PID
\r
1711 fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1712 fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
\r
1713 //end of QA-before pid
\r
1715 // check only for given momentum range
\r
1716 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
\r
1718 //look only at electron nsigma
\r
1719 if(!fElectronOnlyRejection){
\r
1721 //Make the decision based on the n-sigma of electrons
\r
1722 if(nSigma < fElectronRejectionNSigma) continue;
\r
1726 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
\r
1727 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
\r
1728 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
\r
1730 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
\r
1731 if(nSigma < fElectronRejectionNSigma
\r
1732 && nSigmaPions > fElectronRejectionNSigma
\r
1733 && nSigmaKaons > fElectronRejectionNSigma
\r
1734 && nSigmaProtons > fElectronRejectionNSigma ) continue;
\r
1738 //Fill QA after the PID
\r
1739 fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1740 fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
\r
1743 //===========================end of PID (so far only for electron rejection)===============================//
\r
1745 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1746 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1748 // Kinematics cuts from ESD track cuts
\r
1749 if( vPt < fPtMin || vPt > fPtMax) continue;
\r
1750 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
1752 // Extra DCA cuts (for systematic studies [!= -1])
\r
1753 if( fDCAxyCut != -1 && fDCAzCut != -1){
\r
1754 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
\r
1755 continue; // 2D cut
\r
1759 // Extra TPC cuts (for systematic studies [!= -1])
\r
1760 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
\r
1763 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
\r
1767 //Exclude resonances
\r
1768 if(fExcludeResonancesInMC) {
\r
1770 Bool_t kExcludeParticle = kFALSE;
\r
1772 Int_t label = TMath::Abs(aodTrack->GetLabel());
\r
1773 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
\r
1776 //if (AODmcTrack->IsPhysicalPrimary()){
\r
1778 Int_t gMotherIndex = AODmcTrack->GetMother();
\r
1779 if(gMotherIndex != -1) {
\r
1780 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
\r
1782 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
\r
1783 if(pdgCodeOfMother == 113 // rho0
\r
1784 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
\r
1785 // || pdgCodeOfMother == 221 // eta
\r
1786 // || pdgCodeOfMother == 331 // eta'
\r
1787 // || pdgCodeOfMother == 223 // omega
\r
1788 // || pdgCodeOfMother == 333 // phi
\r
1789 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
\r
1790 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
\r
1791 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
\r
1792 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
\r
1793 || pdgCodeOfMother == 111 // pi0 Dalitz
\r
1794 || pdgCodeOfMother == 22 // photon
\r
1796 kExcludeParticle = kTRUE;
\r
1801 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
\r
1802 if(kExcludeParticle) continue;
\r
1805 //Exclude electrons with PDG
\r
1806 if(fExcludeElectronsInMC) {
\r
1808 Int_t label = TMath::Abs(aodTrack->GetLabel());
\r
1809 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
\r
1812 if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
\r
1816 // fill QA histograms
\r
1817 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
\r
1818 fHistDCA->Fill(dcaZ,dcaXY);
\r
1819 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
\r
1820 fHistPt->Fill(vPt,gCentrality);
\r
1821 fHistEta->Fill(vEta,gCentrality);
\r
1822 fHistRapidity->Fill(vY,gCentrality);
\r
1823 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
1824 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
1825 fHistPhi->Fill(vPhi,gCentrality);
\r
1826 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
1827 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
1829 //=======================================correction
\r
1830 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
1831 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
1833 // add the track to the TObjArray
\r
1834 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
1837 //==============================================================================================================
\r
1839 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
\r
1841 AliESDtrack *trackTPC = NULL;
\r
1842 AliMCParticle *mcTrack = NULL;
\r
1844 AliMCEvent* mcEvent = NULL;
\r
1846 // for MC ESDs use also MC information
\r
1847 if(gAnalysisLevel == "MCESD"){
\r
1848 mcEvent = MCEvent();
\r
1850 AliError("mcEvent not available");
\r
1851 return tracksAccepted;
\r
1855 // Loop over tracks in event
\r
1856 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1857 AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));
\r
1859 AliError(Form("Could not receive track %d", iTracks));
\r
1863 // for MC ESDs use also MC information --> MC track not used further???
\r
1864 if(gAnalysisLevel == "MCESD"){
\r
1865 Int_t label = TMath::Abs(track->GetLabel());
\r
1866 if(label > mcEvent->GetNumberOfTracks()) continue;
\r
1867 if(label > mcEvent->GetNumberOfPrimaries()) continue;
\r
1869 mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
\r
1870 if(!mcTrack) continue;
\r
1873 // take only TPC only tracks
\r
1874 trackTPC = new AliESDtrack();
\r
1875 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
\r
1878 if(fESDtrackCuts)
\r
1879 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
\r
1881 // fill QA histograms
\r
1884 trackTPC->GetImpactParameters(b,bCov);
\r
1885 if (bCov[0]<=0 || bCov[2]<=0) {
\r
1886 AliDebug(1, "Estimated b resolution lower or equal zero!");
\r
1887 bCov[0]=0; bCov[2]=0;
\r
1890 Int_t nClustersTPC = -1;
\r
1891 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
\r
1892 //nClustersTPC = track->GetTPCclusters(0); // global track
\r
1893 Float_t chi2PerClusterTPC = -1;
\r
1894 if (nClustersTPC!=0) {
\r
1895 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
\r
1896 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
\r
1899 //===========================PID===============================//
\r
1901 Double_t prob[AliPID::kSPECIES]={0.};
\r
1902 Double_t probTPC[AliPID::kSPECIES]={0.};
\r
1903 Double_t probTOF[AliPID::kSPECIES]={0.};
\r
1904 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
\r
1906 Double_t nSigma = 0.;
\r
1907 UInt_t detUsedTPC = 0;
\r
1908 UInt_t detUsedTOF = 0;
\r
1909 UInt_t detUsedTPCTOF = 0;
\r
1911 //Decide what detector configuration we want to use
\r
1912 switch(fPidDetectorConfig) {
\r
1914 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
\r
1915 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
\r
1916 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
\r
1917 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1918 prob[iSpecies] = probTPC[iSpecies];
\r
1921 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
\r
1922 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
\r
1923 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
\r
1924 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1925 prob[iSpecies] = probTOF[iSpecies];
\r
1928 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
\r
1929 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
\r
1930 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1931 prob[iSpecies] = probTPCTOF[iSpecies];
\r
1935 }//end switch: define detector mask
\r
1937 //Filling the PID QA
\r
1938 Double_t tofTime = -999., length = 999., tof = -999.;
\r
1939 Double_t c = TMath::C()*1.E-9;// m/ns
\r
1940 Double_t beta = -999.;
\r
1941 Double_t nSigmaTOFForParticleOfInterest = -999.;
\r
1942 if ( (track->IsOn(AliESDtrack::kTOFin)) &&
\r
1943 (track->IsOn(AliESDtrack::kTIME)) ) {
\r
1944 tofTime = track->GetTOFsignal();//in ps
\r
1945 length = track->GetIntegratedLength();
\r
1946 tof = tofTime*1E-3; // ns
\r
1949 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
\r
1953 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
\r
1957 length = length*0.01; // in meters
\r
1959 beta = length/tof;
\r
1961 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
\r
1962 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
\r
1963 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
\r
1964 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
\r
1968 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
\r
1969 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
\r
1970 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
\r
1971 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
\r
1972 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
\r
1973 //end of QA-before pid
\r
1975 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
\r
1976 //Make the decision based on the n-sigma
\r
1977 if(fUsePIDnSigma) {
\r
1978 if(nSigma > fPIDNSigma) continue;}
\r
1980 //Make the decision based on the bayesian
\r
1981 else if(fUsePIDPropabilities) {
\r
1982 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
\r
1983 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
\r
1986 //Fill QA after the PID
\r
1987 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
\r
1988 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
\r
1989 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
\r
1991 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
\r
1992 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
\r
1993 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
\r
1994 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
\r
1997 //===========================PID===============================//
\r
1998 vCharge = trackTPC->Charge();
\r
1999 vY = trackTPC->Y();
\r
2000 vEta = trackTPC->Eta();
\r
2001 vPhi = trackTPC->Phi();// * TMath::RadToDeg();
\r
2002 vPt = trackTPC->Pt();
\r
2003 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
\r
2004 fHistDCA->Fill(b[1],b[0]);
\r
2005 fHistChi2->Fill(chi2PerClusterTPC,gCentrality);
\r
2006 fHistPt->Fill(vPt,gCentrality);
\r
2007 fHistEta->Fill(vEta,gCentrality);
\r
2008 fHistPhi->Fill(vPhi,gCentrality);
\r
2009 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
2010 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
2011 fHistRapidity->Fill(vY,gCentrality);
\r
2012 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
2013 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
2015 //=======================================correction
\r
2016 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
2017 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
2019 // add the track to the TObjArray
\r
2020 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
2026 else if(gAnalysisLevel == "MC"){
\r
2028 AliError("mcEvent not available");
\r
2032 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
2034 // Loop over tracks in event
\r
2035 for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
\r
2036 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
\r
2038 AliError(Form("Could not receive particle %d", iTracks));
\r
2042 //exclude non stable particles
\r
2043 if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
\r
2045 vCharge = track->Charge();
\r
2046 vEta = track->Eta();
\r
2047 vPt = track->Pt();
\r
2050 if( vPt < fPtMin || vPt > fPtMax)
\r
2053 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
2055 else if (fUsePID){
\r
2056 if( vY < fEtaMin || vY > fEtaMax) continue;
\r
2059 // Remove neutral tracks
\r
2060 if( vCharge == 0 ) continue;
\r
2062 //analyze one set of particles
\r
2063 if(fUseMCPdgCode) {
\r
2064 TParticle *particle = track->Particle();
\r
2065 if(!particle) continue;
\r
2067 Int_t gPdgCode = particle->GetPdgCode();
\r
2068 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
\r
2072 //Use the acceptance parameterization
\r
2073 if(fAcceptanceParameterization) {
\r
2074 Double_t gRandomNumber = gRandom->Rndm();
\r
2075 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
\r
2079 //Exclude resonances
\r
2080 if(fExcludeResonancesInMC) {
\r
2081 TParticle *particle = track->Particle();
\r
2082 if(!particle) continue;
\r
2084 Bool_t kExcludeParticle = kFALSE;
\r
2085 Int_t gMotherIndex = particle->GetFirstMother();
\r
2086 if(gMotherIndex != -1) {
\r
2087 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
\r
2089 TParticle *motherParticle = motherTrack->Particle();
\r
2090 if(motherParticle) {
\r
2091 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
\r
2092 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
\r
2093 if(pdgCodeOfMother == 113 // rho0
\r
2094 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
\r
2095 // || pdgCodeOfMother == 221 // eta
\r
2096 // || pdgCodeOfMother == 331 // eta'
\r
2097 // || pdgCodeOfMother == 223 // omega
\r
2098 // || pdgCodeOfMother == 333 // phi
\r
2099 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
\r
2100 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
\r
2101 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
\r
2102 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
\r
2103 || pdgCodeOfMother == 111 // pi0 Dalitz
\r
2105 kExcludeParticle = kTRUE;
\r
2111 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
\r
2112 if(kExcludeParticle) continue;
\r
2115 //Exclude electrons with PDG
\r
2116 if(fExcludeElectronsInMC) {
\r
2118 TParticle *particle = track->Particle();
\r
2121 if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
\r
2125 vPhi = track->Phi();
\r
2126 //Printf("phi (before): %lf",vPhi);
\r
2128 fHistPt->Fill(vPt,gCentrality);
\r
2129 fHistEta->Fill(vEta,gCentrality);
\r
2130 fHistPhi->Fill(vPhi,gCentrality);
\r
2131 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
2132 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
2133 //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
\r
2134 fHistRapidity->Fill(vY,gCentrality);
\r
2135 //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
\r
2136 //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
\r
2137 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
2138 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
2140 //Flow after burner
\r
2141 if(fUseFlowAfterBurner) {
\r
2142 Double_t precisionPhi = 0.001;
\r
2143 Int_t maxNumberOfIterations = 100;
\r
2145 Double_t phi0 = vPhi;
\r
2146 Double_t gV2 = fDifferentialV2->Eval(vPt);
\r
2148 for (Int_t j = 0; j < maxNumberOfIterations; j++) {
\r
2149 Double_t phiprev = vPhi;
\r
2150 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
\r
2151 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
\r
2153 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
\r
2155 //Printf("phi (after): %lf\n",vPhi);
\r
2156 Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
\r
2157 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
\r
2158 fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
\r
2160 Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
\r
2161 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
\r
2162 fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
\r
2166 //vPhi *= TMath::RadToDeg();
\r
2168 //=======================================correction
\r
2169 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
2170 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
2172 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
2174 }//MC event object
\r
2177 return tracksAccepted;
\r
2180 //________________________________________________________________________
\r
2181 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
\r
2182 // Clones TObjArray and returns it with tracks after shuffling the charges
\r
2184 TObjArray* tracksShuffled = new TObjArray;
\r
2185 tracksShuffled->SetOwner(kTRUE);
\r
2187 vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks
\r
2189 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
\r
2191 AliVParticle* track = (AliVParticle*) tracks->At(i);
\r
2192 chargeVector->push_back(track->Charge());
\r
2195 random_shuffle(chargeVector->begin(), chargeVector->end());
\r
2197 for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
\r
2198 AliVParticle* track = (AliVParticle*) tracks->At(i);
\r
2199 //==============================correction
\r
2200 Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
\r
2201 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
2202 tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
\r
2205 delete chargeVector;
\r
2207 return tracksShuffled;
\r
2210 //________________________________________________________________________
\r
2211 void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,
\r
2212 const char* lhcPeriod) {
\r
2213 //Function to setup the VZERO gain equalization
\r
2214 //============Get the equilization map============//
\r
2215 TFile *calibrationFile = TFile::Open(filename);
\r
2216 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
\r
2217 Printf("No calibration file found!!!");
\r
2221 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
\r
2223 Printf("Calibration TList not found!!!");
\r
2227 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
\r
2228 if(!fHistVZEROAGainEqualizationMap) {
\r
2229 Printf("VZERO-A calibration object not found!!!");
\r
2232 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
\r
2233 if(!fHistVZEROCGainEqualizationMap) {
\r
2234 Printf("VZERO-C calibration object not found!!!");
\r
2238 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
\r
2239 if(!fHistVZEROChannelGainEqualizationMap) {
\r
2240 Printf("VZERO channel calibration object not found!!!");
\r
2245 //________________________________________________________________________
\r
2246 Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run,
\r
2249 if(!fHistVZEROAGainEqualizationMap) return 1.0;
\r
2251 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
\r
2252 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
\r
2253 if(gRunNumber == run)
\r
2254 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
\r
2260 //________________________________________________________________________
\r
2261 Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run,
\r
2262 const char* side) {
\r
2264 if(!fHistVZEROAGainEqualizationMap) return 1.0;
\r
2266 TString gVZEROSide = side;
\r
2267 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
\r
2268 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
\r
2269 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
\r
2270 if(gRunNumber == run) {
\r
2271 if(gVZEROSide == "A")
\r
2272 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
\r
2273 else if(gVZEROSide == "C")
\r
2274 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
\r
2281 //________________________________________________________________________
\r
2282 void AliAnalysisTaskBFPsi::FinishTaskOutput(){
\r
2283 //Printf("END BF");
\r
2286 AliError("fBalance not available");
\r
2289 if(fRunShuffling) {
\r
2290 if (!fShuffledBalance) {
\r
2291 AliError("fShuffledBalance not available");
\r
2298 //________________________________________________________________________
\r
2299 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
\r
2300 // Draw result to the screen
\r
2301 // Called once at the end of the query
\r
2303 // not implemented ...
\r