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   // fix for temporary bug in ESDfilter
 390   // the AODs with null vertex pointer didn't pass the PhysSel
 391   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 392
 393   // load MC particles and read MC info (for sim only)
 394   TClonesArray* mcArray=0;
 395   if(fOkAODMC){
 396      mcArray =  dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 397      if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 398      ReadAODMCInfo(aod,arrayJPSItoEle);
 399      }
 400
 401   // retrieve AOD primary vertex
 402   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 403
 404   Int_t iOutVerticesHF = 0, iOutJPSItoEle = 0, iOutLikeSign = 0, iOutTracks = 0;
 405
 406   fVerticesHFTClArr->Delete();
 407   iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast();
 408   TClonesArray &verticesHFRef = *fVerticesHFTClArr;
 409
 410   fJpsiToEleTClArr->Delete();
 411   iOutJPSItoEle = fJpsiToEleTClArr->GetEntriesFast();
 412   TClonesArray &arrayJpsiToEleRef = *fJpsiToEleTClArr;
 413
 414   fLikeSignTClArr->Delete();
 415   iOutLikeSign = fLikeSignTClArr->GetEntriesFast();
 416   TClonesArray &arrayLikeSignRef = *fLikeSignTClArr;
 417
 418   fTracksTClArr->Delete();
 419   iOutTracks = fTracksTClArr->GetEntriesFast();
 420   //TClonesArray &arrayTrackRef = *fTracksTClArr;
 421
 422   //
 423   // LOOP OVER LIKE SIGN CANDIDATES
 424   //
 425
 426   // Access to the new AOD container of tracks
 427
 428   Int_t trkIDtoEntry[100000];
 429   AliAODTrack* trackin=0;
 430   for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
 431     trackin = aod->GetTrack(it);
 432     trkIDtoEntry[trackin->GetID()]=it;
 433     }
 434 /*
 435   for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
 436     vtx1->AddDaughter(trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack((*(aod->GetTrack(trkIDtoEntry[it])))));
 437     trackin = (AliAODTrack*)arrayTracks->UncheckedAt(ntracks);
 438     printf("trackin ==== %d \n", trackin);
 439     AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin);
 440     }
 441   const AliAODTrack* trackin;
 442   Int_t numTracks = arrayTracks->GetEntriesFast();
 443   for(Int_t ntracks=0; ntracks<numTracks; ntracks++){
 444   trackin = (AliAODTrack*)arrayTracks->UncheckedAt(ntracks);
 445   printf("trackin ==== %d \n", trackin);
 446     AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin);
 447   }
 448 */
 449
 450   Int_t nPosPairs=0,nNegPairs=0;
 451   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
 452   AliAODRecoDecayHF2Prong* dlsout = 0;
 453   //if(fDebug>1) printf("+++\n+++ Number of total like sign pairs (before cuts)---> %d \n+++\n", nLikeSign);
 454
 455   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
 456     AliAODRecoDecayHF2Prong *dlsin = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
 457     Bool_t unsetvtx=kFALSE;
 458     if(!dlsin->GetOwnPrimaryVtx()) {
 459         dlsin->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 460         unsetvtx=kTRUE;
 461     }
 462     //Int_t okBtoJPSIls=0;
 463     //if(dlsin->Pt() < fPtCuts[1] && dlsin->Pt() > fPtCuts[0]){ // apply pt cut only
 464     //if(dlsin->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSIls)) {
 465        AliAODTrack *trk0 = (AliAODTrack*)dlsin->GetDaughter(0);
 466        fHistMassLS->Fill(dlsin->InvMassJPSIee());
 467        fHistCPtaLS->Fill(dlsin->CosPointingAngle());
 468        fHistd0d0LS->Fill(1e8*dlsin->Prodd0d0());
 469        fHistDCALS->Fill(100*dlsin->GetDCA());
 470        if(!trk0) {
 471           trk0=aod->GetTrack(trkIDtoEntry[dlsin->GetProngID(0)]);
 472           printf("references to standard AOD not available \n");
 473        }
 474        if((trk0->Charge())==1) {
 475           nPosPairs++;
 476           fHistCtsLS->Fill(dlsin->CosThetaStar(0,443,11,11));
 477           fHistCtsLSpos->Fill(dlsin->CosThetaStar(0,443,11,11));
 478         } else {
 479           nNegPairs++;
 480           fHistCtsLS->Fill(dlsin->CosThetaStarJPSI());
 481           fHistCtsLSneg->Fill(dlsin->CosThetaStarJPSI());
 482         }
 483        PostData(1,fOutput);
 484
 485     // Clone like sign candidate for output AOD
 486     dlsout = new(arrayLikeSignRef[iOutLikeSign++]) AliAODRecoDecayHF2Prong(*dlsin);
 487
 488     //}
 489     if(unsetvtx) dlsin->UnsetOwnPrimaryVtx();
 490   }
 491
 492   //if(fDebug>1) printf("+++\n+++ N. of positive pairs passing cuts in Event ----- %d \n+++\n", nPosPairs);
 493   //if(fDebug>1) printf("+++\n+++ N. of negative pairs passing cuts in Event ----- %d \n+++\n", nNegPairs);
 494
 495   fTotPosPairs += nPosPairs;
 496   fTotNegPairs += nNegPairs;
 497
 498   //
 499   // LOOP OVER J/psi->e+e- CANDIDATES
 500   //
 501
 502   Int_t nInJPSItoEle = arrayJPSItoEle->GetEntriesFast();
 503   //if(fDebug>1) printf("+++\n+++ Number of total like JPSI -> ee candidates (before cuts)---> %d \n+++\n", nInJPSItoEle);
 504
 505   //totJPSIin +=  nInJPSItoEle;
 506   Bool_t isOkEle[2] = {kFALSE,kFALSE};
 507
 508   for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
 509
 510     AliAODRecoDecayHF2Prong *dIn = (AliAODRecoDecayHF2Prong*)arrayJPSItoEle->UncheckedAt(iJPSItoEle);
 511     Int_t mcLabel = 0;
 512     if(fOkAODMC) mcLabel = dIn->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
 513
 514     //Secondary vertex
 515     Double_t vtxSec[3] = {0.,0.,0.};
 516     Double_t vtxPrim[3] = {0.,0.,0.};
 517     Double_t vtxSecOut[3] = {0.,0.,0.};
 518     Double_t vtxPrimOut[3] = {0.,0.,0.};
 519     dIn->GetSecondaryVtx(vtxSec);
 520     Bool_t unsetvtx=kFALSE;
 521     if(!dIn->GetOwnPrimaryVtx()) {
 522       dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 523       unsetvtx=kTRUE;
 524     }
 525
 526     //Int_t okBtoJPSI=0;
 527     //if(dIn->Pt() < fPtCuts[1] && dIn->Pt() > fPtCuts[0]){ // apply pt cut only
 528     //if(dIn->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSI)) {
 529       if ( fOkAODMC && mcLabel == -1){AliDebug(2,"No MC particle found");} else {
 530
 531          //fhInvMass->Fill(dIn->InvMassJPSIee());
 532          fhD0->Fill(10000*dIn->ImpParXY());
 533          fhD0D0->Fill(1e8*dIn->Prodd0d0());
 534          fhCosThetaStar->Fill(dIn->CosThetaStarJPSI());
 535          fhCtsVsD0D0->Fill(dIn->CosThetaStarJPSI(),1e8*dIn->Prodd0d0());
 536          fhCosThetaPointing->Fill(dIn->CosPointingAngle());
 537          fhDCA->Fill(100*dIn->GetDCA());
 538
 539          // compute X variable
 540          AliAODVertex* primVertex = dIn->GetOwnPrimaryVtx();
 541          vtxPrim[0] = primVertex->GetX();
 542          vtxPrim[1] = primVertex->GetY();
 543          vtxPrim[2] = primVertex->GetZ();
 544          Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dIn->Px()) + (vtxSec[1]-vtxPrim[1])*(dIn->Py()))/dIn->Pt();
 545          Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dIn->Pt();
 546          fhDecayTime->Fill(10000*pseudoProperDecayTime);
 547
 548          // retrieve daughter tracks
 549          AliAODTrack *trk0 = (AliAODTrack*)dIn->GetDaughter(0);
 550          TBits clusterMap0 = trk0->GetTPCClusterMap();
 551          Int_t npoints0 = clusterMap0.CountBits(0);
 552          AliAODPid *pid0 = trk0->GetDetPid();
 553          if(GetNumberOfSigmas(trk0->P(),pid0->GetTPCsignal(),npoints0,AliPID::kElectron) < 3.) isOkEle[0]=kTRUE;
 554
 555          AliAODTrack *trk1 = (AliAODTrack*)dIn->GetDaughter(1);
 556          TBits clusterMap1 = trk1->GetTPCClusterMap();
 557          Int_t npoints1 = clusterMap1.CountBits(0);
 558          AliAODPid *pid1 = trk1->GetDetPid();
 559          if(GetNumberOfSigmas(trk1->P(),pid1->GetTPCsignal(),npoints1,AliPID::kElectron) < 3.) isOkEle[1]=kTRUE;
 560
 561          if(isOkEle[0] && isOkEle[1]) {
 562            fhInvMass->Fill(dIn->InvMassJPSIee());
 563            fhdEdxTPC->Fill(pid0->GetTPCmomentum(),pid0->GetTPCsignal());
 564            fhdEdxTPC->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal());
 565          }
 566
 567          //printf("TPC momentum %f\n", pid0->GetTPCmomentum());
 568          //printf("TPC signal %f\n", pid0->GetTPCsignal());
 569
 570          // clone candidate for output AOD
 571          AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++])
 572            AliAODVertex(*(dIn->GetSecondaryVtx()));
 573          AliAODRecoDecayHF2Prong* dOut = new(arrayJpsiToEleRef[iOutJPSItoEle++])
 574            AliAODRecoDecayHF2Prong(*dIn); // copy constructor used
 575          dOut->SetSecondaryVtx(v);
 576          dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone()));
 577          v->SetParent(dOut);
 578
 579          // compute X variable from the cloned candidate in the stand-alone AOD file
 580          dOut->GetSecondaryVtx(vtxSecOut);
 581          AliAODVertex* primVertexOut = dOut->GetOwnPrimaryVtx();
 582          vtxPrimOut[0] = primVertexOut->GetX();
 583          vtxPrimOut[1] = primVertexOut->GetY();
 584          vtxPrimOut[2] = primVertexOut->GetZ();
 585          Double_t lxyOut = ((vtxSecOut[0]-vtxPrimOut[0])*(dOut->Px()) + (vtxSecOut[1]-vtxPrimOut[1])*(dOut->Py()))/dOut->Pt();
 586          Double_t pseudoProperDecayTimeOut = lxyOut*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dOut->Pt();
 587          fhDecayTimeOut->Fill(10000*pseudoProperDecayTimeOut);
 588
 589          //AliAODTrack*  trk0 = (AliAODTrack*)dOut->GetDaughter(0);
 590          //printf("+++ Pt first daughter = %d \n +++ \n", trk0->Pt());
 591
 592         }
 593
 594      //} // end of JPSItoEle candidates selection according to cuts
 595
 596     if(unsetvtx) dIn->UnsetOwnPrimaryVtx();
 597
 598     PostData(0,fOutput);
 599
 600  }// end loop on JPSI to ele candidates
 601
 602  //printf("+++\n+++ Number of selected J/psi->e+e-: %d\n+++\n",iOutJPSItoEle);
 603
 604  //totJPSIout += iOutJPSItoEle;
 605
 606  return;
 607
 608 }
 609 //_________________________________________________________________________________________________
 610 void AliAnalysisTaskSEJPSItoEle::Terminate(Option_t */*option*/)
 611 {
 612   //
 613   // Terminate analysis
 614   //
 615
 616   if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle: Terminate() \n");
 617
 618   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 619   if (!fOutput) {
 620     printf("ERROR: fOutput not available\n");
 621     return;
 622   }
 623
 624   fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs);
 625
 626   if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
 627   fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
 628   fhDecayTimeOut = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeOut"));
 629   fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
 630   fhdEdxTPC = dynamic_cast<TH2F*>(fOutput->FindObject("fhdEdxTPC"));
 631   fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
 632   fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
 633   fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
 634   fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
 635   fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
 636   fHistMassLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassLS"));
 637   fHistCtsLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLS"));
 638   fHistCtsLSpos = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSpos"));
 639   fHistCtsLSneg = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSneg"));
 640   fHistCPtaLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaLS"));
 641   fHistd0d0LS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0LS"));
 642   fhDCA = dynamic_cast<TH1F*>(fOutput->FindObject("fhDCA"));
 643   fHistDCALS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistDCALS"));
 644
 645   if(fLsNormalization>0.) {
 646     fHistMassLS->Scale((1/fLsNormalization)*fHistMassLS->GetEntries());
 647     fHistCtsLS->Scale((1/fLsNormalization)*fHistCtsLS->GetEntries());
 648     fHistCtsLSpos->Scale((1/fLsNormalization)*fHistCtsLSpos->GetEntries());
 649     fHistCtsLSneg->Scale((1/fLsNormalization)*fHistCtsLSneg->GetEntries());
 650     fHistCPtaLS->Scale((1/fLsNormalization)*fHistCPtaLS->GetEntries());
 651     fHistd0d0LS->Scale((1/fLsNormalization)*fHistd0d0LS->GetEntries());
 652     fHistDCALS->Scale((1/fLsNormalization)*fHistDCALS->GetEntries());
 653   }
 654
 655   //printf("Tot JPSI in %d\n", totJPSIin);
 656   //printf("Tot JPSI out %d\n", totJPSIout);
 657
 658   return;
 659 }
 660 //_________________________________________________________________________________________________
 661 void AliAnalysisTaskSEJPSItoEle::SetPtCuts(const Double_t ptCuts[2])
 662 {
 663   //
 664   // apply Pt cuts (lower cut, upper cut)
 665   //
 666   for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i];
 667 }
 668 //_________________________________________________________________________________________________
 669 void AliAnalysisTaskSEJPSItoEle::SetCutsJPSI(const Double_t cuts[9])
 670 {
 671   //
 672   // apply JPSI cuts
 673   //
 674   for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
 675 }
 676 //_________________________________________________________________________________________________
 677 void AliAnalysisTaskSEJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray)
 678 {
 679   //
 680   // Reads MC information if needed (for example in case of parametrized PID)
 681   //
 682
 683   // load MC particles
 684   TClonesArray* mcArray =
 685      dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 686   if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 687      AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
 688
 689   // load MC header
 690   AliAODMCHeader* mcHeader =
 691     (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 692   if(!mcHeader){
 693     printf("AliAnalysisTaskSEJPSItoEle::UserExec: MC AOD header branch not found!\n");
 694   }
 695
 696   AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 697
 698   Double_t mcPrimVtx[3];
 699
 700   // get MC primary Vtx
 701   mcHeader->GetVertex(mcPrimVtx);
 702
 703   Int_t nInCandidates = inputArray->GetEntriesFast();
 704   printf("Number of Candidates for MC analysis: %d\n",nInCandidates);
 705
 706   Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother;
 707   Int_t labJPSIdaugh0=0;
 708   Int_t labJPSIdaugh1=0;
 709
 710   for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) {
 711
 712     AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate);
 713     Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
 714
 715     Double_t vtxPrim[3] = {0.,0.,0.};
 716     Double_t vtxSec[3] = {0.,0.,0.};
 717     dd->GetSecondaryVtx(vtxSec);
 718     Bool_t unsetvtx=kFALSE;
 719     if(!dd->GetOwnPrimaryVtx()) {
 720       dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 721       unsetvtx=kTRUE;
 722     }
 723
 724     if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only
 725
 726       AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0);
 727       AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1);
 728       lab0 = trk0->GetLabel();
 729       lab1 = trk1->GetLabel();
 730
 731       AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0);
 732       if(!part0) {
 733         printf("no MC particle\n");
 734         continue;
 735       }
 736
 737       while(part0->GetMother()>=0) {
 738        labMother=part0->GetMother();
 739        part0 = (AliAODMCParticle*)mcArray->At(labMother);
 740        if(!part0) {
 741          printf("no MC mother particle\n");
 742          break;
 743        }
 744        pdgMother = TMath::Abs(part0->GetPdgCode());
 745        if(pdgMother==443) {//this for JPSI
 746          labJPSIdaugh0=labMother;
 747          break;
 748        }
 749       }
 750
 751       AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1);
 752       if(!part1) {
 753         printf("no MC particle\n");
 754         continue;
 755       }
 756
 757       while(part1->GetMother()>=0) {
 758        labMother=part1->GetMother();
 759        part1 = (AliAODMCParticle*)mcArray->At(labMother);
 760        if(!part1) {
 761          printf("no MC mother particle\n");
 762          break;
 763        }
 764        pdgMother = TMath::Abs(part1->GetPdgCode());
 765        if(pdgMother==443) {//this for JPSI
 766          labJPSIdaugh1=labMother;
 767          break;
 768        }
 769       }
 770
 771       if (mcLabelSec == -1)
 772          {
 773           AliDebug(2,"No MC particle found");
 774          }
 775       else {
 776
 777       if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
 778           AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
 779           pdgJPSI = partJPSI->GetPdgCode();
 780           Int_t pdaugh0 = partJPSI->GetDaughter(0);
 781           Int_t pdaugh1 = partJPSI->GetDaughter(1);
 782           AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
 783           AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
 784           pdg0 = partDaugh0->GetPdgCode();
 785           pdg1 = partDaugh1->GetPdgCode();
 786         if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
 787            labJPSIMother = partJPSI->GetMother();
 788            AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
 789            pdgJPSIMother = partJPSIMother->GetPdgCode();
 790
 791         // chech if JPSI comes from B
 792            if( pdgJPSIMother==511   || pdgJPSIMother==521   ||
 793                pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
 794                pdgJPSIMother==513   || pdgJPSIMother==523   ||
 795                pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
 796                pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
 797                pdgJPSIMother==515   || pdgJPSIMother==525   ||
 798                pdgJPSIMother==531   || pdgJPSIMother==10531 ||
 799                pdgJPSIMother==533   || pdgJPSIMother==10533 ||
 800                pdgJPSIMother==20533 || pdgJPSIMother==535   ||
 801                pdgJPSIMother==541   || pdgJPSIMother==10541 ||
 802                pdgJPSIMother==543   || pdgJPSIMother==10543 ||
 803                pdgJPSIMother==20543 || pdgJPSIMother==545)
 804                { //this is for MC JPSI -> ee from B-hadrons
 805
 806                   fhInvMass->Fill(dd->InvMassJPSIee());
 807                   fhD0->Fill(10000*dd->ImpParXY());
 808                   fhD0D0->Fill(1e8*dd->Prodd0d0());
 809                   fhCosThetaStar->Fill(dd->CosThetaStarJPSI());
 810                   fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0());
 811                   fhCosThetaPointing->Fill(dd->CosPointingAngle());
 812
 813                   // compute X variable
 814                   AliAODVertex* primVertex = dd->GetOwnPrimaryVtx();
 815                   vtxPrim[0] = primVertex->GetX();
 816                   vtxPrim[1] = primVertex->GetY();
 817                   vtxPrim[2] = primVertex->GetZ();
 818                   Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt();
 819                   Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt();
 820
 821                   fhDecayTime->Fill(10000*pseudoProperDecayTime);
 822                   // Post the data already here
 823                   PostData(1,fOutput);
 824
 825                   Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt();
 826                   Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt();
 827                   fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1);
 828
 829                   // Post the data already here
 830                   PostData(1,fOutput);
 831
 832             } //this is for MC JPSI -> ee from B-hadrons
 833           } //this is for MC JPSI -> ee
 834         }
 835       } // end of check if JPSI is signal
 836     }
 837   }
 838 }