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