4 #include "TLorentzVector.h"
5 #include "TGraphErrors.h"
12 #include "AliAnalysisTaskSE.h"
13 #include "AliAnalysisManager.h"
15 #include "AliESDVertex.h"
16 #include "AliESDEvent.h"
17 #include "AliESDInputHandler.h"
18 #include "AliAODEvent.h"
19 #include "AliAODTrack.h"
20 #include "AliAODInputHandler.h"
21 #include "AliGenEventHeader.h"
22 #include "AliGenHijingEventHeader.h"
23 #include "AliMCEventHandler.h"
24 #include "AliMCEvent.h"
26 #include "AliESDtrackCuts.h"
31 #include "AliPIDResponse.h"
32 #include "AliPIDCombined.h"
34 #include "AliAnalysisTaskBF.h"
35 #include "AliBalance.h"
38 // Analysis task for the BF code
39 // Authors: Panos.Christakoglou@nikhef.nl
41 ClassImp(AliAnalysisTaskBF)
43 //________________________________________________________________________
44 AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name)
45 : AliAnalysisTaskSE(name),
47 fRunShuffling(kFALSE),
73 fHistdEdxVsPTPCbeforePID(NULL),
74 fHistBetavsPTOFbeforePID(NULL),
75 fHistProbTPCvsPtbeforePID(NULL),
76 fHistProbTOFvsPtbeforePID(NULL),
77 fHistProbTPCTOFvsPtbeforePID(NULL),
78 fHistNSigmaTPCvsPtbeforePID(NULL),
79 fHistNSigmaTOFvsPtbeforePID(NULL),
80 fHistdEdxVsPTPCafterPID(NULL),
81 fHistBetavsPTOFafterPID(NULL),
82 fHistProbTPCvsPtafterPID(NULL),
83 fHistProbTOFvsPtafterPID(NULL),
84 fHistProbTPCTOFvsPtafterPID(NULL),
85 fHistNSigmaTPCvsPtafterPID(NULL),
86 fHistNSigmaTOFvsPtafterPID(NULL),
89 fParticleOfInterest(kPion),
90 fPidDetectorConfig(kTPCTOF),
93 fUsePIDPropabilities(kFALSE),
95 fMinAcceptedPIDProbability(0.8),
97 fCentralityEstimator("V0M"),
98 fUseCentrality(kFALSE),
99 fCentralityPercentileMin(0.),
100 fCentralityPercentileMax(5.),
101 fImpactParameterMin(0.),
102 fImpactParameterMax(20.),
103 fUseMultiplicity(kFALSE),
104 fNumberOfAcceptedTracksMin(0),
105 fNumberOfAcceptedTracksMax(10000),
106 fHistNumberOfAcceptedTracks(0),
107 fUseOfflineTrigger(kFALSE),
111 fAODtrackCutBit(128),
119 fNClustersTPCCut(-1),
120 fAcceptanceParameterization(0),
122 fUseFlowAfterBurner(kFALSE),
123 fExcludeResonancesInMC(kFALSE),
124 fUseMCPdgCode(kFALSE),
125 fPDGCodeToBeAnalyzed(-1) {
127 // Define input and output slots here
128 // Input slot #0 works with a TChain
129 DefineInput(0, TChain::Class());
130 // Output slot #0 writes into a TH1 container
131 DefineOutput(1, TList::Class());
132 DefineOutput(2, TList::Class());
133 DefineOutput(3, TList::Class());
134 DefineOutput(4, TList::Class());
137 //________________________________________________________________________
138 AliAnalysisTaskBF::~AliAnalysisTaskBF() {
141 // delete fShuffledBalance;
146 // delete fHistEventStats;
147 // delete fHistTrackStats;
161 //________________________________________________________________________
162 void AliAnalysisTaskBF::UserCreateOutputObjects() {
166 // global switch disabling the reference
167 // (to avoid "Replacing existing TH1" if several wagons are created in train)
168 Bool_t oldStatus = TH1::AddDirectoryStatus();
169 TH1::AddDirectory(kFALSE);
172 fBalance = new AliBalance();
173 fBalance->SetAnalysisLevel("ESD");
174 //fBalance->SetNumberOfBins(-1,16);
175 fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
178 if(!fShuffledBalance) {
179 fShuffledBalance = new AliBalance();
180 fShuffledBalance->SetAnalysisLevel("ESD");
181 //fShuffledBalance->SetNumberOfBins(-1,16);
182 fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
188 fList->SetName("listQA");
191 //Balance Function list
192 fListBF = new TList();
193 fListBF->SetName("listBF");
197 fListBFS = new TList();
198 fListBFS->SetName("listBFShuffled");
199 fListBFS->SetOwner();
204 fHistListPIDQA = new TList();
205 fHistListPIDQA->SetName("listQAPID");
206 fHistListPIDQA->SetOwner();
210 TString gCutName[4] = {"Total","Offline trigger",
211 "Vertex","Analyzed"};
213 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
215 if ((gAnalysisLevel == "ESD") || (gAnalysisLevel == "AOD") || (gAnalysisLevel == "MCESD")) {
216 fHistEventStats = new TH2D("fHistEventStats",
217 "Event statistics;;Centrality",
218 4,0.5,4.5, 100,0,100);
221 if (gAnalysisLevel == "MC"){
222 fHistEventStats = new TH2D("fHistEventStats",
223 "Event statistics;;Centrality",
224 4,0.5,4.5, 10000,0,15);
228 for(Int_t i = 1; i <= 4; i++)
229 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
230 fList->Add(fHistEventStats);
232 TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
233 fHistCentStats = new TH2F("fHistCentStats",
234 "Centrality statistics;;Cent percentile",
235 9,-0.5,8.5,220,-5,105);
236 for(Int_t i = 1; i <= 9; i++)
237 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
238 fList->Add(fHistCentStats);
240 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);
241 fList->Add(fHistTriggerStats);
243 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);
244 fList->Add(fHistTrackStats);
246 fHistNumberOfAcceptedTracks = new TH2D("fHistNumberOfAcceptedTracks",";N_{acc.};;Centrality",4001,-0.5,4000.5,100,0,100);
247 fList->Add(fHistNumberOfAcceptedTracks);
249 // Vertex distributions
250 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
252 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
254 fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
258 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
259 fList->Add(fHistClus);
260 fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);
261 fList->Add(fHistChi2);
262 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
263 fList->Add(fHistDCA);
264 fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10);
266 fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2);
267 fList->Add(fHistEta);
268 fHistRapidity = new TH1F("fHistRapidity","y distribution",200,-2,2);
269 fList->Add(fHistRapidity);
270 fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380);
271 fList->Add(fHistPhi);
272 fHistPhiBefore = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());
273 fList->Add(fHistPhiBefore);
274 fHistPhiAfter = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());
275 fList->Add(fHistPhiAfter);
276 fHistPhiPos = new TH1F("fHistPhiPos","#phi distribution for positive particles",200,-20,380);
277 fList->Add(fHistPhiPos);
278 fHistPhiNeg = new TH1F("fHistPhiNeg","#phi distribution for negative particles",200,-20,380);
279 fList->Add(fHistPhiNeg);
280 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
281 fList->Add(fHistV0M);
282 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
283 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
284 for(Int_t i = 1; i <= 6; i++)
285 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
286 fList->Add(fHistRefTracks);
288 // QA histograms for HBTinspired and Conversion cuts
289 fList->Add(fBalance->GetQAHistHBTbefore());
290 fList->Add(fBalance->GetQAHistHBTafter());
291 fList->Add(fBalance->GetQAHistConversionbefore());
292 fList->Add(fBalance->GetQAHistConversionafter());
294 // Balance function histograms
295 // Initialize histograms if not done yet
296 if(!fBalance->GetHistNp(0)){
297 AliWarning("Histograms not yet initialized! --> Will be done now");
298 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
299 fBalance->InitHistograms();
303 if(!fShuffledBalance->GetHistNp(0)) {
304 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
305 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
306 fShuffledBalance->InitHistograms();
310 for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
311 fListBF->Add(fBalance->GetHistNp(a));
312 fListBF->Add(fBalance->GetHistNn(a));
313 fListBF->Add(fBalance->GetHistNpn(a));
314 fListBF->Add(fBalance->GetHistNnn(a));
315 fListBF->Add(fBalance->GetHistNpp(a));
316 fListBF->Add(fBalance->GetHistNnp(a));
319 fListBFS->Add(fShuffledBalance->GetHistNp(a));
320 fListBFS->Add(fShuffledBalance->GetHistNn(a));
321 fListBFS->Add(fShuffledBalance->GetHistNpn(a));
322 fListBFS->Add(fShuffledBalance->GetHistNnn(a));
323 fListBFS->Add(fShuffledBalance->GetHistNpp(a));
324 fListBFS->Add(fShuffledBalance->GetHistNnp(a));
328 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
330 //====================PID========================//
332 fPIDCombined = new AliPIDCombined();
333 fPIDCombined->SetDefaultTPCPriors();
335 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
336 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition
338 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
339 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition
341 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
342 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition
344 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
345 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition
347 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
348 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition
350 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
351 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition
353 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
354 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition
356 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
357 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition
359 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
360 fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition
362 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
363 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition
365 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
366 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition
368 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
369 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition
371 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
372 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition
374 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
375 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition
377 //====================PID========================//
381 PostData(2, fListBF);
382 if(fRunShuffling) PostData(3, fListBFS);
383 if(fUsePID) PostData(4, fHistListPIDQA); //PID
385 TH1::AddDirectory(oldStatus);
389 //________________________________________________________________________
390 void AliAnalysisTaskBF::UserExec(Option_t *) {
392 // Called for each event
393 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
395 AliESDtrack *trackTPC = NULL;
397 Int_t gNumberOfAcceptedTracks = 0;
398 Float_t fCentrality = -999.;
400 // for HBT like cuts need magnetic field sign
401 Float_t bSign = 0; // only used in AOD so far
403 // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)
404 vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled
405 vector<Double_t> *chargeVector[9]; // original charge
406 for(Int_t i = 0; i < 9; i++){
407 chargeVectorShuffle[i] = new vector<Double_t>;
408 chargeVector[i] = new vector<Double_t>;
420 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
421 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
425 if(gAnalysisLevel == "ESD") {
426 AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
428 AliError("ERROR: gESD not available");
432 // store offline trigger bits
433 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
435 AliCentrality *centrality = 0x0;
438 centrality = gESD->GetCentrality();
439 fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
442 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
443 fHistEventStats->Fill(1,fCentrality); //all events
444 Bool_t isSelected = kTRUE;
445 if(fUseOfflineTrigger)
446 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
448 fHistEventStats->Fill(2,fCentrality); //triggered events
452 // take only events inside centrality class
453 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
454 fCentralityPercentileMax,
455 fCentralityEstimator.Data()))
457 // centrality QA (V0M)
458 fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());
461 const AliESDVertex *vertex = gESD->GetPrimaryVertex();
463 if(vertex->GetNContributors() > 0) {
464 if(vertex->GetZRes() != 0) {
465 fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
466 if(TMath::Abs(vertex->GetX()) < fVxMax) {
467 if(TMath::Abs(vertex->GetY()) < fVyMax) {
468 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
469 fHistEventStats->Fill(4,fCentrality); //analayzed events
470 fHistVx->Fill(vertex->GetX());
471 fHistVy->Fill(vertex->GetY());
472 fHistVz->Fill(vertex->GetZ());
474 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());
475 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {
476 AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));
478 AliError(Form("ERROR: Could not receive track %d", iTracks));
482 // take only TPC only tracks
483 trackTPC = new AliESDtrack();
484 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
488 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
490 // fill QA histograms
493 trackTPC->GetImpactParameters(b,bCov);
494 if (bCov[0]<=0 || bCov[2]<=0) {
495 AliDebug(1, "Estimated b resolution lower or equal zero!");
496 bCov[0]=0; bCov[2]=0;
499 Int_t nClustersTPC = -1;
500 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
501 //nClustersTPC = track->GetTPCclusters(0); // global track
502 Float_t chi2PerClusterTPC = -1;
503 if (nClustersTPC!=0) {
504 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
505 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
508 //===========================PID===============================//
510 Double_t prob[AliPID::kSPECIES]={0.};
511 Double_t probTPC[AliPID::kSPECIES]={0.};
512 Double_t probTOF[AliPID::kSPECIES]={0.};
513 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
515 Double_t nSigma = 0.;
516 UInt_t detUsedTPC = 0;
517 UInt_t detUsedTOF = 0;
518 UInt_t detUsedTPCTOF = 0;
520 //Decide what detector configuration we want to use
521 switch(fPidDetectorConfig) {
523 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
524 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
525 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
526 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
527 prob[iSpecies] = probTPC[iSpecies];
530 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
531 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
532 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
533 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
534 prob[iSpecies] = probTOF[iSpecies];
537 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
538 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
539 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
540 prob[iSpecies] = probTPCTOF[iSpecies];
544 }//end switch: define detector mask
547 Double_t tofTime = -999., length = 999., tof = -999.;
548 Double_t c = TMath::C()*1.E-9;// m/ns
549 Double_t beta = -999.;
550 Double_t nSigmaTOFForParticleOfInterest = -999.;
551 if ( (track->IsOn(AliESDtrack::kTOFin)) &&
552 (track->IsOn(AliESDtrack::kTIME)) ) {
553 tofTime = track->GetTOFsignal();//in ps
554 length = track->GetIntegratedLength();
555 tof = tofTime*1E-3; // ns
558 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
562 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
566 length = length*0.01; // in meters
570 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
571 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
572 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
573 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
577 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
578 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
579 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
580 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
581 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
582 //end of QA-before pid
584 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
585 //Make the decision based on the n-sigma
587 if(nSigma > fPIDNSigma) continue;}
589 //Make the decision based on the bayesian
590 else if(fUsePIDPropabilities) {
591 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
592 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
595 //Fill QA after the PID
596 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
597 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
598 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
600 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
601 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
602 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
603 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
606 PostData(4, fHistListPIDQA);
608 //===========================PID===============================//
609 vCharge = trackTPC->Charge();
611 vEta = trackTPC->Eta();
612 vPhi = trackTPC->Phi() * TMath::RadToDeg();
614 vPt = trackTPC->Pt();
615 trackTPC->PxPyPz(vP);
616 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
617 fHistDCA->Fill(b[1],b[0]);
618 fHistChi2->Fill(chi2PerClusterTPC);
620 fHistEta->Fill(vEta);
621 fHistPhi->Fill(vPhi);
622 fHistRapidity->Fill(vY);
623 if(vCharge > 0) fHistPhiPos->Fill(vPhi);
624 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
626 // fill charge vector
627 chargeVector[0]->push_back(vCharge);
628 chargeVector[1]->push_back(vY);
629 chargeVector[2]->push_back(vEta);
630 chargeVector[3]->push_back(vPhi);
631 chargeVector[4]->push_back(vP[0]);
632 chargeVector[5]->push_back(vP[1]);
633 chargeVector[6]->push_back(vP[2]);
634 chargeVector[7]->push_back(vPt);
635 chargeVector[8]->push_back(vE);
638 chargeVectorShuffle[0]->push_back(vCharge);
639 chargeVectorShuffle[1]->push_back(vY);
640 chargeVectorShuffle[2]->push_back(vEta);
641 chargeVectorShuffle[3]->push_back(vPhi);
642 chargeVectorShuffle[4]->push_back(vP[0]);
643 chargeVectorShuffle[5]->push_back(vP[1]);
644 chargeVectorShuffle[6]->push_back(vP[2]);
645 chargeVectorShuffle[7]->push_back(vPt);
646 chargeVectorShuffle[8]->push_back(vE);
650 gNumberOfAcceptedTracks += 1;
652 // cout<<"Centrality: "<<fCentrality<<" - Accepted tracks: "<<gNumberOfAcceptedTracks<<endl;
656 }//proper vertex resolution
657 }//proper number of contributors
658 }//vertex object valid
662 //AOD analysis (vertex and track cuts also here!!!!)
663 else if(gAnalysisLevel == "AOD") {
664 AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
666 AliError("ERROR: gAOD not available");
670 // for HBT like cuts need magnetic field sign
671 bSign = (gAOD->GetMagneticField() > 0) ? 1 : -1;
673 AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(gAOD->GetHeader());
674 if(!aodHeader) AliFatal("Not a standard AOD");
676 // store offline trigger bits
677 fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());
680 fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
683 //event selection done in AliAnalysisTaskSE::Exec() --> this is not used
684 fHistEventStats->Fill(1,fCentrality); //all events
685 Bool_t isSelected = kTRUE;
686 if(fUseOfflineTrigger)
687 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
689 fHistEventStats->Fill(2,fCentrality); //triggered events
691 //Centrality stuff (centrality in AOD header)
693 //fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
694 // in OLD AODs (i.e. AOD049) fCentrality can be == 0
698 // QA for centrality estimators
699 fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));
700 fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));
701 fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));
702 fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));
703 fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));
704 fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));
705 fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
706 fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
707 fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
709 // take only events inside centrality class
710 if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax))
713 // centrality QA (V0M)
714 fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());
716 // centrality QA (reference tracks)
717 fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());
718 fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());
719 fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());
720 fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());
721 fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));
722 fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));
723 fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));
724 fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));
725 fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));
728 const AliAODVertex *vertex = gAOD->GetPrimaryVertex();
732 vertex->GetCovarianceMatrix(fCov);
734 if(vertex->GetNContributors() > 0) {
736 fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
737 if(TMath::Abs(vertex->GetX()) < fVxMax) {
738 if(TMath::Abs(vertex->GetY()) < fVyMax) {
739 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
740 fHistEventStats->Fill(4,fCentrality); //analyzed events
741 fHistVx->Fill(vertex->GetX());
742 fHistVy->Fill(vertex->GetY());
743 fHistVz->Fill(vertex->GetZ());
745 //===========================================//
746 TExMap *trackMap = new TExMap();
747 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {
748 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));
750 AliError(Form("ERROR: Could not receive track %d", iTracks));
753 Int_t gID = aodTrack->GetID();
754 //if (!aodTrack->TestFilterBit(fAODtrackCutBit)) trackMap->Add(gID, iTracks);
755 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, iTracks);
757 AliAODTrack* newAodTrack;
758 //===========================================//
760 //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());
761 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {
762 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));
764 AliError(Form("ERROR: Could not receive track %d", iTracks));
769 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
770 //===========================================//
771 // take only TPC only tracks
772 fHistTrackStats->Fill(aodTrack->GetFilterMap());
773 if(!aodTrack->TestFilterBit(fAODtrackCutBit)) continue;
775 Int_t gID = aodTrack->GetID();
776 newAodTrack = gID >= 0 ? aodTrack : dynamic_cast<AliAODTrack*>(gAOD->GetTrack(trackMap->GetValue(-1-gID)));
777 if(!newAodTrack) AliFatal("Not a standard AOD");
778 //Printf("Label: %d - Pt: %lf (old) - %d - Pt: %lf(new)",gID,aodTrack->Pt(), newAodTrack->GetID(), newAodTrack->Pt());
779 //===========================================//
781 //fHistTrackStats->Fill(aodTrack->GetFilterMap());
782 //if(!aodTrack->TestFilterBit(fAODtrackCutBit)) continue;
784 vCharge = aodTrack->Charge();
786 vEta = aodTrack->Eta();
787 vPhi = aodTrack->Phi() * TMath::RadToDeg();
789 vPt = aodTrack->Pt();
790 aodTrack->PxPyPz(vP);
792 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
793 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
796 // Kinematics cuts from ESD track cuts
797 if( vPt < fPtMin || vPt > fPtMax) continue;
800 if( vEta < fEtaMin || vEta > fEtaMax) continue;
804 if( vY < fEtaMin || vY > fEtaMax) continue;
807 // Extra DCA cuts (for systematic studies [!= -1])
808 if( fDCAxyCut != -1 && fDCAzCut != -1){
809 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
814 // Extra TPC cuts (for systematic studies [!= -1])
815 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
818 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
822 //===============================================PID==================================//
824 Double_t prob[AliPID::kSPECIES]={0.};
825 Double_t probTPC[AliPID::kSPECIES]={0.};
826 Double_t probTOF[AliPID::kSPECIES]={0.};
827 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
829 Double_t nSigma = 0.;
830 UInt_t detUsedTPC = 0;
831 UInt_t detUsedTOF = 0;
832 UInt_t detUsedTPCTOF = 0;
834 //Decide what detector configuration we want to use
835 switch(fPidDetectorConfig) {
837 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
838 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest));
839 //detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
840 detUsedTPC = (AliPIDResponse::kDetTPC);
841 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
842 prob[iSpecies] = probTPC[iSpecies];
845 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
846 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest));
847 //detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
848 detUsedTPC = (AliPIDResponse::kDetTPC);
849 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
850 prob[iSpecies] = probTOF[iSpecies];
853 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
854 //detUsedTPCTOF = fPIDCombined->ComputeProbabilities(newAodTrack, fPIDResponse, probTPCTOF);
855 detUsedTPC = (AliPIDResponse::kDetTPC);
856 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
857 prob[iSpecies] = probTPCTOF[iSpecies];
861 }//end switch: define detector mask
864 Double_t tofTime = -999., tof = -999.; //length = 999., tof = -999.;
865 //Double_t c = TMath::C()*1.E-9;// m/ns
866 //Double_t beta = -999.;
867 Double_t nSigmaTOFForParticleOfInterest = -999.;
868 if ( (newAodTrack->IsOn(AliAODTrack::kTOFin)) &&
869 (newAodTrack->IsOn(AliAODTrack::kTIME)) ) {
870 tofTime = newAodTrack->GetTOFsignal();//in ps
871 //length = newAodTrack->GetIntegratedLength();
872 tof = tofTime*1E-3; // ns
875 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
879 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
883 //length = length*0.01; // in meters
887 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest);
888 //fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
889 fHistProbTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]);
890 fHistNSigmaTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest);
894 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest);
895 fHistdEdxVsPTPCbeforePID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal());
896 fHistProbTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]);
897 fHistNSigmaTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest);
898 fHistProbTPCTOFvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
899 //end of QA-before pid
901 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
902 //Make the decision based on the n-sigma
904 if(nSigma > fPIDNSigma) continue;}
906 //Make the decision based on the bayesian
907 else if(fUsePIDPropabilities) {
908 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
909 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
912 //Fill QA after the PID
913 //fHistBetavsPTOFafterPID ->Fill(newAodTrack->P()*newAodTrack->Charge(),beta);
914 fHistProbTOFvsPtafterPID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]);
915 fHistNSigmaTOFvsPtafterPID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest);
917 fHistdEdxVsPTPCafterPID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal());
918 fHistProbTPCvsPtafterPID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]);
919 fHistProbTPCTOFvsPtafterPID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
920 fHistNSigmaTPCvsPtafterPID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest);
923 PostData(4, fHistListPIDQA);
926 //=========================================================PID=================================================================//
928 // fill QA histograms
929 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
930 fHistDCA->Fill(dcaZ,dcaXY);
931 fHistChi2->Fill(aodTrack->Chi2perNDF());
933 fHistEta->Fill(vEta);
934 fHistPhi->Fill(vPhi);
935 fHistRapidity->Fill(vY);
936 if(vCharge > 0) fHistPhiPos->Fill(vPhi);
937 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
939 // fill charge vector
940 chargeVector[0]->push_back(vCharge);
941 chargeVector[1]->push_back(vY);
942 chargeVector[2]->push_back(vEta);
943 chargeVector[3]->push_back(vPhi);
944 chargeVector[4]->push_back(vP[0]);
945 chargeVector[5]->push_back(vP[1]);
946 chargeVector[6]->push_back(vP[2]);
947 chargeVector[7]->push_back(vPt);
948 chargeVector[8]->push_back(vE);
951 chargeVectorShuffle[0]->push_back(vCharge);
952 chargeVectorShuffle[1]->push_back(vY);
953 chargeVectorShuffle[2]->push_back(vEta);
954 chargeVectorShuffle[3]->push_back(vPhi);
955 chargeVectorShuffle[4]->push_back(vP[0]);
956 chargeVectorShuffle[5]->push_back(vP[1]);
957 chargeVectorShuffle[6]->push_back(vP[2]);
958 chargeVectorShuffle[7]->push_back(vPt);
959 chargeVectorShuffle[8]->push_back(vE);
962 gNumberOfAcceptedTracks += 1;
968 }//proper vertex resolution
969 }//proper number of contributors
970 }//vertex object valid
975 if(gAnalysisLevel == "MCESD") {
976 AliMCEvent* mcEvent = MCEvent();
978 AliError("ERROR: mcEvent not available");
982 AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
984 AliError("ERROR: gESD not available");
988 // store offline trigger bits
989 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
991 AliCentrality *centrality = 0x0;
993 centrality = gESD->GetCentrality();
994 fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
997 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
998 fHistEventStats->Fill(1,fCentrality); //all events
999 Bool_t isSelected = kTRUE;
1000 if(fUseOfflineTrigger)
1001 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1003 fHistEventStats->Fill(2,fCentrality); //triggered events
1005 if(fUseCentrality) {
1007 // take only events inside centrality class
1008 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
1009 fCentralityPercentileMax,
1010 fCentralityEstimator.Data()))
1012 // centrality QA (V0M)
1013 fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());
1016 const AliESDVertex *vertex = gESD->GetPrimaryVertex();
1018 if(vertex->GetNContributors() > 0) {
1019 if(vertex->GetZRes() != 0) {
1020 fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
1021 if(TMath::Abs(vertex->GetX()) < fVxMax) {
1022 if(TMath::Abs(vertex->GetY()) < fVyMax) {
1023 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
1024 fHistEventStats->Fill(4,fCentrality); //analayzed events
1025 fHistVx->Fill(vertex->GetX());
1026 fHistVy->Fill(vertex->GetY());
1027 fHistVz->Fill(vertex->GetZ());
1029 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());
1030 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {
1031 AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));
1033 AliError(Form("ERROR: Could not receive track %d", iTracks));
1037 Int_t label = TMath::Abs(track->GetLabel());
1038 if(label > mcEvent->GetNumberOfTracks()) continue;
1039 if(label > mcEvent->GetNumberOfPrimaries()) continue;
1041 AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
1042 if(!mcTrack) continue;
1044 // take only TPC only tracks
1045 trackTPC = new AliESDtrack();
1046 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
1050 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
1052 // fill QA histograms
1055 trackTPC->GetImpactParameters(b,bCov);
1056 if (bCov[0]<=0 || bCov[2]<=0) {
1057 AliDebug(1, "Estimated b resolution lower or equal zero!");
1058 bCov[0]=0; bCov[2]=0;
1061 Int_t nClustersTPC = -1;
1062 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
1063 //nClustersTPC = track->GetTPCclusters(0); // global track
1064 Float_t chi2PerClusterTPC = -1;
1065 if (nClustersTPC!=0) {
1066 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
1067 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
1070 vCharge = trackTPC->Charge();
1072 vEta = trackTPC->Eta();
1073 vPhi = trackTPC->Phi() * TMath::RadToDeg();
1075 vPt = trackTPC->Pt();
1076 trackTPC->PxPyPz(vP);
1078 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
1079 fHistDCA->Fill(b[1],b[0]);
1080 fHistChi2->Fill(chi2PerClusterTPC);
1082 fHistEta->Fill(vEta);
1083 fHistPhi->Fill(vPhi);
1084 fHistRapidity->Fill(vY);
1085 if(vCharge > 0) fHistPhiPos->Fill(vPhi);
1086 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
1088 // fill charge vector
1089 chargeVector[0]->push_back(vCharge);
1090 chargeVector[1]->push_back(vY);
1091 chargeVector[2]->push_back(vEta);
1092 chargeVector[3]->push_back(vPhi);
1093 chargeVector[4]->push_back(vP[0]);
1094 chargeVector[5]->push_back(vP[1]);
1095 chargeVector[6]->push_back(vP[2]);
1096 chargeVector[7]->push_back(vPt);
1097 chargeVector[8]->push_back(vE);
1100 chargeVectorShuffle[0]->push_back(vCharge);
1101 chargeVectorShuffle[1]->push_back(vY);
1102 chargeVectorShuffle[2]->push_back(vEta);
1103 chargeVectorShuffle[3]->push_back(vPhi);
1104 chargeVectorShuffle[4]->push_back(vP[0]);
1105 chargeVectorShuffle[5]->push_back(vP[1]);
1106 chargeVectorShuffle[6]->push_back(vP[2]);
1107 chargeVectorShuffle[7]->push_back(vPt);
1108 chargeVectorShuffle[8]->push_back(vE);
1112 gNumberOfAcceptedTracks += 1;
1115 //cout<<"Centrality: "<<fCentrality<<" - Accepted tracks: "<<gNumberOfAcceptedTracks<<endl;
1119 }//proper vertex resolution
1120 }//proper number of contributors
1121 }//vertex object valid
1126 else if(gAnalysisLevel == "MC") {
1127 AliMCEvent* mcEvent = MCEvent();
1129 AliError("ERROR: mcEvent not available");
1133 //fHistEventStats->Fill(1,fCentrality); //total events
1134 //fHistEventStats->Fill(2,fCentrality); //offline trigger
1136 Double_t gReactionPlane = 0., gImpactParameter = 0.;
1137 if(fUseCentrality) {
1139 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
1141 //Printf("=====================================================");
1142 //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());
1143 //Printf("Impact parameter: %lf",headerH->ImpactParameter());
1144 //Printf("=====================================================");
1145 gReactionPlane = headerH->ReactionPlaneAngle();
1146 gImpactParameter = headerH->ImpactParameter();
1147 fCentrality = gImpactParameter;
1149 fCentrality = gImpactParameter;
1151 // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW)
1152 if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))
1156 fHistEventStats->Fill(1,fCentrality); //total events
1157 fHistEventStats->Fill(2,fCentrality); //offline trigger
1159 AliGenEventHeader *header = mcEvent->GenEventHeader();
1162 TArrayF gVertexArray;
1163 header->PrimaryVertex(gVertexArray);
1164 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
1165 //gVertexArray.At(0),
1166 //gVertexArray.At(1),
1167 //gVertexArray.At(2));
1168 fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
1169 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
1170 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
1171 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
1172 fHistEventStats->Fill(4,fCentrality); //analayzed events
1173 fHistVx->Fill(gVertexArray.At(0));
1174 fHistVy->Fill(gVertexArray.At(1));
1175 fHistVz->Fill(gVertexArray.At(2));
1177 AliInfo(Form("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries()));
1178 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {
1179 AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));
1181 AliError(Form("ERROR: Could not receive particle %d", iTracks));
1185 //exclude non stable particles
1186 if(!mcEvent->IsPhysicalPrimary(iTracks)) continue;
1188 vEta = track->Eta();
1192 if( vPt < fPtMin || vPt > fPtMax)
1195 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1198 if( vY < fEtaMin || vY > fEtaMax) continue;
1201 //analyze one set of particles
1203 TParticle *particle = track->Particle();
1204 if(!particle) continue;
1206 Int_t gPdgCode = particle->GetPdgCode();
1207 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
1211 //Use the acceptance parameterization
1212 if(fAcceptanceParameterization) {
1213 Double_t gRandomNumber = gRandom->Rndm();
1214 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
1218 //Exclude resonances
1219 if(fExcludeResonancesInMC) {
1220 TParticle *particle = track->Particle();
1221 if(!particle) continue;
1223 Bool_t kExcludeParticle = kFALSE;
1224 Int_t gMotherIndex = particle->GetFirstMother();
1225 if(gMotherIndex != -1) {
1226 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1228 TParticle *motherParticle = motherTrack->Particle();
1229 if(motherParticle) {
1230 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
1231 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1232 if(pdgCodeOfMother == 113) {
1233 kExcludeParticle = kTRUE;
1239 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1240 if(kExcludeParticle) continue;
1243 vCharge = track->Charge();
1244 vPhi = track->Phi();
1247 //Printf("phi (before): %lf",vPhi);
1250 fHistEta->Fill(vEta);
1251 fHistPhi->Fill(vPhi);
1252 fHistRapidity->Fill(vY);
1253 if(vCharge > 0) fHistPhiPos->Fill(vPhi);
1254 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
1257 if(fUseFlowAfterBurner) {
1258 Double_t precisionPhi = 0.001;
1259 Int_t maxNumberOfIterations = 100;
1261 Double_t phi0 = vPhi;
1262 Double_t gV2 = fDifferentialV2->Eval(vPt);
1264 for (Int_t j = 0; j < maxNumberOfIterations; j++) {
1265 Double_t phiprev = vPhi;
1266 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane));
1267 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane));
1269 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
1271 //Printf("phi (after): %lf\n",vPhi);
1272 Double_t vDeltaphiBefore = phi0 - gReactionPlane;
1273 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
1274 fHistPhiBefore->Fill(vDeltaphiBefore);
1276 Double_t vDeltaphiAfter = vPhi - gReactionPlane;
1277 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
1278 fHistPhiAfter->Fill(vDeltaphiAfter);
1281 vPhi *= TMath::RadToDeg();
1283 // fill charge vector
1284 chargeVector[0]->push_back(vCharge);
1285 chargeVector[1]->push_back(vY);
1286 chargeVector[2]->push_back(vEta);
1287 chargeVector[3]->push_back(vPhi);
1288 chargeVector[4]->push_back(vP[0]);
1289 chargeVector[5]->push_back(vP[1]);
1290 chargeVector[6]->push_back(vP[2]);
1291 chargeVector[7]->push_back(vPt);
1292 chargeVector[8]->push_back(vE);
1295 chargeVectorShuffle[0]->push_back(vCharge);
1296 chargeVectorShuffle[1]->push_back(vY);
1297 chargeVectorShuffle[2]->push_back(vEta);
1298 chargeVectorShuffle[3]->push_back(vPhi);
1299 chargeVectorShuffle[4]->push_back(vP[0]);
1300 chargeVectorShuffle[5]->push_back(vP[1]);
1301 chargeVectorShuffle[6]->push_back(vP[2]);
1302 chargeVectorShuffle[7]->push_back(vPt);
1303 chargeVectorShuffle[8]->push_back(vE);
1305 gNumberOfAcceptedTracks += 1;
1313 //multiplicity cut (used in pp)
1314 if(fUseMultiplicity) {
1315 if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))
1318 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks, fCentrality);
1320 // calculate balance function
1321 if(fUseMultiplicity)
1322 fBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVector,bSign);
1324 fBalance->CalculateBalance(fCentrality,chargeVector,bSign);
1328 random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );
1329 if(fUseMultiplicity)
1330 fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVectorShuffle,bSign);
1332 fShuffledBalance->CalculateBalance(fCentrality,chargeVectorShuffle,bSign);
1336 //________________________________________________________________________
1337 void AliAnalysisTaskBF::FinishTaskOutput(){
1341 AliError("ERROR: fBalance not available");
1345 if (!fShuffledBalance) {
1346 AliError("ERROR: fShuffledBalance not available");
1353 //________________________________________________________________________
1354 void AliAnalysisTaskBF::Terminate(Option_t *) {
1355 // Draw result to the screen
1356 // Called once at the end of the query
1358 // not implemented ...