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 fHistdEdxVsPTPCbeforePIDelectron(NULL), //+++++++
\r
118 fHistNSigmaTPCvsPtbeforePIDelectron(NULL), //+++++++
\r
119 fHistdEdxVsPTPCafterPIDelectron(NULL), //+++++++
\r
120 fHistNSigmaTPCvsPtafterPIDelectron(NULL), //+++++++
\r
121 fCentralityArrayBinsForCorrections(kCENTRALITY),
\r
122 fCentralityWeights(0x0),
\r
125 fParticleOfInterest(kPion),
\r
126 fPidDetectorConfig(kTPCTOF),
\r
128 fUsePIDnSigma(kTRUE),
\r
129 fUsePIDPropabilities(kFALSE),
\r
131 fMinAcceptedPIDProbability(0.8),
\r
132 fElectronRejection(kFALSE),
\r
133 fElectronOnlyRejection(kFALSE),
\r
134 fElectronRejectionNSigma(-1.),
\r
135 fElectronRejectionMinPt(0.),
\r
136 fElectronRejectionMaxPt(1000.),
\r
138 fCentralityEstimator("V0M"),
\r
139 fUseCentrality(kFALSE),
\r
140 fCentralityPercentileMin(0.),
\r
141 fCentralityPercentileMax(5.),
\r
142 fImpactParameterMin(0.),
\r
143 fImpactParameterMax(20.),
\r
144 fMultiplicityEstimator("V0A"),
\r
145 fUseMultiplicity(kFALSE),
\r
146 fNumberOfAcceptedTracksMin(0),
\r
147 fNumberOfAcceptedTracksMax(10000),
\r
148 fHistNumberOfAcceptedTracks(0),
\r
149 fHistMultiplicity(0),
\r
150 fUseOfflineTrigger(kFALSE),
\r
151 fCheckFirstEventInChunk(kFALSE),
\r
152 fCheckPileUp(kFALSE),
\r
153 fCheckPrimaryFlagAOD(kFALSE),
\r
154 fUseMCforKinematics(kFALSE),
\r
158 fnAODtrackCutBit(128),
\r
161 fPtMinForCorrections(0.3),//=================================correction
\r
162 fPtMaxForCorrections(1.5),//=================================correction
\r
163 fPtBinForCorrections(36), //=================================correction
\r
166 fEtaMinForCorrections(-0.8),//=================================correction
\r
167 fEtaMaxForCorrections(0.8),//=================================correction
\r
168 fEtaBinForCorrections(16), //=================================correction
\r
171 fPhiMinForCorrections(0.),//=================================correction
\r
172 fPhiMaxForCorrections(360.),//=================================correction
\r
173 fPhiBinForCorrections(100), //=================================correction
\r
177 fNClustersTPCCut(-1),
\r
178 fAcceptanceParameterization(0),
\r
179 fDifferentialV2(0),
\r
180 fUseFlowAfterBurner(kFALSE),
\r
181 fExcludeResonancesInMC(kFALSE),
\r
182 fExcludeElectronsInMC(kFALSE),
\r
183 fUseMCPdgCode(kFALSE),
\r
184 fPDGCodeToBeAnalyzed(-1),
\r
185 fEventClass("EventPlane"),
\r
186 fCustomBinning(""),
\r
187 fHistVZEROAGainEqualizationMap(0),
\r
188 fHistVZEROCGainEqualizationMap(0),
\r
189 fHistVZEROChannelGainEqualizationMap(0) {
\r
191 // Define input and output slots here
\r
192 // Input slot #0 works with a TChain
\r
194 //======================================================correction
\r
195 for (Int_t i=0; i<kCENTRALITY; i++){
\r
196 fHistCorrectionPlus[i] = NULL;
\r
197 fHistCorrectionMinus[i] = NULL;
\r
198 fCentralityArrayForCorrections[i] = -1.;
\r
200 //=====================================================correction
\r
202 DefineInput(0, TChain::Class());
\r
203 // Output slot #0 writes into a TH1 container
\r
204 DefineOutput(1, TList::Class());
\r
205 DefineOutput(2, TList::Class());
\r
206 DefineOutput(3, TList::Class());
\r
207 DefineOutput(4, TList::Class());
\r
208 DefineOutput(5, TList::Class());
\r
211 //________________________________________________________________________
\r
212 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
\r
214 // delete fBalance;
\r
215 // delete fShuffledBalance;
\r
217 // delete fListBF;
\r
218 // delete fListBFS;
\r
220 // delete fHistEventStats;
\r
221 // delete fHistTrackStats;
\r
222 // delete fHistVx;
\r
223 // delete fHistVy;
\r
224 // delete fHistVz;
\r
226 // delete fHistClus;
\r
227 // delete fHistDCA;
\r
228 // delete fHistChi2;
\r
230 // delete fHistEta;
\r
231 // delete fHistPhi;
\r
232 // delete fHistEtaPhiPos;
\r
233 // delete fHistEtaPhiNeg;
\r
234 // delete fHistV0M;
\r
237 //________________________________________________________________________
\r
238 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
\r
239 // Create histograms
\r
242 // global switch disabling the reference
\r
243 // (to avoid "Replacing existing TH1" if several wagons are created in train)
\r
244 Bool_t oldStatus = TH1::AddDirectoryStatus();
\r
245 TH1::AddDirectory(kFALSE);
\r
248 fBalance = new AliBalancePsi();
\r
249 fBalance->SetAnalysisLevel("ESD");
\r
250 fBalance->SetEventClass(fEventClass);
\r
251 //fBalance->SetNumberOfBins(-1,16);
\r
252 //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
\r
254 if(fRunShuffling) {
\r
255 if(!fShuffledBalance) {
\r
256 fShuffledBalance = new AliBalancePsi();
\r
257 fShuffledBalance->SetAnalysisLevel("ESD");
\r
258 fShuffledBalance->SetEventClass(fEventClass);
\r
259 //fShuffledBalance->SetNumberOfBins(-1,16);
\r
260 //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
\r
264 if(!fMixedBalance) {
\r
265 fMixedBalance = new AliBalancePsi();
\r
266 fMixedBalance->SetAnalysisLevel("ESD");
\r
267 fMixedBalance->SetEventClass(fEventClass);
\r
272 fList = new TList();
\r
273 fList->SetName("listQA");
\r
276 //Balance Function list
\r
277 fListBF = new TList();
\r
278 fListBF->SetName("listBF");
\r
279 fListBF->SetOwner();
\r
281 if(fRunShuffling) {
\r
282 fListBFS = new TList();
\r
283 fListBFS->SetName("listBFShuffled");
\r
284 fListBFS->SetOwner();
\r
288 fListBFM = new TList();
\r
289 fListBFM->SetName("listTriggeredBFMixed");
\r
290 fListBFM->SetOwner();
\r
294 if(fUsePID || fElectronRejection) {
\r
295 fHistListPIDQA = new TList();
\r
296 fHistListPIDQA->SetName("listQAPID");
\r
297 fHistListPIDQA->SetOwner();
\r
301 TString gCutName[7] = {"Total","Offline trigger",
\r
302 "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};
\r
303 fHistEventStats = new TH2F("fHistEventStats",
\r
304 "Event statistics;;Centrality percentile;N_{events}",
\r
305 7,0.5,7.5,220,-5,105);
\r
306 for(Int_t i = 1; i <= 7; i++)
\r
307 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
\r
308 fList->Add(fHistEventStats);
\r
310 TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
\r
311 fHistCentStats = new TH2F("fHistCentStats",
\r
312 "Centrality statistics;;Cent percentile",
\r
313 13,-0.5,12.5,220,-5,105);
\r
314 for(Int_t i = 1; i <= 13; i++){
\r
315 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
\r
316 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
\r
318 fList->Add(fHistCentStats);
\r
320 fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);
\r
321 fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());
\r
322 fList->Add(fHistCentStatsUsed);
\r
324 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
\r
325 fList->Add(fHistTriggerStats);
\r
327 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
\r
328 fList->Add(fHistTrackStats);
\r
330 fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
\r
331 fList->Add(fHistNumberOfAcceptedTracks);
\r
333 fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);
\r
334 fList->Add(fHistMultiplicity);
\r
336 // Vertex distributions
\r
337 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
\r
338 fList->Add(fHistVx);
\r
339 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
\r
340 fList->Add(fHistVy);
\r
341 fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
\r
342 fList->Add(fHistVz);
\r
344 //TPC vs VZERO multiplicity
\r
345 fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
\r
346 if(fMultiplicityEstimator == "V0A")
\r
347 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");
\r
348 else if(fMultiplicityEstimator == "V0C")
\r
349 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");
\r
351 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");
\r
352 fList->Add(fHistTPCvsVZEROMultiplicity);
\r
354 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
355 fList->Add(fHistVZEROSignal);
\r
358 fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
\r
359 fList->Add(fHistEventPlane);
\r
362 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
\r
363 fList->Add(fHistClus);
\r
364 fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);
\r
365 fList->Add(fHistChi2);
\r
366 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
\r
367 fList->Add(fHistDCA);
\r
368 fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);
\r
369 fList->Add(fHistPt);
\r
370 fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);
\r
371 fList->Add(fHistEta);
\r
372 fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);
\r
373 fList->Add(fHistRapidity);
\r
374 fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);
\r
375 fList->Add(fHistPhi);
\r
376 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
377 fList->Add(fHistEtaPhiPos);
\r
378 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
379 fList->Add(fHistEtaPhiNeg);
\r
380 fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
\r
381 fList->Add(fHistPhiBefore);
\r
382 fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
\r
383 fList->Add(fHistPhiAfter);
\r
384 fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
\r
385 fList->Add(fHistPhiPos);
\r
386 fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);
\r
387 fList->Add(fHistPhiNeg);
\r
388 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
\r
389 fList->Add(fHistV0M);
\r
390 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
\r
391 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
\r
392 for(Int_t i = 1; i <= 6; i++)
\r
393 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
\r
394 fList->Add(fHistRefTracks);
\r
396 // Balance function histograms
\r
397 // Initialize histograms if not done yet (including the custom binning)
\r
398 if(!fBalance->GetHistNp()){
\r
399 AliInfo("Histograms not yet initialized! --> Will be done now");
\r
400 fBalance->SetCustomBinning(fCustomBinning);
\r
401 fBalance->InitHistograms();
\r
404 if(fRunShuffling) {
\r
405 if(!fShuffledBalance->GetHistNp()) {
\r
406 AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");
\r
407 fShuffledBalance->SetCustomBinning(fCustomBinning);
\r
408 fShuffledBalance->InitHistograms();
\r
413 if(!fMixedBalance->GetHistNp()) {
\r
414 AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");
\r
415 fMixedBalance->SetCustomBinning(fCustomBinning);
\r
416 fMixedBalance->InitHistograms();
\r
420 // QA histograms for different cuts
\r
421 fList->Add(fBalance->GetQAHistHBTbefore());
\r
422 fList->Add(fBalance->GetQAHistHBTafter());
\r
423 fList->Add(fBalance->GetQAHistConversionbefore());
\r
424 fList->Add(fBalance->GetQAHistConversionafter());
\r
425 fList->Add(fBalance->GetQAHistPsiMinusPhi());
\r
426 fList->Add(fBalance->GetQAHistResonancesBefore());
\r
427 fList->Add(fBalance->GetQAHistResonancesRho());
\r
428 fList->Add(fBalance->GetQAHistResonancesK0());
\r
429 fList->Add(fBalance->GetQAHistResonancesLambda());
\r
430 fList->Add(fBalance->GetQAHistQbefore());
\r
431 fList->Add(fBalance->GetQAHistQafter());
\r
433 //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
\r
434 fListBF->Add(fBalance->GetHistNp());
\r
435 fListBF->Add(fBalance->GetHistNn());
\r
436 fListBF->Add(fBalance->GetHistNpn());
\r
437 fListBF->Add(fBalance->GetHistNnn());
\r
438 fListBF->Add(fBalance->GetHistNpp());
\r
439 fListBF->Add(fBalance->GetHistNnp());
\r
441 if(fRunShuffling) {
\r
442 fListBFS->Add(fShuffledBalance->GetHistNp());
\r
443 fListBFS->Add(fShuffledBalance->GetHistNn());
\r
444 fListBFS->Add(fShuffledBalance->GetHistNpn());
\r
445 fListBFS->Add(fShuffledBalance->GetHistNnn());
\r
446 fListBFS->Add(fShuffledBalance->GetHistNpp());
\r
447 fListBFS->Add(fShuffledBalance->GetHistNnp());
\r
451 fListBFM->Add(fMixedBalance->GetHistNp());
\r
452 fListBFM->Add(fMixedBalance->GetHistNn());
\r
453 fListBFM->Add(fMixedBalance->GetHistNpn());
\r
454 fListBFM->Add(fMixedBalance->GetHistNnn());
\r
455 fListBFM->Add(fMixedBalance->GetHistNpp());
\r
456 fListBFM->Add(fMixedBalance->GetHistNnp());
\r
463 Int_t trackDepth = fMixingTracks;
\r
464 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
\r
467 Double_t* centbins = NULL;
\r
468 Int_t nCentralityBins;
\r
469 if(fBalance->IsUseVertexBinning()){
\r
470 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);
\r
473 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);
\r
476 // multiplicity bins
\r
477 Double_t* multbins = NULL;
\r
478 Int_t nMultiplicityBins;
\r
479 multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);
\r
482 Double_t* vtxbins = NULL;
\r
484 if(fBalance->IsUseVertexBinning()){
\r
485 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);
\r
488 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);
\r
491 // Event plane angle (Psi) bins
\r
492 Double_t* psibins = NULL;
\r
494 psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);
\r
497 // run the event mixing also in bins of event plane (statistics!)
\r
498 if(fRunMixingEventPlane){
\r
499 if(fEventClass=="Multiplicity"){
\r
500 if(multbins && vtxbins && psibins){
\r
501 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
\r
505 if(centbins && vtxbins && psibins){
\r
506 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
\r
511 if(fEventClass=="Multiplicity"){
\r
512 if(multbins && vtxbins){
\r
513 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
\r
517 if(centbins && vtxbins){
\r
518 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
\r
523 // check pool manager
\r
525 AliError("Event Mixing required, but Pool Manager not initialized...");
\r
530 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
\r
532 //====================PID========================//
\r
534 fPIDCombined = new AliPIDCombined();
\r
535 fPIDCombined->SetDefaultTPCPriors();
\r
537 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
\r
538 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
\r
540 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
\r
541 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID);
\r
543 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
\r
544 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID);
\r
546 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
\r
547 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);
\r
549 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
\r
550 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);
\r
552 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
\r
553 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);
\r
555 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
\r
556 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID);
\r
558 fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2);
\r
559 fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++
\r
561 fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500);
\r
562 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++
\r
564 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
\r
565 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
\r
567 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
\r
568 fHistListPIDQA->Add(fHistBetavsPTOFafterPID);
\r
570 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
\r
571 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);
\r
573 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
\r
574 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID);
\r
576 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
\r
577 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);
\r
579 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
\r
580 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
\r
582 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
\r
583 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);
\r
585 fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2);
\r
586 fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++
\r
588 fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500);
\r
589 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++
\r
592 // for electron rejection only TPC nsigma histograms
\r
593 if(fElectronRejection) {
\r
595 fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000);
\r
596 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);
\r
598 fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500);
\r
599 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);
\r
601 fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000);
\r
602 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);
\r
604 fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500);
\r
605 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron);
\r
607 //====================PID========================//
\r
609 // Post output data.
\r
610 PostData(1, fList);
\r
611 PostData(2, fListBF);
\r
612 if(fRunShuffling) PostData(3, fListBFS);
\r
613 if(fRunMixing) PostData(4, fListBFM);
\r
614 if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID
\r
616 AliInfo("Finished setting up the Output");
\r
618 TH1::AddDirectory(oldStatus);
\r
622 //________________________________________________________________________
\r
623 void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename,
\r
624 Int_t nCentralityBins,
\r
625 Double_t *centralityArrayForCorrections) {
\r
626 //Open files that will be used for correction
\r
627 fCentralityArrayBinsForCorrections = nCentralityBins;
\r
628 for (Int_t i=0; i<nCentralityBins; i++)
\r
629 fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];
\r
631 // No file specified -> run without corrections
\r
632 if(!filename.Contains(".root")) {
\r
633 AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));
\r
637 //Open the input file
\r
638 TFile *f = TFile::Open(filename);
\r
640 AliInfo(Form("File %s not found --> run without corrections",filename.Data()));
\r
644 //TString listEffName = "";
\r
645 for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {
\r
646 //Printf("iCent %d:",iCent);
\r
647 TString histoName = "fHistCorrectionPlus";
\r
648 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
\r
649 fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));
\r
650 if(!fHistCorrectionPlus[iCent]) {
\r
651 AliError(Form("fHist %s not found!!!",histoName.Data()));
\r
655 histoName = "fHistCorrectionMinus";
\r
656 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
\r
657 fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data()));
\r
658 if(!fHistCorrectionMinus[iCent]) {
\r
659 AliError(Form("fHist %s not found!!!",histoName.Data()));
\r
662 }//loop over centralities: ONLY the PbPb case is covered
\r
664 if(fHistCorrectionPlus[0]){
\r
665 fEtaMinForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmin();
\r
666 fEtaMaxForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmax();
\r
667 fEtaBinForCorrections = fHistCorrectionPlus[0]->GetNbinsX();
\r
669 fPtMinForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmin();
\r
670 fPtMaxForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmax();
\r
671 fPtBinForCorrections = fHistCorrectionPlus[0]->GetNbinsY();
\r
673 fPhiMinForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmin();
\r
674 fPhiMaxForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmax();
\r
675 fPhiBinForCorrections = fHistCorrectionPlus[0]->GetNbinsZ();
\r
679 //________________________________________________________________________
\r
680 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
\r
682 // Called for each event
\r
684 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
685 Int_t gNumberOfAcceptedTracks = 0;
\r
686 Double_t lMultiplicityVar = -999.; //-1
\r
687 Double_t gReactionPlane = -1.;
\r
688 Float_t bSign = 0.;
\r
690 // get the event (for generator level: MCEvent())
\r
691 AliVEvent* eventMain = NULL;
\r
692 if(gAnalysisLevel == "MC") {
\r
693 eventMain = dynamic_cast<AliVEvent*>(MCEvent());
\r
696 eventMain = dynamic_cast<AliVEvent*>(InputEvent());
\r
697 // for HBT like cuts need magnetic field sign
\r
698 bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
\r
701 AliError("eventMain not available");
\r
705 // PID Response task active?
\r
706 if(fUsePID || fElectronRejection) {
\r
707 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
\r
708 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
\r
711 // check event cuts and fill event histograms
\r
712 if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){
\r
715 // get the reaction plane
\r
716 if(fEventClass != "Multiplicity") {
\r
717 gReactionPlane = GetEventPlane(eventMain);
\r
718 fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
\r
719 if(gReactionPlane < 0){
\r
724 // get the accepted tracks in main event
\r
725 TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
\r
726 gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
\r
728 //multiplicity cut (used in pp)
\r
729 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
\r
731 // store charges of all accepted tracks,shuffle and reassign(two extra loops)
\r
732 TObjArray* tracksShuffled = NULL;
\r
734 tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
\r
740 // 1. First get an event pool corresponding in mult (cent) and
\r
741 // zvertex to the current event. Once initialized, the pool
\r
742 // should contain nMix (reduced) events. This routine does not
\r
743 // pre-scan the chain. The first several events of every chain
\r
744 // will be skipped until the needed pools are filled to the
\r
745 // specified depth. If the pool categories are not too rare, this
\r
746 // should not be a problem. If they are rare, you could lose`
\r
749 // 2. Collect the whole pool's content of tracks into one TObjArray
\r
750 // (bgTracks), which is effectively a single background super-event.
\r
752 // 3. The reduced and bgTracks arrays must both be passed into
\r
753 // FillCorrelations(). Also nMix should be passed in, so a weight
\r
754 // of 1./nMix can be applied.
\r
756 AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
\r
759 AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
\r
763 //pool->SetDebug(1);
\r
765 if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){
\r
768 Int_t nMix = pool->GetCurrentNEvents();
\r
769 //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
\r
771 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
\r
772 //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
\r
773 //if (pool->IsReady())
\r
774 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
\r
776 // Fill mixed-event histos here
\r
777 for (Int_t jMix=0; jMix<nMix; jMix++)
\r
779 TObjArray* tracksMixed = pool->GetEvent(jMix);
\r
780 fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
\r
784 // Update the Event pool
\r
785 pool->UpdatePool(tracksMain);
\r
786 //pool->PrintInfo();
\r
788 }//pool NULL check
\r
791 // calculate balance function
\r
792 fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
\r
794 // calculate shuffled balance function
\r
795 if(fRunShuffling && tracksShuffled != NULL) {
\r
796 fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
\r
800 //________________________________________________________________________
\r
801 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
\r
802 // Checks the Event cuts
\r
803 // Fills Event statistics histograms
\r
805 Bool_t isSelectedMain = kTRUE;
\r
806 Float_t gRefMultiplicity = -1.;
\r
807 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
809 AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
\r
811 fHistEventStats->Fill(1,gRefMultiplicity); //all events
\r
813 // check first event in chunk (is not needed for new reconstructions)
\r
814 if(fCheckFirstEventInChunk){
\r
815 AliAnalysisUtils ut;
\r
816 if(ut.IsFirstEventInChunk(event))
\r
818 fHistEventStats->Fill(6,gRefMultiplicity);
\r
820 // check for pile-up event
\r
822 AliAnalysisUtils ut;
\r
823 ut.SetUseMVPlpSelection(kTRUE);
\r
824 ut.SetUseOutOfBunchPileUp(kTRUE);
\r
825 if(ut.IsPileUpEvent(event))
\r
827 fHistEventStats->Fill(7,gRefMultiplicity);
\r
830 // Event trigger bits
\r
831 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
\r
832 if(fUseOfflineTrigger)
\r
833 isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
\r
835 if(isSelectedMain) {
\r
836 fHistEventStats->Fill(2,gRefMultiplicity); //triggered events
\r
839 if(gAnalysisLevel == "MC") {
\r
841 AliError("mcEvent not available");
\r
846 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
\r
848 TArrayF gVertexArray;
\r
849 header->PrimaryVertex(gVertexArray);
\r
850 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
\r
851 //gVertexArray.At(0),
\r
852 //gVertexArray.At(1),
\r
853 //gVertexArray.At(2));
\r
854 fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex
\r
855 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
\r
856 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
\r
857 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
\r
858 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
\r
860 // get the reference multiplicty or centrality
\r
861 gRefMultiplicity = GetRefMultiOrCentrality(event);
\r
863 fHistVx->Fill(gVertexArray.At(0));
\r
864 fHistVy->Fill(gVertexArray.At(1));
\r
865 fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);
\r
867 // take only events inside centrality class
\r
868 if(fUseCentrality) {
\r
869 if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){
\r
870 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
\r
871 return gRefMultiplicity;
\r
872 }//centrality class
\r
874 // take events only within the same multiplicity class
\r
875 else if(fUseMultiplicity) {
\r
876 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
\r
877 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
\r
878 return gRefMultiplicity;
\r
880 }//multiplicity range
\r
888 // Event Vertex AOD, ESD, ESDMC
\r
890 const AliVVertex *vertex = event->GetPrimaryVertex();
\r
893 Double32_t fCov[6];
\r
894 vertex->GetCovarianceMatrix(fCov);
\r
895 if(vertex->GetNContributors() > 0) {
\r
897 fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex
\r
898 if(TMath::Abs(vertex->GetX()) < fVxMax) {
\r
899 if(TMath::Abs(vertex->GetY()) < fVyMax) {
\r
900 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
\r
901 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
\r
903 // get the reference multiplicty or centrality
\r
904 gRefMultiplicity = GetRefMultiOrCentrality(event);
\r
906 fHistVx->Fill(vertex->GetX());
\r
907 fHistVy->Fill(vertex->GetY());
\r
908 fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);
\r
910 // take only events inside centrality class
\r
911 if(fUseCentrality) {
\r
912 if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){
\r
914 // centrality weighting (optional for 2011 if central and semicentral triggers are used)
\r
915 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity)){
\r
916 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
\r
920 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
\r
921 return gRefMultiplicity;
\r
922 }//centrality class
\r
924 // take events only within the same multiplicity class
\r
925 else if(fUseMultiplicity) {
\r
927 //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,
\r
928 //fNumberOfAcceptedTracksMax,gRefMultiplicity);
\r
930 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
\r
931 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
\r
932 return gRefMultiplicity;
\r
934 }//multiplicity range
\r
938 }//proper vertex resolution
\r
939 }//proper number of contributors
\r
940 }//vertex object valid
\r
941 }//triggered event
\r
944 // in all other cases return -1 (event not accepted)
\r
949 //________________________________________________________________________
\r
950 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
\r
951 // Checks the Event cuts
\r
952 // Fills Event statistics histograms
\r
954 Float_t gCentrality = -1.;
\r
955 Double_t gMultiplicity = -1.;
\r
956 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
959 // calculate centrality always (not only in centrality mode)
\r
960 if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++
\r
961 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
\r
963 gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
\r
965 // QA for centrality estimators
\r
966 fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
\r
967 fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));
\r
968 fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));
\r
969 fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
\r
970 fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
\r
971 fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL"));
\r
972 fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
\r
973 fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
\r
974 fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));
\r
975 fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));
\r
976 fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
\r
977 fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
\r
978 fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
\r
980 // Centrality estimator USED ++++++++++++++++++++++++++++++
\r
981 fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));
\r
983 // centrality QA (V0M)
\r
984 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
\r
986 // centrality QA (reference tracks)
\r
987 fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
\r
988 fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
\r
989 fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
\r
990 fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
\r
991 fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
\r
992 fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
\r
993 fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
\r
994 fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
\r
995 fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
\r
1000 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
\r
1001 AliCentrality *centrality = event->GetCentrality();
\r
1002 gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
\r
1004 // QA for centrality estimators
\r
1005 fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));
\r
1006 fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));
\r
1007 fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));
\r
1008 fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));
\r
1009 fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));
\r
1010 fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));
\r
1011 fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));
\r
1012 fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));
\r
1013 fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));
\r
1014 fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));
\r
1015 fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));
\r
1016 fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));
\r
1017 fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));
\r
1019 // Centrality estimator USED ++++++++++++++++++++++++++++++
\r
1020 fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
\r
1022 // centrality QA (V0M)
\r
1023 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
\r
1026 else if(gAnalysisLevel == "MC"){
\r
1027 Double_t gImpactParameter = 0.;
\r
1028 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
1030 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
\r
1032 gImpactParameter = headerH->ImpactParameter();
\r
1033 gCentrality = gImpactParameter;
\r
1039 gCentrality = -1.;
\r
1042 // calculate reference multiplicity always (not only in multiplicity mode)
\r
1043 if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){
\r
1044 AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
\r
1046 gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
\r
1047 fHistMultiplicity->Fill(gMultiplicity);
\r
1048 }//AliESDevent cast
\r
1051 else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){
\r
1052 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
\r
1053 if ((fMultiplicityEstimator == "V0M")||
\r
1054 (fMultiplicityEstimator == "V0A")||
\r
1055 (fMultiplicityEstimator == "V0C") ||
\r
1056 (fMultiplicityEstimator == "TPC")) {
\r
1057 gMultiplicity = GetReferenceMultiplicityFromAOD(event);
\r
1058 if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);
\r
1062 gMultiplicity = header->GetRefMultiplicity();
\r
1063 if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);
\r
1065 fHistMultiplicity->Fill(gMultiplicity);
\r
1067 else if(gAnalysisLevel == "MC") {
\r
1068 AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
1069 //Calculating the multiplicity as the number of charged primaries
\r
1070 //within \pm 0.8 in eta and pT > 0.1 GeV/c
\r
1071 for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
\r
1072 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
\r
1074 AliError(Form("Could not receive particle %d", iParticle));
\r
1078 //exclude non stable particles
\r
1079 if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
\r
1081 //++++++++++++++++
\r
1082 if (fMultiplicityEstimator == "V0M") {
\r
1083 if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7))
\r
1085 else if (fMultiplicityEstimator == "V0A") {
\r
1086 if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}
\r
1087 else if (fMultiplicityEstimator == "V0C") {
\r
1088 if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}
\r
1089 else if (fMultiplicityEstimator == "TPC") {
\r
1090 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
\r
1091 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
\r
1094 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
\r
1095 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
\r
1097 //++++++++++++++++
\r
1099 if(track->Charge() == 0) continue;
\r
1101 gMultiplicity += 1;
\r
1102 }//loop over primaries
\r
1103 fHistMultiplicity->Fill(gMultiplicity);
\r
1106 gMultiplicity = -1;
\r
1110 // decide what should be returned only here
\r
1111 Double_t lReturnVal = -100;
\r
1112 if(fEventClass=="Multiplicity"){
\r
1113 lReturnVal = gMultiplicity;
\r
1114 }else if(fEventClass=="Centrality"){
\r
1115 lReturnVal = gCentrality;
\r
1117 return lReturnVal;
\r
1120 //________________________________________________________________________
\r
1121 Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){
\r
1122 //Function that returns the reference multiplicity from AODs (data or reco MC)
\r
1123 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
\r
1124 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
\r
1125 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
\r
1127 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
\r
1129 Printf("ERROR: AOD header not available");
\r
1132 Int_t gRunNumber = header->GetRunNumber();
\r
1134 // Loop over tracks in event
\r
1135 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1136 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
\r
1138 AliError(Form("Could not receive track %d", iTracks));
\r
1143 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
\r
1145 if(aodTrack->Charge() == 0) continue;
\r
1146 // Kinematics cuts from ESD track cuts
\r
1147 if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;
\r
1148 if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;
\r
1150 //=================PID (so far only for electron rejection)==========================//
\r
1151 if(fElectronRejection) {
\r
1152 // get the electron nsigma
\r
1153 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
\r
1155 // check only for given momentum range
\r
1156 if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){
\r
1157 //look only at electron nsigma
\r
1158 if(!fElectronOnlyRejection) {
\r
1159 //Make the decision based on the n-sigma of electrons
\r
1160 if(nSigma < fElectronRejectionNSigma) continue;
\r
1163 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
\r
1164 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
\r
1165 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
\r
1167 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
\r
1168 if(nSigma < fElectronRejectionNSigma
\r
1169 && nSigmaPions > fElectronRejectionNSigma
\r
1170 && nSigmaKaons > fElectronRejectionNSigma
\r
1171 && nSigmaProtons > fElectronRejectionNSigma ) continue;
\r
1174 }//electron rejection
\r
1176 gRefMultiplicityTPC += 1.0;
\r
1179 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
\r
1180 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
\r
1181 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
\r
1183 if(iChannel < 32)
\r
1184 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
\r
1185 else if(iChannel >= 32)
\r
1186 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
\r
1189 //Equalization of gain
\r
1190 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
\r
1192 gRefMultiplicityVZEROA /= gFactorA;
\r
1193 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
\r
1195 gRefMultiplicityVZEROC /= gFactorC;
\r
1196 if((gFactorA != 0)&&(gFactorC != 0))
\r
1197 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
\r
1200 Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);
\r
1202 fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
\r
1204 if(fMultiplicityEstimator == "TPC")
\r
1205 gRefMultiplicity = gRefMultiplicityTPC;
\r
1206 else if(fMultiplicityEstimator == "V0M")
\r
1207 gRefMultiplicity = gRefMultiplicityVZERO;
\r
1208 else if(fMultiplicityEstimator == "V0A")
\r
1209 gRefMultiplicity = gRefMultiplicityVZEROA;
\r
1210 else if(fMultiplicityEstimator == "V0C")
\r
1211 gRefMultiplicity = gRefMultiplicityVZEROC;
\r
1213 return gRefMultiplicity;
\r
1216 //________________________________________________________________________
\r
1217 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
\r
1218 // Get the event plane
\r
1220 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
1222 Float_t gVZEROEventPlane = -10.;
\r
1223 Float_t gReactionPlane = -10.;
\r
1224 Double_t qxTot = 0.0, qyTot = 0.0;
\r
1226 //MC: from reaction plane
\r
1227 if(gAnalysisLevel == "MC"){
\r
1229 AliError("mcEvent not available");
\r
1233 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
1235 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
\r
1237 gReactionPlane = headerH->ReactionPlaneAngle();
\r
1238 //gReactionPlane *= TMath::RadToDeg();
\r
1243 // AOD,ESD,ESDMC: from VZERO Event Plane
\r
1246 AliEventplane *ep = event->GetEventplane();
\r
1248 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
\r
1249 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
\r
1250 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
\r
1251 gReactionPlane = gVZEROEventPlane;
\r
1255 return gReactionPlane;
\r
1258 //________________________________________________________________________
\r
1259 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta,
\r
1263 Double_t gCentrality) {
\r
1264 // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality)
\r
1266 Double_t correction = 1.;
\r
1267 Int_t binEta = 0, binPt = 0, binPhi = 0;
\r
1269 //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);
\r
1270 if(fEtaBinForCorrections != 0) {
\r
1271 Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;
\r
1272 if(fEtaMaxForCorrections != fEtaMinForCorrections)
\r
1273 binEta = (Int_t)((vEta-fEtaMinForCorrections)/widthEta)+1;
\r
1276 if(fPtBinForCorrections != 0) {
\r
1277 Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;
\r
1278 if(fPtMaxForCorrections != fPtMinForCorrections)
\r
1279 binPt = (Int_t)((vPt-fPtMinForCorrections)/widthPt) + 1;
\r
1282 if(fPhiBinForCorrections != 0) {
\r
1283 Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;
\r
1284 if(fPhiMaxForCorrections != fPhiMinForCorrections)
\r
1285 binPhi = (Int_t)((vPhi-fPhiMinForCorrections)/widthPhi)+ 1;
\r
1288 Int_t gCentralityInt = -1;
\r
1289 for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){
\r
1290 if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){
\r
1291 gCentralityInt = i;
\r
1296 // centrality not in array --> no correction
\r
1297 if(gCentralityInt < 0){
\r
1302 //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);
\r
1304 if(fHistCorrectionPlus[gCentralityInt]){
\r
1305 if (vCharge > 0) {
\r
1306 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));
\r
1307 //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);
\r
1309 if (vCharge < 0) {
\r
1310 correction = fHistCorrectionMinus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));
\r
1311 //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);
\r
1317 }//centrality in array
\r
1319 if (correction == 0.) {
\r
1320 AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt));
\r
1324 return correction;
\r
1327 //________________________________________________________________________
\r
1328 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){
\r
1329 // Returns TObjArray with tracks after all track cuts (only for AOD!)
\r
1330 // Fills QA histograms
\r
1332 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
1334 //output TObjArray holding all good tracks
\r
1335 TObjArray* tracksAccepted = new TObjArray;
\r
1336 tracksAccepted->SetOwner(kTRUE);
\r
1345 if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
\r
1346 // Loop over tracks in event
\r
1348 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1349 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
\r
1351 AliError(Form("Could not receive track %d", iTracks));
\r
1357 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
\r
1358 // take only TPC only tracks
\r
1359 //fHistTrackStats->Fill(aodTrack->GetFilterMap());
\r
1360 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
\r
1361 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
\r
1364 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
\r
1367 // additional check on kPrimary flag
\r
1368 if(fCheckPrimaryFlagAOD){
\r
1369 if(aodTrack->GetType() != AliAODTrack::kPrimary)
\r
1374 vCharge = aodTrack->Charge();
\r
1375 vEta = aodTrack->Eta();
\r
1376 vY = aodTrack->Y();
\r
1377 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
\r
1378 vPt = aodTrack->Pt();
\r
1380 //===========================PID (so far only for electron rejection)===============================//
\r
1381 if(fElectronRejection) {
\r
1383 // get the electron nsigma
\r
1384 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
\r
1386 //Fill QA before the PID
\r
1387 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1388 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
\r
1389 //end of QA-before pid
\r
1391 // check only for given momentum range
\r
1392 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
\r
1394 //look only at electron nsigma
\r
1395 if(!fElectronOnlyRejection){
\r
1397 //Make the decision based on the n-sigma of electrons
\r
1398 if(nSigma < fElectronRejectionNSigma) continue;
\r
1402 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
\r
1403 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
\r
1404 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
\r
1406 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
\r
1407 if(nSigma < fElectronRejectionNSigma
\r
1408 && nSigmaPions > fElectronRejectionNSigma
\r
1409 && nSigmaKaons > fElectronRejectionNSigma
\r
1410 && nSigmaProtons > fElectronRejectionNSigma ) continue;
\r
1414 //Fill QA after the PID
\r
1415 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1416 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
\r
1419 //===========================end of PID (so far only for electron rejection)===============================//
\r
1421 //+++++++++++++++++++++++++++++//
\r
1422 //===========================PID===============================//
\r
1424 Double_t prob[AliPID::kSPECIES]={0.};
\r
1425 Double_t probTPC[AliPID::kSPECIES]={0.};
\r
1426 Double_t probTOF[AliPID::kSPECIES]={0.};
\r
1427 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
\r
1429 Double_t nSigma = 0.;
\r
1430 Double_t nSigmaTPC = 0.; //++++
\r
1431 Double_t nSigmaTOF = 0.; //++++
\r
1432 Double_t nSigmaTPCTOF = 0.; //++++
\r
1433 UInt_t detUsedTPC = 0;
\r
1434 UInt_t detUsedTOF = 0;
\r
1435 UInt_t detUsedTPCTOF = 0;
\r
1437 nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
\r
1438 nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
\r
1439 nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);
\r
1441 //Decide what detector configuration we want to use
\r
1442 switch(fPidDetectorConfig) {
\r
1444 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
\r
1445 nSigma = nSigmaTPC;
\r
1446 detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
\r
1447 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1448 prob[iSpecies] = probTPC[iSpecies];
\r
1451 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
\r
1452 nSigma = nSigmaTOF;
\r
1453 detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
\r
1454 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1455 prob[iSpecies] = probTOF[iSpecies];
\r
1458 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
\r
1459 nSigma = nSigmaTPCTOF;
\r
1460 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);
\r
1461 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1462 prob[iSpecies] = probTPCTOF[iSpecies];
\r
1466 }//end switch: define detector mask
\r
1468 //Filling the PID QA
\r
1469 Double_t tofTime = -999., length = 999., tof = -999.;
\r
1470 Double_t c = TMath::C()*1.E-9;// m/ns
\r
1471 Double_t beta = -999.;
\r
1472 if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) &&
\r
1473 (aodTrack->IsOn(AliAODTrack::kTIME)) ) {
\r
1474 tofTime = aodTrack->GetTOFsignal();//in ps
\r
1475 length = aodTrack->GetIntegratedLength();
\r
1476 tof = tofTime*1E-3; // ns
\r
1479 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
\r
1483 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
\r
1487 length = length*0.01; // in meters
\r
1489 beta = length/tof;
\r
1491 fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
\r
1492 fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
\r
1493 fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);
\r
1496 fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1497 fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
\r
1498 fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC);
\r
1500 fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
\r
1502 //combined TPC&TOF
\r
1503 fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
\r
1504 fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF);
\r
1505 Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF);
\r
1507 //end of QA-before pid
\r
1509 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
\r
1510 //Make the decision based on the n-sigma
\r
1511 if(fUsePIDnSigma) {
\r
1512 if(nSigma > fPIDNSigma) continue;
\r
1514 fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF);
\r
1515 fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC);
\r
1516 fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF);
\r
1518 //Make the decision based on the bayesian
\r
1519 else if(fUsePIDPropabilities) {
\r
1520 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
\r
1521 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
\r
1523 fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
\r
1524 fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
\r
1525 fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
\r
1529 //Fill QA after the PID
\r
1530 fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
\r
1531 fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1532 fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
\r
1535 //===========================PID===============================//
\r
1536 //+++++++++++++++++++++++++++++//
\r
1539 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1540 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1543 // Kinematics cuts from ESD track cuts
\r
1544 if( vPt < fPtMin || vPt > fPtMax) continue;
\r
1545 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
1547 // Extra DCA cuts (for systematic studies [!= -1])
\r
1548 if( fDCAxyCut != -1 && fDCAzCut != -1){
\r
1549 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
\r
1550 continue; // 2D cut
\r
1554 // Extra TPC cuts (for systematic studies [!= -1])
\r
1555 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
\r
1558 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
\r
1562 // fill QA histograms
\r
1563 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
\r
1564 fHistDCA->Fill(dcaZ,dcaXY);
\r
1565 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
\r
1566 fHistPt->Fill(vPt,gCentrality);
\r
1567 fHistEta->Fill(vEta,gCentrality);
\r
1568 fHistRapidity->Fill(vY,gCentrality);
\r
1569 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
1570 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
1571 fHistPhi->Fill(vPhi,gCentrality);
\r
1572 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
1573 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
1575 //=======================================correction
\r
1576 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
1577 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
1579 // add the track to the TObjArray
\r
1580 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
1584 //==============================================================================================================
\r
1585 else if(gAnalysisLevel == "MCAOD") {
\r
1587 AliMCEvent* mcEvent = MCEvent();
\r
1589 AliError("ERROR: Could not retrieve MC event");
\r
1593 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
\r
1594 AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
\r
1596 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
\r
1600 if(!aodTrack->IsPhysicalPrimary()) continue;
\r
1602 vCharge = aodTrack->Charge();
\r
1603 vEta = aodTrack->Eta();
\r
1604 vY = aodTrack->Y();
\r
1605 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
\r
1606 vPt = aodTrack->Pt();
\r
1608 // Kinematics cuts from ESD track cuts
\r
1609 if( vPt < fPtMin || vPt > fPtMax) continue;
\r
1610 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
1612 // Remove neutral tracks
\r
1613 if( vCharge == 0 ) continue;
\r
1615 //Exclude resonances
\r
1616 if(fExcludeResonancesInMC) {
\r
1618 Bool_t kExcludeParticle = kFALSE;
\r
1619 Int_t gMotherIndex = aodTrack->GetMother();
\r
1620 if(gMotherIndex != -1) {
\r
1621 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
\r
1623 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
\r
1624 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
\r
1625 //if(pdgCodeOfMother == 113) {
\r
1626 if(pdgCodeOfMother == 113 // rho0
\r
1627 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
\r
1628 // || pdgCodeOfMother == 221 // eta
\r
1629 // || pdgCodeOfMother == 331 // eta'
\r
1630 // || pdgCodeOfMother == 223 // omega
\r
1631 // || pdgCodeOfMother == 333 // phi
\r
1632 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
\r
1633 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
\r
1634 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
\r
1635 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
\r
1636 || pdgCodeOfMother == 111 // pi0 Dalitz
\r
1637 || pdgCodeOfMother == 22 // photon
\r
1639 kExcludeParticle = kTRUE;
\r
1644 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
\r
1645 if(kExcludeParticle) continue;
\r
1648 //Exclude electrons with PDG
\r
1649 if(fExcludeElectronsInMC) {
\r
1651 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
\r
1655 // fill QA histograms
\r
1656 fHistPt->Fill(vPt,gCentrality);
\r
1657 fHistEta->Fill(vEta,gCentrality);
\r
1658 fHistRapidity->Fill(vY,gCentrality);
\r
1659 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
1660 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
1661 fHistPhi->Fill(vPhi,gCentrality);
\r
1662 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
1663 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
1665 //=======================================correction
\r
1666 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
1667 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
1669 // add the track to the TObjArray
\r
1670 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
1674 //==============================================================================================================
\r
1676 //==============================================================================================================
\r
1677 else if(gAnalysisLevel == "MCAODrec") {
\r
1679 /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
\r
1681 printf("ERROR: fAOD not available\n");
\r
1685 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
\r
1687 AliError("No array of MC particles found !!!");
\r
1690 AliMCEvent* mcEvent = MCEvent();
\r
1692 AliError("ERROR: Could not retrieve MC event");
\r
1693 return tracksAccepted;
\r
1696 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1697 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
\r
1699 AliError(Form("Could not receive track %d", iTracks));
\r
1703 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
\r
1704 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
\r
1706 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
\r
1708 vCharge = aodTrack->Charge();
\r
1709 vEta = aodTrack->Eta();
\r
1710 vY = aodTrack->Y();
\r
1711 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
\r
1712 vPt = aodTrack->Pt();
\r
1714 //===========================use MC information for Kinematics===============================//
\r
1715 if(fUseMCforKinematics){
\r
1717 Int_t label = TMath::Abs(aodTrack->GetLabel());
\r
1718 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
\r
1721 vCharge = AODmcTrack->Charge();
\r
1722 vEta = AODmcTrack->Eta();
\r
1723 vY = AODmcTrack->Y();
\r
1724 vPhi = AODmcTrack->Phi();// * TMath::RadToDeg();
\r
1725 vPt = AODmcTrack->Pt();
\r
1728 AliDebug(1, "no MC particle for this track");
\r
1732 //===========================end of use MC information for Kinematics========================//
\r
1735 //===========================PID (so far only for electron rejection)===============================//
\r
1736 if(fElectronRejection) {
\r
1738 // get the electron nsigma
\r
1739 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
\r
1741 //Fill QA before the PID
\r
1742 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1743 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
\r
1744 //end of QA-before pid
\r
1746 // check only for given momentum range
\r
1747 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
\r
1749 //look only at electron nsigma
\r
1750 if(!fElectronOnlyRejection){
\r
1752 //Make the decision based on the n-sigma of electrons
\r
1753 if(nSigma < fElectronRejectionNSigma) continue;
\r
1757 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
\r
1758 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
\r
1759 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
\r
1761 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
\r
1762 if(nSigma < fElectronRejectionNSigma
\r
1763 && nSigmaPions > fElectronRejectionNSigma
\r
1764 && nSigmaKaons > fElectronRejectionNSigma
\r
1765 && nSigmaProtons > fElectronRejectionNSigma ) continue;
\r
1769 //Fill QA after the PID
\r
1770 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
\r
1771 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
\r
1774 //===========================end of PID (so far only for electron rejection)===============================//
\r
1776 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1777 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
\r
1779 // Kinematics cuts from ESD track cuts
\r
1780 if( vPt < fPtMin || vPt > fPtMax) continue;
\r
1781 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
1783 // Extra DCA cuts (for systematic studies [!= -1])
\r
1784 if( fDCAxyCut != -1 && fDCAzCut != -1){
\r
1785 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
\r
1786 continue; // 2D cut
\r
1790 // Extra TPC cuts (for systematic studies [!= -1])
\r
1791 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
\r
1794 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
\r
1798 //Exclude resonances
\r
1799 if(fExcludeResonancesInMC) {
\r
1801 Bool_t kExcludeParticle = kFALSE;
\r
1803 Int_t label = TMath::Abs(aodTrack->GetLabel());
\r
1804 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
\r
1807 //if (AODmcTrack->IsPhysicalPrimary()){
\r
1809 Int_t gMotherIndex = AODmcTrack->GetMother();
\r
1810 if(gMotherIndex != -1) {
\r
1811 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
\r
1813 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
\r
1814 if(pdgCodeOfMother == 113 // rho0
\r
1815 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
\r
1816 // || pdgCodeOfMother == 221 // eta
\r
1817 // || pdgCodeOfMother == 331 // eta'
\r
1818 // || pdgCodeOfMother == 223 // omega
\r
1819 // || pdgCodeOfMother == 333 // phi
\r
1820 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
\r
1821 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
\r
1822 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
\r
1823 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
\r
1824 || pdgCodeOfMother == 111 // pi0 Dalitz
\r
1825 || pdgCodeOfMother == 22 // photon
\r
1827 kExcludeParticle = kTRUE;
\r
1832 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
\r
1833 if(kExcludeParticle) continue;
\r
1836 //Exclude electrons with PDG
\r
1837 if(fExcludeElectronsInMC) {
\r
1839 Int_t label = TMath::Abs(aodTrack->GetLabel());
\r
1840 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
\r
1843 if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
\r
1847 // fill QA histograms
\r
1848 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
\r
1849 fHistDCA->Fill(dcaZ,dcaXY);
\r
1850 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
\r
1851 fHistPt->Fill(vPt,gCentrality);
\r
1852 fHistEta->Fill(vEta,gCentrality);
\r
1853 fHistRapidity->Fill(vY,gCentrality);
\r
1854 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
1855 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
1856 fHistPhi->Fill(vPhi,gCentrality);
\r
1857 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
1858 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
1860 //=======================================correction
\r
1861 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
1862 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
1864 // add the track to the TObjArray
\r
1865 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
1868 //==============================================================================================================
\r
1870 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
\r
1872 AliESDtrack *trackTPC = NULL;
\r
1873 AliMCParticle *mcTrack = NULL;
\r
1875 AliMCEvent* mcEvent = NULL;
\r
1877 // for MC ESDs use also MC information
\r
1878 if(gAnalysisLevel == "MCESD"){
\r
1879 mcEvent = MCEvent();
\r
1881 AliError("mcEvent not available");
\r
1882 return tracksAccepted;
\r
1886 // Loop over tracks in event
\r
1887 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
\r
1888 AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));
\r
1890 AliError(Form("Could not receive track %d", iTracks));
\r
1894 // for MC ESDs use also MC information --> MC track not used further???
\r
1895 if(gAnalysisLevel == "MCESD"){
\r
1896 Int_t label = TMath::Abs(track->GetLabel());
\r
1897 if(label > mcEvent->GetNumberOfTracks()) continue;
\r
1898 if(label > mcEvent->GetNumberOfPrimaries()) continue;
\r
1900 mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
\r
1901 if(!mcTrack) continue;
\r
1904 // take only TPC only tracks
\r
1905 trackTPC = new AliESDtrack();
\r
1906 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
\r
1909 if(fESDtrackCuts)
\r
1910 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
\r
1912 // fill QA histograms
\r
1915 trackTPC->GetImpactParameters(b,bCov);
\r
1916 if (bCov[0]<=0 || bCov[2]<=0) {
\r
1917 AliDebug(1, "Estimated b resolution lower or equal zero!");
\r
1918 bCov[0]=0; bCov[2]=0;
\r
1921 Int_t nClustersTPC = -1;
\r
1922 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
\r
1923 //nClustersTPC = track->GetTPCclusters(0); // global track
\r
1924 Float_t chi2PerClusterTPC = -1;
\r
1925 if (nClustersTPC!=0) {
\r
1926 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
\r
1927 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
\r
1930 //===========================PID===============================//
\r
1932 Double_t prob[AliPID::kSPECIES]={0.};
\r
1933 Double_t probTPC[AliPID::kSPECIES]={0.};
\r
1934 Double_t probTOF[AliPID::kSPECIES]={0.};
\r
1935 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
\r
1937 Double_t nSigma = 0.;
\r
1938 UInt_t detUsedTPC = 0;
\r
1939 UInt_t detUsedTOF = 0;
\r
1940 UInt_t detUsedTPCTOF = 0;
\r
1942 //Decide what detector configuration we want to use
\r
1943 switch(fPidDetectorConfig) {
\r
1945 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
\r
1946 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
\r
1947 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
\r
1948 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1949 prob[iSpecies] = probTPC[iSpecies];
\r
1952 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
\r
1953 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
\r
1954 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
\r
1955 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1956 prob[iSpecies] = probTOF[iSpecies];
\r
1959 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
\r
1960 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
\r
1961 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
\r
1962 prob[iSpecies] = probTPCTOF[iSpecies];
\r
1966 }//end switch: define detector mask
\r
1968 //Filling the PID QA
\r
1969 Double_t tofTime = -999., length = 999., tof = -999.;
\r
1970 Double_t c = TMath::C()*1.E-9;// m/ns
\r
1971 Double_t beta = -999.;
\r
1972 Double_t nSigmaTOFForParticleOfInterest = -999.;
\r
1973 if ( (track->IsOn(AliESDtrack::kTOFin)) &&
\r
1974 (track->IsOn(AliESDtrack::kTIME)) ) {
\r
1975 tofTime = track->GetTOFsignal();//in ps
\r
1976 length = track->GetIntegratedLength();
\r
1977 tof = tofTime*1E-3; // ns
\r
1980 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
\r
1984 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
\r
1988 length = length*0.01; // in meters
\r
1990 beta = length/tof;
\r
1992 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
\r
1993 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
\r
1994 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
\r
1995 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
\r
1999 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
\r
2000 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
\r
2001 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
\r
2002 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
\r
2003 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
\r
2004 //end of QA-before pid
\r
2006 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
\r
2007 //Make the decision based on the n-sigma
\r
2008 if(fUsePIDnSigma) {
\r
2009 if(nSigma > fPIDNSigma) continue;}
\r
2011 //Make the decision based on the bayesian
\r
2012 else if(fUsePIDPropabilities) {
\r
2013 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
\r
2014 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
\r
2017 //Fill QA after the PID
\r
2018 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
\r
2019 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
\r
2020 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
\r
2022 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
\r
2023 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
\r
2024 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
\r
2025 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
\r
2028 //===========================PID===============================//
\r
2029 vCharge = trackTPC->Charge();
\r
2030 vY = trackTPC->Y();
\r
2031 vEta = trackTPC->Eta();
\r
2032 vPhi = trackTPC->Phi();// * TMath::RadToDeg();
\r
2033 vPt = trackTPC->Pt();
\r
2034 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
\r
2035 fHistDCA->Fill(b[1],b[0]);
\r
2036 fHistChi2->Fill(chi2PerClusterTPC,gCentrality);
\r
2037 fHistPt->Fill(vPt,gCentrality);
\r
2038 fHistEta->Fill(vEta,gCentrality);
\r
2039 fHistPhi->Fill(vPhi,gCentrality);
\r
2040 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
2041 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
2042 fHistRapidity->Fill(vY,gCentrality);
\r
2043 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
2044 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
2046 //=======================================correction
\r
2047 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
2048 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
2050 // add the track to the TObjArray
\r
2051 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
2057 else if(gAnalysisLevel == "MC"){
\r
2059 AliError("mcEvent not available");
\r
2063 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
\r
2065 // Loop over tracks in event
\r
2066 for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
\r
2067 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
\r
2069 AliError(Form("Could not receive particle %d", iTracks));
\r
2073 //exclude non stable particles
\r
2074 if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
\r
2076 vCharge = track->Charge();
\r
2077 vEta = track->Eta();
\r
2078 vPt = track->Pt();
\r
2081 if( vPt < fPtMin || vPt > fPtMax)
\r
2084 if( vEta < fEtaMin || vEta > fEtaMax) continue;
\r
2086 else if (fUsePID){
\r
2087 if( vY < fEtaMin || vY > fEtaMax) continue;
\r
2090 // Remove neutral tracks
\r
2091 if( vCharge == 0 ) continue;
\r
2093 //analyze one set of particles
\r
2094 if(fUseMCPdgCode) {
\r
2095 TParticle *particle = track->Particle();
\r
2096 if(!particle) continue;
\r
2098 Int_t gPdgCode = particle->GetPdgCode();
\r
2099 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
\r
2103 //Use the acceptance parameterization
\r
2104 if(fAcceptanceParameterization) {
\r
2105 Double_t gRandomNumber = gRandom->Rndm();
\r
2106 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
\r
2110 //Exclude resonances
\r
2111 if(fExcludeResonancesInMC) {
\r
2112 TParticle *particle = track->Particle();
\r
2113 if(!particle) continue;
\r
2115 Bool_t kExcludeParticle = kFALSE;
\r
2116 Int_t gMotherIndex = particle->GetFirstMother();
\r
2117 if(gMotherIndex != -1) {
\r
2118 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
\r
2120 TParticle *motherParticle = motherTrack->Particle();
\r
2121 if(motherParticle) {
\r
2122 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
\r
2123 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
\r
2124 if(pdgCodeOfMother == 113 // rho0
\r
2125 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
\r
2126 // || pdgCodeOfMother == 221 // eta
\r
2127 // || pdgCodeOfMother == 331 // eta'
\r
2128 // || pdgCodeOfMother == 223 // omega
\r
2129 // || pdgCodeOfMother == 333 // phi
\r
2130 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
\r
2131 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
\r
2132 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
\r
2133 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
\r
2134 || pdgCodeOfMother == 111 // pi0 Dalitz
\r
2136 kExcludeParticle = kTRUE;
\r
2142 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
\r
2143 if(kExcludeParticle) continue;
\r
2146 //Exclude electrons with PDG
\r
2147 if(fExcludeElectronsInMC) {
\r
2149 TParticle *particle = track->Particle();
\r
2152 if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
\r
2156 vPhi = track->Phi();
\r
2157 //Printf("phi (before): %lf",vPhi);
\r
2159 fHistPt->Fill(vPt,gCentrality);
\r
2160 fHistEta->Fill(vEta,gCentrality);
\r
2161 fHistPhi->Fill(vPhi,gCentrality);
\r
2162 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
\r
2163 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
\r
2164 //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
\r
2165 fHistRapidity->Fill(vY,gCentrality);
\r
2166 //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
\r
2167 //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
\r
2168 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
\r
2169 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
\r
2171 //Flow after burner
\r
2172 if(fUseFlowAfterBurner) {
\r
2173 Double_t precisionPhi = 0.001;
\r
2174 Int_t maxNumberOfIterations = 100;
\r
2176 Double_t phi0 = vPhi;
\r
2177 Double_t gV2 = fDifferentialV2->Eval(vPt);
\r
2179 for (Int_t j = 0; j < maxNumberOfIterations; j++) {
\r
2180 Double_t phiprev = vPhi;
\r
2181 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
\r
2182 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
\r
2184 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
\r
2186 //Printf("phi (after): %lf\n",vPhi);
\r
2187 Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
\r
2188 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
\r
2189 fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
\r
2191 Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
\r
2192 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
\r
2193 fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
\r
2197 //vPhi *= TMath::RadToDeg();
\r
2199 //=======================================correction
\r
2200 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
\r
2201 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
2203 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
\r
2205 }//MC event object
\r
2208 return tracksAccepted;
\r
2211 //________________________________________________________________________
\r
2212 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
\r
2213 // Clones TObjArray and returns it with tracks after shuffling the charges
\r
2215 TObjArray* tracksShuffled = new TObjArray;
\r
2216 tracksShuffled->SetOwner(kTRUE);
\r
2218 vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks
\r
2220 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
\r
2222 AliVParticle* track = (AliVParticle*) tracks->At(i);
\r
2223 chargeVector->push_back(track->Charge());
\r
2226 random_shuffle(chargeVector->begin(), chargeVector->end());
\r
2228 for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
\r
2229 AliVParticle* track = (AliVParticle*) tracks->At(i);
\r
2230 //==============================correction
\r
2231 Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
\r
2232 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
\r
2233 tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
\r
2236 delete chargeVector;
\r
2238 return tracksShuffled;
\r
2241 //________________________________________________________________________
\r
2242 void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,
\r
2243 const char* lhcPeriod) {
\r
2244 //Function to setup the VZERO gain equalization
\r
2245 //============Get the equilization map============//
\r
2246 TFile *calibrationFile = TFile::Open(filename);
\r
2247 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
\r
2248 Printf("No calibration file found!!!");
\r
2252 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
\r
2254 Printf("Calibration TList not found!!!");
\r
2258 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
\r
2259 if(!fHistVZEROAGainEqualizationMap) {
\r
2260 Printf("VZERO-A calibration object not found!!!");
\r
2263 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
\r
2264 if(!fHistVZEROCGainEqualizationMap) {
\r
2265 Printf("VZERO-C calibration object not found!!!");
\r
2269 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
\r
2270 if(!fHistVZEROChannelGainEqualizationMap) {
\r
2271 Printf("VZERO channel calibration object not found!!!");
\r
2276 //________________________________________________________________________
\r
2277 Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run,
\r
2280 if(!fHistVZEROAGainEqualizationMap) return 1.0;
\r
2282 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
\r
2283 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
\r
2284 if(gRunNumber == run)
\r
2285 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
\r
2291 //________________________________________________________________________
\r
2292 Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run,
\r
2293 const char* side) {
\r
2295 if(!fHistVZEROAGainEqualizationMap) return 1.0;
\r
2297 TString gVZEROSide = side;
\r
2298 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
\r
2299 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
\r
2300 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
\r
2301 if(gRunNumber == run) {
\r
2302 if(gVZEROSide == "A")
\r
2303 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
\r
2304 else if(gVZEROSide == "C")
\r
2305 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
\r
2312 //____________________________________________________________________
\r
2313 Bool_t AliAnalysisTaskBFPsi::AcceptEventCentralityWeight(Double_t centrality)
\r
2315 // copied from AliAnalysisTaskPhiCorrelations
\r
2317 // rejects "randomly" events such that the centrality gets flat
\r
2318 // uses fCentralityWeights histogram
\r
2320 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
\r
2322 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
\r
2323 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
\r
2325 Bool_t result = kFALSE;
\r
2326 if (centralityDigits < weight)
\r
2329 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
\r
2334 //________________________________________________________________________
\r
2335 void AliAnalysisTaskBFPsi::FinishTaskOutput(){
\r
2336 //Printf("END BF");
\r
2339 AliError("fBalance not available");
\r
2342 if(fRunShuffling) {
\r
2343 if (!fShuffledBalance) {
\r
2344 AliError("fShuffledBalance not available");
\r
2351 //________________________________________________________________________
\r
2352 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
\r
2353 // Draw result to the screen
\r
2354 // Called once at the end of the query
\r
2356 // not implemented ...
\r
git://git.uio.no / u / mrichter / AliRoot.git / blob