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 }