PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.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 AliAnalysisTaskSEBtoJPSItoEle
  17 //                AliAnalysisTaskSE for the reconstruction
  18 //                   of heavy-flavour decay candidates
  19 //            Author: C.Di Giglio, carmelo.digiglio@ba.infn.it
  20 //*************************************************************************
  21 class AliAnalysisTaskSE;
  22 class AliESDEvent;
  23 class AliVEvent;
  24 class iostream;
  25 class TSystem;
  26 class TROOT;
  27 #include <TFile.h>
  28 #include <TClonesArray.h>
  29 #include <TDatabasePDG.h>
  30
  31 #include "AliAODEvent.h"
  32 #include "AliAODMCParticle.h"
  33 #include "AliAODMCHeader.h"
  34 #include "AliAODRecoDecayHF2Prong.h"
  35
  36 #include "AliAnalysisBtoJPSItoEle.h"
  37 #include "AliAnalysisTaskSEBtoJPSItoEle.h"
  38
  39 ClassImp(AliAnalysisTaskSEBtoJPSItoEle)
  40
  41 //______________________________________________________________________________
  42 AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle() :
  43 AliAnalysisTaskSE(),
  44 fOutput(0),
  45 fCdfFit(0),
  46 fNtupleJPSI(0),
  47 fhDecayTimeMCjpsifromB(0),
  48 fhDecayTime(0),
  49 fhInvMass(0),
  50 fhD0(0),
  51 fhD0D0(0),
  52 fhCosThetaStar(0),
  53 fhCosThetaPointing(0),
  54 fhCtsVsD0D0(0),
  55 fOkAODMC(kFALSE),
  56 fNameMCfile(""),
  57 fOkMinimize(kFALSE)
  58 {
  59   // default constructor
  60
  61   // Output slot #1 writes into a TList container
  62   DefineOutput(1,TList::Class());  //My private output
  63 }
  64 //_________________________________________________________________________________________________
  65 AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle(const char* name) :
  66 AliAnalysisTaskSE(name),
  67 fOutput(0),
  68 fCdfFit(0),
  69 fNtupleJPSI(0),
  70 fhDecayTimeMCjpsifromB(0),
  71 fhDecayTime(0),
  72 fhInvMass(0),
  73 fhD0(0),
  74 fhD0D0(0),
  75 fhCosThetaStar(0),
  76 fhCosThetaPointing(0),
  77 fhCtsVsD0D0(0),
  78 fOkAODMC(kFALSE),
  79 fNameMCfile(""),
  80 fOkMinimize(kFALSE)
  81 {
  82   // default constructor
  83
  84   // Output slot #1 writes into a TList container
  85   DefineOutput(1,TList::Class());  //My private output
  86 }
  87 //_________________________________________________________________________________________________
  88 AliAnalysisTaskSEBtoJPSItoEle::~AliAnalysisTaskSEBtoJPSItoEle()
  89 {
  90   // destructor
  91
  92   if (fOutput) {
  93     delete fOutput;
  94     fOutput = 0;
  95   }
  96 }
  97 //_________________________________________________________________________________________________
  98 void AliAnalysisTaskSEBtoJPSItoEle::Init()
  99 {
 100   // Initialization
 101   //Double_t ptCuts[2] = {1.,1.5}; // 1 GeV < pt < 1.5 GeV
 102   Double_t ptCuts[2] = {1.,100}; //
 103   SetPtCuts(ptCuts);
 104   SetMinimize(kTRUE);
 105   SetAODMCInfo(kTRUE);
 106
 107   if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::Init() \n");
 108
 109   return;
 110
 111 }
 112 //_________________________________________________________________________________________________
 113 void AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects()
 114 {
 115   // Create the output container
 116
 117   if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() \n");
 118
 119   // Several histograms are more conveniently managed in a TList
 120   fOutput = new TList();
 121   fOutput->SetOwner();
 122
 123   if(fOkAODMC){
 124   fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.);
 125   fhDecayTimeMCjpsifromB->Sumw2();
 126   fhDecayTimeMCjpsifromB->SetMinimum(0);
 127   fOutput->Add(fhDecayTimeMCjpsifromB);
 128   }
 129
 130   fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.);
 131   fhDecayTime->Sumw2();
 132   fhDecayTime->SetMinimum(0);
 133   fOutput->Add(fhDecayTime);
 134
 135   fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; M [GeV]; Entries",100,0.,3.2);
 136   fhInvMass->Sumw2();
 137   fhInvMass->SetMinimum(0);
 138   fOutput->Add(fhInvMass);
 139
 140   fhD0 = new TH1F("fhD0", "Impact parameter; D0 [#mu m]; Entries",100,-5000.,5000.);
 141   fhD0->Sumw2();
 142   fhD0->SetMinimum(0);
 143   fOutput->Add(fhD0);
 144
 145   fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.);
 146   fhD0D0->Sumw2();
 147   fhD0D0->SetMinimum(0);
 148   fOutput->Add(fhD0D0);
 149
 150   fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2);
 151   fhCosThetaStar->Sumw2();
 152   fhCosThetaStar->SetMinimum(0);
 153   fOutput->Add(fhCosThetaStar);
 154
 155   fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1);
 156   fhCosThetaPointing->Sumw2();
 157   fhCosThetaPointing->SetMinimum(0);
 158   fOutput->Add(fhCosThetaPointing);
 159
 160   fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000);
 161   fhCtsVsD0D0->Sumw2();
 162   fhCtsVsD0D0->SetMinimum(0);
 163   fOutput->Add(fhCtsVsD0D0);
 164
 165   fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#Psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass");
 166   fOutput->Add(fNtupleJPSI);
 167
 168   return;
 169 }
 170 //_________________________________________________________________________________________________
 171 void AliAnalysisTaskSEBtoJPSItoEle::UserExec(Option_t */*option*/)
 172 {
 173   // Execute analysis for current event:
 174
 175   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 176
 177   // AOD primary vertex
 178   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 179
 180   // load JPSI candidates
 181   TClonesArray *inputArrayJPSItoEle =
 182     (TClonesArray*)aod->GetList()->FindObject("JPSItoEle");
 183   if(!inputArrayJPSItoEle) {
 184     printf("AliAnalysisTaskSECompareHF::UserExec: JPSItoEle branch not found!\n");
 185     return;
 186   }
 187
 188   // Read AOD Monte Carlo information
 189   if(fOkAODMC) ReadAODMCInfo(aod,inputArrayJPSItoEle);
 190
 191   // loop over J/Psi->ee candidates
 192   Int_t nInJPSItoEle = inputArrayJPSItoEle->GetEntriesFast();
 193   printf("Number of B->JPSI->e+e-: %d\n",nInJPSItoEle);
 194
 195   for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
 196     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayJPSItoEle->UncheckedAt(iJPSItoEle);
 197     //Secondary vertex
 198     Double_t vtxSec[3] = {0.,0.,0.};
 199     Double_t vtxPrim[3] = {0.,0.,0.};
 200     d->GetSecondaryVtx(vtxSec);
 201     Bool_t unsetvtx=kFALSE;
 202     if(!d->GetOwnPrimaryVtx()) {
 203       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 204       unsetvtx=kTRUE;
 205     }
 206     //Int_t okJPSI=0;
 207     // here analyze only if cuts are passed
 208     //if(d->SelectBtoJPSI(fCuts,okJPSI)) {
 209     if(d->Pt() < fPtCuts[1] && d->Pt() > fPtCuts[0]){ // apply pt cut only
 210       fhInvMass->Fill(d->InvMassJPSIee());
 211       fhD0->Fill(10000*d->ImpParXY());
 212       fhD0D0->Fill(1e8*d->Prodd0d0());
 213       fhCosThetaStar->Fill(d->CosThetaStarJPSI());
 214       fhCtsVsD0D0->Fill(d->CosThetaStarJPSI(),1e8*d->Prodd0d0());
 215       fhCosThetaPointing->Fill(d->CosPointingAngle());
 216
 217       // compute X variable
 218       AliAODVertex* primVertex = d->GetOwnPrimaryVtx();
 219       vtxPrim[0] = primVertex->GetX();
 220       vtxPrim[1] = primVertex->GetY();
 221       vtxPrim[2] = primVertex->GetZ();
 222       Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(d->Px()) + (vtxSec[1]-vtxPrim[1])*(d->Py()))/d->Pt();
 223       Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/d->Pt();
 224
 225       fhDecayTime->Fill(10000*pseudoProperDecayTime);
 226       // Post the data already here
 227       PostData(1,fOutput);
 228
 229       fNtupleJPSI->Fill(pseudoProperDecayTime,d->InvMassJPSIee()); // fill N-tuple with X,M values
 230
 231      } // end of JPSItoEle candidates selection according to cuts
 232
 233     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 234
 235    }// end loop on JPSI to ele candidates
 236
 237 }
 238 //_________________________________________________________________________________________________
 239 void AliAnalysisTaskSEBtoJPSItoEle::Terminate(Option_t */*option*/)
 240 {
 241   //
 242   // Terminate analysis
 243   //
 244
 245   if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle: Terminate() \n");
 246
 247   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 248   if (!fOutput) {
 249     printf("ERROR: fOutput not available\n");
 250     return;
 251   }
 252
 253   if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
 254   fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
 255   fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
 256   fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
 257   fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
 258   fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
 259   fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
 260   fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
 261   fNtupleJPSI = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleJPSI"));
 262
 263   if(fOkMinimize){
 264      Double_t* pseudoProperValues=0x0;
 265      Double_t* invMassValues=0x0;
 266      Int_t ncand=0;
 267      fCdfFit = new AliAnalysisBtoJPSItoEle();
 268
 269      printf("+++\n+++  AliAnalysisBtoJPSItoEle object instantiated ---> OK\n+++\n");
 270      fCdfFit->ReadCandidates(fNtupleJPSI,pseudoProperValues,invMassValues,ncand);
 271      printf("+++\n+++  X and M vectors filled starting from N-tuple ---> OK\n+++\n");
 272      printf("+++\n+++  Number of candidates ---> %d J/psi \n+++\n",ncand);
 273      fCdfFit->SetPtBin(0);
 274      printf("+++\n+++  Pt bin setted n. ---> %d  \n+++\n",fCdfFit->GetPtBin());
 275
 276      if(fOkAODMC) {fCdfFit->CloneMCtemplate(fhDecayTimeMCjpsifromB);}
 277         else {
 278           //SetPathMCfile(".");
 279           TH1F* hLocal = OpenLocalMCFile(fNameMCfile,fCdfFit->GetPtBin());
 280           fCdfFit->CloneMCtemplate(hLocal);
 281         }
 282
 283      printf("+++\n+++  MC template histo copied ---> OK\n+++\n");
 284      fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand);
 285    }
 286
 287   return;
 288 }
 289 //_________________________________________________________________________________________________
 290 void AliAnalysisTaskSEBtoJPSItoEle::SetPtCuts(const Double_t ptCuts[2])
 291 {
 292   //
 293   // apply Pt cuts (lower cut, upper cut)
 294   //
 295   for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i];
 296 }
 297 //_________________________________________________________________________________________________
 298 void AliAnalysisTaskSEBtoJPSItoEle::SetCutsD0(const Double_t cuts[9])
 299 {
 300   //
 301   // apply D0 and JPSI cuts
 302   //
 303   for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
 304 }
 305 //_________________________________________________________________________________________________
 306 void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(const AliAODEvent* aodEvent, const TClonesArray* inputArray)
 307 {
 308   //
 309   // Reads MC information if needed (for example in case of parametrized PID)
 310   //
 311
 312   // load MC particles
 313   TClonesArray *mcArray =
 314     (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 315   if(!mcArray) {
 316     printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC particles branch not found!\n");
 317     return;
 318   }
 319
 320   // load MC header
 321   AliAODMCHeader* mcHeader =
 322     (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 323   if(!mcHeader){
 324     printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC AOD header branch not found!\n");
 325   }
 326
 327   Double_t rmax = 0.000005;
 328   Double_t mcPrimVtx[3];
 329
 330   // get MC primary Vtx
 331   mcHeader->GetVertex(mcPrimVtx);
 332
 333   Int_t nInCandidates = inputArray->GetEntriesFast();
 334   printf("Number of Candidates for MC analysis: %d\n",nInCandidates);
 335
 336   Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother;
 337   Int_t labJPSIdaugh0=0;
 338   Int_t labJPSIdaugh1=0;
 339
 340   for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) {
 341     AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate);
 342     if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only
 343
 344       AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0);
 345       AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1);
 346       lab0 = trk0->GetLabel();
 347       lab1 = trk1->GetLabel();
 348
 349       AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0);
 350       if(!part0) {
 351         printf("no MC particle\n");
 352         continue;
 353       }
 354
 355       while(part0->GetMother()>=0) {
 356        labMother=part0->GetMother();
 357        part0 = (AliAODMCParticle*)mcArray->At(labMother);
 358        if(!part0) {
 359          printf("no MC mother particle\n");
 360          break;
 361        }
 362        pdgMother = TMath::Abs(part0->GetPdgCode());
 363        if(pdgMother==443) {//this for JPSI
 364          labJPSIdaugh0=labMother;
 365          break;
 366        }
 367       }
 368
 369       AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1);
 370       if(!part1) {
 371         printf("no MC particle\n");
 372         continue;
 373       }
 374
 375       while(part1->GetMother()>=0) {
 376        labMother=part1->GetMother();
 377        part1 = (AliAODMCParticle*)mcArray->At(labMother);
 378        if(!part1) {
 379          printf("no MC mother particle\n");
 380          break;
 381        }
 382        pdgMother = TMath::Abs(part1->GetPdgCode());
 383        if(pdgMother==443) {//this for JPSI
 384          labJPSIdaugh1=labMother;
 385          break;
 386        }
 387       }
 388
 389       if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
 390         AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
 391         pdgJPSI = partJPSI->GetPdgCode();
 392         if(pdgJPSI == 443){//this is for MC JPSI
 393            if(TMath::Sqrt((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Xv()-mcPrimVtx[0])+
 394                           (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Yv()-mcPrimVtx[1])+
 395                           (partJPSI->Zv()-mcPrimVtx[2])*(partJPSI->Zv()-mcPrimVtx[2])>rmax)){ //this is for MC displaced JPSI
 396               Int_t pdaugh0 = partJPSI->GetDaughter(0);
 397               Int_t pdaugh1 = partJPSI->GetDaughter(1);
 398               AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
 399               AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
 400               pdg0 = partDaugh0->GetPdgCode();
 401               pdg1 = partDaugh1->GetPdgCode();
 402               if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
 403                 labJPSIMother = partJPSI->GetMother();
 404                 AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
 405                 pdgJPSIMother = partJPSIMother->GetPdgCode();
 406                if(pdgJPSIMother==511   || pdgJPSIMother==521   ||
 407                   pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
 408                   pdgJPSIMother==513   || pdgJPSIMother==523   ||
 409                   pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
 410                   pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
 411                   pdgJPSIMother==515   || pdgJPSIMother==525   ||
 412                   pdgJPSIMother==531   || pdgJPSIMother==10531 ||
 413                   pdgJPSIMother==533   || pdgJPSIMother==10533 ||
 414                   pdgJPSIMother==20533 || pdgJPSIMother==535   ||
 415                   pdgJPSIMother==541   || pdgJPSIMother==10541 ||
 416                   pdgJPSIMother==543   || pdgJPSIMother==10543 ||
 417                   pdgJPSIMother==20543 || pdgJPSIMother==545)
 418                   { //this is for MC JPSI -> ee from B-hadrons
 419
 420                   Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt();
 421                   Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt();
 422                   fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1);
 423
 424                   // Post the data already here
 425                   PostData(1,fOutput);
 426
 427                 } //this is for MC JPSI -> ee from B-hadrons
 428               } //this is for MC JPSI -> ee
 429             } //this is for MC displaced JPSI
 430           } //this is for MC JPSI
 431         } // end of check if JPSI is signal
 432       }
 433     }
 434
 435 }
 436 //_________________________________________________________________________________________________
 437 TH1F* AliAnalysisTaskSEBtoJPSItoEle::OpenLocalMCFile(TString datadir, Int_t binNum)
 438 {
 439   //
 440   // Open a local file with MC x distribution for JPSI from B-hadron
 441   //
 442
 443   TH1F* hCsiMC = new TH1F();
 444   TString useFile = datadir.Data();
 445   useFile.Append("/CsiMCfromKine_10PtBins.root");
 446   TFile *f = new TFile(useFile);
 447   Double_t scale = 0.;
 448   char processCsiMCXHisto[1024];
 449   if(binNum == 0) sprintf(processCsiMCXHisto,"hPseudoProper");
 450   else sprintf(processCsiMCXHisto,"hPseudoProper%d",binNum);
 451   hCsiMC = (TH1F*)f->Get(processCsiMCXHisto);
 452   scale=1/hCsiMC->Integral();
 453   hCsiMC->Scale(scale);
 454   printf ("Opening Histo with Csi_MC(x) template with n. %f entries", hCsiMC->GetEntries());
 455
 456   return hCsiMC;
 457
 458 }
 459