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