PWG3/vertexingHF/AliAnalysisTaskCharmFraction.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 /////////////////////////////////////////////////////////////
  17 //
  18 // Class AliAnalysisTaskCharmFraction
  19 // AliAnalysisTask for the extraction of the fraction of prompt charm
  20 // using the charm hadron impact parameter to the primary vertex
  21 //
  22 // Author: Andrea Rossi, andrea.rossi@ts.infn.it
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TChain.h>
  26 #include <TTree.h>
  27 #include <TH1F.h>
  28 #include <TH2F.h>
  29 #include <TCanvas.h>
  30 #include <TDatabasePDG.h>
  31 #include <TMath.h>
  32
  33 #include "AliAnalysisTask.h"
  34 #include "AliAnalysisManager.h"
  35 #include "AliAODEvent.h"
  36 #include "AliAODInputHandler.h"
  37 #include "AliAnalysisVertexingHF.h"
  38 #include "AliAODRecoDecayHF2Prong.h"
  39 #include "AliAODRecoDecayHF.h"
  40 #include "AliAODRecoDecay.h"
  41 #include "AliAODTrack.h"
  42 #include "AliAODVertex.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliAODMCHeader.h"
  45 #include "AliAnalysisTaskCharmFraction.h"
  46
  47 ClassImp(AliAnalysisTaskCharmFraction)
  48
  49   //________________________________________________________________________
  50   AliAnalysisTaskCharmFraction::AliAnalysisTaskCharmFraction(const char *name)
  51     : AliAnalysisTask(name, ""), fAOD(0), fVHF(0),
  52       fArrayD0toKpi(0),
  53       fArrayMC(0),
  54       fAODmcHeader(0),
  55       fhMass(0),
  56       fhMassTrue(0),
  57       fhCPtaVSd0d0(0),
  58       fhd0xd0(0),
  59       fhCPta(0),
  60       fhSecVtxZ(0),
  61       fhSecVtxX(0),
  62       fhSecVtxY(0),
  63       fhSecVtxXY(0),
  64       fhSecVtxPhi(0),
  65       fhd0D0(0),
  66       fhd0D0pt(0x0),
  67       fhd0D0VtxTrue(0),
  68       fhd0D0VtxTruept(0x0),
  69       fhMCd0D0(0),
  70       fhMCd0D0pt(0x0),
  71       fnbins(),
  72       fD0usecuts(kTRUE),
  73       fcheckMC(kFALSE),
  74       fcheckMCD0(kFALSE),
  75       fcheckMC2prongs(kFALSE),
  76       fcheckMCprompt(kFALSE),
  77       fcheckMCfromB(kFALSE),
  78       fSkipD0star(kFALSE),
  79       fStudyd0fromBTrue(kTRUE),
  80       fStudyPureBackground(kFALSE),
  81       fSideBands(0)
  82 {
  83   // Constructor
  84   //Side bands selection (see cxx): 999 -> no inv mass selection at all,
  85   //                                >0 ->  both D0 and D0bar hypothesis are required to be kFALSE,
  86   //                                <0 -> only 1 inv mass hypothesis is required to be KFALSE
  87
  88   SetNPtBins();
  89   // Define input and output slots here
  90   // Input slot #0 works with a TChain
  91   DefineInput(0, TChain::Class());
  92   //  DefineInput(1, TChain::Class());
  93   // Output slot #0 writes into a TH1 container
  94   DefineOutput(0, TH2F::Class());// hCPtaVsd0d0
  95   DefineOutput(1, TH2F::Class());// hVtxXY
  96   for(Int_t j=2;j<3*fnbins+11;j++){ // hCPta, hd0xd0, hVtZ, hVtX hVtY hVtPhi + nptbin hd0D0 +hd0D0 integreated +nptbin hd0VtxTrueD0 +hd0VtxTrueD0 integrated +nptbin hMCd0D0 +hMCd0D0 integrated
  97     DefineOutput(j, TH1F::Class());
  98   }
  99   //  DefineOutput(fnbins+8, TH1F::Class());// hd0D0 pt integrated
 100 }
 101
 102 //________________________________________________________________________
 103 AliAnalysisTaskCharmFraction::AliAnalysisTaskCharmFraction(const char *name,Int_t nptbins)
 104   : AliAnalysisTask(name, ""), fAOD(0), fVHF(0),
 105     fArrayD0toKpi(0),
 106     fArrayMC(0),
 107     fAODmcHeader(0),
 108     fhMass(0),
 109     fhMassTrue(0),
 110     fhCPtaVSd0d0(0),
 111     fhd0xd0(0),
 112     fhCPta(0),
 113     fhSecVtxZ(0),
 114     fhSecVtxX(0),
 115     fhSecVtxY(0),
 116     fhSecVtxXY(0),
 117     fhSecVtxPhi(0),
 118     fhd0D0(0),
 119     fhd0D0pt(0x0),
 120     fhd0D0VtxTrue(0),
 121     fhd0D0VtxTruept(0x0),
 122     fhMCd0D0(0),
 123     fhMCd0D0pt(0x0),
 124     fnbins(),
 125     fD0usecuts(kTRUE),
 126     fcheckMC(kFALSE),
 127     fcheckMCD0(kFALSE),
 128     fcheckMC2prongs(kFALSE),
 129     fcheckMCprompt(kFALSE),
 130     fcheckMCfromB(kFALSE),
 131     fSkipD0star(kFALSE),
 132     fStudyd0fromBTrue(kTRUE),
 133     fStudyPureBackground(kFALSE),
 134     fSideBands(0)
 135 {
 136   // Constructor
 137   //Side bands selection (see cxx): 999 -> no inv mass selection at all,
 138   //                                >0 ->  both D0 and D0bar hypothesis are required to be kFALSE,
 139   //                                <0 -> only 1 inv mass hypothesis is required to be KFALSE
 140   SetNPtBins(nptbins);
 141   // Define input and output slots here
 142   // Input slot #0 works with a TChain
 143   DefineInput(0, TChain::Class());
 144   //  DefineInput(1, TChain::Class());
 145   // Output slot #0 writes into a TH1 container
 146
 147   DefineOutput(0, TH2F::Class());// hCPtaVsd0d0
 148   DefineOutput(1, TH2F::Class());// hVtxXY
 149   for(Int_t j=2;j<3*fnbins+13;j++){ // hCPta, hd0xd0, hVtZ, hVtX hVtY hVtPhi + nptbin hd0D0 +hd0D0 integreated +nptbin hd0VtxTrueD0 +hd0VtxTrueD0 integrated +nptbin hMCd0D0 +hMCd0D0 integrated
 150     DefineOutput(j, TH1F::Class());
 151   }
 152   //  DefineOutput(fnbins+8, TH1F::Class());// hd0D0 pt integrated
 153 }
 154
 155 //________________________________________________________________________
 156 void AliAnalysisTaskCharmFraction::ConnectInputData(Option_t *)
 157 {
 158   // Connect ESD or AOD here
 159   // Called once
 160
 161   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 162   //  TTree* treeFriend = dynamic_cast<TTree*> (GetInputData(1));
 163   if (!tree) {
 164     Printf("ERROR: Could not read chain from input slot 0");
 165   } else {
 166     //  tree->AddFriend(treeFriend);
 167
 168     // Disable all branches and enable only the needed ones
 169     // The next two lines are different when data produced as AliESDEvent is read
 170
 171     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 172     fAOD=new AliAODEvent();
 173     fAOD->ReadFromTree(tree);
 174
 175     fArrayD0toKpi = (TClonesArray*)fAOD->GetList()->FindObject("D0toKpi");
 176     fArrayMC =      (TClonesArray*)fAOD->GetList()->FindObject("mcparticles");
 177
 178     fAODmcHeader=(AliAODMCHeader*)fAOD->GetList()->FindObject("mcHeader");
 179
 180     if (!aodH) {
 181       Printf("ERROR: Could not get AODInputHandler");
 182     } else
 183       fAOD = aodH->GetEvent();
 184
 185   }
 186 }
 187
 188 //________________________________________________________________________
 189 void AliAnalysisTaskCharmFraction::CreateOutputObjects()
 190 {
 191   // Create histograms
 192   // Called once
 193   fhd0D0 = new TH1F("hd0D0","D^{0} impact par. plot  (All momenta)",1000,-1000.,1000.);
 194   fhd0D0->SetXTitle("Impact parameter [#mum]");
 195   fhd0D0->SetYTitle("Entries");
 196
 197   fhd0D0VtxTrue = new TH1F("hd0D0VtxTrue","D^{0} impact par. w.r.t. True Vtx (All momenta)",1000,-1000.,1000.);
 198   fhd0D0VtxTrue->SetXTitle("Impact parameter [#mum]");
 199   fhd0D0VtxTrue->SetYTitle("Entries");
 200
 201   fhMCd0D0 = new TH1F("hMCd0D0","D^{0} impact par. plot  (All momenta)",1000,-1000.,1000.);
 202   fhMCd0D0->SetXTitle("MC Impact parameter [#mum]");
 203   fhMCd0D0->SetYTitle("Entries");
 204
 205   fhCPtaVSd0d0=new TH2F("hCPtaVSd0d0","hCPtaVSd0d0",1000,-100000.,100000.,100,0.,1.);
 206   fhSecVtxZ=new TH1F("hSecVtxZ","hSecVtxZ",1000,-8.,8.);
 207   fhSecVtxX=new TH1F("hSecVtxX","hSecVtxX",1000,-3000.,3000.);
 208   fhSecVtxY=new TH1F("hSecVtxY","hSecVtxY",1000,-3000.,3000.);
 209   fhSecVtxXY=new TH2F("hSecVtxXY","hSecVtxXY",1000,-3000.,3000.,1000,-3000.,3000.);
 210   fhSecVtxPhi=new TH1F("hSecVtxPhi","hSecVtxPhi",180,-180.1,180.1);
 211   fhCPta=new TH1F("hCPta","hCPta",100,0.,1.);
 212   fhd0xd0=new TH1F("hd0xd0","hd0xd0",1000,-100000.,100000.);
 213   fhMassTrue=new TH1F("hMassTrue","D^{0} MC inv. Mass (All momenta)",600,1.600,2.200);
 214   fhMass=new TH1F("hMass","D^{0} inv. Mass (All momenta)",600,1.600,2.200);
 215   Double_t standardptbins[14]={0.,0.5,1.,2.,3.,5.,7.,10.,15.,20.,30.,40.,50.,100.};
 216
 217   fhd0D0pt=new TH1F*[fnbins];
 218   fhMCd0D0pt=new TH1F*[fnbins];
 219   fhd0D0VtxTruept=new TH1F*[fnbins];
 220   TString namehist="hd0D0pt";
 221   TString titlehist="D^{0} impact par. plot ";
 222   TString strnamept,strtitlept;
 223   for(Int_t i=0;i<fnbins;i++){
 224     strnamept=namehist;
 225     strnamept+=i;// strnamept+=standardptbins[i+1];
 226     //    strnamept.Append("_");
 227     //strnamept+=standardptbins[i+1];
 228
 229     strtitlept=namehist;
 230     strtitlept+=standardptbins[i];
 231     strtitlept.Append("<= pt <");
 232     strtitlept+=standardptbins[i+1];
 233     strtitlept.Append(" [GeV/c]");
 234
 235     fhd0D0pt[i] = new TH1F(strnamept.Data(),strtitlept.Data(),1000,-1000.,1000.);
 236     fhd0D0pt[i]->SetXTitle("Impact parameter [#mum] ");
 237     fhd0D0pt[i]->SetYTitle("Entries");
 238
 239     strnamept.ReplaceAll("hd0D0","hMCd0D0");
 240     fhMCd0D0pt[i] = new TH1F(strnamept.Data(),strtitlept.Data(),1000,-1000.,1000.);
 241     fhMCd0D0pt[i]->SetXTitle("MC Impact parameter [#mum] ");
 242     fhMCd0D0pt[i]->SetYTitle("Entries");
 243
 244     strnamept.ReplaceAll("hMCd0D0","hd0D0VtxTrue");
 245     fhd0D0VtxTruept[i] = new TH1F(strnamept.Data(),strtitlept.Data(),1000,-1000.,1000.);
 246     fhd0D0VtxTruept[i]->SetXTitle("Impact parameter w.r.t. True Vtx [#mum] ");
 247     fhd0D0VtxTruept[i]->SetYTitle("Entries");
 248
 249   }
 250
 251
 252   /*  fHistPt = new TH1F("fHistPt", "P_{T} distribution", 15, 0.1, 3.1);
 253   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 254   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 255   fHistPt->SetMarkerStyle(kFullCircle);*/
 256 }
 257
 258 //________________________________________________________________________
 259 void AliAnalysisTaskCharmFraction::Exec(Option_t *)
 260 {
 261   // Main loop
 262   // Called for each event
 263
 264   if (!fAOD) {
 265     Printf("ERROR: fAOD not available");
 266     return;
 267   }
 268
 269   //Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
 270   //  Int_t nTotHF=0,nTotDstar=0,nTot3Prong=0;
 271   Int_t nTotD0toKpi=0;
 272   //  const  Int_t nptbins=10;
 273   Double_t standardptbins[14]={0.,0.5,1.,2.,3.,5.,7.,10.,15.,20.,30.,40.,50.,100.};
 274   Double_t pD[3],xD[3],pXtrTrue[2],pYtrTrue[2],pZtrTrue[2],xtr1[3],xtr2[3];
 275   Double_t cutsD0[9]=
 276         // cutsD0[0] = inv. mass half width [GeV]
 277         // cutsD0[1] = dca [cm]
 278         // cutsD0[2] = cosThetaStar
 279         // cutsD0[3] = pTK [GeV/c]
 280         // cutsD0[4] = pTPi [GeV/c]
 281         // cutsD0[5] = d0K [cm]   upper limit!
 282         // cutsD0[6] = d0Pi [cm]  upper limit!
 283         // cutsD0[7] = d0d0 [cm^2]
 284         // cutsD0[8] = cosThetaPoint
 285     {1000.,
 286      100000.,
 287      1.1,
 288      0.,
 289      0.,
 290      100000.,
 291      100000.,
 292      100000000.,
 293      -1.1};
 294   if(fD0usecuts){
 295     cutsD0[0]=0.033;
 296     cutsD0[1]=0.03;
 297     cutsD0[2]=0.8;
 298     cutsD0[3]=0.;
 299     cutsD0[4]=0.;
 300     cutsD0[5]=100000.;
 301     cutsD0[6]=100000.;
 302     cutsD0[7]=-0.00005;
 303     cutsD0[8]=0.8;
 304   }
 305   const Double_t fixedInvMassCut=cutsD0[0];
 306   AliAODVertex *vtx1=0;
 307   // primary vertex
 308   vtx1 = (AliAODVertex*)fAOD->GetPrimaryVertex();
 309   Double_t vtxTrue[3];
 310   fAODmcHeader->GetVertex(vtxTrue);
 311   //    vtx1->Print();
 312   AliAODRecoDecayHF *aodDMC=0x0;
 313   Bool_t nomum=kFALSE,background=kFALSE;
 314   AliAODMCParticle *mum1=0x0,*mum2=0x0;
 315   AliAODMCParticle *b1=0x0,*b2=0x0;
 316   AliAODMCParticle *sonpart=0x0,*grandmoth1=0x0;
 317
 318   // make trkIDtoEntry register (temporary)
 319   Int_t trkIDtoEntry[100000];
 320   for(Int_t it=0;it<fAOD->GetNumberOfTracks();it++) {
 321     AliAODTrack *track = fAOD->GetTrack(it);
 322     if(track->GetID()<0) {
 323       //printf("Track ID <0, id= %d\n",track->GetID());
 324       return;
 325     }
 326     trkIDtoEntry[track->GetID()]=it;
 327   }
 328
 329
 330   // loop over D0->Kpi candidates
 331   Int_t nD0toKpi = fArrayD0toKpi->GetEntriesFast();
 332   nTotD0toKpi += nD0toKpi;
 333   //    cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
 334
 335   for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
 336     if(aodDMC!=0x0)delete aodDMC;
 337     mum1=0x0;
 338     mum2=0x0;
 339     b1=0x0;
 340     b2=0x0;
 341     sonpart=0x0;
 342     grandmoth1=0x0;
 343
 344     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)fArrayD0toKpi->UncheckedAt(iD0toKpi);
 345     Bool_t unsetvtx=kFALSE;
 346     if(!d->GetOwnPrimaryVtx()) {
 347       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 348       unsetvtx=kTRUE;
 349     }
 350
 351     //    d->SetOwnPrimaryVtx(vtx1);
 352     //unsetvtx=kTRUE;
 353
 354     Int_t okD0=0,okD0bar=0;
 355     Double_t invMassD0,invMassD0bar,cutmassvalue=-1.;
 356     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 357     nomum=kFALSE;
 358     background=kFALSE;
 359
 360     cutsD0[0]=fixedInvMassCut;
 361     if(fSideBands){
 362       cutmassvalue=cutsD0[0];
 363       if(fSideBands==999.)cutsD0[0]=1000.;
 364       else cutsD0[0]=2.*TMath::Abs(fSideBands)*cutsD0[0];
 365     }
 366     if(d->SelectD0(cutsD0,okD0,okD0bar)) {
 367       if(fSideBands){
 368         if(fSideBands==999.)cutmassvalue=-1.;
 369         d->InvMassD0(invMassD0,invMassD0bar);
 370
 371         if(TMath::Abs(invMassD0-mD0PDG)    > TMath::Abs(fSideBands)*cutmassvalue) okD0 = okD0 && kTRUE;
 372         else okD0=kFALSE;
 373         if(TMath::Abs(invMassD0bar-mD0PDG) > TMath::Abs(fSideBands)*cutmassvalue) okD0bar = okD0bar && kTRUE;    // SIDE BANDS SELECTION
 374         else okD0bar=kFALSE;
 375
 376         if(fSideBands<0.){
 377           if(!okD0 && !okD0bar) continue;
 378         }
 379         else if (!okD0 || !okD0bar)continue;
 380       }
 381
 382
 383       // get daughter AOD tracks
 384       AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
 385       AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
 386
 387       //Check if it's a True D0
 388       if(fcheckMC){
 389         while(!background){
 390
 391           if(trk0->GetLabel()==-1||trk1->GetLabel()==-1){
 392             //fake tracks!!
 393             background=kTRUE;  //FAKE TRACK
 394             break;
 395           }
 396           //      printf("Before entering the MC checks \n");
 397
 398           b1=(AliAODMCParticle*)fArrayMC->At(trk0->GetLabel());
 399           b2=(AliAODMCParticle*)fArrayMC->At(trk1->GetLabel());
 400
 401
 402           if(b1->GetMother()==-1||b2->GetMother()==-1){
 403             nomum=kTRUE; //Tracks with no mother
 404             //printf("Inside problem mothering \n");// FAKE DECAY VERTEX
 405             background=kTRUE;
 406             break;
 407             // if(!fStudyPureBackground)continue;
 408           }
 409
 410           if(b1->GetMother()!=b2->GetMother()){
 411             //Check the label of the mother is the same
 412             background=kTRUE; // NOT SAME MOTHER
 413             break;
 414           }
 415           //      printf("After first mothering \n");
 416           mum1=(AliAODMCParticle*)fArrayMC->At(b1->GetMother());
 417           mum2=(AliAODMCParticle*)fArrayMC->At(b2->GetMother());
 418           //    if((mum1->GetPdgCode()!=mum2->GetPdgCode()))continue; //Check the mother is the same particle
 419           // printf("After second mothering \n");
 420           //      printf("Particle codes: tr1: %d, tr2: %d, mum1: %d, mum 2: %d \n",b1->GetPdgCode(),b2->GetPdgCode(),mum1->GetPdgCode(),mum2->GetPdgCode());
 421           if(fcheckMCD0)if(TMath::Abs(mum1->GetPdgCode())!=421){
 422             //^fStudyPureBackground)continue;//NPRONG
 423             background=kTRUE;// NOT A D0
 424             break;
 425           }
 426           if(fcheckMC2prongs)if(mum1->GetDaughter(1)-mum1->GetDaughter(0)+1!=2){
 427             //^fStudyPureBackground)continue;// NPRONG
 428             background=kTRUE; // NOT A 2 PRONG DECAY
 429             break;
 430           }
 431
 432           if(mum1->GetMother()==-1){
 433             background=kTRUE; // A particle coming from nothing
 434             break;
 435           }
 436
 437           // ^fStudyPureBackground)continue;
 438           //      printf("After second mothering \n");
 439           grandmoth1=(AliAODMCParticle*)fArrayMC->At(mum1->GetMother());
 440           if(fSkipD0star)if(TMath::Abs(grandmoth1->GetPdgCode())==413||TMath::Abs(grandmoth1->GetPdgCode())==423){
 441             //^fStudyPureBackground)continue;
 442             background=kTRUE;// D0 COMING FROM A D0*
 443             break;
 444           }
 445           if(fcheckMCD0){//THIS CHECK IS NEEDE TO AVOID POSSIBLE FAILURE IN THE SECOND WHILE
 446             while(TMath::Abs(grandmoth1->GetPdgCode())!=4&&TMath::Abs(grandmoth1->GetPdgCode())!=5){
 447               if(grandmoth1->GetMother()==-1){
 448                 //            printf("Inside grandmothering \n");
 449                 //printf("mother=-1, pdgcode: %d \n",grandmoth1->GetPdgCode());
 450                 Int_t son=grandmoth1->GetDaughter(0);
 451                 sonpart=(AliAODMCParticle*)fArrayMC->At(son);
 452                 while(TMath::Abs(sonpart->GetPdgCode())!=421){
 453                   //printf("mother=-1, pdgcode: %d \n",sonpart->GetPdgCode());
 454                   son++;
 455                   sonpart=(AliAODMCParticle*)fArrayMC->At(son);
 456                 }
 457                 break;
 458               }
 459               grandmoth1=(AliAODMCParticle*)fArrayMC->At(grandmoth1->GetMother());
 460             }
 461           }
 462           if(fcheckMCprompt)if(TMath::Abs(grandmoth1->GetPdgCode())!=4){
 463             //^fStudyPureBackground)continue;
 464             background=kTRUE;
 465             break;
 466           }
 467           if(fcheckMCfromB){
 468             if(TMath::Abs(grandmoth1->GetPdgCode())!=5){
 469               //^fStudyPureBackground)continue;
 470               background=kTRUE;
 471               break;
 472             }
 473           }
 474           //else if(!fStudyPureBackground)continue;
 475           break;
 476         }
 477         if(background^fStudyPureBackground)continue;
 478         if(fStudyd0fromBTrue){
 479           if(b1!=0x0&&b2!=0x0){
 480             Int_t charge[2]={0,0};
 481             if(b1->Charge()==-1)charge[0]=1;
 482             else {
 483               if(b2->Charge()==-1){
 484                 //printf("Same charges for prongs \n");
 485                 continue;
 486               }
 487               charge[1]=1;
 488             }
 489             /* if(!pXtrTrue[charge[0]]=b1->Px()){printf("!first MCtrack:Get momentum X failed \n");continue;}
 490                if(!b1->Py(pYtrTrue[charge[0]])){printf("!first MCtrack:Get momentum Y failed \n");continue;}
 491                if(!b1->Pz(pZtrTrue[charge[0]])){printf("!first MCtrack:Get momentum Z failed \n");continue;}
 492             */
 493
 494             pXtrTrue[charge[0]]=b1->Px();
 495             pYtrTrue[charge[0]]=b1->Py();
 496             pZtrTrue[charge[0]]=b1->Pz();
 497             if(!b1->XvYvZv(xtr1)){
 498               //printf("!first MCtrack:Get position failed \n");
 499               continue;
 500             }
 501             /*if(!b2->Px(pXtrTrue[charge[1]])){printf("!second MCtrack:Get momentum X failed \n");continue;}
 502               if(!b2->Py(pYtrTrue[charge[1]])){printf("!second MCtrack:Get momentum Y failed \n");continue;}
 503               if(!b2->Pz(pZtrTrue[charge[1]])){printf("!second MCtrack:Get momentum Z failed \n");continue;}*/
 504             pXtrTrue[charge[1]]=b2->Px();
 505             pYtrTrue[charge[1]]=b2->Py();
 506             pZtrTrue[charge[1]]=b2->Pz();
 507             if(!nomum){
 508               if(!(mum1==0x0||mum2==0x0)){
 509                 if(!mum1->PxPyPz(pD)){
 510                   //printf("!D from B:Get momentum failed \n");
 511                   continue;
 512                 }
 513                 if(!mum1->XvYvZv(xD)){
 514                   //printf("!D from B:Get position failed \n");
 515                   continue;
 516                 }
 517                 if(pXtrTrue[0]+pXtrTrue[1]!=pD[0]){
 518                   //printf("Warning! MCinfo: Pt of the mother doesn't match the sum of the daughters pt: X component  %f + %f = %f Vs %f \n",pXtrTrue[0],pXtrTrue[1],pXtrTrue[0]+pXtrTrue[1],pD[0]);
 519                   //printf("Warning! MCinfo: Pt of the mother doesn't match the sum of the daughters pt: Y componente %f + %f = %f Vs %f \n",pYtrTrue[0],pYtrTrue[1],pYtrTrue[0]+pYtrTrue[1],pD[1]);
 520                   //printf("Warning!: Pt comparison from the sum of the daughters:%f Vs mother pT: %f \n",TMath::Sqrt((pXtrTrue[0]+pXtrTrue[1])*(pXtrTrue[0]+pXtrTrue[1])+(pYtrTrue[0]+pYtrTrue[1])*(pYtrTrue[0]+pYtrTrue[1])),TMath::Sqrt(pD[0]*pD[0]+pD[1]*pD[1]));
 521                   //          return;
 522                 }
 523
 524                 //          if(!b2->PxPyPz(ptr2True)){printf("!second MCtrack:Get momentum failed \n");continue;}
 525                 if(!b2->XvYvZv(xtr2)){
 526                   //printf("!second MCtrack:Get position failed \n");
 527                   continue;
 528                 }
 529                 Double_t d0dummy[2]={0.,0.};//TEMPORARY : dummy d0 for AliAODRecoDeay constructor
 530                 aodDMC=new AliAODRecoDecayHF(vtxTrue,xD,2,0,pXtrTrue,pYtrTrue,pZtrTrue,d0dummy);
 531                 fhMassTrue->Fill(mum1->GetCalcMass());
 532               }
 533               //else printf("Warning! The candidate comes from two track with a different mother, \n -> the true mother d0 doesn't exist ! -> There will be different entries between observed and MC d0distributions \n");
 534             }
 535           }
 536         }
 537         //printf("End of MC checks \n");
 538       }
 539
 540       //printf("A real D0 was found!! \n");
 541
 542
 543       //cout<<1e8*d->Prodd0d0()<<endl;
 544       /*      if(fSideBands){   // The integral of the Hists will in general differ
 545                         // from the number of candidates which passed the selection
 546         if(okD0)fhMass->Fill(d->InvMassD0(),1.);
 547         if(okD0bar)fhMass->Fill(d->InvMassD0bar(),1.);
 548       }
 549       else { // The integral of the Hists will in general differ
 550         // from the number of candidates which passed the selection
 551         if(okD0)fhMass->Fill(d->InvMassD0(),1.);
 552         if(okD0bar)fhMass->Fill(d->InvMassD0bar(),1.);
 553       }
 554       */
 555
 556       if(okD0)fhMass->Fill(d->InvMassD0(),1.);
 557       if(okD0bar)fhMass->Fill(d->InvMassD0bar(),1.);
 558
 559       fhCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
 560       fhSecVtxZ->Fill(d->GetSecVtxZ());
 561       fhSecVtxX->Fill(d->GetSecVtxX()*10000.);
 562       fhSecVtxY->Fill(d->GetSecVtxY()*10000.);
 563       fhSecVtxXY->Fill(d->GetSecVtxX()*10000.,d->GetSecVtxY()*10000.);
 564       fhSecVtxPhi->Fill(TMath::ATan2(d->GetSecVtxY()*10000.,d->GetSecVtxX()*10000.)*TMath::RadToDeg());
 565       fhd0xd0->Fill(1e8*d->Prodd0d0());
 566       fhCPta->Fill(d->CosPointingAngle());
 567
 568
 569       //cout<<d->GetSecVtxX() <<endl;
 570
 571       if(!trk0 || !trk1) {
 572         /*      if(d->GetProngID(0)<0||d->GetProngID(1)<0){
 573                 printf("Wrong ProngID,%d or %d\n",d->GetProngID(0),d->GetProngID(1));continue;
 574                 }
 575         */
 576         trk0=fAOD->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
 577         trk1=fAOD->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
 578         //            cout<<"references to standard AOD not available"<<endl;
 579       }
 580       //            cout<<"pt of positive track: "<<trk0->Pt()<<endl;
 581
 582       // make a AliExternalTrackParam from the D0
 583       // and calculate impact parameters to primary vertex
 584
 585
 586       //        AliExternalTrackParam trackD0(d);
 587       //cout<<"pt of D0: "<<d->Pt()<<" (AliAODRecoDecay); "<<trackD0.Pt()<<" (track)"<<endl;
 588       //trackD0.Print();
 589       Double_t ptD0=d->Pt();
 590       for(Int_t j=0;j<fnbins;j++){
 591         if(ptD0<standardptbins[j+1]){
 592           fhd0D0VtxTruept[j]->Fill(d->AliAODRecoDecay::ImpParXY(vtxTrue)*10000.);
 593           fhd0D0pt[j]->Fill(d->ImpParXY()*10000.);
 594           // printf("HERE in tha class \n");
 595           if(aodDMC!=0x0)fhMCd0D0pt[j]->Fill(aodDMC->ImpParXY()*10000.);
 596           // printf("HERE2 in tha class \n");
 597           break;
 598         }
 599       }
 600       fhd0D0->Fill(d->ImpParXY()*10000.);
 601       fhd0D0VtxTrue->Fill(d->AliAODRecoDecay::ImpParXY(vtxTrue)*10000.);
 602       if(aodDMC!=0x0)fhMCd0D0->Fill(aodDMC->ImpParXY()*10000.);
 603       //  if((1.864-3.*0.011.<d->InvMassD0()&&d->InvMassD0()<1.864+3.*0.011)||(1.864-3.*0.011<d->InvMassD0bar()&&d->InvMassD0bar()<1.864+3.*0.011)){
 604       if(aodDMC!=0x0){
 605         delete aodDMC;
 606         aodDMC=0x0;
 607       }
 608       mum1=0x0;
 609       mum2=0x0;
 610       b1=0x0;
 611       b2=0x0;
 612       sonpart=0x0;
 613       grandmoth1=0x0;
 614
 615       //            }
 616       //        Double_t dz[2],covdz[3];
 617       //trackD0.PropagateToDCA(vtx1,fAOD->GetMagneticField(),1000.,dz,covdz);
 618       //        cout<<"D0 impact parameter rphi: "<<dz[0]<<" +- "<<TMath::Sqrt(covdz[0])<<endl;
 619     }
 620     // printf("Before end of the event \n");
 621     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 622
 623   }
 624
 625
 626
 627   // Post output data.
 628
 629   PostData(0,fhCPtaVSd0d0);
 630   PostData(1,fhSecVtxXY);
 631   PostData(2,fhd0xd0);
 632   PostData(3,fhCPta);
 633   PostData(4,fhSecVtxZ);
 634   PostData(5,fhSecVtxX);
 635   PostData(6,fhSecVtxY);
 636   PostData(7,fhSecVtxPhi);
 637
 638
 639   for(Int_t j=0;j<fnbins;j++){
 640
 641     PostData(j+8, fhd0D0pt[j]);
 642   }
 643
 644   PostData(fnbins+8, fhd0D0);
 645
 646   for(Int_t j=0;j<fnbins;j++){
 647
 648     PostData(j+fnbins+9,fhMCd0D0pt[j]);
 649   }
 650
 651   PostData(2*fnbins+9, fhMCd0D0);
 652
 653   for(Int_t j=0;j<fnbins;j++){
 654
 655     PostData(j+2*fnbins+10,fhd0D0VtxTruept[j]);
 656   }
 657
 658   PostData(3*fnbins+10,fhd0D0VtxTrue);
 659   PostData(3*fnbins+11,fhMassTrue);
 660   PostData(3*fnbins+12,fhMass);
 661 }
 662
 663 //________________________________________________________________________
 664 void AliAnalysisTaskCharmFraction::Terminate(Option_t *)
 665 {
 666   // Draw result to the screen
 667   // Called once at the end of the query
 668
 669   fhd0D0 = dynamic_cast<TH1F*> (GetOutputData(fnbins+8));
 670   if (!fhd0D0) {
 671     Printf("ERROR: fhd0D0 not available");
 672     return;
 673   }
 674
 675   TCanvas *c1 = new TCanvas("AliAnalysisTaskCharmFraction","D0 impact par",10,10,510,510);
 676   c1->cd(1)->SetLogy();
 677   fhd0D0->DrawCopy("E");
 678
 679   TCanvas *cd0D0pt=new TCanvas("cd0D0pt","cd0D0pt",0,0,800,700);
 680   cd0D0pt->Divide(fnbins/2+fnbins%2,2);
 681   for(Int_t j=0;j<fnbins;j++){
 682     fhd0D0pt[j] = dynamic_cast<TH1F*> (GetOutputData(8+j));
 683     if (!fhd0D0pt[j]) {
 684       Printf("ERROR: fhd0D0pt not available");
 685       return;
 686     }
 687     cd0D0pt->cd(j+1);
 688     fhd0D0pt[j]->Draw();
 689   }
 690
 691 }