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