AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskEventMixingBF.cxx
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
41 ClassImp(AliAnalysisTaskEventMixingBF)
42
43 //________________________________________________________________________
44 AliAnalysisTaskEventMixingBF::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 //________________________________________________________________________
138 AliAnalysisTaskEventMixingBF::~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 //________________________________________________________________________
162 void 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 //________________________________________________________________________
356 void AliAnalysisTaskEventMixingBF::UserExec(Option_t *) {
357   // Main loop
358   // Called for each event
359   // NOTHING TO DO for event mixing!
360 }      
361
362 //________________________________________________________________________
363 void  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 //________________________________________________________________________
380 void 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
388 void 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       
437      AliAODHeader *aodHeaderMain = dynamic_cast<AliAODHeader*>(aodEventMain->GetHeader());
438      if(!aodHeaderMain) AliFatal("Not a standard AOD");  
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
668 AliMixInputEventHandler *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