adding electron rejection to eff cont task (same way as for default BF analysis)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskEffContBF.cxx
1 #include "TChain.h"
2 #include "TList.h"
3 #include "TTree.h"
4 #include "TH1F.h"
5 #include "TH2F.h"
6 #include "TH3D.h"
7 #include "TCanvas.h"
8 #include "TParticle.h"
9 #include "TObjArray.h"
10
11 #include "AliAnalysisTask.h"
12 #include "AliAnalysisManager.h"
13
14 #include "AliStack.h"
15 #include "AliMCEvent.h"
16 #include "AliAODEvent.h" 
17 #include "AliAODInputHandler.h"
18 #include "AliAODMCParticle.h"
19 #include "AliAODMCHeader.h"
20 #include "AliCentrality.h"
21 #include "AliGenEventHeader.h"
22
23 #include "AliLog.h"
24 #include "AliAnalysisTaskEffContBF.h"
25
26 // ---------------------------------------------------------------------
27 //
28 // Task for calculating the efficiency of the Balance Function 
29 // for single particles and pairs
30 // 
31 // ---------------------------------------------------------------------
32
33 ClassImp(AliAnalysisTaskEffContBF)
34
35 //________________________________________________________________________
36 AliAnalysisTaskEffContBF::AliAnalysisTaskEffContBF(const char *name) 
37   : AliAnalysisTaskSE(name), 
38     fAOD(0),
39     fArrayMC(0), 
40     fQAList(0), 
41     fOutputList(0), 
42     fHistEventStats(0), 
43     fHistCentrality(0),
44     fHistNMult(0), 
45     fHistVz(0), 
46     fHistNSigmaTPCvsPtbeforePID(0),
47     fHistNSigmaTPCvsPtafterPID(0),  
48     fHistContaminationSecondariesPlus(0),
49     fHistContaminationSecondariesMinus(0), //
50     fHistContaminationPrimariesPlus(0),
51     fHistContaminationPrimariesMinus(0), //
52     fHistGeneratedEtaPtPhiPlus(0), 
53     fHistSurvivedEtaPtPhiPlus(0),
54     fHistGeneratedEtaPtPhiMinus(0),
55     fHistSurvivedEtaPtPhiMinus(0),
56     fHistGeneratedEtaPtPlusControl(0),
57     fHistSurvivedEtaPtPlusControl(0),
58     fHistGeneratedEtaPtMinusControl(0),
59     fHistSurvivedEtaPtMinusControl(0),
60     fHistGeneratedEtaPtPlusPlus(0),
61     fHistSurvivedEtaPtPlusPlus(0),
62     fHistGeneratedEtaPtMinusMinus(0),
63     fHistSurvivedEtaPtMinusMinus(0),
64     fHistGeneratedEtaPtPlusMinus(0),
65     fHistSurvivedEtaPtPlusMinus(0),
66     fHistGeneratedPhiEtaPlusPlus(0),
67     fHistSurvivedPhiEtaPlusPlus(0),
68     fHistGeneratedPhiEtaMinusMinus(0),
69     fHistSurvivedPhiEtaMinusMinus(0),
70     fHistGeneratedPhiEtaPlusMinus(0),
71     fHistSurvivedPhiEtaPlusMinus(0),
72     fUseCentrality(kFALSE),
73     fCentralityEstimator("V0M"), 
74     fCentralityPercentileMin(0.0), 
75     fCentralityPercentileMax(5.0), 
76     fInjectedSignals(kFALSE),
77     fPIDResponse(0),
78     fElectronRejection(kFALSE),
79     fElectronOnlyRejection(kFALSE),
80     fElectronRejectionNSigma(-1.),
81     fElectronRejectionMinPt(0.),
82     fElectronRejectionMaxPt(1000.),
83     fVxMax(3.0), 
84     fVyMax(3.0),
85     fVzMax(10.), 
86     fAODTrackCutBit(128),
87     fMinNumberOfTPCClusters(80),
88     fMaxChi2PerTPCCluster(4.0),
89     fMaxDCAxy(3.0),
90     fMaxDCAz(3.0),
91     fMinPt(0.0),
92     fMaxPt(20.0),
93     fMinEta(-0.8), 
94     fMaxEta(0.8),
95     fEtaRangeMin(0.0), 
96     fEtaRangeMax(1.6), 
97     fPtRangeMin(0.0), 
98     fPtRangeMax(20.0), 
99     fEtaBin(100), //=100 (BF) 16
100     fdEtaBin(64), //=64 (BF)  16
101     fPtBin(100), //=100 (BF)  36
102     fHistSurvived4EtaPtPhiPlus(0),
103     fHistSurvived8EtaPtPhiPlus(0)
104   
105 {   
106   // Define input and output slots here
107   // Input slot #0 works with a TChain
108   DefineInput(0, TChain::Class());
109   // Output slot #0 id reserved by the base class for AOD
110   // Output slot #1 writes into a TH1 container
111   DefineOutput(1, TList::Class());
112   DefineOutput(2, TList::Class());
113 }
114
115 //________________________________________________________________________
116 void AliAnalysisTaskEffContBF::UserCreateOutputObjects() {
117   // Create histograms
118   // Called once
119
120   fQAList = new TList();
121   fQAList->SetName("QAList");
122   fQAList->SetOwner();
123   
124   fOutputList = new TList();
125   fOutputList->SetName("OutputList");
126   fOutputList->SetOwner();
127   
128   //Event stats.
129   TString gCutName[4] = {"Total","Offline trigger",
130                          "Vertex","Analyzed"};
131   fHistEventStats = new TH1F("fHistEventStats",
132                              "Event statistics;;N_{events}",
133                              4,0.5,4.5);
134   for(Int_t i = 1; i <= 4; i++)
135     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
136   fQAList->Add(fHistEventStats);
137
138   //====================================================//
139   Int_t ptBin = 36;
140   Int_t etaBin = 16;
141   Int_t phiBin = 100;
142
143   Double_t nArrayPt[37]={0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10.0, 15.0, 20.0};
144   Double_t nArrayEta[17]={-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}; 
145
146   Double_t nArrayPhi[phiBin+1];
147   for(Int_t iBin = 0; iBin <= phiBin; iBin++) 
148     nArrayPhi[iBin] = iBin*TMath::TwoPi()/phiBin;
149
150   Int_t detaBin = 16;
151   Int_t dphiBin = 100;
152   Double_t nArrayDPhi[dphiBin+1];
153   for(Int_t iBin = 0; iBin <= phiBin; iBin++) 
154     nArrayDPhi[iBin] = iBin*TMath::TwoPi()/dphiBin;
155   Double_t nArrayDEta[17]={0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6}; 
156   //====================================================//
157
158   //AOD analysis
159   fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
160                              1001,-0.5,100.5);
161   fQAList->Add(fHistCentrality);
162   
163   //multiplicity (good MC tracks)
164   TString histName;
165   histName = "fHistNMult";
166   fHistNMult = new TH1F(histName.Data(), 
167                         ";N_{mult.}",
168                         200,0,20000);
169   fQAList->Add(fHistNMult);
170
171   //Vz addition+++++++++++++++++++++++++++++
172   fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
173   fQAList->Add(fHistVz);
174
175   //Electron cuts -> PID QA
176   fHistNSigmaTPCvsPtbeforePID = new TH2F ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore",200, 0, 20, 200, -10, 10); 
177   fQAList->Add(fHistNSigmaTPCvsPtbeforePID);
178
179   fHistNSigmaTPCvsPtafterPID = new TH2F ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter",200, 0, 20, 200, -10, 10); 
180   fQAList->Add(fHistNSigmaTPCvsPtafterPID);
181
182   //Contamination for Secondaries 
183   fHistContaminationSecondariesPlus = new TH3D("fHistContaminationSecondariesPlus","Secondaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
184   fOutputList->Add(fHistContaminationSecondariesPlus);
185
186   fHistContaminationSecondariesMinus = new TH3D("fHistContaminationSecondariesMinus","Secondaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
187   fOutputList->Add(fHistContaminationSecondariesMinus);
188
189   //Contamination for Primaries
190   fHistContaminationPrimariesPlus = new TH3D("fHistContaminationPrimariesPlus","Primaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
191   fOutputList->Add(fHistContaminationPrimariesPlus);
192
193   fHistContaminationPrimariesMinus = new TH3D("fHistContaminationPrimariesMinus","Primaries;#eta;p_{T} (GeV/c);#varphi",etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
194   fOutputList->Add(fHistContaminationPrimariesMinus);
195   
196   //eta vs pt for MC positives
197   fHistGeneratedEtaPtPhiPlus = new TH3D("fHistGeneratedEtaPtPhiPlus",
198                                         "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
199                                         etaBin,nArrayEta, ptBin, nArrayPt,phiBin, nArrayPhi);
200   // fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,nArrayPhi);
201   fOutputList->Add(fHistGeneratedEtaPtPhiPlus);
202   fHistSurvivedEtaPtPhiPlus = new TH3D("fHistSurvivedEtaPtPhiPlus",
203                                        "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
204                                        etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
205   fOutputList->Add(fHistSurvivedEtaPtPhiPlus);
206   
207   //eta vs pt for MC negatives
208   fHistGeneratedEtaPtPhiMinus = new TH3D("fHistGeneratedEtaPtPhiMinus",
209                                          "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
210                                          etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
211   fOutputList->Add(fHistGeneratedEtaPtPhiMinus);
212   fHistSurvivedEtaPtPhiMinus = new TH3D("fHistSurvivedEtaPtPhiMinus",
213                                         "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
214                                         etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
215   fOutputList->Add(fHistSurvivedEtaPtPhiMinus);
216  
217   //eta vs pt for MC positives (control)
218   fHistGeneratedEtaPtPlusControl = new TH2F("fHistGeneratedEtaPtPlusControl",
219                                             "Generated positive primaries;#eta;p_{T} (GeV/c)",
220                                             etaBin,nArrayEta,ptBin,nArrayPt);
221   fOutputList->Add(fHistGeneratedEtaPtPlusControl);
222   fHistSurvivedEtaPtPlusControl = new TH2F("fHistSurvivedEtaPtPlusControl",
223                                            "Survived positive primaries;#eta;p_{T} (GeV/c)",
224                                            etaBin,nArrayEta,ptBin,nArrayPt);
225   fOutputList->Add(fHistSurvivedEtaPtPlusControl);
226   
227   //eta vs pt for MC negatives (control)
228   fHistGeneratedEtaPtMinusControl = new TH2F("fHistGeneratedEtaPtMinusControl",
229                                              "Generated positive primaries;#eta;p_{T} (GeV/c)",
230                                              etaBin,nArrayEta,ptBin,nArrayPt);
231   fOutputList->Add(fHistGeneratedEtaPtMinusControl);
232   fHistSurvivedEtaPtMinusControl = new TH2F("fHistSurvivedEtaPtMinusControl",
233                                             "Survived positive primaries;#eta;p_{T} (GeV/c)",
234                                             etaBin,nArrayEta,ptBin,nArrayPt);
235   fOutputList->Add(fHistSurvivedEtaPtMinusControl);
236   
237   //eta vs pt for MC ++
238   fHistGeneratedEtaPtPlusPlus = new TH2F("fHistGeneratedEtaPtPlusPlus",
239                                          "Generated ++ primaries;#Delta#eta;p_{T} (GeV/c)",
240                                          detaBin,nArrayDEta,ptBin,nArrayPt);
241   fOutputList->Add(fHistGeneratedEtaPtPlusPlus);
242   fHistSurvivedEtaPtPlusPlus = new TH2F("fHistSurvivedEtaPtPlusPlus",
243                                         "Survived ++ primaries;#Delta#eta;p_{T} (GeV/c)",
244                                         detaBin,nArrayDEta,ptBin,nArrayPt);
245   fOutputList->Add(fHistSurvivedEtaPtPlusPlus);
246   
247   //eta vs pt for MC --
248   fHistGeneratedEtaPtMinusMinus = new TH2F("fHistGeneratedEtaPtMinusMinus",
249                                            "Generated -- primaries;#Delta#eta;p_{T} (GeV/c)",
250                                            detaBin,nArrayDEta,ptBin,nArrayPt);
251   fOutputList->Add(fHistGeneratedEtaPtMinusMinus);
252   fHistSurvivedEtaPtMinusMinus = new TH2F("fHistSurvivedEtaPtMinusMinus",
253                                           "Survived -- primaries;#Delta#eta;p_{T} (GeV/c)",
254                                           detaBin,nArrayDEta,ptBin,nArrayPt);
255   fOutputList->Add(fHistSurvivedEtaPtMinusMinus);
256  
257   //eta vs pt for MC +-
258   fHistGeneratedEtaPtPlusMinus = new TH2F("fHistGeneratedEtaPtPlusMinus",
259                                           "Generated +- primaries;#Delta#eta;p_{T} (GeV/c)",
260                                           detaBin,nArrayDEta,ptBin,nArrayPt);
261   fOutputList->Add(fHistGeneratedEtaPtPlusMinus);
262   fHistSurvivedEtaPtPlusMinus = new TH2F("fHistSurvivedEtaPtPlusMinus",
263                                          "Survived +- primaries;#Delta#eta;p_{T} (GeV/c)",
264                                          detaBin,nArrayDEta,ptBin,nArrayPt);
265   fOutputList->Add(fHistSurvivedEtaPtPlusMinus);
266  
267   //=============================//
268   //phi vs eta for MC ++
269   fHistGeneratedPhiEtaPlusPlus = new TH2F("fHistGeneratedPhiEtaPlusPlus",
270                                           "Generated ++ primaries;#Delta#phi",
271                                           dphiBin,nArrayDPhi,detaBin,nArrayDEta);
272   fOutputList->Add(fHistGeneratedPhiEtaPlusPlus);
273   fHistSurvivedPhiEtaPlusPlus = new TH2F("fHistSurvivedPhiEtaPlusPlus",
274                                          "Survived ++ primaries;#Delta#phi;#Delta#eta",
275                                          dphiBin,nArrayDPhi,detaBin,nArrayDEta);
276   fOutputList->Add(fHistSurvivedPhiEtaPlusPlus);
277   
278   //phi vs eta for MC --
279   fHistGeneratedPhiEtaMinusMinus = new TH2F("fHistGeneratedPhiEtaMinusMinus",
280                                             "Generated -- primaries;#Delta#phi;#Delta#eta",
281                                             dphiBin,nArrayDPhi,detaBin,nArrayDEta);
282   fOutputList->Add(fHistGeneratedPhiEtaMinusMinus);
283   fHistSurvivedPhiEtaMinusMinus = new TH2F("fHistSurvivedPhiEtaMinusMinus",
284                                            "Survived -- primaries;#Delta#phi;#Delta#eta",
285                                            dphiBin,nArrayDPhi,detaBin,nArrayDEta);
286   fOutputList->Add(fHistSurvivedPhiEtaMinusMinus);
287   
288   //phi vs eta for MC +-
289   fHistGeneratedPhiEtaPlusMinus = new TH2F("fHistGeneratedPhiEtaPlusMinus",
290                                            "Generated +- primaries;#Delta#phi;#Delta#eta",
291                                            dphiBin,nArrayDPhi,detaBin,nArrayDEta);
292   fOutputList->Add(fHistGeneratedPhiEtaPlusMinus);
293   fHistSurvivedPhiEtaPlusMinus = new TH2F("fHistSurvivedPhiEtaPlusMinus",
294                                           "Survived +- primaries;#Delta#phi;#Delta#eta",
295                                           dphiBin,nArrayDPhi,detaBin,nArrayDEta);
296   fOutputList->Add(fHistSurvivedPhiEtaPlusMinus);
297   
298   fHistSurvived4EtaPtPhiPlus = new TH3F("fHistSurvived4EtaPtPhiPlus",
299                                         "Survived4 + primaries;#eta;p_{T} (GeV/c);#phi",
300                                         etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
301   fOutputList->Add(fHistSurvived4EtaPtPhiPlus);
302   fHistSurvived8EtaPtPhiPlus = new TH3F("fHistSurvived8EtaPtPhiPlus",
303                                         "Survived8 + primaries;#eta;p_{T} (GeV/c);#phi",
304                                         etaBin,nArrayEta,ptBin,nArrayPt,phiBin,nArrayPhi);
305   fOutputList->Add(fHistSurvived8EtaPtPhiPlus);
306   
307   //fQAList->Print();
308   //fOutputList->Print(); 
309   PostData(1, fQAList);
310   PostData(2, fOutputList);
311 }
312
313 //________________________________________________________________________
314 void AliAnalysisTaskEffContBF::UserExec(Option_t *) {
315   // Main loop
316   // Called for each event
317   
318   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
319   if (!fAOD) {
320     printf("ERROR: fAOD not available\n");
321     return;
322   }
323
324   fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));   
325   if (!fArrayMC)  
326     AliFatal("No array of MC particles found !!!"); // MW  no AliFatal use return values   
327   
328   AliMCEvent* mcEvent = MCEvent();
329   if (!mcEvent) {
330     AliError("ERROR: Could not retrieve MC event");
331     return;
332   }
333
334   // PID Response task active?
335   if(fElectronRejection) {
336     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
337     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
338   }
339
340   // ==============================================================================================
341   // Copy from AliAnalysisTaskPhiCorrelations:
342   // For productions with injected signals, figure out above which label to skip particles/tracks
343   Int_t skipParticlesAbove = 0;
344   if (fInjectedSignals)
345   {
346     AliGenEventHeader* eventHeader = 0;
347     Int_t headers = 0;
348     
349     // AOD only
350     AliAODMCHeader* header = (AliAODMCHeader*) fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
351     if (!header)
352       AliFatal("fInjectedSignals set but no MC header found");
353     
354     headers = header->GetNCocktailHeaders();
355     eventHeader = header->GetCocktailHeader(0);
356     
357     
358     if (!eventHeader)
359       {
360         // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing 
361         // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
362         AliError("First event header not found. Skipping this event.");
363         return;
364       }
365     
366     skipParticlesAbove = eventHeader->NProduced();
367     AliInfo(Form("Injected signals in this event (%d headers). Keeping particles/tracks of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove)); 
368   }
369   // ==============================================================================================
370
371   
372   // arrays for 2 particle histograms
373   Int_t nMCLabelCounter         = 0;
374   const Int_t maxMCLabelCounter = 20000;
375   
376   Double_t eta[maxMCLabelCounter];
377   Double_t pt[maxMCLabelCounter];
378   Double_t phi[maxMCLabelCounter];
379   Int_t level[maxMCLabelCounter];
380   Int_t charge[maxMCLabelCounter];
381   
382   //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fAOD->GetNumberOfTracks()));
383
384   fHistEventStats->Fill(1); //all events
385   
386   //Centrality stuff
387   AliAODHeader *header = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
388   Double_t nCentrality = 0;
389   
390   if(fUseCentrality) {
391     AliCentrality *centrality = header->GetCentralityP();
392     nCentrality =centrality->GetCentralityPercentile(fCentralityEstimator.Data());
393     
394
395     if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
396                                              fCentralityPercentileMax,
397                                              fCentralityEstimator.Data()))
398       return;
399     else {    
400       fHistEventStats->Fill(2); //triggered + centrality
401       fHistCentrality->Fill(nCentrality);
402     }
403   }
404   //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
405
406   const AliAODVertex *vertex = fAOD->GetPrimaryVertex(); 
407   if(vertex) {
408     if(vertex->GetNContributors() > 0) {
409       Double32_t fCov[6];    
410       vertex->GetCovarianceMatrix(fCov);   
411       if(fCov[5] != 0) {
412         fHistEventStats->Fill(3); //events with a proper vertex
413         if(TMath::Abs(vertex->GetX()) < fVxMax) {    // antes Xv
414           //Printf("X Vertex: %lf", vertex->GetX());
415           //Printf("Y Vertex: %lf", vertex->GetY());
416           if(TMath::Abs(vertex->GetY()) < fVyMax) {  // antes Yv
417             if(TMath::Abs(vertex->GetZ()) < fVzMax) {  // antes Zv
418               //Printf("Z Vertex: %lf", vertex->GetZ());
419               
420               fHistEventStats->Fill(4); //analyzed events
421               fHistVz->Fill(vertex->GetZ()); 
422               
423               //++++++++++++++++++CONTAMINATION++++++++++++++++++//
424               Int_t nGoodAODTracks = fAOD->GetNumberOfTracks();
425               Int_t nMCParticles = mcEvent->GetNumberOfTracks();
426               TArrayI labelMCArray(nMCParticles);
427               
428               for(Int_t jTracks = 0; jTracks < nGoodAODTracks; jTracks++) {
429                 AliAODTrack* track = fAOD->GetTrack(jTracks);
430                 if(!track) continue;
431                 
432                 if (!track->TestFilterBit(fAODTrackCutBit))
433                   continue;
434                 
435                 //acceptance
436                 if(TMath::Abs(track->Eta()) > fMaxEta) 
437                   continue;
438                 if((track->Pt() > fMaxPt)||(track->Pt() <  fMinPt)) 
439                   continue;
440                 
441                 Double_t phiRad = track->Phi(); 
442                 
443                 Int_t label = TMath::Abs(track->GetLabel());
444                 if(label > nMCParticles) continue;
445                 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) mcEvent->GetTrack(label); 
446                 Short_t gAODmcCharge = AODmcTrack->Charge();////
447                 //fHistContaminationPrimaries->Fill(track->Eta(),track->Pt(),phiDeg);
448                 //if (!(AODmcTrack->IsPhysicalPrimary())) {
449                 //fHistContaminationSecondaries->Fill(track->Eta(),track->Pt(),phiDeg);
450                 //}
451
452                 // ==============================================================================================
453                 // Partial copy from AliAnalyseLeadingTrackUE::RemoveInjectedSignals:
454                 // Skip tracks that come from injected signals
455                 if (fInjectedSignals)
456                   {    
457      
458                     AliAODMCParticle* mother = AODmcTrack;
459                     
460                     // find the primary mother (if not already physical primary)
461                     while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
462                       {
463                         if (((AliAODMCParticle*)mother)->GetMother() < 0)
464                           {
465                             mother = 0;
466                             break;
467                           }
468                         
469                         mother = (AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
470                         if (!mother)
471                           break;
472                       }
473                     
474                     
475                     if (!mother)
476                       {
477                         AliError(Form("WARNING: No mother found for particle %d:", AODmcTrack->GetLabel()));
478                         continue;
479                       }
480
481                     if (mother->GetLabel() >= skipParticlesAbove)
482                       {
483                         //AliInfo(Form("Remove particle %d (>= %d)",mother->GetLabel(),skipParticlesAbove));
484                         continue;
485                       }
486                   }
487                 // ==============================================================================================
488
489                 if (AODmcTrack->IsPhysicalPrimary()) {
490                   if(gAODmcCharge > 0){
491                     fHistContaminationPrimariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
492                   }
493                   if(gAODmcCharge < 0){
494                     fHistContaminationPrimariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
495                   }
496                 }
497                 else{
498                   if(gAODmcCharge > 0){
499                     fHistContaminationSecondariesPlus->Fill(track->Eta(),track->Pt(),phiRad);
500                   }
501                   if(gAODmcCharge < 0){
502                     fHistContaminationSecondariesMinus->Fill(track->Eta(),track->Pt(),phiRad);
503                   }
504                 }
505               }//loop over tracks
506               //++++++++++++++++++CONTAMINATION++++++++++++++++++//
507               
508               //++++++++++++++++++EFFICIENCY+++++++++++++++++++++//
509               for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
510                 AliAODMCParticle *mcTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); 
511                 if (!mcTrack) {
512                   AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
513                   continue;
514                 }
515                 
516                 //exclude particles generated out of the acceptance
517                 Double_t vz = mcTrack->Zv();
518                 if (TMath::Abs(vz) > 50.) continue;
519                 //acceptance
520                 if(TMath::Abs(mcTrack->Eta()) > fMaxEta) 
521                   continue;
522                 if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt)) 
523                   continue;
524                 
525                 if(!mcTrack->IsPhysicalPrimary()) continue;   
526                 
527                 Short_t gMCCharge = mcTrack->Charge();
528                 Double_t phiRad = mcTrack->Phi(); 
529                 
530                 if(gMCCharge > 0)
531                   fHistGeneratedEtaPtPhiPlus->Fill(mcTrack->Eta(),
532                                                    mcTrack->Pt(),
533                                                    phiRad);
534                 else if(gMCCharge < 0)
535                   fHistGeneratedEtaPtPhiMinus->Fill(mcTrack->Eta(),
536                                                     mcTrack->Pt(),
537                                                     phiRad);
538                 
539                 Bool_t labelTPC = kTRUE;
540                 if(labelTPC) {
541                   labelMCArray.AddAt(iTracks,nMCLabelCounter);
542                   if(nMCLabelCounter >= maxMCLabelCounter){
543                     AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
544                     break;
545                   }
546                   //fill the arrays for 2 particle analysis
547                   eta[nMCLabelCounter]    = mcTrack->Eta();
548                   pt[nMCLabelCounter]     = mcTrack->Pt();
549                   phi[nMCLabelCounter]    = mcTrack->Phi();
550                   charge[nMCLabelCounter] = gMCCharge;
551                   
552                   level[nMCLabelCounter]  = 1;
553                   nMCLabelCounter += 1;
554                 }  
555               }//loop over MC particles
556               
557               fHistNMult->Fill(nMCLabelCounter);
558               
559               //AOD track loop
560               Int_t nGoodTracks = fAOD->GetNumberOfTracks();   
561               TArrayI labelArray(nGoodTracks);
562               Int_t labelCounter = 0;
563               
564               for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
565                 AliAODTrack *trackAOD = static_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));    
566                 if(!trackAOD) continue;
567                 
568                 //track cuts
569                 if (!trackAOD->TestFilterBit(fAODTrackCutBit)) 
570                   continue;
571
572                 Int_t label = TMath::Abs(trackAOD->GetLabel()); 
573                 if(IsLabelUsed(labelArray,label)) continue;
574                 labelArray.AddAt(label,labelCounter);
575                 labelCounter += 1;
576                 
577                 Bool_t iFound = kFALSE;
578                 Int_t mcGoods = nMCLabelCounter;
579                 for (Int_t k = 0; k < mcGoods; k++) {
580                   Int_t mcLabel = labelMCArray.At(k);
581                   iFound = kFALSE;
582                   
583                   if (mcLabel != TMath::Abs(label)) continue;
584                   if(mcLabel != label) continue;                    
585                   if(label > trackAOD->GetLabel()) continue; 
586                   
587                   //acceptance
588                   if(TMath::Abs(trackAOD->Eta()) > fMaxEta) 
589                     continue;
590                   if((trackAOD->Pt() > fMaxPt)||(trackAOD->Pt() <  fMinPt)) 
591                     continue;
592                   
593                   Short_t gCharge = trackAOD->Charge();
594                   Double_t phiRad = trackAOD->Phi();
595
596                   //===========================PID (so far only for electron rejection)===============================//                    
597                   if(fElectronRejection) {
598                     
599                     // get the electron nsigma
600                     Double_t nSigma = fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kElectron);
601                     fHistNSigmaTPCvsPtbeforePID->Fill(trackAOD->Pt(),nSigma);               
602                     
603                     // check only for given momentum range
604                     if( trackAOD->Pt() > fElectronRejectionMinPt && trackAOD->Pt() < fElectronRejectionMaxPt ){
605                       
606                       //look only at electron nsigma
607                       if(!fElectronOnlyRejection){
608                         
609                         //Make the decision based on the n-sigma of electrons
610                         if(TMath::Abs(nSigma) < fElectronRejectionNSigma) continue;
611                       }
612                       else{
613                         
614                         Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kPion));
615                         Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kKaon));
616                         Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackAOD,(AliPID::EParticleType)AliPID::kProton));
617                         
618                         //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
619                         if(TMath::Abs(nSigma) < fElectronRejectionNSigma
620                            && nSigmaPions   > fElectronRejectionNSigma
621                            && nSigmaKaons   > fElectronRejectionNSigma
622                            && nSigmaProtons > fElectronRejectionNSigma ) continue;
623                       }
624                     }
625                     
626                     fHistNSigmaTPCvsPtafterPID->Fill(trackAOD->Pt(),nSigma);                
627
628                   }
629                   //===========================end of PID (so far only for electron rejection)===============================//           
630                   
631                   if(TMath::Abs(trackAOD->Eta()) < fMaxEta && trackAOD->Pt() > fMinPt&&trackAOD->Pt() < fMaxPt){ 
632                     level[k]  = 2;
633                     
634                     if(gCharge > 0)
635                       fHistSurvivedEtaPtPhiPlus->Fill(trackAOD->Eta(),
636                                                       trackAOD->Pt(),
637                                                       phiRad);
638                     else if(gCharge < 0)
639                       fHistSurvivedEtaPtPhiMinus->Fill(trackAOD->Eta(),
640                                                        trackAOD->Pt(),
641                                                        phiRad);
642                   }//tracks                
643                 }//end of mcGoods
644               }//AOD track loop
645               
646               labelMCArray.Reset();
647               labelArray.Reset();              
648               
649             }//Vz cut
650           }//Vy cut
651         }//Vx cut
652       }//Vz resolution
653     }//number of contributors
654   }//valid vertex  
655   
656   // Here comes the 2 particle analysis
657   // loop over all good MC particles
658   for (Int_t i = 0; i < nMCLabelCounter ; i++) {
659     // control 1D histograms (charge might be different?)
660     if(charge[i] > 0){
661       if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
662       if(level[i] > 1) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
663     }
664     else if(charge[i] < 0){
665       if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
666       if(level[i] > 1) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
667     }
668     
669     
670     for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
671       
672       if(charge[i] > 0 && charge[j] > 0 ){
673         if(level[i] > 0 && level[j] > 0) {  
674           fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
675           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
676             fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
677         }
678         if(level[i] > 1 && level[j] > 1) {
679           fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
680           if (TMath::Abs(phi[i]-phi[j]) < TMath::Pi())
681             fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
682         }
683       }
684       
685       else if(charge[i] < 0 && charge[j] < 0 ){
686         if(level[i] > 0 && level[j] > 0) {
687           fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);     
688           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
689             fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));              
690         }
691         if(level[i] > 2 && level[j] > 1) {
692           fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
693           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
694             fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
695         }
696       }
697       
698       else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
699         if(level[i] > 0 && level[j] > 0) {
700           fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
701           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
702             fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));   
703         }
704         if(level[i] > 2 && level[j] > 1) {
705           fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
706           if (TMath::Abs(phi[i]-phi[j]) <  TMath::Pi())
707             fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
708         }       
709       }
710     }
711   }
712 }
713
714 //________________________________________________________________________
715 void AliAnalysisTaskEffContBF::Terminate(Option_t *) {
716   // Draw result to the screen
717   // Called once at the end of the query
718
719   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
720   if (!fOutputList) {
721     printf("ERROR: Output list not available\n");
722     return;
723   }
724 }
725
726 //____________________________________________________________________//
727 Bool_t AliAnalysisTaskEffContBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
728   //Checks if the label is used already
729   Bool_t status = kFALSE;
730   for(Int_t i = 0; i < labelArray.GetSize(); i++) {
731     if(labelArray.At(i) == label)
732       status = kTRUE;
733   }
734
735   return status;
736 }