]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskEventMixingBF.cxx
CommitLineData
c683985a 1#include "TChain.h"
2#include "TList.h"
3#include "TCanvas.h"
4#include "TLorentzVector.h"
5#include "TGraphErrors.h"
6#include "TH1F.h"
7#include "TH2F.h"
8#include "TArrayF.h"
9#include "TF1.h"
10#include "TRandom.h"
11
12#include "AliAnalysisTaskSE.h"
13#include "AliAnalysisManager.h"
14
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"
25#include "AliMixInputEventHandler.h"
26#include "AliStack.h"
27#include "AliESDtrackCuts.h"
28
29#include "TH2D.h"
30#include "AliPID.h"
31#include "AliPIDResponse.h"
32#include "AliPIDCombined.h"
33
34#include "AliAnalysisTaskEventMixingBF.h"
35#include "AliBalanceEventMixing.h"
36
37
38// Analysis task for the EventMixingBF code
39// Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch
40
41ClassImp(AliAnalysisTaskEventMixingBF)
42
43//________________________________________________________________________
44AliAnalysisTaskEventMixingBF::AliAnalysisTaskEventMixingBF(const char *name)
45: AliAnalysisTaskSE(name),
46 fBalance(0),
47 fRunShuffling(kFALSE),
48 fShuffledBalance(0),
49 fList(0),
50 fListEventMixingBF(0),
51 fListEventMixingBFS(0),
52 fHistListPIDQA(0),
53 fHistEventStats(0),
54 fHistCentStats(0),
55 fHistTriggerStats(0),
56 fHistTrackStats(0),
57 fHistVx(0),
58 fHistVy(0),
59 fHistVz(0),
60 fHistClus(0),
61 fHistDCA(0),
62 fHistChi2(0),
63 fHistPt(0),
64 fHistEta(0),
65 fHistPhi(0),
66 fHistPhiBefore(0),
67 fHistPhiAfter(0),
68 fHistV0M(0),
69 fHistRefTracks(0),
70 fHistdEdxVsPTPCbeforePID(NULL),
71 fHistBetavsPTOFbeforePID(NULL),
72 fHistProbTPCvsPtbeforePID(NULL),
73 fHistProbTOFvsPtbeforePID(NULL),
74 fHistProbTPCTOFvsPtbeforePID(NULL),
75 fHistNSigmaTPCvsPtbeforePID(NULL),
76 fHistNSigmaTOFvsPtbeforePID(NULL),
77 fHistdEdxVsPTPCafterPID(NULL),
78 fHistBetavsPTOFafterPID(NULL),
79 fHistProbTPCvsPtafterPID(NULL),
80 fHistProbTOFvsPtafterPID(NULL),
81 fHistProbTPCTOFvsPtafterPID(NULL),
82 fHistNSigmaTPCvsPtafterPID(NULL),
83 fHistNSigmaTOFvsPtafterPID(NULL),
84 fPIDResponse(0x0),
85 fPIDCombined(0x0),
86 fParticleOfInterest(kPion),
87 fPidDetectorConfig(kTPCTOF),
88 fUsePID(kFALSE),
89 fUsePIDnSigma(kTRUE),
90 fUsePIDPropabilities(kFALSE),
91 fPIDNSigma(3.),
92 fMinAcceptedPIDProbability(0.8),
93 fESDtrackCuts(0),
94 fCentralityEstimator("V0M"),
95 fUseCentrality(kFALSE),
96 fCentralityPercentileMin(0.),
97 fCentralityPercentileMax(5.),
98 fImpactParameterMin(0.),
99 fImpactParameterMax(20.),
100 fUseMultiplicity(kFALSE),
101 fNumberOfAcceptedTracksMin(0),
102 fNumberOfAcceptedTracksMax(10000),
103 fHistNumberOfAcceptedTracks(0),
104 fUseOfflineTrigger(kFALSE),
105 fVxMax(0.3),
106 fVyMax(0.3),
107 fVzMax(10.),
108 nAODtrackCutBit(128),
109 fPtMin(0.3),
110 fPtMax(1.5),
111 fEtaMin(-0.8),
112 fEtaMax(-0.8),
113 fDCAxyCut(-1),
114 fDCAzCut(-1),
115 fTPCchi2Cut(-1),
116 fNClustersTPCCut(-1),
117 fAcceptanceParameterization(0),
118 fDifferentialV2(0),
119 fUseFlowAfterBurner(kFALSE),
120 fExcludeResonancesInMC(kFALSE),
121 fUseMCPdgCode(kFALSE),
122 fPDGCodeToBeAnalyzed(-1),
123 fMainEvent(0x0),
124 fMixEvent(0x0)
125{
126 // Constructor
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());
135}
136
137//________________________________________________________________________
138AliAnalysisTaskEventMixingBF::~AliAnalysisTaskEventMixingBF() {
139
140 // delete fBalance;
141 // delete fShuffledBalance;
142 // delete fList;
143 // delete fListEventMixingBF;
144 // delete fListEventMixingBFS;
145
146 // delete fHistEventStats;
147 // delete fHistTrackStats;
148 // delete fHistVx;
149 // delete fHistVy;
150 // delete fHistVz;
151
152 // delete fHistClus;
153 // delete fHistDCA;
154 // delete fHistChi2;
155 // delete fHistPt;
156 // delete fHistEta;
157 // delete fHistPhi;
158 // delete fHistV0M;
159}
160
161//________________________________________________________________________
162void AliAnalysisTaskEventMixingBF::UserCreateOutputObjects() {
163 // Create histograms
164 // Called once
165 if(!fBalance) {
166 fBalance = new AliBalanceEventMixing();
167 fBalance->SetAnalysisLevel("ESD");
168 //fBalance->SetNumberOfBins(-1,16);
169 fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
170 }
171 if(fRunShuffling) {
172 if(!fShuffledBalance) {
173 fShuffledBalance = new AliBalanceEventMixing();
174 fShuffledBalance->SetAnalysisLevel("ESD");
175 //fShuffledBalance->SetNumberOfBins(-1,16);
176 fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
177 }
178 }
179
180 //QA list
181 fList = new TList();
182 fList->SetName("listQA");
183 fList->SetOwner();
184
185 //Balance Function list
186 fListEventMixingBF = new TList();
187 fListEventMixingBF->SetName("listEventMixingBF");
188 fListEventMixingBF->SetOwner();
189
190 if(fRunShuffling) {
191 fListEventMixingBFS = new TList();
192 fListEventMixingBFS->SetName("listEventMixingBFShuffled");
193 fListEventMixingBFS->SetOwner();
194 }
195
196 //PID QA list
197 if(fUsePID) {
198 fHistListPIDQA = new TList();
199 fHistListPIDQA->SetName("listQAPID");
200 fHistListPIDQA->SetOwner();
201 }
202
203 //Event stats.
204 TString gCutName[4] = {"Total","Offline trigger",
205 "Vertex","Analyzed"};
206 fHistEventStats = new TH1F("fHistEventStats",
207 "Event statistics;;N_{events}",
208 4,0.5,4.5);
209 for(Int_t i = 1; i <= 4; i++)
210 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
211 fList->Add(fHistEventStats);
212
213 TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
214 fHistCentStats = new TH2F("fHistCentStats",
215 "Centrality statistics;;Cent percentile",
216 9,-0.5,8.5,220,-5,105);
217 for(Int_t i = 1; i <= 9; i++)
218 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
219 fList->Add(fHistCentStats);
220
221 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);
222 fList->Add(fHistTriggerStats);
223
224 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);
225 fList->Add(fHistTrackStats);
226
227 fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);
228 fList->Add(fHistNumberOfAcceptedTracks);
229
230 // Vertex distributions
231 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
232 fList->Add(fHistVx);
233 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
234 fList->Add(fHistVy);
235 fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
236 fList->Add(fHistVz);
237
238 // QA histograms
239 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
240 fList->Add(fHistClus);
241 fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);
242 fList->Add(fHistChi2);
243 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
244 fList->Add(fHistDCA);
245 fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10);
246 fList->Add(fHistPt);
247 fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2);
248 fList->Add(fHistEta);
249 fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380);
250 fList->Add(fHistPhi);
251 fHistPhiBefore = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());
252 fList->Add(fHistPhiBefore);
253 fHistPhiAfter = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());
254 fList->Add(fHistPhiAfter);
255 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
256 fList->Add(fHistV0M);
257 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
258 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
259 for(Int_t i = 1; i <= 6; i++)
260 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
261 fList->Add(fHistRefTracks);
262
263 // Balance function histograms
264 // Initialize histograms if not done yet
265 if(!fBalance->GetHistNp(0)){
266 AliWarning("Histograms not yet initialized! --> Will be done now");
267 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
268 fBalance->InitHistograms();
269 }
270
271 if(fRunShuffling) {
272 if(!fShuffledBalance->GetHistNp(0)) {
273 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
274 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
275 fShuffledBalance->InitHistograms();
276 }
277 }
278
279 for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
280 fListEventMixingBF->Add(fBalance->GetHistNp(a));
281 fListEventMixingBF->Add(fBalance->GetHistNn(a));
282 fListEventMixingBF->Add(fBalance->GetHistNpn(a));
283 fListEventMixingBF->Add(fBalance->GetHistNnn(a));
284 fListEventMixingBF->Add(fBalance->GetHistNpp(a));
285 fListEventMixingBF->Add(fBalance->GetHistNnp(a));
286
287 if(fRunShuffling) {
288 fListEventMixingBFS->Add(fShuffledBalance->GetHistNp(a));
289 fListEventMixingBFS->Add(fShuffledBalance->GetHistNn(a));
290 fListEventMixingBFS->Add(fShuffledBalance->GetHistNpn(a));
291 fListEventMixingBFS->Add(fShuffledBalance->GetHistNnn(a));
292 fListEventMixingBFS->Add(fShuffledBalance->GetHistNpp(a));
293 fListEventMixingBFS->Add(fShuffledBalance->GetHistNnp(a));
294 }
295 }
296
297 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
298
299 //====================PID========================//
300 if(fUsePID) {
301 fPIDCombined = new AliPIDCombined();
302 fPIDCombined->SetDefaultTPCPriors();
303
304 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
305 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition
306
307 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
308 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition
309
310 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
311 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition
312
313 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
314 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition
315
316 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
317 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition
318
319 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
320 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition
321
322 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
323 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition
324
325 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
326 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition
327
328 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
329 fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition
330
331 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
332 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition
333
334 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
335 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition
336
337 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
338 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition
339
340 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
341 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition
342
343 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
344 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition
345 }
346 //====================PID========================//
347
348 // Post output data.
349 PostData(1, fList);
350 PostData(2, fListEventMixingBF);
351 if(fRunShuffling) PostData(3, fListEventMixingBFS);
352 if(fUsePID) PostData(4, fHistListPIDQA); //PID
353}
354
355//________________________________________________________________________
356void AliAnalysisTaskEventMixingBF::UserExec(Option_t *) {
357 // Main loop
358 // Called for each event
359 // NOTHING TO DO for event mixing!
360}
361
362//________________________________________________________________________
363void AliAnalysisTaskEventMixingBF::FinishTaskOutput(){
364 //Printf("END EventMixingBF");
365
366 if (!fBalance) {
367 AliError("ERROR: fBalance not available");
368 return;
369 }
370 if(fRunShuffling) {
371 if (!fShuffledBalance) {
372 AliError("ERROR: fShuffledBalance not available");
373 return;
374 }
375 }
376
377}
378
379//________________________________________________________________________
380void AliAnalysisTaskEventMixingBF::Terminate(Option_t *) {
381 // Draw result to the screen
382 // Called once at the end of the query
383
384 // not implemented ...
385
386}
387
388void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
389{
390 // Main loop for event mixing
391
392 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
393
394 AliMixInputEventHandler *mixIEH = SetupEventsForMixing();
395
396 Float_t fCentrality = 0.;
397
398 // for HBT like cuts need magnetic field sign
399 Float_t bSign = 0; // only used in AOD so far
400
401 // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)
402 vector<Double_t> *chargeVector[9]; // original charge
403 for(Int_t i = 0; i < 9; i++){
404 chargeVector[i] = new vector<Double_t>;
405 }
406
407 Double_t vCharge;
408 Double_t vY;
409 Double_t vEta;
410 Double_t vPhi;
411 Double_t vP[3];
412 Double_t vPt;
413 Double_t vE;
414
415 Int_t iMainTrackUsed = -1;
416
417 // -------------------------------------------------------------
418 // At the moment MIXING only for AODs
419 if(mixIEH){
420
421 //AOD analysis (vertex and track cuts also here!!!!)
422 if(gAnalysisLevel == "AOD") {
423 AliAODEvent* aodEventMain = dynamic_cast<AliAODEvent*>(fMainEvent);
424 if(!aodEventMain) {
425 AliError("ERROR: aodEventMain not available");
426 return;
427 }
428 AliAODEvent *aodEventMix = dynamic_cast<AliAODEvent *>(fMixEvent);
429 if(!aodEventMix) {
430 AliError("ERROR: aodEventMix not available");
431 return;
432 }
433
434 // for HBT like cuts need magnetic field sign
435 bSign = (aodEventMain->GetMagneticField() > 0) ? 1 : -1;
436
0a918d8d 437 AliAODHeader *aodHeaderMain = dynamic_cast<AliAODHeader*>(aodEventMain->GetHeader());
438 if(!aodHeaderMain) AliFatal("Not a standard AOD");
c683985a 439
440 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
441 fHistEventStats->Fill(1); //all events
442
443 // this is not needed (checked in mixing handler!)
444 Bool_t isSelectedMain = kTRUE;
445 Bool_t isSelectedMix = kTRUE;
446
447 if(fUseOfflineTrigger){
448 isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
449 isSelectedMix = ((AliInputEventHandler*)((AliMultiInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetFirstMultiInputHandler())->IsEventSelected();
450 }
451
452 if(isSelectedMain && isSelectedMix) {
453 fHistEventStats->Fill(2); //triggered events
454
455 //Centrality stuff (centrality in AOD header)
456 if(fUseCentrality) {
457 fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
458
459 // QA for centrality estimators
460 fHistCentStats->Fill(0.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));
461 fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));
462 fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));
463 fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));
464 fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));
465 fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));
466 fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
467 fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
468 fHistCentStats->Fill(8.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
469
470 // take only events inside centrality class
471 if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax))
472 return;
473
474 // centrality QA (V0M)
475 fHistV0M->Fill(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());
476
477 // centrality QA (reference tracks)
478 fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());
479 fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());
480 fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());
481 fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());
482 fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));
483 fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));
484 fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));
485 fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));
486 fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));
487 }
488
489 // // this is crashing (Bug in ROOT (to be solved)) but not needed (checked in mixing handler)
490 // const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();
491 // const AliAODVertex *vertexMix = aodEventMix->GetPrimaryVertex();
492
493 // if(vertexMain && vertexMix) {
494 // Double32_t fCovMain[6];
495 // Double32_t fCovMix[6];
496 // vertexMain->GetCovarianceMatrix(fCovMain);
497 // vertexMix->GetCovarianceMatrix(fCovMix);
498
499 // if(vertexMain->GetNContributors() > 0 && vertexMix->GetNContributors() > 0) {
500 // if(fCovMain[5] != 0 && fCovMix[5] != 0) {
501 // fHistEventStats->Fill(3); //events with a proper vertex
502 // if(TMath::Abs(vertexMain->GetX()) < fVxMax && TMath::Abs(vertexMix->GetX()) < fVxMax ) {
503 // if(TMath::Abs(vertexMain->GetY()) < fVyMax && TMath::Abs(vertexMix->GetY()) < fVyMax) {
504 // if(TMath::Abs(vertexMain->GetZ()) < fVzMax && TMath::Abs(vertexMix->GetZ()) < fVzMax) {
505 // fHistEventStats->Fill(4); //analyzed events
506 // fHistVx->Fill(vertexMain->GetX());
507 // fHistVy->Fill(vertexMain->GetY());
508 // fHistVz->Fill(vertexMain->GetZ());
509
510 // Loop over tracks in main event
511 for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {
512 AliAODTrack* aodTrackMain = dynamic_cast<AliAODTrack *>(aodEventMain->GetTrack(iTracksMain));
513 if (!aodTrackMain) {
514 AliError(Form("ERROR: Could not receive track %d", iTracksMain));
515 continue;
516 }
517
518 // AOD track cuts
519
520 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
521 // take only TPC only tracks
522 fHistTrackStats->Fill(aodTrackMain->GetFilterMap());
523 if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;
524
525 vCharge = aodTrackMain->Charge();
526 vY = aodTrackMain->Y();
527 vEta = aodTrackMain->Eta();
528 vPhi = aodTrackMain->Phi() * TMath::RadToDeg();
529 vE = aodTrackMain->E();
530 vPt = aodTrackMain->Pt();
531 aodTrackMain->PxPyPz(vP);
532
533 Float_t dcaXYMain = aodTrackMain->DCA(); // this is the DCA from global track (not exactly what is cut on)
534 Float_t dcaZMain = aodTrackMain->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
535
536
537 // Kinematics cuts from ESD track cuts
538 if( vPt < fPtMin || vPt > fPtMax) continue;
539 if( vEta < fEtaMin || vEta > fEtaMax) continue;
540
541 // Extra DCA cuts (for systematic studies [!= -1])
542 if( fDCAxyCut != -1 && fDCAzCut != -1){
543 if(TMath::Sqrt((dcaXYMain*dcaXYMain)/(fDCAxyCut*fDCAxyCut)+(dcaZMain*dcaZMain)/(fDCAzCut*fDCAzCut)) > 1 ){
544 continue; // 2D cut
545 }
546 }
547
548 // Extra TPC cuts (for systematic studies [!= -1])
549 if( fTPCchi2Cut != -1 && aodTrackMain->Chi2perNDF() > fTPCchi2Cut){
550 continue;
551 }
552 if( fNClustersTPCCut != -1 && aodTrackMain->GetTPCNcls() < fNClustersTPCCut){
553 continue;
554 }
555
556 // fill QA histograms
557 fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());
558 fHistDCA->Fill(dcaZMain,dcaXYMain);
559 fHistChi2->Fill(aodTrackMain->Chi2perNDF());
560 fHistPt->Fill(vPt);
561 fHistEta->Fill(vEta);
562 fHistPhi->Fill(vPhi);
563
564 // fill charge vector
565 chargeVector[0]->push_back(vCharge);
566 chargeVector[1]->push_back(vY);
567 chargeVector[2]->push_back(vEta);
568 chargeVector[3]->push_back(vPhi);
569 chargeVector[4]->push_back(vP[0]);
570 chargeVector[5]->push_back(vP[1]);
571 chargeVector[6]->push_back(vP[2]);
572 chargeVector[7]->push_back(vPt);
573 chargeVector[8]->push_back(vE);
574
575 // -------------------------------------------------------------
576 // for each track in main event loop over all tracks in mix event
577 for (Int_t iTracksMix = 0; iTracksMix < aodEventMix->GetNumberOfTracks(); iTracksMix++) {
578
579 AliAODTrack* aodTrackMix = dynamic_cast<AliAODTrack *>(aodEventMix->GetTrack(iTracksMix));
580 if (!aodTrackMix) {
581 AliError(Form("ERROR: Could not receive track %d", iTracksMix));
582 continue;
583 }
584
585 // AOD track cuts
586
587 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
588 // take only TPC only tracks
589 fHistTrackStats->Fill(aodTrackMix->GetFilterMap());
590 if(!aodTrackMix->TestFilterBit(nAODtrackCutBit)) continue;
591
592 vCharge = aodTrackMix->Charge();
593 vY = aodTrackMix->Y();
594 vEta = aodTrackMix->Eta();
595 vPhi = aodTrackMix->Phi() * TMath::RadToDeg();
596 vE = aodTrackMix->E();
597 vPt = aodTrackMix->Pt();
598 aodTrackMix->PxPyPz(vP);
599
600 Float_t dcaXYMix = aodTrackMix->DCA(); // this is the DCA from global track (not exactly what is cut on)
601 Float_t dcaZMix = aodTrackMix->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
602
603
604 // Kinematics cuts from ESD track cuts
605 if( vPt < fPtMin || vPt > fPtMax) continue;
606 if( vEta < fEtaMin || vEta > fEtaMax) continue;
607
608 // Extra DCA cuts (for systematic studies [!= -1])
609 if( fDCAxyCut != -1 && fDCAxyCut != -1){
610 if(TMath::Sqrt((dcaXYMix*dcaXYMix)/(fDCAxyCut*fDCAxyCut)+(dcaZMix*dcaZMix)/(fDCAzCut*fDCAzCut)) > 1 ){
611 continue; // 2D cut
612 }
613 }
614
615 // Extra TPC cuts (for systematic studies [!= -1])
616 if( fTPCchi2Cut != -1 && aodTrackMix->Chi2perNDF() > fTPCchi2Cut){
617 continue;
618 }
619 if( fNClustersTPCCut != -1 && aodTrackMix->GetTPCNcls() < fNClustersTPCCut){
620 continue;
621 }
622
623 // fill QA histograms
624 fHistClus->Fill(aodTrackMix->GetITSNcls(),aodTrackMix->GetTPCNcls());
625 fHistDCA->Fill(dcaZMix,dcaXYMix);
626 fHistChi2->Fill(aodTrackMix->Chi2perNDF());
627 fHistPt->Fill(vPt);
628 fHistEta->Fill(vEta);
629 fHistPhi->Fill(vPhi);
630
631 // fill charge vector
632 chargeVector[0]->push_back(vCharge);
633 chargeVector[1]->push_back(vY);
634 chargeVector[2]->push_back(vEta);
635 chargeVector[3]->push_back(vPhi);
636 chargeVector[4]->push_back(vP[0]);
637 chargeVector[5]->push_back(vP[1]);
638 chargeVector[6]->push_back(vP[2]);
639 chargeVector[7]->push_back(vPt);
640 chargeVector[8]->push_back(vE);
641
642
643
644 } //mix track loop
645
646 // calculate balance function for each track in main event
647 iMainTrackUsed++; // is needed to do no double counting in Balance Function calculation
648 if(iMainTrackUsed >= (Int_t)chargeVector[0]->size()) break; //do not allow more tracks than in mixed event!
649 fBalance->CalculateBalance(fCentrality,chargeVector,iMainTrackUsed,bSign);
650 // clean charge vector afterwards
651 for(Int_t i = 0; i < 9; i++){
652 chargeVector[i]->clear();
653 }
654
655
656 } //main track loop
657 // }//Vz cut
658 // }//Vy cut
659 // }//Vx cut
660 // }//proper vertexresolution
661 // }//proper number of contributors
662 // }//vertex object valid
663 }//triggered event
664 }//AOD analysis
665 }
666}
667
668AliMixInputEventHandler *AliAnalysisTaskEventMixingBF::SetupEventsForMixing() {
669 //sets the input handlers for event mixing
670
671 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
672 AliMultiInputEventHandler *inEvHMain = dynamic_cast<AliMultiInputEventHandler *>(mgr->GetInputEventHandler());
673 if (inEvHMain) {
674
675 AliMixInputEventHandler *mixEH = dynamic_cast<AliMixInputEventHandler *>(inEvHMain->GetFirstMultiInputHandler());
676 if (!mixEH) return nullptr;
677 if (mixEH->CurrentBinIndex() < 0) {
678 AliDebug(AliLog::kDebug + 1, "Current event mixEH->CurrentEntry() == -1");
679 return nullptr;
680 }
681 AliDebug(AliLog::kDebug, Form("Mixing %lld %d [%lld,%lld] %d", mixEH->CurrentEntry(), mixEH->CurrentBinIndex(), mixEH->CurrentEntryMain(), mixEH->CurrentEntryMix(), mixEH->NumberMixed()));
682
683 AliInputEventHandler *ihMainCurrent = inEvHMain->GetFirstInputEventHandler();
684 fMainEvent = ihMainCurrent->GetEvent();
685
686 AliMultiInputEventHandler *inEvHMixedCurrent = mixEH->GetFirstMultiInputHandler(); // for buffer = 1
687 AliInputEventHandler *ihMixedCurrent = inEvHMixedCurrent->GetFirstInputEventHandler();
688 fMixEvent = ihMixedCurrent->GetEvent();
689
690 return mixEH;
691 }
692 return NULL;
693}