4ea17f10b8d329614baf1493c73ee5e06074c7e3
[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 = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
331     if(!track) AliFatal("Not a standard AOD");
332     if (!fTrackCuts->IsSelected(track,kTRUE)) continue;
333     if(track->Charge()==0) continue;
334     if(track->Pt()<fPtCut) continue;
335     if(track->Eta()>fTrackCuts->GetEtaMax() || track->Eta()<fTrackCuts->GetEtaMin()) continue;
336
337     nReco++;
338
339     Int_t dataID=fHelperPID->GetParticleSpecies(track,kTRUE);
340     if(fDoPID && (dataID>3 || dataID < -3)) continue;
341
342     // PID ID, pt, y, eta, phi
343     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()};
344     hTrackReco->Fill(varfill3);
345
346     if(!fDoPID || dataID==1)
347       {
348         if(track->Charge() > 0) kaonsPosData->Add(track);
349         else if(track->Charge() < 0) kaonsNegData->Add(track);
350       }
351
352     if (fIsMC)
353       {
354         Int_t mcLabel = track->GetLabel();
355         AliAODMCParticle* tempMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(mcLabel));
356         if(!tempMC)
357           {
358             AliError("Cannot get MC particle");
359             continue; 
360           }
361         
362         Int_t genID = fHelperPID->GetParticleSpecies(tempMC);
363         if(fDoPID && (genID != dataID)) continue;
364         if(fDoPID && (genID>3 || genID < -3)) continue;
365
366         // PID ID, pt, y, eta, phi
367         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()};
368         hTrackMatch->Fill(varfill4);
369
370         if(genID==1)
371           {
372             Int_t motherId = TMath::Abs(tempMC->GetMother());
373             
374             AliAODMCParticle* genMother = (AliAODMCParticle*)arrayMC->At(motherId);
375             if(!genMother) continue;
376             
377             if(genMother->GetPdgCode() != 333) continue;
378             
379             if(track->Charge() > 0)
380               {
381                 for(Int_t k = 0; k < kaonsNegGen->GetEntriesFast(); k++)
382                   {
383                     AliVParticle* tempGenMatch = (AliVParticle*)kaonsNegGen->UncheckedAt(k);
384                     AliAODMCParticle* tempMatchMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(tempGenMatch->GetLabel()));
385                     if(!tempMatchMC) continue;
386                     if(tempMatchMC->GetMother() != motherId) continue;
387                     
388                     TLorentzVector*c = (TLorentzVector*)makePhi(track,tempGenMatch);
389                     Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
390                     if(invMass < 1.11)
391                       {
392                         Double_t fillgen[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 0};
393                         hPairMatch->Fill(fillgen);
394                       }
395                     delete c;
396                   } 
397               }
398             else if(track->Charge() < 0)
399               {
400                 for(Int_t k = 0; k < kaonsPosGen->GetEntriesFast(); k++)
401                   {
402                     AliVParticle* tempGenMatch = (AliVParticle*)kaonsPosGen->UncheckedAt(k);
403                     AliAODMCParticle* tempMatchMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(tempGenMatch->GetLabel()));
404                     if(!tempMatchMC) continue;
405                     if(tempMatchMC->GetMother() != motherId) continue;
406                     
407                     TLorentzVector*c = (TLorentzVector*)makePhi(track,tempGenMatch);
408                     Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
409                     if(invMass < 1.11)
410                       {
411                         Double_t fillgen2[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 0};
412                         hPairMatch->Fill(fillgen2);
413                       }
414                     delete c;
415                   }
416               }
417             
418             if(track->Charge() > 0) kaonsPosGen->Add(track);
419             else if(track->Charge() < 0) kaonsNegGen->Add(track);
420           }
421       }
422         
423   } // end loop on tracks
424   
425   UnlikeSign(kaonsPosData,kaonsNegData,hPairReco,cent);
426   LikeSign(kaonsPosData,kaonsNegData,hPairReco,cent);
427   
428   Double_t varfill5[4] = {cent, fAOD->GetPrimaryVertex()->GetZ(), static_cast<Double_t>(nMc), static_cast<Double_t>(nReco)};
429   hEvent->Fill(varfill5);
430
431   kaonsPosMc->Clear();
432   kaonsNegMc->Clear();
433   kaonsPosData->Clear();
434   kaonsNegData->Clear();
435   kaonsPosGen->Clear();
436   kaonsNegGen->Clear();
437
438   if(kaonsPosMc) delete kaonsPosMc;
439   if(kaonsNegMc) delete kaonsNegMc;
440   if(kaonsPosData) delete kaonsPosData;
441   if(kaonsNegData) delete kaonsNegData;
442   if(kaonsPosGen) delete kaonsPosGen;
443   if(kaonsNegGen) delete kaonsNegGen;
444
445   PostData(1,fOutput);
446   PostData(2, fHelperPID);
447   PostData(3, fEventCuts);
448   PostData(4, fTrackCuts);
449   //Printf("............. end of Exec");
450   
451 }
452
453 //_________________________________________________________________
454 void   AliAnalysisTaskPhiEffMc::Terminate(Option_t *)
455 {
456   // Terminate analysis
457   //
458   fOutput = dynamic_cast<TList*>(GetOutputData(1));
459   if (!fOutput) {
460     printf("ERROR: fOutput not available\n");
461     return;
462   } 
463   
464     printf("AliAnalysisTaskPhiEffMc: Terminate() \n");
465   
466 }
467
468 //_________________________________________________________________
469 void AliAnalysisTaskPhiEffMc::UnlikeSign(TObjArray* kaonsPos, TObjArray* kaonsNeg, THnSparseF* h, Double_t cent)
470 {
471   Int_t countPos = kaonsPos->GetEntriesFast();
472   Int_t countNeg = kaonsNeg->GetEntriesFast();
473
474   if(countPos+countNeg < 2) return;
475
476   for(Int_t cp = 0; cp < countPos; cp++) // form k+k- pairs and compute invariant mass
477     {
478       AliVParticle* temp1 = (AliVParticle*)kaonsPos->UncheckedAt(cp);
479       for(Int_t cn = 0; cn < countNeg; cn++)
480         {
481           AliVParticle* temp2 = (AliVParticle*)kaonsNeg->UncheckedAt(cn);
482           TLorentzVector*c = (TLorentzVector*)makePhi(temp1,temp2);
483           // InvMass, pt, y, eta, phi, pair ID
484           Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
485           if(invMass < 1.11)
486             {
487               Double_t varfill[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 0};
488               h->Fill(varfill);
489             }
490           delete c;
491         }
492     }
493 }
494
495 void AliAnalysisTaskPhiEffMc::LikeSign(TObjArray* kaonsPos, TObjArray* kaonsNeg, THnSparseF* h, Double_t cent)
496 {
497   Int_t countPos = kaonsPos->GetEntriesFast();
498   Int_t countNeg = kaonsNeg->GetEntriesFast();
499
500   if(countPos < 2 && countNeg < 2) return;
501
502   for(Int_t cp = 0; cp < countPos; cp++) // form k+k+ pairs and compute invariant mass
503     {
504       AliVParticle* temp1 = (AliVParticle*)kaonsPos->UncheckedAt(cp);
505       for(Int_t cn = cp+1; cn < countPos; cn++)
506         {
507           AliVParticle* temp2 = (AliVParticle*)kaonsPos->UncheckedAt(cn);
508           TLorentzVector*c = (TLorentzVector*)makePhi(temp1,temp2);
509           // InvMass, pt, y, eta, phi, pair ID
510           Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
511           if(invMass < 1.11)
512             {
513               //cout << "y = " << 0.5*log((c->E()+c->Pz())/(c->E()-c->Pz())) << "   " << c->Rapidity() << endl;
514               //cout << "eta = " << -1.*log(tan(c->Theta()/2.)) << "   " << c->Eta() << endl;
515               Double_t varfill1[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 1};
516               h->Fill(varfill1);
517             }
518           delete c;
519         }
520     }
521   for(Int_t cp = 0; cp < countNeg; cp++) // form k-k- pairs and compute invariant mass
522     {
523       AliVParticle* temp1 = (AliVParticle*)kaonsNeg->UncheckedAt(cp);
524       for(Int_t cn = cp+1; cn < countNeg; cn++)
525         {
526           AliVParticle* temp2 = (AliVParticle*)kaonsNeg->UncheckedAt(cn);
527           TLorentzVector*c = (TLorentzVector*)makePhi(temp1,temp2);
528           Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0);
529           if(invMass < 1.11)
530             {
531               Double_t varfill2[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 2};
532               h->Fill(varfill2);
533             }
534           delete c;
535         } 
536     }
537 }
538
539 TLorentzVector* AliAnalysisTaskPhiEffMc::makePhi(AliVParticle* p1, AliVParticle* p2)
540 {
541   TLorentzVector* a = new TLorentzVector(p1->Px(),p1->Py(),p1->Pz(),sqrt(pow(p1->P(),2)+pow(4.93676999999999977e-01, 2)));
542   TLorentzVector* b = new TLorentzVector(p2->Px(),p2->Py(),p2->Pz(),sqrt(pow(p2->P(),2)+pow(4.93676999999999977e-01, 2)));
543   TLorentzVector* c = new TLorentzVector((*a)+(*b));
544   delete a;
545   delete b;
546   return c;
547 }
548