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->GetXv()) < fVxMax) {
467 if(TMath::Abs(vertex->GetYv()) < fVyMax) {
468 if(TMath::Abs(vertex->GetZv()) < fVzMax) {
469 fHistEventStats->Fill(4,fCentrality); //analayzed events
470 fHistVx->Fill(vertex->GetXv());
471 fHistVy->Fill(vertex->GetYv());
472 fHistVz->Fill(vertex->GetZv());
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 = gAOD->GetHeader();
675 // store offline trigger bits
676 fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());
679 fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
682 //event selection done in AliAnalysisTaskSE::Exec() --> this is not used
683 fHistEventStats->Fill(1,fCentrality); //all events
684 Bool_t isSelected = kTRUE;
685 if(fUseOfflineTrigger)
686 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
688 fHistEventStats->Fill(2,fCentrality); //triggered events
690 //Centrality stuff (centrality in AOD header)
692 //fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
693 // in OLD AODs (i.e. AOD049) fCentrality can be == 0
697 // QA for centrality estimators
698 fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));
699 fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));
700 fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));
701 fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));
702 fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));
703 fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));
704 fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
705 fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
706 fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
708 // take only events inside centrality class
709 if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax))
712 // centrality QA (V0M)
713 fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());
715 // centrality QA (reference tracks)
716 fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());
717 fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());
718 fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());
719 fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());
720 fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));
721 fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));
722 fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));
723 fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));
724 fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));
727 const AliAODVertex *vertex = gAOD->GetPrimaryVertex();
731 vertex->GetCovarianceMatrix(fCov);
733 if(vertex->GetNContributors() > 0) {
735 fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
736 if(TMath::Abs(vertex->GetX()) < fVxMax) {
737 if(TMath::Abs(vertex->GetY()) < fVyMax) {
738 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
739 fHistEventStats->Fill(4,fCentrality); //analyzed events
740 fHistVx->Fill(vertex->GetX());
741 fHistVy->Fill(vertex->GetY());
742 fHistVz->Fill(vertex->GetZ());
744 //===========================================//
745 TExMap *trackMap = new TExMap();
746 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {
747 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));
749 AliError(Form("ERROR: Could not receive track %d", iTracks));
752 Int_t gID = aodTrack->GetID();
753 //if (!aodTrack->TestFilterBit(fAODtrackCutBit)) trackMap->Add(gID, iTracks);
754 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, iTracks);
756 AliAODTrack* newAodTrack;
757 //===========================================//
759 //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());
760 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {
761 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));
763 AliError(Form("ERROR: Could not receive track %d", iTracks));
768 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
769 //===========================================//
770 // take only TPC only tracks
771 fHistTrackStats->Fill(aodTrack->GetFilterMap());
772 if(!aodTrack->TestFilterBit(fAODtrackCutBit)) continue;
774 Int_t gID = aodTrack->GetID();
775 newAodTrack = gID >= 0 ? aodTrack : gAOD->GetTrack(trackMap->GetValue(-1-gID));
776 //Printf("Label: %d - Pt: %lf (old) - %d - Pt: %lf(new)",gID,aodTrack->Pt(), newAodTrack->GetID(), newAodTrack->Pt());
777 //===========================================//
779 //fHistTrackStats->Fill(aodTrack->GetFilterMap());
780 //if(!aodTrack->TestFilterBit(fAODtrackCutBit)) continue;
782 vCharge = aodTrack->Charge();
784 vEta = aodTrack->Eta();
785 vPhi = aodTrack->Phi() * TMath::RadToDeg();
787 vPt = aodTrack->Pt();
788 aodTrack->PxPyPz(vP);
790 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
791 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
794 // Kinematics cuts from ESD track cuts
795 if( vPt < fPtMin || vPt > fPtMax) continue;
798 if( vEta < fEtaMin || vEta > fEtaMax) continue;
802 if( vY < fEtaMin || vY > fEtaMax) continue;
805 // Extra DCA cuts (for systematic studies [!= -1])
806 if( fDCAxyCut != -1 && fDCAzCut != -1){
807 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
812 // Extra TPC cuts (for systematic studies [!= -1])
813 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
816 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
820 //===============================================PID==================================//
822 Double_t prob[AliPID::kSPECIES]={0.};
823 Double_t probTPC[AliPID::kSPECIES]={0.};
824 Double_t probTOF[AliPID::kSPECIES]={0.};
825 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
827 Double_t nSigma = 0.;
828 UInt_t detUsedTPC = 0;
829 UInt_t detUsedTOF = 0;
830 UInt_t detUsedTPCTOF = 0;
832 //Decide what detector configuration we want to use
833 switch(fPidDetectorConfig) {
835 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
836 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest));
837 //detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
838 detUsedTPC = (AliPIDResponse::kDetTPC);
839 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
840 prob[iSpecies] = probTPC[iSpecies];
843 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
844 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest));
845 //detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
846 detUsedTPC = (AliPIDResponse::kDetTPC);
847 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
848 prob[iSpecies] = probTOF[iSpecies];
851 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
852 //detUsedTPCTOF = fPIDCombined->ComputeProbabilities(newAodTrack, fPIDResponse, probTPCTOF);
853 detUsedTPC = (AliPIDResponse::kDetTPC);
854 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
855 prob[iSpecies] = probTPCTOF[iSpecies];
859 }//end switch: define detector mask
862 Double_t tofTime = -999., tof = -999.; //length = 999., tof = -999.;
863 //Double_t c = TMath::C()*1.E-9;// m/ns
864 //Double_t beta = -999.;
865 Double_t nSigmaTOFForParticleOfInterest = -999.;
866 if ( (newAodTrack->IsOn(AliAODTrack::kTOFin)) &&
867 (newAodTrack->IsOn(AliAODTrack::kTIME)) ) {
868 tofTime = newAodTrack->GetTOFsignal();//in ps
869 //length = newAodTrack->GetIntegratedLength();
870 tof = tofTime*1E-3; // ns
873 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
877 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
881 //length = length*0.01; // in meters
885 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest);
886 //fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
887 fHistProbTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]);
888 fHistNSigmaTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest);
892 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest);
893 fHistdEdxVsPTPCbeforePID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal());
894 fHistProbTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]);
895 fHistNSigmaTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest);
896 fHistProbTPCTOFvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
897 //end of QA-before pid
899 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
900 //Make the decision based on the n-sigma
902 if(nSigma > fPIDNSigma) continue;}
904 //Make the decision based on the bayesian
905 else if(fUsePIDPropabilities) {
906 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
907 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
910 //Fill QA after the PID
911 //fHistBetavsPTOFafterPID ->Fill(newAodTrack->P()*newAodTrack->Charge(),beta);
912 fHistProbTOFvsPtafterPID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]);
913 fHistNSigmaTOFvsPtafterPID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest);
915 fHistdEdxVsPTPCafterPID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal());
916 fHistProbTPCvsPtafterPID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]);
917 fHistProbTPCTOFvsPtafterPID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
918 fHistNSigmaTPCvsPtafterPID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest);
921 PostData(4, fHistListPIDQA);
924 //=========================================================PID=================================================================//
926 // fill QA histograms
927 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
928 fHistDCA->Fill(dcaZ,dcaXY);
929 fHistChi2->Fill(aodTrack->Chi2perNDF());
931 fHistEta->Fill(vEta);
932 fHistPhi->Fill(vPhi);
933 fHistRapidity->Fill(vY);
934 if(vCharge > 0) fHistPhiPos->Fill(vPhi);
935 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
937 // fill charge vector
938 chargeVector[0]->push_back(vCharge);
939 chargeVector[1]->push_back(vY);
940 chargeVector[2]->push_back(vEta);
941 chargeVector[3]->push_back(vPhi);
942 chargeVector[4]->push_back(vP[0]);
943 chargeVector[5]->push_back(vP[1]);
944 chargeVector[6]->push_back(vP[2]);
945 chargeVector[7]->push_back(vPt);
946 chargeVector[8]->push_back(vE);
949 chargeVectorShuffle[0]->push_back(vCharge);
950 chargeVectorShuffle[1]->push_back(vY);
951 chargeVectorShuffle[2]->push_back(vEta);
952 chargeVectorShuffle[3]->push_back(vPhi);
953 chargeVectorShuffle[4]->push_back(vP[0]);
954 chargeVectorShuffle[5]->push_back(vP[1]);
955 chargeVectorShuffle[6]->push_back(vP[2]);
956 chargeVectorShuffle[7]->push_back(vPt);
957 chargeVectorShuffle[8]->push_back(vE);
960 gNumberOfAcceptedTracks += 1;
966 }//proper vertex resolution
967 }//proper number of contributors
968 }//vertex object valid
973 if(gAnalysisLevel == "MCESD") {
974 AliMCEvent* mcEvent = MCEvent();
976 AliError("ERROR: mcEvent not available");
980 AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
982 AliError("ERROR: gESD not available");
986 // store offline trigger bits
987 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
989 AliCentrality *centrality = 0x0;
991 centrality = gESD->GetCentrality();
992 fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
995 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
996 fHistEventStats->Fill(1,fCentrality); //all events
997 Bool_t isSelected = kTRUE;
998 if(fUseOfflineTrigger)
999 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1001 fHistEventStats->Fill(2,fCentrality); //triggered events
1003 if(fUseCentrality) {
1005 // take only events inside centrality class
1006 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
1007 fCentralityPercentileMax,
1008 fCentralityEstimator.Data()))
1010 // centrality QA (V0M)
1011 fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());
1014 const AliESDVertex *vertex = gESD->GetPrimaryVertex();
1016 if(vertex->GetNContributors() > 0) {
1017 if(vertex->GetZRes() != 0) {
1018 fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
1019 if(TMath::Abs(vertex->GetXv()) < fVxMax) {
1020 if(TMath::Abs(vertex->GetYv()) < fVyMax) {
1021 if(TMath::Abs(vertex->GetZv()) < fVzMax) {
1022 fHistEventStats->Fill(4,fCentrality); //analayzed events
1023 fHistVx->Fill(vertex->GetXv());
1024 fHistVy->Fill(vertex->GetYv());
1025 fHistVz->Fill(vertex->GetZv());
1027 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());
1028 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {
1029 AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));
1031 AliError(Form("ERROR: Could not receive track %d", iTracks));
1035 Int_t label = TMath::Abs(track->GetLabel());
1036 if(label > mcEvent->GetNumberOfTracks()) continue;
1037 if(label > mcEvent->GetNumberOfPrimaries()) continue;
1039 AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
1040 if(!mcTrack) continue;
1042 // take only TPC only tracks
1043 trackTPC = new AliESDtrack();
1044 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
1048 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
1050 // fill QA histograms
1053 trackTPC->GetImpactParameters(b,bCov);
1054 if (bCov[0]<=0 || bCov[2]<=0) {
1055 AliDebug(1, "Estimated b resolution lower or equal zero!");
1056 bCov[0]=0; bCov[2]=0;
1059 Int_t nClustersTPC = -1;
1060 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
1061 //nClustersTPC = track->GetTPCclusters(0); // global track
1062 Float_t chi2PerClusterTPC = -1;
1063 if (nClustersTPC!=0) {
1064 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
1065 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
1068 vCharge = trackTPC->Charge();
1070 vEta = trackTPC->Eta();
1071 vPhi = trackTPC->Phi() * TMath::RadToDeg();
1073 vPt = trackTPC->Pt();
1074 trackTPC->PxPyPz(vP);
1076 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
1077 fHistDCA->Fill(b[1],b[0]);
1078 fHistChi2->Fill(chi2PerClusterTPC);
1080 fHistEta->Fill(vEta);
1081 fHistPhi->Fill(vPhi);
1082 fHistRapidity->Fill(vY);
1083 if(vCharge > 0) fHistPhiPos->Fill(vPhi);
1084 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
1086 // fill charge vector
1087 chargeVector[0]->push_back(vCharge);
1088 chargeVector[1]->push_back(vY);
1089 chargeVector[2]->push_back(vEta);
1090 chargeVector[3]->push_back(vPhi);
1091 chargeVector[4]->push_back(vP[0]);
1092 chargeVector[5]->push_back(vP[1]);
1093 chargeVector[6]->push_back(vP[2]);
1094 chargeVector[7]->push_back(vPt);
1095 chargeVector[8]->push_back(vE);
1098 chargeVectorShuffle[0]->push_back(vCharge);
1099 chargeVectorShuffle[1]->push_back(vY);
1100 chargeVectorShuffle[2]->push_back(vEta);
1101 chargeVectorShuffle[3]->push_back(vPhi);
1102 chargeVectorShuffle[4]->push_back(vP[0]);
1103 chargeVectorShuffle[5]->push_back(vP[1]);
1104 chargeVectorShuffle[6]->push_back(vP[2]);
1105 chargeVectorShuffle[7]->push_back(vPt);
1106 chargeVectorShuffle[8]->push_back(vE);
1110 gNumberOfAcceptedTracks += 1;
1113 //cout<<"Centrality: "<<fCentrality<<" - Accepted tracks: "<<gNumberOfAcceptedTracks<<endl;
1117 }//proper vertex resolution
1118 }//proper number of contributors
1119 }//vertex object valid
1124 else if(gAnalysisLevel == "MC") {
1125 AliMCEvent* mcEvent = MCEvent();
1127 AliError("ERROR: mcEvent not available");
1131 //fHistEventStats->Fill(1,fCentrality); //total events
1132 //fHistEventStats->Fill(2,fCentrality); //offline trigger
1134 Double_t gReactionPlane = 0., gImpactParameter = 0.;
1135 if(fUseCentrality) {
1137 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
1139 //Printf("=====================================================");
1140 //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());
1141 //Printf("Impact parameter: %lf",headerH->ImpactParameter());
1142 //Printf("=====================================================");
1143 gReactionPlane = headerH->ReactionPlaneAngle();
1144 gImpactParameter = headerH->ImpactParameter();
1145 fCentrality = gImpactParameter;
1147 fCentrality = gImpactParameter;
1149 // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW)
1150 if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))
1154 fHistEventStats->Fill(1,fCentrality); //total events
1155 fHistEventStats->Fill(2,fCentrality); //offline trigger
1157 AliGenEventHeader *header = mcEvent->GenEventHeader();
1160 TArrayF gVertexArray;
1161 header->PrimaryVertex(gVertexArray);
1162 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
1163 //gVertexArray.At(0),
1164 //gVertexArray.At(1),
1165 //gVertexArray.At(2));
1166 fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
1167 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
1168 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
1169 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
1170 fHistEventStats->Fill(4,fCentrality); //analayzed events
1171 fHistVx->Fill(gVertexArray.At(0));
1172 fHistVy->Fill(gVertexArray.At(1));
1173 fHistVz->Fill(gVertexArray.At(2));
1175 AliInfo(Form("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries()));
1176 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {
1177 AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));
1179 AliError(Form("ERROR: Could not receive particle %d", iTracks));
1183 //exclude non stable particles
1184 if(!mcEvent->IsPhysicalPrimary(iTracks)) continue;
1186 vEta = track->Eta();
1190 if( vPt < fPtMin || vPt > fPtMax)
1193 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1196 if( vY < fEtaMin || vY > fEtaMax) continue;
1199 //analyze one set of particles
1201 TParticle *particle = track->Particle();
1202 if(!particle) continue;
1204 Int_t gPdgCode = particle->GetPdgCode();
1205 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
1209 //Use the acceptance parameterization
1210 if(fAcceptanceParameterization) {
1211 Double_t gRandomNumber = gRandom->Rndm();
1212 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
1216 //Exclude resonances
1217 if(fExcludeResonancesInMC) {
1218 TParticle *particle = track->Particle();
1219 if(!particle) continue;
1221 Bool_t kExcludeParticle = kFALSE;
1222 Int_t gMotherIndex = particle->GetFirstMother();
1223 if(gMotherIndex != -1) {
1224 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1226 TParticle *motherParticle = motherTrack->Particle();
1227 if(motherParticle) {
1228 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
1229 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1230 if(pdgCodeOfMother == 113) {
1231 kExcludeParticle = kTRUE;
1237 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1238 if(kExcludeParticle) continue;
1241 vCharge = track->Charge();
1242 vPhi = track->Phi();
1245 //Printf("phi (before): %lf",vPhi);
1248 fHistEta->Fill(vEta);
1249 fHistPhi->Fill(vPhi);
1250 fHistRapidity->Fill(vY);
1251 if(vCharge > 0) fHistPhiPos->Fill(vPhi);
1252 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
1255 if(fUseFlowAfterBurner) {
1256 Double_t precisionPhi = 0.001;
1257 Int_t maxNumberOfIterations = 100;
1259 Double_t phi0 = vPhi;
1260 Double_t gV2 = fDifferentialV2->Eval(vPt);
1262 for (Int_t j = 0; j < maxNumberOfIterations; j++) {
1263 Double_t phiprev = vPhi;
1264 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane));
1265 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane));
1267 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
1269 //Printf("phi (after): %lf\n",vPhi);
1270 Double_t vDeltaphiBefore = phi0 - gReactionPlane;
1271 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
1272 fHistPhiBefore->Fill(vDeltaphiBefore);
1274 Double_t vDeltaphiAfter = vPhi - gReactionPlane;
1275 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
1276 fHistPhiAfter->Fill(vDeltaphiAfter);
1279 vPhi *= TMath::RadToDeg();
1281 // fill charge vector
1282 chargeVector[0]->push_back(vCharge);
1283 chargeVector[1]->push_back(vY);
1284 chargeVector[2]->push_back(vEta);
1285 chargeVector[3]->push_back(vPhi);
1286 chargeVector[4]->push_back(vP[0]);
1287 chargeVector[5]->push_back(vP[1]);
1288 chargeVector[6]->push_back(vP[2]);
1289 chargeVector[7]->push_back(vPt);
1290 chargeVector[8]->push_back(vE);
1293 chargeVectorShuffle[0]->push_back(vCharge);
1294 chargeVectorShuffle[1]->push_back(vY);
1295 chargeVectorShuffle[2]->push_back(vEta);
1296 chargeVectorShuffle[3]->push_back(vPhi);
1297 chargeVectorShuffle[4]->push_back(vP[0]);
1298 chargeVectorShuffle[5]->push_back(vP[1]);
1299 chargeVectorShuffle[6]->push_back(vP[2]);
1300 chargeVectorShuffle[7]->push_back(vPt);
1301 chargeVectorShuffle[8]->push_back(vE);
1303 gNumberOfAcceptedTracks += 1;
1311 //multiplicity cut (used in pp)
1312 if(fUseMultiplicity) {
1313 if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))
1316 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks, fCentrality);
1318 // calculate balance function
1319 if(fUseMultiplicity)
1320 fBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVector,bSign);
1322 fBalance->CalculateBalance(fCentrality,chargeVector,bSign);
1326 random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );
1327 if(fUseMultiplicity)
1328 fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVectorShuffle,bSign);
1330 fShuffledBalance->CalculateBalance(fCentrality,chargeVectorShuffle,bSign);
1334 //________________________________________________________________________
1335 void AliAnalysisTaskBF::FinishTaskOutput(){
1339 AliError("ERROR: fBalance not available");
1343 if (!fShuffledBalance) {
1344 AliError("ERROR: fShuffledBalance not available");
1351 //________________________________________________________________________
1352 void AliAnalysisTaskBF::Terminate(Option_t *) {
1353 // Draw result to the screen
1354 // Called once at the end of the query
1356 // not implemented ...