Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / PhiEffMC / AliAnalysisTaskPhiEffMc.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //-----------------------------------------------------------------
17 //         AliAnalysisTaskPhiEffMc class
18 //-----------------------------------------------------------------
19
20 #include "TChain.h"
21 #include "TTree.h"
22 #include "TLegend.h"
23 #include "TH1F.h"
24 #include "TH2F.h"
25 #include "TH3F.h"
26 #include "THnSparse.h"
27 #include "TCanvas.h"
28 #include "TObjArray.h"
29 #include "AliAnalysisTask.h"
30 #include "AliAnalysisManager.h"
31 #include "AliVTrack.h"
32 #include "AliAODMCParticle.h"
33 #include "AliVParticle.h"
34 #include "AliAODEvent.h"
35 #include "AliAODInputHandler.h"
36 #include "AliAnalysisTaskPhiEffMc.h"
37 #include "AliAnalysisTaskESDfilter.h"
38 #include "AliAnalysisDataContainer.h"
39 #include "AliCentrality.h"
40 #include "TProof.h"
41 #include "AliPID.h"
42 #include "AliVEvent.h"
43 #include "AliPIDResponse.h"
44 #include "PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODTrackCuts.h"
45 #include "PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.h"
46 #include "AliStack.h"
47 #include <TMCProcess.h>
48 #include "AliAnalyseLeadingTrackUE.h"
49
50 #include <iostream>
51
52 using namespace std;
53
54 ClassImp(AliAnalysisTaskPhiEffMc) 
55
56 //________________________________________________________________________
57 AliAnalysisTaskPhiEffMc::AliAnalysisTaskPhiEffMc(const char *name) : AliAnalysisTaskSE(name), 
58   fAOD(0x0), 
59   fIsMC(0),
60   fOutput(0x0),
61   fHelperPID(0x0),
62   fTrackCuts(0x0),
63   fEventCuts(0x0),
64   fPtCut(0.),
65   fDoPID(kTRUE)
66 {
67   // Default constructor
68   
69
70   DefineInput(0, TChain::Class());
71   DefineOutput(1, TList::Class());
72   DefineOutput(2, AliHelperPID::Class());
73   DefineOutput(3, AliSpectraAODEventCuts::Class());
74   DefineOutput(4, AliSpectraAODTrackCuts::Class());
75 }
76 //________________________________________________________________________
77 //________________________________________________________________________
78 void AliAnalysisTaskPhiEffMc::UserCreateOutputObjects()
79 {
80   Printf("\n\n\n\n\n\n In CreateOutput Object:");
81   
82   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
83   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
84   if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");
85
86   fOutput = new TList();
87   fOutput->SetOwner();
88   fOutput->SetName("BlaBla");
89
90   // event histogram
91   const Int_t nEventDim = 4;
92   //                      cent  vz MCmult RecoMult
93   Int_t nEventBins[4] =   { 10,  41,  50,   50};
94   Double_t nEventMin[4] = {  0, -10,   0,    0};
95   Double_t nEventMax[4] = {100,  10, 100,  100};
96
97   THnSparseF* hEvent = new THnSparseF("hEvent", "hEvent", nEventDim, nEventBins, nEventMin, nEventMax);
98   hEvent->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
99   hEvent->GetAxis(1)->SetTitle("v_{z} (cm)");
100   hEvent->GetAxis(2)->SetTitle("# MC particles");
101   hEvent->GetAxis(3)->SetTitle("# reconstructed tracks");
102   fOutput->Add(hEvent);
103
104   TH3F* hKcorr = new TH3F("hKcorr","hKcorr",100,0,10,100,0,10,100,0,10);
105   hKcorr->GetXaxis()->SetTitle("#phi p_{T} (GeV/c)");
106   hKcorr->GetYaxis()->SetTitle("k+ p_{T} (GeV/c)");
107   hKcorr->GetZaxis()->SetTitle("k- p_{T} (GeV/c)");
108   fOutput->Add(hKcorr);
109
110   // single charged particles -- bins for THnSparse histograms
111   const Int_t nTrackDim = 6;
112   //                                PID ID, pt, y, cent, eta, phi
113   Int_t    nTrackBins[nTrackDim] =   {   7, 50, 20,  10,  20, 30};
114   Double_t nTrackMin[nTrackDim] =    {-3.5,  0, -1,   0,  -1,  0};
115   Double_t nTrackMax[nTrackDim] =    { 3.5,  5,  1, 100,   1, 2*TMath::Pi()};
116
117   // single charged particles -- Monte Carlo particles
118   THnSparseF* hTrackMc = new THnSparseF("hTrackMc", "hTrackMc", nTrackDim, nTrackBins, nTrackMin, nTrackMax);
119   hTrackMc->GetAxis(0)->SetTitle("PID ID * Charge");
120   hTrackMc->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
121   hTrackMc->GetAxis(2)->SetTitle("y");
122   hTrackMc->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
123   hTrackMc->GetAxis(4)->SetTitle("#eta");
124   hTrackMc->GetAxis(5)->SetTitle("#phi (rad)");
125   fOutput->Add(hTrackMc);
126
127   // single charged particles -- reconstructed tracks
128   THnSparseF* hTrackReco = new THnSparseF("hTrackReco", "hTrackReco", nTrackDim, nTrackBins, nTrackMin, nTrackMax);
129   hTrackReco->GetAxis(0)->SetTitle("PID ID * Charge");
130   hTrackReco->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
131   hTrackReco->GetAxis(2)->SetTitle("y");
132   hTrackReco->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
133   hTrackReco->GetAxis(4)->SetTitle("#eta");
134   hTrackReco->GetAxis(5)->SetTitle("#phi (rad)");
135   fOutput->Add(hTrackReco);
136
137   // single charged particles -- Monte Carlo particles matched to reconstructed tracks
138   THnSparseF* hTrackMatch = new THnSparseF("hTrackMatch", "hTrackMatch", nTrackDim, nTrackBins, nTrackMin, nTrackMax);
139   hTrackMatch->GetAxis(0)->SetTitle("PID ID * Charge");
140   hTrackMatch->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
141   hTrackMatch->GetAxis(2)->SetTitle("y");
142   hTrackMatch->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
143   hTrackMatch->GetAxis(4)->SetTitle("#eta");
144   hTrackMatch->GetAxis(5)->SetTitle("#phi (rad)");
145   fOutput->Add(hTrackMatch);
146
147   // kaon pairs -- bins for THnSparse histograms
148   const Int_t nPairDim = 7;
149   //                             InvMass, pt, y, cent, eta, phi,          pair ID
150   Int_t    nPairBins[nPairDim] =   { 100, 50, 20,  10,  20, 30,              3};
151   Double_t nPairMin[nPairDim] =    {0.98,  0, -1,   0,  -1,  0,           -0.5};
152   Double_t nPairMax[nPairDim] =    { 1.1,  5,  1, 100,   1, 2*TMath::Pi(), 2.5};
153   // pair ID = 0 unlike sign k+k-
154   //           1 like sign pos k+k+
155   //           2 like sign neg k-k-
156
157   // kaon pairs -- Monte Carlo particles -- real phi's
158   THnSparseF* hPairMcPhi = new THnSparseF("hPairMcPhi", "hPairMcPhi", nPairDim, nPairBins, nPairMin, nPairMax);
159   hPairMcPhi->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)");
160   hPairMcPhi->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
161   hPairMcPhi->GetAxis(2)->SetTitle("y");
162   hPairMcPhi->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
163   hPairMcPhi->GetAxis(4)->SetTitle("#eta");
164   hPairMcPhi->GetAxis(5)->SetTitle("#phi (rad)");
165   hPairMcPhi->GetAxis(6)->SetTitle("pair ID");
166   fOutput->Add(hPairMcPhi);
167
168   // kaon pairs -- Monte Carlo particles -- real phi's with kaon daughters that pass pt and eta cuts
169   THnSparseF* hPairMcPhiCuts = new THnSparseF("hPairMcPhiCuts", "hPairMcPhiCuts", nPairDim, nPairBins, nPairMin, nPairMax);
170   hPairMcPhiCuts->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)");
171   hPairMcPhiCuts->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
172   hPairMcPhiCuts->GetAxis(2)->SetTitle("y");
173   hPairMcPhiCuts->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
174   hPairMcPhiCuts->GetAxis(4)->SetTitle("#eta");
175   hPairMcPhiCuts->GetAxis(5)->SetTitle("#phi (rad)");
176   hPairMcPhiCuts->GetAxis(6)->SetTitle("pair ID");
177   fOutput->Add(hPairMcPhiCuts);
178
179   // kaon pairs -- Monte Carlo particles -- phi's reconstructed from MC kaons
180   THnSparseF* hPairMc = new THnSparseF("hPairMc", "hPairMc", nPairDim, nPairBins, nPairMin, nPairMax);
181   hPairMc->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)");
182   hPairMc->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
183   hPairMc->GetAxis(2)->SetTitle("y");
184   hPairMc->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
185   hPairMc->GetAxis(4)->SetTitle("#eta");
186   hPairMc->GetAxis(5)->SetTitle("#phi (rad)");
187   hPairMc->GetAxis(6)->SetTitle("pair ID");
188   fOutput->Add(hPairMc);
189
190   // kaon pairs -- reconstructed tracks
191   THnSparseF* hPairReco = new THnSparseF("hPairReco", "hPairReco", nPairDim, nPairBins, nPairMin, nPairMax);
192   hPairReco->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)");
193   hPairReco->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
194   hPairReco->GetAxis(2)->SetTitle("y");
195   hPairReco->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
196   hPairReco->GetAxis(4)->SetTitle("#eta");
197   hPairReco->GetAxis(5)->SetTitle("#phi (rad)");
198   hPairReco->GetAxis(6)->SetTitle("pair ID");
199   fOutput->Add(hPairReco);
200
201   // kaon pairs -- Monte Carlo particles matched to reconstructed tracks
202   THnSparseF* hPairMatch = new THnSparseF("hPairMatch", "hPairMatch", nPairDim, nPairBins, nPairMin, nPairMax);
203   hPairMatch->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)");
204   hPairMatch->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
205   hPairMatch->GetAxis(2)->SetTitle("y");
206   hPairMatch->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
207   hPairMatch->GetAxis(4)->SetTitle("#eta");
208   hPairMatch->GetAxis(5)->SetTitle("#phi (rad)");
209   hPairMatch->GetAxis(6)->SetTitle("pair ID");
210   fOutput->Add(hPairMatch);
211
212   PostData(1, fOutput);
213   PostData(2, fHelperPID);
214   PostData(3, fEventCuts);
215   PostData(4, fTrackCuts);
216 }
217
218 //________________________________________________________________________
219 void AliAnalysisTaskPhiEffMc::UserExec(Option_t *)
220 {
221   // main event loop
222   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
223   if (!fAOD) {
224     AliWarning("ERROR: AliAODEvent not available \n");
225     return;
226   }
227   
228   if (strcmp(fAOD->ClassName(), "AliAODEvent"))
229     {
230       AliFatal("Not processing AODs");
231     }
232
233   if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
234
235   Double_t cent = fEventCuts->GetCent();
236
237   THnSparseF* hEvent = (THnSparseF*)fOutput->FindObject("hEvent");
238   THnSparseF* hTrackMc = (THnSparseF*)fOutput->FindObject("hTrackMc");
239   THnSparseF* hTrackReco = (THnSparseF*)fOutput->FindObject("hTrackReco");
240   THnSparseF* hTrackMatch = (THnSparseF*)fOutput->FindObject("hTrackMatch");
241   THnSparseF* hPairMcPhi = (THnSparseF*)fOutput->FindObject("hPairMcPhi");
242   THnSparseF* hPairMcPhiCuts = (THnSparseF*)fOutput->FindObject("hPairMcPhiCuts");
243   THnSparseF* hPairMc = (THnSparseF*)fOutput->FindObject("hPairMc");
244   THnSparseF* hPairReco = (THnSparseF*)fOutput->FindObject("hPairReco");
245   THnSparseF* hPairMatch = (THnSparseF*)fOutput->FindObject("hPairMatch");
246   TH3F* hKcorr = (TH3F*)fOutput->FindObject("hKcorr");
247
248   TObjArray* kaonsPosMc = new TObjArray();
249   TObjArray* kaonsNegMc = new TObjArray();
250   TObjArray* kaonsPosData = new TObjArray();
251   TObjArray* kaonsNegData = new TObjArray();
252   TObjArray* kaonsPosGen = new TObjArray();
253   TObjArray* kaonsNegGen = new TObjArray();
254
255   Int_t nMc = 0; // number of charged MC particles which satisfy cuts (except PID)
256   Int_t nReco = 0; // number of reco tracks which satisfy cuts (except PID)
257
258   //MC Loop
259   TClonesArray *arrayMC = 0;
260   if (fIsMC)
261     {
262       arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
263       if (!arrayMC) {
264         AliFatal("Error: MC particles branch not found!\n");
265       }
266
267
268       Int_t nMC = arrayMC->GetEntries();
269       for (Int_t iMC = 0; iMC < nMC; iMC++)
270         {
271           AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
272
273           if(partMC->Pt()<fPtCut) continue;
274           if(partMC->Eta()>fTrackCuts->GetEtaMax() || partMC->Eta()<fTrackCuts->GetEtaMin() ) continue;
275
276           // PDG codes: pion(211), kaon(+/-321), proton(2212), phi(333)
277
278           // make pt spectrum of only MC phis that decay to k+k-
279           if(partMC->GetPdgCode()==333)
280             {
281               if(partMC->GetNDaughters()==2)
282                 {
283                   if(TMath::Abs(((AliAODMCParticle*) arrayMC->At(TMath::Abs(partMC->GetDaughter(0))))->GetPdgCode()) == 321 && TMath::Abs(((AliAODMCParticle*) arrayMC->At(TMath::Abs(partMC->GetDaughter(1))))->GetPdgCode()) == 321)
284                     {
285                       AliAODMCParticle *d1 = (AliAODMCParticle*) arrayMC->At(TMath::Abs(partMC->GetDaughter(0)));
286                       AliAODMCParticle *d2 = (AliAODMCParticle*) arrayMC->At(TMath::Abs(partMC->GetDaughter(1)));
287                       if(d1->Charge() > 0) hKcorr->Fill(partMC->Pt(),d1->Pt(),d2->Pt());
288                       else hKcorr->Fill(partMC->Pt(),d2->Pt(),d1->Pt());
289                       TLorentzVector *d3 = (TLorentzVector*)makePhi(d1,d2);
290                       // InvMass, pt, y, eta, phi, pair ID
291                       Double_t invMass = (d3->M2() > 0 ? sqrt(d3->M2()) : 0);
292                       if(invMass < 1.11)
293                         {
294                           Double_t varfill1[7] = {invMass, d3->Pt(), d3->Rapidity(), cent, d3->Eta(), (d3->Phi() > 0 ? d3->Phi() : d3->Phi()+2*TMath::Pi()), 0};
295                           hPairMcPhi->Fill(varfill1);
296                           if(d1->Pt()>fPtCut && d2->Pt()>fPtCut && d1->Eta()<fTrackCuts->GetEtaMax() && d1->Eta()>fTrackCuts->GetEtaMin() && d2->Eta()<fTrackCuts->GetEtaMax() && d2->Eta()>fTrackCuts->GetEtaMin())
297                             hPairMcPhiCuts->Fill(varfill1);
298                         }
299                       delete d3;
300                     }
301                 }
302             }
303
304           if(partMC->Charge()==0) continue; // remove neutrals
305           if(!partMC->IsPhysicalPrimary()) continue;  // need to remove all non-final state particles
306
307           nMc++;
308
309           Int_t mcID = fHelperPID->GetParticleSpecies(partMC);
310           if(fDoPID && (mcID>3 || mcID < -3)) continue;
311           // PID ID, pt, y, eta, phi
312           Double_t varfill2[6] = {static_cast<Double_t>((partMC->Charge() > 0 ? mcID+1 : -1*(mcID+1))), partMC->Pt(), partMC->Y(), cent, partMC->Eta(), partMC->Phi()};
313           hTrackMc->Fill(varfill2);
314
315           if(!fDoPID || mcID == 1)
316             {
317               if(partMC->Charge() > 0) kaonsPosMc->Add(partMC);
318               else if(partMC->Charge() < 0) kaonsNegMc->Add(partMC);
319             }
320           
321         } // end MC track loop
322       
323       UnlikeSign(kaonsPosMc,kaonsNegMc,hPairMc,cent);
324       LikeSign(kaonsPosMc,kaonsNegMc,hPairMc,cent);
325     }
326
327
328   //track loop
329   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
330     AliAODTrack* track = fAOD->GetTrack(iTracks);
331     if (!fTrackCuts->IsSelected(track,kTRUE)) continue;
332     if(track->Charge()==0) continue;
333     if(track->Pt()<fPtCut) continue;
334     if(track->Eta()>fTrackCuts->GetEtaMax() || track->Eta()<fTrackCuts->GetEtaMin()) continue;
335
336     nReco++;
337
338     Int_t dataID=fHelperPID->GetParticleSpecies(track,kTRUE);
339     if(fDoPID && (dataID>3 || dataID < -3)) continue;
340
341     // PID ID, pt, y, eta, phi
342     Double_t varfill3[6] = {static_cast<Double_t>((track->Charge() > 0 ? dataID+1 : -1*(dataID+1))), track->Pt(), track->Y(), cent, track->Eta(), track->Phi()};
343     hTrackReco->Fill(varfill3);
344
345     if(!fDoPID || dataID==1)
346       {
347         if(track->Charge() > 0) kaonsPosData->Add(track);
348         else if(track->Charge() < 0) kaonsNegData->Add(track);
349       }
350
351     if (fIsMC)
352       {
353         Int_t mcLabel = track->GetLabel();
354         AliAODMCParticle* tempMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(mcLabel));
355         if(!tempMC)
356           {
357             AliError("Cannot get MC particle");
358             continue; 
359           }
360         
361         Int_t genID = fHelperPID->GetParticleSpecies(tempMC);
362         if(fDoPID && (genID != dataID)) continue;
363         if(fDoPID && (genID>3 || genID < -3)) continue;
364
365         // PID ID, pt, y, eta, phi
366         Double_t varfill4[6] = {static_cast<Double_t>((track->Charge() > 0 ? genID+1 : -1*(genID+1))), track->Pt(), track->Y(), cent, track->Eta(), track->Phi()};
367         hTrackMatch->Fill(varfill4);
368
369         if(genID==1)
370           {
371             Int_t motherId = TMath::Abs(tempMC->GetMother());
372             
373             AliAODMCParticle* genMother = (AliAODMCParticle*)arrayMC->At(motherId);
374             if(!genMother) continue;
375             
376             if(genMother->GetPdgCode() != 333) continue;
377             
378             if(track->Charge() > 0)
379               {
380                 for(Int_t k = 0; k < kaonsNegGen->GetEntriesFast(); k++)
381                   {
382                     AliVParticle* tempGenMatch = (AliVParticle*)kaonsNegGen->UncheckedAt(k);
383                     AliAODMCParticle* tempMatchMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(tempGenMatch->GetLabel()));
384                     if(!tempMatchMC) continue;
385                     if(tempMatchMC->GetMother() != motherId) continue;
386                     
387                     TLorentzVector*c = (TLorentzVector*)makePhi(track,tempGenMatch);
388                     Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
389                     if(invMass < 1.11)
390                       {
391                         Double_t fillgen[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 0};
392                         hPairMatch->Fill(fillgen);
393                       }
394                     delete c;
395                   } 
396               }
397             else if(track->Charge() < 0)
398               {
399                 for(Int_t k = 0; k < kaonsPosGen->GetEntriesFast(); k++)
400                   {
401                     AliVParticle* tempGenMatch = (AliVParticle*)kaonsPosGen->UncheckedAt(k);
402                     AliAODMCParticle* tempMatchMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(tempGenMatch->GetLabel()));
403                     if(!tempMatchMC) continue;
404                     if(tempMatchMC->GetMother() != motherId) continue;
405                     
406                     TLorentzVector*c = (TLorentzVector*)makePhi(track,tempGenMatch);
407                     Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
408                     if(invMass < 1.11)
409                       {
410                         Double_t fillgen2[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 0};
411                         hPairMatch->Fill(fillgen2);
412                       }
413                     delete c;
414                   }
415               }
416             
417             if(track->Charge() > 0) kaonsPosGen->Add(track);
418             else if(track->Charge() < 0) kaonsNegGen->Add(track);
419           }
420       }
421         
422   } // end loop on tracks
423   
424   UnlikeSign(kaonsPosData,kaonsNegData,hPairReco,cent);
425   LikeSign(kaonsPosData,kaonsNegData,hPairReco,cent);
426   
427   Double_t varfill5[4] = {cent, fAOD->GetPrimaryVertex()->GetZ(), static_cast<Double_t>(nMc), static_cast<Double_t>(nReco)};
428   hEvent->Fill(varfill5);
429
430   kaonsPosMc->Clear();
431   kaonsNegMc->Clear();
432   kaonsPosData->Clear();
433   kaonsNegData->Clear();
434   kaonsPosGen->Clear();
435   kaonsNegGen->Clear();
436
437   if(kaonsPosMc) delete kaonsPosMc;
438   if(kaonsNegMc) delete kaonsNegMc;
439   if(kaonsPosData) delete kaonsPosData;
440   if(kaonsNegData) delete kaonsNegData;
441   if(kaonsPosGen) delete kaonsPosGen;
442   if(kaonsNegGen) delete kaonsNegGen;
443
444   PostData(1,fOutput);
445   PostData(2, fHelperPID);
446   PostData(3, fEventCuts);
447   PostData(4, fTrackCuts);
448   //Printf("............. end of Exec");
449   
450 }
451
452 //_________________________________________________________________
453 void   AliAnalysisTaskPhiEffMc::Terminate(Option_t *)
454 {
455   // Terminate analysis
456   //
457   fOutput = dynamic_cast<TList*>(GetOutputData(1));
458   if (!fOutput) {
459     printf("ERROR: fOutput not available\n");
460     return;
461   } 
462   
463     printf("AliAnalysisTaskPhiEffMc: Terminate() \n");
464   
465 }
466
467 //_________________________________________________________________
468 void AliAnalysisTaskPhiEffMc::UnlikeSign(TObjArray* kaonsPos, TObjArray* kaonsNeg, THnSparseF* h, Double_t cent)
469 {
470   Int_t countPos = kaonsPos->GetEntriesFast();
471   Int_t countNeg = kaonsNeg->GetEntriesFast();
472
473   if(countPos+countNeg < 2) return;
474
475   for(Int_t cp = 0; cp < countPos; cp++) // form k+k- pairs and compute invariant mass
476     {
477       AliVParticle* temp1 = (AliVParticle*)kaonsPos->UncheckedAt(cp);
478       for(Int_t cn = 0; cn < countNeg; cn++)
479         {
480           AliVParticle* temp2 = (AliVParticle*)kaonsNeg->UncheckedAt(cn);
481           TLorentzVector*c = (TLorentzVector*)makePhi(temp1,temp2);
482           // InvMass, pt, y, eta, phi, pair ID
483           Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
484           if(invMass < 1.11)
485             {
486               Double_t varfill[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 0};
487               h->Fill(varfill);
488             }
489           delete c;
490         }
491     }
492 }
493
494 void AliAnalysisTaskPhiEffMc::LikeSign(TObjArray* kaonsPos, TObjArray* kaonsNeg, THnSparseF* h, Double_t cent)
495 {
496   Int_t countPos = kaonsPos->GetEntriesFast();
497   Int_t countNeg = kaonsNeg->GetEntriesFast();
498
499   if(countPos < 2 && countNeg < 2) return;
500
501   for(Int_t cp = 0; cp < countPos; cp++) // form k+k+ pairs and compute invariant mass
502     {
503       AliVParticle* temp1 = (AliVParticle*)kaonsPos->UncheckedAt(cp);
504       for(Int_t cn = cp+1; cn < countPos; cn++)
505         {
506           AliVParticle* temp2 = (AliVParticle*)kaonsPos->UncheckedAt(cn);
507           TLorentzVector*c = (TLorentzVector*)makePhi(temp1,temp2);
508           // InvMass, pt, y, eta, phi, pair ID
509           Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
510           if(invMass < 1.11)
511             {
512               //cout << "y = " << 0.5*log((c->E()+c->Pz())/(c->E()-c->Pz())) << "   " << c->Rapidity() << endl;
513               //cout << "eta = " << -1.*log(tan(c->Theta()/2.)) << "   " << c->Eta() << endl;
514               Double_t varfill1[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 1};
515               h->Fill(varfill1);
516             }
517           delete c;
518         }
519     }
520   for(Int_t cp = 0; cp < countNeg; cp++) // form k-k- pairs and compute invariant mass
521     {
522       AliVParticle* temp1 = (AliVParticle*)kaonsNeg->UncheckedAt(cp);
523       for(Int_t cn = cp+1; cn < countNeg; cn++)
524         {
525           AliVParticle* temp2 = (AliVParticle*)kaonsNeg->UncheckedAt(cn);
526           TLorentzVector*c = (TLorentzVector*)makePhi(temp1,temp2);
527           Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
528           if(invMass < 1.11)
529             {
530               Double_t varfill2[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 2};
531               h->Fill(varfill2);
532             }
533           delete c;
534         } 
535     }
536 }
537
538 TLorentzVector* AliAnalysisTaskPhiEffMc::makePhi(AliVParticle* p1, AliVParticle* p2)
539 {
540   TLorentzVector* a = new TLorentzVector(p1->Px(),p1->Py(),p1->Pz(),sqrt(pow(p1->P(),2)+pow(4.93676999999999977e-01, 2)));
541   TLorentzVector* b = new TLorentzVector(p2->Px(),p2->Py(),p2->Pz(),sqrt(pow(p2->P(),2)+pow(4.93676999999999977e-01, 2)));
542   TLorentzVector* c = new TLorentzVector((*a)+(*b));
543   delete a;
544   delete b;
545   return c;
546 }
547