PWG3/vertexingHF/AliAnalysisTaskSECompareHF.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, 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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // AliAnalysisTaskSE for the comparison of heavy flavor
  21 // decay candidates with the MC truth.
  22 //
  23 // Author: A.Dainese, andrea.dainese@lnl.infn.it
  24 /////////////////////////////////////////////////////////////
  25
  26 #include <TClonesArray.h>
  27 #include <TNtuple.h>
  28 #include <TList.h>
  29 #include <TH1F.h>
  30
  31 #include "AliAnalysisManager.h"
  32 #include "AliAODHandler.h"
  33 #include "AliAODEvent.h"
  34 #include "AliAODVertex.h"
  35 #include "AliESDtrack.h"
  36 #include "AliAODTrack.h"
  37 #include "AliAODMCHeader.h"
  38 #include "AliAODMCParticle.h"
  39 #include "AliAODRecoDecayHF2Prong.h"
  40 #include "AliAODRecoDecayHF3Prong.h"
  41 #include "AliAODRecoDecayHF4Prong.h"
  42 #include "AliAODRecoCascadeHF.h"
  43 #include "AliAnalysisVertexingHF.h"
  44 #include "AliAnalysisTaskSE.h"
  45 #include "AliAnalysisTaskSECompareHF.h"
  46
  47 ClassImp(AliAnalysisTaskSECompareHF)
  48
  49
  50 //________________________________________________________________________
  51 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
  52 AliAnalysisTaskSE(),
  53 fOutput(0),
  54 fNtupleCmp(0),
  55 fHistMass(0),
  56 fHistNEvents(0)
  57 {
  58   // Default constructor
  59
  60   // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
  61 }
  62
  63 //________________________________________________________________________
  64 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
  65 AliAnalysisTaskSE(name),
  66 fOutput(0),
  67 fNtupleCmp(0),
  68 fHistMass(0),
  69 fHistNEvents(0)
  70 {
  71   // Standard constructor
  72
  73   // Output slot #1 writes into a TList container
  74   DefineOutput(1,TList::Class());  //My private output
  75   // Output slot #2 writes into a TNtuple container
  76   DefineOutput(2,TNtuple::Class());  //My private output
  77 }
  78
  79 //________________________________________________________________________
  80 AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
  81 {
  82   // Destructor
  83   if (fOutput) {
  84     delete fOutput;
  85     fOutput = 0;
  86   }
  87 }
  88
  89 //________________________________________________________________________
  90 void AliAnalysisTaskSECompareHF::Init()
  91 {
  92   // Initialization
  93
  94   if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
  95
  96   return;
  97 }
  98
  99 //________________________________________________________________________
 100 void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
 101 {
 102   // Create the output container
 103   //
 104   if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
 105
 106   // Several histograms are more conveniently managed in a TList
 107   fOutput = new TList();
 108   fOutput->SetOwner();
 109
 110   fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
 111   fHistMass->Sumw2();
 112   fHistMass->SetMinimum(0);
 113   fOutput->Add(fHistMass);
 114
 115   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
 116   fHistNEvents->Sumw2();
 117   fHistNEvents->SetMinimum(0);
 118   fOutput->Add(fHistNEvents);
 119
 120   OpenFile(2); // 2 is the slot number of the ntuple
 121   fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0");
 122
 123   return;
 124 }
 125
 126 //________________________________________________________________________
 127 void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
 128 {
 129   // Execute analysis for current event:
 130   // heavy flavor candidates association to MC truth
 131
 132
 133   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 134
 135   TClonesArray *inputArrayVertices = 0;
 136   TClonesArray *inputArrayD0toKpi = 0;
 137   TClonesArray *inputArrayDstar = 0;
 138
 139   if(!aod && AODEvent() && IsStandardAOD()) {
 140     // In case there is an AOD handler writing a standard AOD, use the AOD
 141     // event in memory rather than the input (ESD) event.
 142     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 143     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 144     // have to taken from the AOD event hold by the AliAODExtension
 145     AliAODHandler* aodHandler = (AliAODHandler*)
 146       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 147     if(aodHandler->GetExtensions()) {
 148       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 149       AliAODEvent *aodFromExt = ext->GetAOD();
 150       // load HF vertices
 151       inputArrayVertices = (TClonesArray*)aodFromExt->GetList()->FindObject("VerticesHF");
 152       // load D0->Kpi candidates
 153       inputArrayD0toKpi = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 154       // load D*+ candidates
 155       inputArrayDstar = (TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 156     }
 157   } else if(aod) {
 158     // load HF vertices
 159     inputArrayVertices = (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
 160     // load D0->Kpi candidates
 161     inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 162     // load D*+ candidates
 163     inputArrayDstar = (TClonesArray*)aod->GetList()->FindObject("Dstar");
 164   }
 165
 166
 167   if(!inputArrayVertices || !aod) {
 168     printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
 169     return;
 170   }
 171   if(!inputArrayD0toKpi) {
 172     printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
 173     return;
 174   }
 175   if(!inputArrayDstar) {
 176     printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
 177     return;
 178   }
 179
 180
 181   fHistNEvents->Fill(0); // count event
 182   // Post the data already here
 183   PostData(1,fOutput);
 184
 185   // AOD primary vertex
 186   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 187   //vtx1->Print();
 188
 189   // load MC particles
 190   TClonesArray *mcArray =
 191     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 192   if(!mcArray) {
 193     printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
 194     return;
 195   }
 196
 197   // load MC header
 198   AliAODMCHeader *mcHeader =
 199     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 200   if(!mcHeader) {
 201     printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
 202     return;
 203   }
 204
 205   Int_t nprongs,lab,okD0,okD0bar,pdg;
 206   Bool_t unsetvtx;
 207   Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz;
 208   AliAODRecoDecayHF2Prong *d2=0;
 209   AliAODRecoDecayHF3Prong *d3=0;
 210   AliAODRecoDecayHF4Prong *d4=0;
 211
 212   Int_t pdgDgD0toKpi[2]={321,211};
 213   Int_t pdgDgDplustoKpipi[3]={321,211,211};
 214   Int_t pdgDgD0toKpipipi[4]={321,211,211,211};
 215
 216   // loop over vertices
 217   Int_t nVertices = inputArrayVertices->GetEntriesFast();
 218   if(fDebug>1) printf("Number of vertices: %d\n",nVertices);
 219
 220   for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
 221     AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);
 222
 223     vtx->GetXYZ(posRec);
 224     vtx->GetCovarianceMatrix(covRec);
 225     errx=1.; erry=1.; errz=1.;
 226     if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]);
 227     if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]);
 228     if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]);
 229
 230
 231     // get parent AliAODRecoDecayHF
 232     nprongs= vtx->GetNDaughters();
 233     // check that the daughters have kITSrefit and kTPCrefit
 234     Bool_t allDgOK=kTRUE;
 235     for(Int_t idg=0; idg<nprongs; idg++) {
 236       AliAODTrack *track = (AliAODTrack*)vtx->GetDaughter(idg);
 237       if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE;
 238       if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE;
 239     }
 240     if(!allDgOK) continue;
 241
 242
 243     switch(nprongs) {
 244     case 2: // look for D0->Kpi
 245       d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
 246       if(d2->IsLikeSign()) continue;
 247       if(d2->Charge() != 0) continue; // these are D*
 248       lab = d2->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
 249       if(lab>=0) {
 250         unsetvtx=kFALSE;
 251         if(!d2->GetOwnPrimaryVtx()) {
 252           d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 253           unsetvtx=kTRUE;
 254         }
 255         okD0=1; okD0bar=1;
 256         AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
 257         pdg = dMC->GetPdgCode();
 258         invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
 259         // get a daughter for true pos of decay vertex
 260         AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
 261         dg0MC->XvYvZv(posTrue);
 262         fHistMass->Fill(invmass);
 263         // Post the data already here
 264         PostData(1,fOutput);
 265         Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
 266                          (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
 267                          (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
 268                          (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
 269                          (Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass,
 270                          (Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()};
 271         fNtupleCmp->Fill(tmp);
 272         PostData(2,fNtupleCmp);
 273
 274         if(unsetvtx) d2->UnsetOwnPrimaryVtx();
 275       }
 276       break;
 277     case 3: // look for D+
 278       d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
 279       if(d3->IsLikeSign()) continue;
 280       lab = d3->MatchToMC(411,mcArray,3,pdgDgDplustoKpipi);
 281       if(lab>=0) {
 282         unsetvtx=kFALSE;
 283         if(!d3->GetOwnPrimaryVtx()) {
 284           d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 285           unsetvtx=kTRUE;
 286         }
 287         AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
 288         pdg = dMC->GetPdgCode();
 289         invmass = d3->InvMassDplus();
 290         // get a daughter for true pos of decay vertex
 291         AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
 292         dg0MC->XvYvZv(posTrue);
 293         Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
 294                          (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
 295                          (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
 296                          (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
 297                          (Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass,
 298                          (Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)*d3->Getd0Prong(1)*d3->Getd0Prong(2))};
 299         fNtupleCmp->Fill(tmp);
 300         PostData(2,fNtupleCmp);
 301         if(unsetvtx) d3->UnsetOwnPrimaryVtx();
 302       }
 303       break;
 304     case 4: // look for D0->Kpipipi
 305       d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
 306       if(d4->IsLikeSign()) continue;
 307       lab = d4->MatchToMC(421,mcArray,4,pdgDgD0toKpipipi);
 308       if(lab>=0) {
 309         unsetvtx=kFALSE;
 310         if(!d4->GetOwnPrimaryVtx()) {
 311           d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 312           unsetvtx=kTRUE;
 313         }
 314         okD0=0; okD0bar=0;
 315         AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
 316         pdg = dMC->GetPdgCode();
 317         //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
 318         invmass = 10.;    // dummy
 319         // get a daughter for true pos of decay vertex
 320         AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
 321         dg0MC->XvYvZv(posTrue);
 322         Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
 323                          (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
 324                          (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
 325                          (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
 326                          (Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass,
 327                          (Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)*d4->Getd0Prong(1)*d4->Getd0Prong(2)*d4->Getd0Prong(3))};
 328         fNtupleCmp->Fill(tmp);
 329         PostData(2,fNtupleCmp);
 330         if(unsetvtx) d4->UnsetOwnPrimaryVtx();
 331       }
 332       break;
 333     }
 334
 335   } // end loop on vertices
 336
 337
 338   // loop over D*+ candidates
 339   /*
 340   for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
 341     AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
 342     Int_t labDstar = c->MatchToMC(413,421,mcArray);
 343     if(labDstar>=0) printf("GOOD MATCH FOR D*+\n");
 344   }
 345   */
 346
 347   // Post the data already here
 348   PostData(1,fOutput);
 349   PostData(2,fNtupleCmp);
 350
 351   return;
 352 }
 353 //________________________________________________________________________
 354 void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
 355 {
 356   // Terminate analysis
 357   //
 358   if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
 359
 360   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 361   if (!fOutput) {
 362     printf("ERROR: fOutput not available\n");
 363     return;
 364   }
 365
 366   fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
 367   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 368
 369   //fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
 370
 371   return;
 372 }
 373