New task for filtering of AODs with Jpsi candidates (Carmelo)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEJPSItoEle.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 //                    Class AliAnalysisTaskSEJPSItoEle
17 // AliAnalysisTaskSE class to read both J/psi -> e+e- candidates and 
18 // like-sign pair candidates and to store them into a specific
19 // stand-alone AOD file for J/psi into dieletrons analysis. 
20 // The current task:
21 //
22 //  1) produces several histograms (including invariant mass distributions)
23 //     for both unlike sign (US) and like sign (LS) pairs.
24 //
25 //  2) selects only J/Psi to dielectrons candidates and copies it to a 
26 //     specific AOD file, namely AliAODjpsi.root. The references
27 //     to AliAODTrack objects are copied as well. 
28 //     The specific AliAODjpsi.root is thought  as the input file 
29 //     to the AliAnalysisBtoJPSItoEle.h class, which performs (locally)
30 //     the analysis of J/Psi from beauty hadrons decays. 
31 //     
32 //     Author: C.Di Giglio, carmelo.digiglio@ba.infn.it
33 //*************************************************************************
34 class AliAnalysisTaskSE;
35 class AliESDEvent;
36 class AliVEvent;
37 class iostream;
38 class TSystem;
39 class TROOT;
40 #include "TFile.h"
41 #include "TClonesArray.h"
42 #include "TDatabasePDG.h"
43 #include "TROOT.h"
44 #include "TCanvas.h"
45
46 #include "AliAODEvent.h"
47 #include "AliAODMCParticle.h"
48 #include "AliAODMCHeader.h"
49 #include "AliAODRecoDecayHF2Prong.h"
50 #include "AliAODHandler.h"
51 #include "AliAnalysisManager.h"
52 #include "AliAnalysisTaskSEJPSItoEle.h"
53 #include "AliAODTrack.h"
54 #include "AliExternalTrackParam.h"
55
56 ClassImp(AliAnalysisTaskSEJPSItoEle)
57
58 //______________________________________________________________________________
59 AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle() :
60 AliAnalysisTaskSE(),
61 fOutput(0),
62 fhDecayTimeMCjpsifromB(0),
63 fhDecayTime(0),                         
64 fhDecayTimeOut(0),
65 fhInvMass(0),                           
66 fhD0(0),                                
67 fhD0D0(0),                              
68 fhCosThetaStar(0),                      
69 fhCosThetaPointing(0),                  
70 fhDCA(0),
71 fhCtsVsD0D0(0),
72 fHistMassLS(0),
73 fHistCtsLS(0),
74 fHistCtsLSpos(0),
75 fHistCtsLSneg(0),
76 fHistCPtaLS(0),
77 fHistd0d0LS(0),
78 fHistDCALS(0),
79 fVHF(0),
80 fTotPosPairs(0),
81 fTotNegPairs(0),
82 fLsNormalization(1.),
83 fOkAODMC(kFALSE),
84 fOkLikeSign(kFALSE),
85 fVerticesHFTClArr(0),
86 fJpsiToEleTClArr(0),
87 fLikeSignTClArr(0),
88 fTracksTClArr(0),
89 fChain(0),
90 fOrigAOD(0),
91 fNewAOD(0)
92 {
93   // default constructor
94 }                 
95 //_________________________________________________________________________________________________
96 AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle(const char* name) :
97 AliAnalysisTaskSE(name),
98 fOutput(0),
99 fhDecayTimeMCjpsifromB(0),
100 fhDecayTime(0),
101 fhDecayTimeOut(0),
102 fhInvMass(0),                           
103 fhD0(0),                                
104 fhD0D0(0),                              
105 fhCosThetaStar(0),                      
106 fhCosThetaPointing(0),                  
107 fhDCA(0),
108 fhCtsVsD0D0(0),
109 fHistMassLS(0),
110 fHistCtsLS(0),
111 fHistCtsLSpos(0),
112 fHistCtsLSneg(0),
113 fHistCPtaLS(0),
114 fHistd0d0LS(0),
115 fHistDCALS(0),
116 fVHF(0),
117 fTotPosPairs(0),
118 fTotNegPairs(0),
119 fLsNormalization(1.),
120 fOkAODMC(kFALSE),
121 fOkLikeSign(kFALSE),
122 fVerticesHFTClArr(0),
123 fJpsiToEleTClArr(0),
124 fLikeSignTClArr(0),
125 fTracksTClArr(0),
126 fChain(0),
127 fOrigAOD(0),
128 fNewAOD(0)
129 {
130   // standard constructor
131
132   // Input slot #0 works with an Ntuple
133   DefineInput(0, TChain::Class());
134
135   // Output slot #0 writes into a TTree container
136   // Output slot #1 writes into a TList container
137   DefineOutput(0, TTree::Class());
138   DefineOutput(1,TList::Class());  //My private output
139 }
140 //_________________________________________________________________________________________________
141 AliAnalysisTaskSEJPSItoEle::~AliAnalysisTaskSEJPSItoEle()
142 {
143   // destructor
144
145   if (fOutput) {
146     delete fOutput;
147     fOutput = 0;
148   } 
149
150   if (fVHF) {
151     delete fVHF;
152     fVHF = 0;
153   }
154
155 }
156 //_________________________________________________________________________________________________
157 void AliAnalysisTaskSEJPSItoEle::Init()
158 {
159   // Initialization
160
161   if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle::Init() \n");
162
163   gROOT->LoadMacro("ConfigVertexingHF.C");
164
165   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
166   fVHF->PrintStatus();
167
168   return;
169 }
170 //_________________________________________________________________________________________________
171 void AliAnalysisTaskSEJPSItoEle::UserCreateOutputObjects() 
172 {
173   // Create the output container
174
175   if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle::UserCreateOutputObjects() \n");
176
177   if (!AODEvent()) {
178      Fatal("UserCreateOutputObjects", "This task needs an AOD handler");
179      return;
180   }   
181
182   if(!fNewAOD) fNewAOD = new AliAODEvent();
183   fNewAOD->CreateStdContent();
184
185   TString filename = "AliAOD.Jpsitoele.root";
186   if (!IsStandardAOD()) filename = "";
187
188   // Array of HF vertices  
189   fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0);
190   fVerticesHFTClArr->SetName("VerticesHF");
191   AddAODBranch("TClonesArray", &fVerticesHFTClArr);
192
193   // Array of J/psi -> e+e- candidates 
194   fJpsiToEleTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0);
195   fJpsiToEleTClArr->SetName("JpsiToEleEvents");
196   AddAODBranch("TClonesArray", &fJpsiToEleTClArr, filename);
197
198   // Array of like sign candidates 
199   fLikeSignTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0);
200   fLikeSignTClArr->SetName("LikeSignEvents");
201   AddAODBranch("TClonesArray", &fLikeSignTClArr, filename);
202
203   // Array of tracks from J/psi -> e+e- candidates
204   fTracksTClArr = new TClonesArray("AliAODTrack", 0);
205   fTracksTClArr->SetName("Tracks");
206   AddAODBranch("TClonesArray", &fTracksTClArr, filename);
207
208   fOutput = new TList();
209   fOutput->SetOwner();
210
211   if(fOkAODMC){
212   fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.);
213   fhDecayTimeMCjpsifromB->Sumw2();
214   fhDecayTimeMCjpsifromB->SetMinimum(0);
215   fOutput->Add(fhDecayTimeMCjpsifromB);
216   }
217
218   // Invariant mass
219   fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; M [GeV]; Entries",100,0.,3.2);
220   fhInvMass->Sumw2();
221   fhInvMass->SetMinimum(0);
222   fOutput->Add(fhInvMass);
223
224   fHistMassLS = new TH1F("fHistMassLS", "Like sign pairs invariant mass; M [GeV]; Entries",200,2.8,3.25);
225   fHistMassLS->Sumw2();
226   fHistMassLS->SetMinimum(0);
227   fOutput->Add(fHistMassLS);
228
229   // Pseudo proper decay time
230   fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.);
231   fhDecayTime->Sumw2();
232   fhDecayTime->SetMinimum(0);
233   fOutput->Add(fhDecayTime);
234  
235   // Pseudo proper decay time
236   fhDecayTimeOut = new TH1F("fhDecayTimeOut", "J/#Psi pseudo proper decay time (output standalone AOD); X [#mu m]; Entries",200,-2000.,2000.);
237   fhDecayTimeOut->Sumw2();
238   fhDecayTimeOut->SetMinimum(0);
239   fOutput->Add(fhDecayTimeOut);
240
241   // Dictance of closest approach
242   fhDCA = new TH1F("fhDCA", "J/#Psi distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
243   fhDCA->Sumw2();
244   fhDCA->SetMinimum(0);
245   fOutput->Add(fhDCA);
246
247   fHistDCALS = new TH1F("fHistDCALS", "Like sign pairs distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
248   fHistDCALS->Sumw2();
249   fHistDCALS->SetMinimum(0);
250   fOutput->Add(fHistDCALS);
251
252   // Impact parameter
253   fhD0 = new TH1F("fhD0", "Impact parameter; d0 [#mu m]; Entries",100,-5000.,5000.);
254   fhD0->Sumw2();
255   fhD0->SetMinimum(0);
256   fOutput->Add(fhD0);
257
258   // Product of impact parameters
259   fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.);
260   fhD0D0->Sumw2();
261   fhD0D0->SetMinimum(0);
262   fOutput->Add(fhD0D0);
263
264   fHistd0d0LS = new TH1F("fHistd0d0LS", "Like sign pairs product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.);
265   fHistd0d0LS->Sumw2();
266   fHistd0d0LS->SetMinimum(0);
267   fOutput->Add(fHistd0d0LS);
268
269   // Cosine of the decay angle
270   fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2);
271   fhCosThetaStar->Sumw2();
272   fhCosThetaStar->SetMinimum(0);
273   fOutput->Add(fhCosThetaStar);
274
275   fHistCtsLS = new TH1F("fHistCtsLS", "Like sign pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
276   fHistCtsLS->Sumw2();
277   fHistCtsLS->SetMinimum(0);
278   fOutput->Add(fHistCtsLS);
279
280   fHistCtsLSpos = new TH1F("fHistCtsLSpos", "Like sign ++ pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
281   fHistCtsLSpos->Sumw2();
282   fHistCtsLSpos->SetMinimum(0);
283   fOutput->Add(fHistCtsLSpos);
284
285   fHistCtsLSneg = new TH1F("fHistCtsLSneg", "Like sign -- pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
286   fHistCtsLSneg->Sumw2();
287   fHistCtsLSneg->SetMinimum(0);
288   fOutput->Add(fHistCtsLSneg);
289
290   // Cosine of pointing angle
291   fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1);
292   fhCosThetaPointing->Sumw2();
293   fhCosThetaPointing->SetMinimum(0);
294   fOutput->Add(fhCosThetaPointing);
295
296   fHistCPtaLS = new TH1F("fHistCPtaLS", "Like sign pairs cosine of pointing angle; Cos#Theta_{point}; Entries",200,-1.,1.);
297   fHistCPtaLS->Sumw2();
298   fHistCPtaLS->SetMinimum(0);
299   fOutput->Add(fHistCPtaLS);
300
301   // CosThetaStar vs. d0xd0
302   fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000);
303   fhCtsVsD0D0->Sumw2();
304   fhCtsVsD0D0->SetMinimum(0);
305   fOutput->Add(fhCtsVsD0D0);
306
307   return;
308 }
309 //_________________________________________________________________________________________________
310 void AliAnalysisTaskSEJPSItoEle::UserExec(Option_t */*option*/)
311 {
312   // Execute analysis for current event:
313
314   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
315
316   TClonesArray *arrayJPSItoEle = 0;
317   TClonesArray *arrayLikeSign = 0;
318   TClonesArray *arrayTracks = 0;
319
320   if(!aod && AODEvent() && IsStandardAOD()) {
321     // In case there is an AOD handler writing a standard AOD, use the AOD 
322     // event in memory rather than the input (ESD) event.    
323     aod = dynamic_cast<AliAODEvent*> (AODEvent());
324     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
325     // have to taken from the AOD event hold by the AliAODExtension
326     AliAODHandler* aodHandler = (AliAODHandler*)
327       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
328
329     if(aodHandler->GetExtensions()) {
330       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
331       AliAODEvent *aodFromExt = ext->GetAOD();
332
333       // load tracks from AOD event   
334       arrayTracks=(TClonesArray*)aod->GetList()->FindObject("tracks");
335       //Int_t totTracks = arrayTracks->GetEntriesFast();
336
337       // load Jpsi candidates from AOD friend  
338       arrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle");
339       //Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast();
340
341       // load like sign candidates from AOD friend
342       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign2Prong");
343       //Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast();
344
345     }
346
347   } else {
348     // load Jpsi candidates                                                   
349     arrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle");
350     //Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast();
351     // load like sign candidates
352     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong");
353     //Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast();
354   }
355
356   fOrigAOD = aod; // copy pointer to the current AliAODEvent in the data member fOrigAOD
357   if (!aod) return;
358   Int_t nTracks = fOrigAOD->GetNumberOfTracks();
359   printf("+++\n+++ Number of tracks in Event---> %d\n+++\n",nTracks);
360
361   if(!arrayJPSItoEle) {
362     printf("AliAnalysisTaskSEJPSItoEle::UserExec: JPSItoEle branch not found!\n");
363     return;
364   }
365   if(!arrayLikeSign) {
366     printf("AliAnalysisTaskSEJPSItoEle::UserExec: LikeSign2Prong branch not found!\n");
367     return;
368   }
369
370   // load MC particles and read MC info (for sim only)
371   TClonesArray* mcArray =  dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
372   if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
373      //AliInfo(Form("+++\n+++ MC particles found in mcArray ---> %d \n+++\n",mcArray->GetEntriesFast()));
374   if(fOkAODMC) ReadAODMCInfo(aod,arrayJPSItoEle);
375
376   // retrieve AOD primary vertex
377   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
378
379   Int_t iOutVerticesHF = 0, iOutJPSItoEle = 0, iOutLikeSign = 0, iOutTracks = 0; 
380
381   fVerticesHFTClArr->Delete();
382   iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast();
383   TClonesArray &verticesHFRef = *fVerticesHFTClArr;
384
385   fJpsiToEleTClArr->Delete();
386   iOutJPSItoEle = fJpsiToEleTClArr->GetEntriesFast();
387   TClonesArray &arrayJpsiToEleRef = *fJpsiToEleTClArr;
388
389   fLikeSignTClArr->Delete();
390   iOutLikeSign = fLikeSignTClArr->GetEntriesFast();
391   TClonesArray &arrayLikeSignRef = *fLikeSignTClArr;
392
393   fTracksTClArr->Delete();
394   iOutTracks = fTracksTClArr->GetEntriesFast();
395   //TClonesArray &arrayTrackRef = *fTracksTClArr;
396
397   //
398   // LOOP OVER LIKE SIGN CANDIDATES
399   //
400
401   // Access to the new AOD container of tracks
402   Int_t trkIDtoEntry[100000];
403   for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
404     AliAODTrack *track = aod->GetTrack(it);
405     trkIDtoEntry[track->GetID()]=it;
406   }
407
408   Int_t nPosPairs=0,nNegPairs=0;
409   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
410   if(fDebug>1) printf("+++\n+++ Number of total like sign pairs (before cuts)---> %d \n+++\n", nLikeSign);
411
412   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
413     AliAODRecoDecayHF2Prong *dlsin = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
414     Bool_t unsetvtx=kFALSE;
415     if(!dlsin->GetOwnPrimaryVtx()) {
416         dlsin->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
417         unsetvtx=kTRUE;
418     }
419     Int_t okBtoJPSIls=0;
420     if(dlsin->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSIls)) {
421        AliAODTrack *trk0 = (AliAODTrack*)dlsin->GetDaughter(0);
422        fHistMassLS->Fill(dlsin->InvMassJPSIee());
423        fHistCPtaLS->Fill(dlsin->CosPointingAngle());
424        fHistd0d0LS->Fill(1e8*dlsin->Prodd0d0());
425        fHistDCALS->Fill(100*dlsin->GetDCA());
426        if(!trk0) {
427           trk0=aod->GetTrack(trkIDtoEntry[dlsin->GetProngID(0)]);
428           printf("references to standard AOD not available \n");
429        }
430        if((trk0->Charge())==1) {
431           nPosPairs++;
432           fHistCtsLS->Fill(dlsin->CosThetaStar(0,443,11,11));
433           fHistCtsLSpos->Fill(dlsin->CosThetaStar(0,443,11,11));
434         } else {
435           nNegPairs++;
436           fHistCtsLS->Fill(dlsin->CosThetaStarJPSI());
437           fHistCtsLSneg->Fill(dlsin->CosThetaStarJPSI());
438         }
439        PostData(1,fOutput);
440      
441        // Clone like sign candidate for output AOD
442        new(arrayLikeSignRef[iOutLikeSign++]) AliAODRecoDecayHF2Prong(*dlsin);
443
444     }
445     if(unsetvtx) dlsin->UnsetOwnPrimaryVtx();
446   }
447
448   if(fDebug>1) printf("+++\n+++ N. of positive pairs passing cuts in Event ----- %d \n+++\n", nPosPairs);
449   if(fDebug>1) printf("+++\n+++ N. of negative pairs passing cuts in Event ----- %d \n+++\n", nNegPairs);
450
451   fTotPosPairs += nPosPairs;
452   fTotNegPairs += nNegPairs;
453
454   //
455   // LOOP OVER J/psi->e+e- CANDIDATES
456   //
457
458   Int_t nInJPSItoEle = arrayJPSItoEle->GetEntriesFast();
459   if(fDebug>1) printf("+++\n+++ Number of total like JPSI -> ee candidates (before cuts)---> %d \n+++\n", nInJPSItoEle);
460
461   //totJPSIin +=  nInJPSItoEle;
462
463   for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
464
465     AliAODRecoDecayHF2Prong *dIn = (AliAODRecoDecayHF2Prong*)arrayJPSItoEle->UncheckedAt(iJPSItoEle);
466     Int_t mcLabel = 0;
467     if(fOkAODMC) mcLabel = dIn->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
468
469     //Secondary vertex
470     Double_t vtxSec[3] = {0.,0.,0.};
471     Double_t vtxPrim[3] = {0.,0.,0.};
472     Double_t vtxSecOut[3] = {0.,0.,0.};
473     Double_t vtxPrimOut[3] = {0.,0.,0.};
474     dIn->GetSecondaryVtx(vtxSec); 
475     Bool_t unsetvtx=kFALSE;
476     if(!dIn->GetOwnPrimaryVtx()) {
477       dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
478       unsetvtx=kTRUE;
479     }
480
481     Int_t okBtoJPSI=0;
482     if(dIn->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSI)) {
483       if ( fOkAODMC && mcLabel == -1){AliDebug(2,"No MC particle found");} else {
484
485          fhInvMass->Fill(dIn->InvMassJPSIee()); 
486          fhD0->Fill(10000*dIn->ImpParXY());
487          fhD0D0->Fill(1e8*dIn->Prodd0d0());
488          fhCosThetaStar->Fill(dIn->CosThetaStarJPSI());      
489          fhCtsVsD0D0->Fill(dIn->CosThetaStarJPSI(),1e8*dIn->Prodd0d0());
490          fhCosThetaPointing->Fill(dIn->CosPointingAngle());
491          fhDCA->Fill(100*dIn->GetDCA());
492
493          // compute X variable
494          AliAODVertex* primVertex = dIn->GetOwnPrimaryVtx();
495          vtxPrim[0] = primVertex->GetX();
496          vtxPrim[1] = primVertex->GetY();
497          vtxPrim[2] = primVertex->GetZ();
498          Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dIn->Px()) + (vtxSec[1]-vtxPrim[1])*(dIn->Py()))/dIn->Pt();
499          Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dIn->Pt();
500          fhDecayTime->Fill(10000*pseudoProperDecayTime);
501   
502          // clone candidate for output AOD
503          AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++]) 
504            AliAODVertex(*(dIn->GetSecondaryVtx()));
505          AliAODRecoDecayHF2Prong* dOut = new(arrayJpsiToEleRef[iOutJPSItoEle++]) 
506            AliAODRecoDecayHF2Prong(*dIn); // copy constructor used
507          dOut->SetSecondaryVtx(v);
508          dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone()));
509          v->SetParent(dOut);
510
511          // compute X variable from the cloned candidate in the stand-alone AOD file
512          dOut->GetSecondaryVtx(vtxSecOut);
513          AliAODVertex* primVertexOut = dOut->GetOwnPrimaryVtx();
514          vtxPrimOut[0] = primVertexOut->GetX();
515          vtxPrimOut[1] = primVertexOut->GetY();
516          vtxPrimOut[2] = primVertexOut->GetZ();
517          Double_t lxyOut = ((vtxSecOut[0]-vtxPrimOut[0])*(dOut->Px()) + (vtxSecOut[1]-vtxPrimOut[1])*(dOut->Py()))/dOut->Pt();
518          Double_t pseudoProperDecayTimeOut = lxyOut*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dOut->Pt();
519          fhDecayTimeOut->Fill(10000*pseudoProperDecayTimeOut);
520
521          // retrieve tracks using the clone J/psi-->e+e- candidate stored in the output AOD
522          //AliAODTrack* daugh0Out = (AliAODTrack*)dOut->GetSecondaryVtx()->GetDaughter(0);
523          //AliAODTrack* daugh1Out= (AliAODTrack*)dOut->GetSecondaryVtx()->GetDaughter(1);
524          //printf("pt of positive track: %f\n",daugh0Out->Pt());
525          //printf("pt of negative track: %f\n",daugh1Out->Pt());
526
527         }
528
529      } // end of JPSItoEle candidates selection according to cuts
530
531     if(unsetvtx) dIn->UnsetOwnPrimaryVtx();
532
533     PostData(0,fOutput);
534
535  }// end loop on JPSI to ele candidates
536
537  printf("+++\n+++ Number of selected J/psi->e+e-: %d\n+++\n",iOutJPSItoEle);
538  
539  //totJPSIout += iOutJPSItoEle;
540
541  return;
542
543 }
544 //_________________________________________________________________________________________________
545 void AliAnalysisTaskSEJPSItoEle::Terminate(Option_t */*option*/)
546 {
547   //
548   // Terminate analysis
549   //
550
551   if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle: Terminate() \n");
552
553   fOutput = dynamic_cast<TList*> (GetOutputData(1));
554   if (!fOutput) {
555     printf("ERROR: fOutput not available\n");
556     return;
557   }
558
559   fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs);
560
561   if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
562   fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
563   fhDecayTimeOut = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeOut"));
564   fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
565   fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
566   fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
567   fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
568   fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
569   fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
570   fHistMassLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassLS"));
571   fHistCtsLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLS"));
572   fHistCtsLSpos = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSpos"));
573   fHistCtsLSneg = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSneg"));
574   fHistCPtaLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaLS"));
575   fHistd0d0LS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0LS"));
576   fhDCA = dynamic_cast<TH1F*>(fOutput->FindObject("fhDCA"));
577   fHistDCALS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistDCALS"));
578
579   if(fLsNormalization>0.) {
580     fHistMassLS->Scale((1/fLsNormalization)*fHistMassLS->GetEntries());
581     fHistCtsLS->Scale((1/fLsNormalization)*fHistCtsLS->GetEntries());
582     fHistCtsLSpos->Scale((1/fLsNormalization)*fHistCtsLSpos->GetEntries());
583     fHistCtsLSneg->Scale((1/fLsNormalization)*fHistCtsLSneg->GetEntries());
584     fHistCPtaLS->Scale((1/fLsNormalization)*fHistCPtaLS->GetEntries());
585     fHistd0d0LS->Scale((1/fLsNormalization)*fHistd0d0LS->GetEntries());
586     fHistDCALS->Scale((1/fLsNormalization)*fHistDCALS->GetEntries());
587   }
588
589   //printf("Tot JPSI in %d\n", totJPSIin);
590   //printf("Tot JPSI out %d\n", totJPSIout);
591
592   return;
593 }
594 //_________________________________________________________________________________________________
595 void AliAnalysisTaskSEJPSItoEle::SetPtCuts(const Double_t ptCuts[2])
596 {
597   //
598   // apply Pt cuts (lower cut, upper cut)
599   //
600   for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i];
601 }
602 //_________________________________________________________________________________________________
603 void AliAnalysisTaskSEJPSItoEle::SetCutsJPSI(const Double_t cuts[9])
604 {
605   //
606   // apply JPSI cuts 
607   //
608   for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
609 }
610 //_________________________________________________________________________________________________
611 void AliAnalysisTaskSEJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray) 
612 {
613   //
614   // Reads MC information if needed (for example in case of parametrized PID)
615   //
616
617   // load MC particles
618   TClonesArray* mcArray =
619      dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
620   if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
621      AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
622
623   // load MC header 
624   AliAODMCHeader* mcHeader =
625     (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
626   if(!mcHeader){
627     printf("AliAnalysisTaskSEJPSItoEle::UserExec: MC AOD header branch not found!\n");
628   }
629
630   AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
631
632   Double_t mcPrimVtx[3];
633
634   // get MC primary Vtx
635   mcHeader->GetVertex(mcPrimVtx);
636
637   Int_t nInCandidates = inputArray->GetEntriesFast();
638   printf("Number of Candidates for MC analysis: %d\n",nInCandidates);
639
640   Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother;
641   Int_t labJPSIdaugh0=0;
642   Int_t labJPSIdaugh1=0;
643
644   for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) {
645
646     AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate);
647     Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
648
649     Double_t vtxPrim[3] = {0.,0.,0.};
650     Double_t vtxSec[3] = {0.,0.,0.};
651     dd->GetSecondaryVtx(vtxSec);
652     Bool_t unsetvtx=kFALSE;
653     if(!dd->GetOwnPrimaryVtx()) {
654       dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
655       unsetvtx=kTRUE;
656     }
657
658     if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only
659
660       AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0);
661       AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1);
662       lab0 = trk0->GetLabel();
663       lab1 = trk1->GetLabel();
664
665       AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0);
666       if(!part0) {
667         printf("no MC particle\n");
668         continue;
669       }
670
671       while(part0->GetMother()>=0) {
672        labMother=part0->GetMother();
673        part0 = (AliAODMCParticle*)mcArray->At(labMother);
674        if(!part0) {
675          printf("no MC mother particle\n");
676          break;
677        }
678        pdgMother = TMath::Abs(part0->GetPdgCode());
679        if(pdgMother==443) {//this for JPSI
680          labJPSIdaugh0=labMother;
681          break;
682        }
683       }
684
685       AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1);
686       if(!part1) {
687         printf("no MC particle\n");
688         continue;
689       }
690
691       while(part1->GetMother()>=0) {
692        labMother=part1->GetMother();
693        part1 = (AliAODMCParticle*)mcArray->At(labMother);
694        if(!part1) {
695          printf("no MC mother particle\n");
696          break;
697        }
698        pdgMother = TMath::Abs(part1->GetPdgCode());
699        if(pdgMother==443) {//this for JPSI
700          labJPSIdaugh1=labMother;
701          break;
702        }
703       }
704
705       if (mcLabelSec == -1)
706          {
707           AliDebug(2,"No MC particle found");
708          }
709       else {
710
711       if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
712           AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
713           pdgJPSI = partJPSI->GetPdgCode();
714           Int_t pdaugh0 = partJPSI->GetDaughter(0);
715           Int_t pdaugh1 = partJPSI->GetDaughter(1);
716           AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
717           AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
718           pdg0 = partDaugh0->GetPdgCode();
719           pdg1 = partDaugh1->GetPdgCode();
720         if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
721            labJPSIMother = partJPSI->GetMother();
722            AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
723            pdgJPSIMother = partJPSIMother->GetPdgCode();
724
725         // chech if JPSI comes from B 
726            if( pdgJPSIMother==511   || pdgJPSIMother==521   ||
727                pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
728                pdgJPSIMother==513   || pdgJPSIMother==523   ||
729                pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
730                pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
731                pdgJPSIMother==515   || pdgJPSIMother==525   ||
732                pdgJPSIMother==531   || pdgJPSIMother==10531 ||
733                pdgJPSIMother==533   || pdgJPSIMother==10533 ||
734                pdgJPSIMother==20533 || pdgJPSIMother==535   ||
735                pdgJPSIMother==541   || pdgJPSIMother==10541 ||
736                pdgJPSIMother==543   || pdgJPSIMother==10543 ||
737                pdgJPSIMother==20543 || pdgJPSIMother==545)
738                { //this is for MC JPSI -> ee from B-hadrons
739
740                   fhInvMass->Fill(dd->InvMassJPSIee()); 
741                   fhD0->Fill(10000*dd->ImpParXY());
742                   fhD0D0->Fill(1e8*dd->Prodd0d0());
743                   fhCosThetaStar->Fill(dd->CosThetaStarJPSI());      
744                   fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0());
745                   fhCosThetaPointing->Fill(dd->CosPointingAngle());
746   
747                   // compute X variable
748                   AliAODVertex* primVertex = dd->GetOwnPrimaryVtx();
749                   vtxPrim[0] = primVertex->GetX();
750                   vtxPrim[1] = primVertex->GetY();
751                   vtxPrim[2] = primVertex->GetZ();
752                   Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt();
753                   Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt();
754   
755                   fhDecayTime->Fill(10000*pseudoProperDecayTime);
756                   // Post the data already here
757                   PostData(1,fOutput);
758
759                   Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt();
760                   Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt();
761                   fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1);
762
763                   // Post the data already here
764                   PostData(1,fOutput);
765
766             } //this is for MC JPSI -> ee from B-hadrons
767           } //this is for MC JPSI -> ee
768         }
769       } // end of check if JPSI is signal
770     }
771   }
772 }