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