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