PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // AliAnalysisTaskSE for D0 candidates invariant mass histogram
  21 // and comparison with the MC truth and cut variables distributions.
  22 //
  23 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
  24 // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
  25 // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
  26 /////////////////////////////////////////////////////////////
  27
  28 #include <Riostream.h>
  29 #include <TClonesArray.h>
  30 #include <TCanvas.h>
  31 #include <TNtuple.h>
  32 #include <TTree.h>
  33 #include <TList.h>
  34 #include <TH1F.h>
  35 #include <TH2F.h>
  36 #include <TDatabasePDG.h>
  37
  38 #include <AliAnalysisDataSlot.h>
  39 #include <AliAnalysisDataContainer.h>
  40 #include "AliAnalysisManager.h"
  41 #include "AliESDtrack.h"
  42 #include "AliVertexerTracks.h"
  43 #include "AliAODHandler.h"
  44 #include "AliAODEvent.h"
  45 #include "AliAODVertex.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODMCHeader.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliAODRecoDecayHF2Prong.h"
  50 #include "AliAODRecoCascadeHF.h"
  51 #include "AliAnalysisVertexingHF.h"
  52 #include "AliAnalysisTaskSE.h"
  53 #include "AliAnalysisTaskSED0Mass.h"
  54 #include "AliNormalizationCounter.h"
  55
  56 using std::cout;
  57 using std::endl;
  58
  59 ClassImp(AliAnalysisTaskSED0Mass)
  60
  61
  62 //________________________________________________________________________
  63 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
  64 AliAnalysisTaskSE(),
  65   fOutputMass(0),
  66   fOutputMassPt(0),
  67   fDistr(0),
  68   fNentries(0),
  69   fCuts(0),
  70   fArray(0),
  71   fReadMC(0),
  72   fCutOnDistr(0),
  73   fUsePid4Distr(0),
  74   fCounter(0),
  75   fNPtBins(1),
  76   fLsNormalization(1.),
  77   fFillOnlyD0D0bar(0),
  78   fDaughterTracks(),
  79   fIsSelectedCandidate(0),
  80   fFillVarHists(kTRUE),
  81   fSys(0),
  82   fIsRejectSDDClusters(0),
  83   fFillPtHist(kTRUE),
  84   fFillImpParHist(kFALSE),
  85   fWriteVariableTree(kFALSE),
  86   fVariablesTree(0),
  87   fCandidateVariables()
  88 {
  89   // Default constructor
  90   for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
  91
  92 }
  93
  94 //________________________________________________________________________
  95 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
  96   AliAnalysisTaskSE(name),
  97   fOutputMass(0),
  98   fOutputMassPt(0),
  99   fDistr(0),
 100   fNentries(0),
 101   fCuts(0),
 102   fArray(0),
 103   fReadMC(0),
 104   fCutOnDistr(0),
 105   fUsePid4Distr(0),
 106   fCounter(0),
 107   fNPtBins(1),
 108   fLsNormalization(1.),
 109   fFillOnlyD0D0bar(0),
 110   fDaughterTracks(),
 111   fIsSelectedCandidate(0),
 112   fFillVarHists(kTRUE),
 113   fSys(0),
 114   fIsRejectSDDClusters(0),
 115   fFillPtHist(kTRUE),
 116   fFillImpParHist(kFALSE),
 117   fWriteVariableTree(kFALSE),
 118   fVariablesTree(0),
 119   fCandidateVariables()
 120 {
 121   // Default constructor
 122
 123   fNPtBins=cuts->GetNPtBins();
 124
 125   fCuts=cuts;
 126   for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
 127
 128   // Output slot #1 writes into a TList container (mass with cuts)
 129   DefineOutput(1,TList::Class());  //My private output
 130   // Output slot #2 writes into a TList container (distributions)
 131   DefineOutput(2,TList::Class());  //My private output
 132   // Output slot #3 writes into a TH1F container (number of events)
 133   DefineOutput(3,TH1F::Class());  //My private output
 134   // Output slot #4 writes into a TList container (cuts)
 135   DefineOutput(4,AliRDHFCutsD0toKpi::Class());  //My private output
 136   // Output slot #5 writes Normalization Counter
 137   DefineOutput(5,AliNormalizationCounter::Class());
 138   // Output slot #6 stores the mass vs pt and impact parameter distributions
 139   DefineOutput(6,TList::Class());  //My private output
 140   // Output slot #7 keeps a tree of the candidate variables after track selection
 141   DefineOutput(7,TTree::Class());  //My private outpu
 142 }
 143
 144 //________________________________________________________________________
 145 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
 146 {
 147   if (fOutputMass) {
 148     delete fOutputMass;
 149     fOutputMass = 0;
 150   }
 151   if (fOutputMassPt) {
 152     delete fOutputMassPt;
 153     fOutputMassPt = 0;
 154   }
 155   if (fDistr) {
 156     delete fDistr;
 157     fDistr = 0;
 158   }
 159   if (fCuts) {
 160     delete fCuts;
 161     fCuts = 0;
 162   }
 163   for(Int_t i=0; i<5; i++){
 164     if(fHistMassPtImpParTC[i]) delete fHistMassPtImpParTC[i];
 165   }
 166   if (fNentries){
 167     delete fNentries;
 168     fNentries = 0;
 169   }
 170   if(fCounter){
 171     delete fCounter;
 172     fCounter=0;
 173   }
 174   if(fVariablesTree){
 175     delete fVariablesTree;
 176     fVariablesTree = 0;
 177   }
 178
 179 }
 180
 181 //________________________________________________________________________
 182 void AliAnalysisTaskSED0Mass::Init()
 183 {
 184   // Initialization
 185
 186   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
 187
 188
 189   AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
 190   const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 191   copyfCuts->SetName(nameoutput);
 192   // Post the data
 193   PostData(4,copyfCuts);
 194
 195
 196   return;
 197 }
 198
 199 //________________________________________________________________________
 200 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
 201 {
 202
 203   // Create the output container
 204   //
 205   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
 206
 207   // Several histograms are more conveniently managed in a TList
 208   fOutputMass = new TList();
 209   fOutputMass->SetOwner();
 210   fOutputMass->SetName("listMass");
 211
 212   fOutputMassPt = new TList();
 213   fOutputMassPt->SetOwner();
 214   fOutputMassPt->SetName("listMassPt");
 215
 216   fDistr = new TList();
 217   fDistr->SetOwner();
 218   fDistr->SetName("distributionslist");
 219
 220   TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
 221   TString nameMassPt="", nameSgnPt="", nameBkgPt="", nameRflPt="";
 222   Int_t nbins2dPt=60; Float_t binInPt=0., binFinPt=30.;
 223
 224   for(Int_t i=0;i<fCuts->GetNPtBins();i++){
 225
 226     nameMass="histMass_";
 227     nameMass+=i;
 228     nameSgn27="histSgn27_";
 229     nameSgn27+=i;
 230     nameSgn="histSgn_";
 231     nameSgn+=i;
 232     nameBkg="histBkg_";
 233     nameBkg+=i;
 234     nameRfl="histRfl_";
 235     nameRfl+=i;
 236     nameMassNocutsS="hMassS_";
 237     nameMassNocutsS+=i;
 238     nameMassNocutsB="hMassB_";
 239     nameMassNocutsB+=i;
 240
 241     //histograms of cut variable distributions
 242     if(fFillVarHists){
 243       if(fReadMC){
 244         //  dca
 245         namedistr="hdcaS_";
 246         namedistr+=i;
 247         TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 248         // impact parameter
 249         namedistr="hd0piS_";
 250         namedistr+=i;
 251         TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
 252
 253         namedistr="hd0KS_";
 254         namedistr+=i;
 255         TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
 256         namedistr="hd0d0S_";
 257         namedistr+=i;
 258         TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 259
 260         //decay lenght
 261         namedistr="hdeclS_";
 262         namedistr+=i;
 263         TH1F *hdeclengthS = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm]",200,0,0.015);
 264
 265         namedistr="hnormdeclS_";
 266         namedistr+=i;
 267         TH1F *hnormdeclengthS = new TH1F(namedistr.Data(), "Normalized Decay Length^{2} distribution;(Decay Length/Err)^{2} ",200,0,12.);
 268
 269         namedistr="hdeclxyS_";
 270         namedistr+=i;
 271         TH1F* hdeclxyS=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
 272         namedistr="hnormdeclxyS_";
 273         namedistr+=i;
 274         TH1F* hnormdeclxyS=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
 275
 276         namedistr="hdeclxyd0d0S_";
 277         namedistr+=i;
 278         TH2F* hdeclxyd0d0S=new TH2F(namedistr.Data(),"Correlation decay Length XY - d_{0}#times d_{0};Decay Length XY [cm];d_{0}#times d_{0}[cm^{2}]",200,0,0.15,200,-0.001,0.001);
 279
 280         namedistr="hnormdeclxyd0d0S_";
 281         namedistr+=i;
 282         TH2F* hnormdeclxyd0d0S=new TH2F(namedistr.Data(),"Correlation normalized decay Length XY - d_{0}#times d_{0};Decay Length XY/Err;d_{0}#times d_{0}[cm^{2}]",200,0,6,200,-0.001,0.001);
 283
 284         //  costhetapoint
 285         namedistr="hcosthetapointS_";
 286         namedistr+=i;
 287         TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 288
 289         namedistr="hcosthetapointxyS_";
 290         namedistr+=i;
 291         TH1F *hcosthetapointxyS = new TH1F(namedistr.Data(), "cos#theta_{Point} XYdistribution;cos#theta_{Point}",300,0.,1.);
 292
 293         TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
 294         tmpMS->Sumw2();
 295
 296
 297         fDistr->Add(hdcaS);
 298
 299         fDistr->Add(hd0piS);
 300         fDistr->Add(hd0KS);
 301
 302         fDistr->Add(hd0d0S);
 303
 304         fDistr->Add(hcosthetapointS);
 305
 306         fDistr->Add(hcosthetapointxyS);
 307
 308         fDistr->Add(hdeclengthS);
 309
 310         fDistr->Add(hnormdeclengthS);
 311
 312         fDistr->Add(hdeclxyS);
 313
 314         fDistr->Add(hnormdeclxyS);
 315
 316         fDistr->Add(hdeclxyd0d0S);
 317         fDistr->Add(hnormdeclxyd0d0S);
 318
 319         fDistr->Add(tmpMS);
 320       }
 321
 322       //  dca
 323       namedistr="hdcaB_";
 324       namedistr+=i;
 325       TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 326
 327       // impact parameter
 328       namedistr="hd0B_";
 329       namedistr+=i;
 330       TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
 331
 332       namedistr="hd0d0B_";
 333       namedistr+=i;
 334       TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 335
 336       //decay lenght
 337       namedistr="hdeclB_";
 338       namedistr+=i;
 339       TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm^{2}]",200,0,0.015);
 340
 341       namedistr="hnormdeclB_";
 342       namedistr+=i;
 343       TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;(Decay Length/Err)^{2} ",200,0,12.);
 344
 345       namedistr="hdeclxyB_";
 346       namedistr+=i;
 347       TH1F* hdeclxyB=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
 348       namedistr="hnormdeclxyB_";
 349       namedistr+=i;
 350       TH1F* hnormdeclxyB=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
 351
 352       namedistr="hdeclxyd0d0B_";
 353       namedistr+=i;
 354       TH2F* hdeclxyd0d0B=new TH2F(namedistr.Data(),"Correlation decay Length XY - d_{0}#times d_{0};Decay Length XY [cm];d_{0}#times d_{0}[cm^{2}]",200,0,0.15,200,-0.001,0.001);
 355
 356       namedistr="hnormdeclxyd0d0B_";
 357       namedistr+=i;
 358       TH2F* hnormdeclxyd0d0B=new TH2F(namedistr.Data(),"Correlation normalized decay Length XY - d_{0}#times d_{0};Decay Length XY/Err;d_{0}#times d_{0}[cm^{2}]",200,0,6,200,-0.001,0.001);
 359
 360       //  costhetapoint
 361       namedistr="hcosthetapointB_";
 362       namedistr+=i;
 363       TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 364
 365       namedistr="hcosthetapointxyB_";
 366       namedistr+=i;
 367       TH1F *hcosthetapointxyB = new TH1F(namedistr.Data(), "cos#theta_{Point} XY distribution;cos#theta_{Point} XY",300,0.,1.);
 368
 369       TH1F* tmpMB = new TH1F(nameMassNocutsB.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
 370       tmpMB->Sumw2();
 371
 372
 373
 374       fDistr->Add(hdcaB);
 375
 376       fDistr->Add(hd0B);
 377
 378       fDistr->Add(hd0d0B);
 379
 380       fDistr->Add(hcosthetapointB);
 381
 382       fDistr->Add(hcosthetapointxyB);
 383
 384       fDistr->Add(hdeclengthB);
 385
 386       fDistr->Add(hnormdeclengthB);
 387
 388       fDistr->Add(hdeclxyB);
 389
 390       fDistr->Add(hnormdeclxyB);
 391
 392       fDistr->Add(hdeclxyd0d0B);
 393       fDistr->Add(hnormdeclxyd0d0B);
 394
 395       fDistr->Add(tmpMB);
 396
 397       //histograms filled only when the secondary vertex is recalculated w/o the daughter tracks (as requested in the cut object)
 398
 399       if(fCuts->GetIsPrimaryWithoutDaughters()){
 400         if(fReadMC){
 401           namedistr="hd0vpiS_";
 402           namedistr+=i;
 403           TH1F *hd0vpiS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions)(vtx w/o these tracks);d0(#pi) [cm]",200,-0.1,0.1);
 404
 405           namedistr="hd0vKS_";
 406           namedistr+=i;
 407           TH1F *hd0vKS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons) (vtx w/o these tracks);d0(K) [cm]",200,-0.1,0.1);
 408
 409           namedistr="hd0d0vS_";
 410           namedistr+=i;
 411           TH1F *hd0d0vS = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution (vtx w/o these tracks);d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 412
 413           namedistr="hdeclvS_";
 414           namedistr+=i;
 415           TH1F *hdeclengthvS = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
 416
 417           namedistr="hnormdeclvS_";
 418           namedistr+=i;
 419           TH1F *hnormdeclengthvS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
 420
 421           fDistr->Add(hd0vpiS);
 422           fDistr->Add(hd0vKS);
 423           fDistr->Add(hd0d0vS);
 424           fDistr->Add(hdeclengthvS);
 425           fDistr->Add(hnormdeclengthvS);
 426
 427         }
 428
 429         namedistr="hd0vmoresB_";
 430         namedistr+=i;
 431         TH1F *hd0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
 432
 433         namedistr="hd0d0vmoresB_";
 434         namedistr+=i;
 435         TH1F *hd0d0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
 436
 437
 438         namedistr="hd0vB_";
 439         namedistr+=i;
 440         TH1F *hd0vB = new TH1F(namedistr.Data(), "Impact parameter distribution (vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
 441
 442         namedistr="hd0vp0B_";
 443         namedistr+=i;
 444         TH1F *hd0vp0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong + ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
 445
 446         namedistr="hd0vp1B_";
 447         namedistr+=i;
 448         TH1F *hd0vp1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong - ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
 449
 450         namedistr="hd0d0vB_";
 451         namedistr+=i;
 452         TH1F *hd0d0vB = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution (vtx w/o these tracks);d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 453
 454         namedistr="hdeclvB_";
 455         namedistr+=i;
 456         TH1F *hdeclengthvB = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
 457
 458         namedistr="hnormdeclvB_";
 459         namedistr+=i;
 460         TH1F *hnormdeclengthvB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
 461
 462         fDistr->Add(hd0vB);
 463         fDistr->Add(hd0vp0B);
 464         fDistr->Add(hd0vp1B);
 465         fDistr->Add(hd0vmoresB);
 466
 467         fDistr->Add(hd0d0vB);
 468         fDistr->Add(hd0d0vmoresB);
 469
 470         fDistr->Add(hdeclengthvB);
 471
 472         fDistr->Add(hnormdeclengthvB);
 473       }
 474
 475
 476     }
 477     //histograms of invariant mass distributions
 478
 479
 480     //MC signal
 481     if(fReadMC){
 482       TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.5648,2.1648);
 483
 484       TH1F *tmpSl=(TH1F*)tmpSt->Clone();
 485       tmpSt->Sumw2();
 486       tmpSl->Sumw2();
 487
 488       //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 489       TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.5648,2.1648);
 490       //TH1F *tmpRl=(TH1F*)tmpRt->Clone();
 491       TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.5648,2.1648);
 492       //TH1F *tmpBl=(TH1F*)tmpBt->Clone();
 493       tmpBt->Sumw2();
 494       //tmpBl->Sumw2();
 495       tmpRt->Sumw2();
 496       //tmpRl->Sumw2();
 497
 498       fOutputMass->Add(tmpSt);
 499       fOutputMass->Add(tmpRt);
 500       fOutputMass->Add(tmpBt);
 501
 502     }
 503
 504     //mass
 505     TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.5648,2.1648);
 506     //TH1F *tmpMl=(TH1F*)tmpMt->Clone();
 507     tmpMt->Sumw2();
 508     //tmpMl->Sumw2();
 509     //distribution w/o cuts large range
 510     // TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
 511
 512     fOutputMass->Add(tmpMt);
 513
 514     if(fSys==0){ //histograms filled only in pp to save time in PbPb
 515       if(fFillVarHists){
 516
 517         if(fReadMC){
 518           //  pT
 519           namedistr="hptpiS_";
 520           namedistr+=i;
 521           TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
 522
 523           namedistr="hptKS_";
 524           namedistr+=i;
 525           TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
 526
 527           //  costhetastar
 528           namedistr="hcosthetastarS_";
 529           namedistr+=i;
 530           TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 531
 532           //pT no mass cut
 533
 534           namedistr="hptpiSnoMcut_";
 535           namedistr+=i;
 536           TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
 537
 538           namedistr="hptKSnoMcut_";
 539           namedistr+=i;
 540           TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
 541
 542           fDistr->Add(hptpiS);
 543           fDistr->Add(hptKS);
 544           fDistr->Add(hcosthetastarS);
 545
 546           fDistr->Add(hptpiSnoMcut);
 547           fDistr->Add(hptKSnoMcut);
 548
 549           //  costhetapoint vs d0 or d0d0
 550           namedistr="hcosthpointd0S_";
 551           namedistr+=i;
 552           TH2F *hcosthpointd0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
 553           namedistr="hcosthpointd0d0S_";
 554           namedistr+=i;
 555           TH2F *hcosthpointd0d0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
 556
 557           fDistr->Add(hcosthpointd0S);
 558           fDistr->Add(hcosthpointd0d0S);
 559
 560           //to compare with AliAnalysisTaskCharmFraction
 561           TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",200,1.5648,2.1648);
 562           TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
 563           tmpS27t->Sumw2();
 564           tmpS27l->Sumw2();
 565
 566           fOutputMass->Add(tmpS27t);
 567           fOutputMass->Add(tmpS27l);
 568
 569         }
 570
 571         //  pT
 572         namedistr="hptB_";
 573         namedistr+=i;
 574         TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 575
 576         //  costhetastar
 577         namedistr="hcosthetastarB_";
 578         namedistr+=i;
 579         TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 580
 581         //pT no mass cut
 582         namedistr="hptB1prongnoMcut_";
 583         namedistr+=i;
 584         TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 585
 586         namedistr="hptB2prongsnoMcut_";
 587         namedistr+=i;
 588         TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 589
 590         fDistr->Add(hptB);
 591         fDistr->Add(hcosthetastarB);
 592
 593         fDistr->Add(hptB1pnoMcut);
 594         fDistr->Add(hptB2pnoMcut);
 595
 596         //impact parameter of negative/positive track
 597         namedistr="hd0p0B_";
 598         namedistr+=i;
 599         TH1F *hd0p0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.1,0.1);
 600
 601         namedistr="hd0p1B_";
 602         namedistr+=i;
 603         TH1F *hd0p1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong -);d0 [cm]",200,-0.1,0.1);
 604
 605         //impact parameter corrected for strangeness
 606         namedistr="hd0moresB_";
 607         namedistr+=i;
 608         TH1F *hd0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
 609
 610         namedistr="hd0d0moresB_";
 611         namedistr+=i;
 612         TH1F *hd0d0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
 613
 614
 615         namedistr="hcosthetapointmoresB_";
 616         namedistr+=i;
 617         TH1F *hcosthetapointmoresB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 618
 619         //  costhetapoint vs d0 or d0d0
 620         namedistr="hcosthpointd0B_";
 621         namedistr+=i;
 622         TH2F *hcosthpointd0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
 623
 624         namedistr="hcosthpointd0d0B_";
 625         namedistr+=i;
 626         TH2F *hcosthpointd0d0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
 627
 628         fDistr->Add(hd0p0B);
 629         fDistr->Add(hd0p1B);
 630
 631         fDistr->Add(hd0moresB);
 632         fDistr->Add(hd0d0moresB);
 633         fDistr->Add(hcosthetapointmoresB);
 634
 635
 636         fDistr->Add(hcosthpointd0B);
 637
 638
 639         fDistr->Add(hcosthpointd0d0B);
 640       }
 641
 642     } //end pp histo
 643   }
 644
 645
 646   //for Like sign analysis
 647
 648   if(fArray==1){
 649     namedistr="hpospair";
 650     TH1F* hpospair=new TH1F(namedistr.Data(),"Number of positive pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
 651     namedistr="hnegpair";
 652     TH1F* hnegpair=new TH1F(namedistr.Data(),"Number of negative pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
 653     fOutputMass->Add(hpospair);
 654     fOutputMass->Add(hnegpair);
 655   }
 656
 657
 658   // 2D Pt distributions and impact parameter histograms
 659   if(fFillPtHist) {
 660
 661     nameMassPt="histMassPt";
 662     nameSgnPt="histSgnPt";
 663     nameBkgPt="histBkgPt";
 664     nameRflPt="histRflPt";
 665
 666     //MC signal
 667     if(fReadMC){
 668       TH2F* tmpStPt = new TH2F(nameSgnPt.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",200,1.5648,2.16484,nbins2dPt,binInPt,binFinPt);
 669       TH2F *tmpSlPt=(TH2F*)tmpStPt->Clone();
 670       tmpStPt->Sumw2();
 671       tmpSlPt->Sumw2();
 672
 673       //Reflection: histo filled with D0MassV1 which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 674       TH2F* tmpRtPt = new TH2F(nameRflPt.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",200,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
 675       TH2F* tmpBtPt = new TH2F(nameBkgPt.Data(), "Background invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",200,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
 676       tmpBtPt->Sumw2();
 677       tmpRtPt->Sumw2();
 678
 679       fOutputMassPt->Add(tmpStPt);
 680       fOutputMassPt->Add(tmpRtPt);
 681       fOutputMassPt->Add(tmpBtPt);
 682
 683       //       cout<<endl<<endl<<"***************************************"<<endl;
 684       //       cout << " created and added to the list "<< nameSgnPt.Data() <<" "<< tmpStPt <<
 685       //        ", "<<nameRflPt.Data() <<" " << tmpRtPt<<", "<<nameBkgPt.Data()<< " " << tmpBtPt <<endl;
 686       //       cout<<"***************************************"<<endl<<endl;
 687     }
 688
 689     TH2F* tmpMtPt = new TH2F(nameMassPt.Data(),"D^{0} invariant mass; M [GeV]; Entries; Pt[GeV/c]",200,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
 690     tmpMtPt->Sumw2();
 691
 692     fOutputMassPt->Add(tmpMtPt);
 693   }
 694
 695   if(fFillImpParHist) CreateImpactParameterHistos();
 696
 697   const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
 698
 699   fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex ***  Integral(5,6) = pt out of bounds", 18,-0.5,17.5);
 700
 701   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 702   fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
 703   if(fReadMC) fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
 704   else fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
 705   fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
 706   fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
 707   fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
 708   if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
 709   if(fFillVarHists){
 710     fNentries->GetXaxis()->SetBinLabel(8,"PID=0");
 711     fNentries->GetXaxis()->SetBinLabel(9,"PID=1");
 712     fNentries->GetXaxis()->SetBinLabel(10,"PID=2");
 713     fNentries->GetXaxis()->SetBinLabel(11,"PID=3");
 714   }
 715   if(fReadMC && fSys==0){
 716     fNentries->GetXaxis()->SetBinLabel(12,"K");
 717     fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
 718   }
 719   fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
 720   fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
 721   if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
 722   if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
 723   fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
 724   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
 725
 726   fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
 727   fCounter->Init();
 728
 729   //
 730   // Output slot 7 : tree of the candidate variables after track selection
 731   //
 732   nameoutput = GetOutputSlot(7)->GetContainer()->GetName();
 733   fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
 734   Int_t nVar = 15;
 735   fCandidateVariables = new Double_t [nVar];
 736   TString * fCandidateVariableNames = new TString[nVar];
 737   fCandidateVariableNames[0]="massD0"; fCandidateVariableNames[1]="massD0bar"; fCandidateVariableNames[2]="pt";
 738   fCandidateVariableNames[3]="dca"; fCandidateVariableNames[4]="costhsD0"; fCandidateVariableNames[5]="costhsD0bar";
 739   fCandidateVariableNames[6]="ptk"; fCandidateVariableNames[7]="ptpi";
 740   fCandidateVariableNames[8]="d0k"; fCandidateVariableNames[9]="d0pi"; fCandidateVariableNames[10]="d0xd0";
 741   fCandidateVariableNames[11]="costhp"; fCandidateVariableNames[12]="costhpxy";
 742   fCandidateVariableNames[13]="lxy";  fCandidateVariableNames[14]="specialcuts";
 743   for(Int_t ivar=0; ivar<nVar; ivar++){
 744     fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/d",fCandidateVariableNames[ivar].Data()));
 745   }
 746
 747   // Post the data
 748   PostData(1,fOutputMass);
 749   PostData(2,fDistr);
 750   PostData(3,fNentries);
 751   PostData(5,fCounter);
 752   PostData(6,fOutputMassPt);
 753   PostData(7,fVariablesTree);
 754
 755   return;
 756 }
 757
 758 //________________________________________________________________________
 759 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
 760 {
 761   // Execute analysis for current event:
 762   // heavy flavor candidates association to MC truth
 763   //cout<<"I'm in UserExec"<<endl;
 764
 765
 766   //cuts order
 767   //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 768   //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 769   //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 770   //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 771   //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 772   //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 773   //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 774   //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 775   //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 776
 777
 778   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 779
 780   TString bname;
 781   if(fArray==0){ //D0 candidates
 782     // load D0->Kpi candidates
 783     //cout<<"D0 candidates"<<endl;
 784     bname="D0toKpi";
 785   } else { //LikeSign candidates
 786     // load like sign candidates
 787     //cout<<"LS candidates"<<endl;
 788     bname="LikeSign2Prong";
 789   }
 790
 791   TClonesArray *inputArray=0;
 792
 793   if(!aod && AODEvent() && IsStandardAOD()) {
 794     // In case there is an AOD handler writing a standard AOD, use the AOD
 795     // event in memory rather than the input (ESD) event.
 796     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 797     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 798     // have to taken from the AOD event hold by the AliAODExtension
 799     AliAODHandler* aodHandler = (AliAODHandler*)
 800       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 801
 802     if(aodHandler->GetExtensions()) {
 803       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 804       AliAODEvent* aodFromExt = ext->GetAOD();
 805       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
 806     }
 807   } else if(aod) {
 808     inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
 809   }
 810
 811
 812   if(!inputArray || !aod) {
 813     printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
 814     return;
 815   }
 816
 817   // fix for temporary bug in ESDfilter
 818   // the AODs with null vertex pointer didn't pass the PhysSel
 819   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 820
 821
 822   TClonesArray *mcArray = 0;
 823   AliAODMCHeader *mcHeader = 0;
 824
 825   if(fReadMC) {
 826     // load MC particles
 827     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 828     if(!mcArray) {
 829       printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
 830       return;
 831     }
 832
 833     // load MC header
 834     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 835     if(!mcHeader) {
 836       printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
 837       return;
 838     }
 839   }
 840
 841   //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
 842
 843   //histogram filled with 1 for every AOD
 844   fNentries->Fill(0);
 845   fCounter->StoreEvent(aod,fCuts,fReadMC);
 846   //fCounter->StoreEvent(aod,fReadMC);
 847
 848   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
 849   TString trigclass=aod->GetFiredTriggerClasses();
 850   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
 851
 852   if(!fCuts->IsEventSelected(aod)) {
 853     if(fCuts->GetWhyRejection()==1) // rejected for pileup
 854       fNentries->Fill(13);
 855     if(fSys==1 && (fCuts->GetWhyRejection()==2 || fCuts->GetWhyRejection()==3)) fNentries->Fill(15);
 856     if(fCuts->GetWhyRejection()==7) fNentries->Fill(17);
 857     return;
 858   }
 859
 860   // Check the Nb of SDD clusters
 861   if (fIsRejectSDDClusters) {
 862     Bool_t skipEvent = kFALSE;
 863     Int_t ntracks = 0;
 864     if (aod) ntracks = aod->GetNTracks();
 865     for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
 866       //    ... get the track
 867       AliAODTrack * track = aod->GetTrack(itrack);
 868       if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
 869         skipEvent=kTRUE;
 870         fNentries->Fill(16);
 871         break;
 872       }
 873     }
 874     if (skipEvent) return;
 875   }
 876
 877   // AOD primary vertex
 878   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 879
 880   Bool_t isGoodVtx=kFALSE;
 881
 882   //vtx1->Print();
 883   TString primTitle = vtx1->GetTitle();
 884   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
 885     isGoodVtx=kTRUE;
 886     fNentries->Fill(3);
 887   }
 888
 889   // loop over candidates
 890   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 891   if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 892
 893   // FILE *f=fopen("4display.txt","a");
 894   // fprintf(f,"Number of D0->Kpi: %d\n",nInD0toKpi);
 895   Int_t nSelectedloose=0,nSelectedtight=0;
 896   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 897     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
 898
 899     if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
 900         fNentries->Fill(2);
 901         continue; //skip the D0 from Dstar
 902       }
 903
 904     // Bool_t unsetvtx=kFALSE;
 905     // if(!d->GetOwnPrimaryVtx()) {
 906     //   d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 907     //   unsetvtx=kTRUE;
 908     // }
 909
 910
 911     if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
 912       nSelectedloose++;
 913       nSelectedtight++;
 914       if(fSys==0){
 915         if(fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
 916       }
 917       Int_t ptbin=fCuts->PtBin(d->Pt());
 918       if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 919       fIsSelectedCandidate=fCuts->IsSelected(d,AliRDHFCuts::kAll,aod); //selected
 920       if(fFillVarHists) {
 921         //if(!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate)) {
 922         fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(0),0);
 923         fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(1),1);
 924         //check daughters
 925         if(!fDaughterTracks.UncheckedAt(0) || !fDaughterTracks.UncheckedAt(1)) {
 926           AliDebug(1,"at least one daughter not found!");
 927           fNentries->Fill(5);
 928           fDaughterTracks.Clear();
 929           continue;
 930         }
 931         //}
 932         FillVarHists(aod,d,mcArray,fCuts,fDistr);
 933       }
 934
 935       FillMassHists(d,mcArray,mcHeader,fCuts,fOutputMass);
 936     }
 937
 938     fDaughterTracks.Clear();
 939     //if(unsetvtx) d->UnsetOwnPrimaryVtx();
 940   } //end for prongs
 941   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
 942   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
 943   // Post the data
 944   PostData(1,fOutputMass);
 945   PostData(2,fDistr);
 946   PostData(3,fNentries);
 947   PostData(5,fCounter);
 948   PostData(6,fOutputMassPt);
 949   PostData(7,fVariablesTree);
 950
 951   return;
 952 }
 953
 954 //____________________________________________________________________________
 955 void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
 956   //
 957   // function used in UserExec to fill variable histograms:
 958   //
 959
 960
 961   Int_t pdgDgD0toKpi[2]={321,211};
 962   Int_t lab=-9999;
 963   if(fReadMC) lab=part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
 964   //Double_t pt = d->Pt(); //mother pt
 965   Int_t isSelectedPID=3;
 966   if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
 967   if (isSelectedPID==0)fNentries->Fill(7);
 968   if (isSelectedPID==1)fNentries->Fill(8);
 969   if (isSelectedPID==2)fNentries->Fill(9);
 970   if (isSelectedPID==3)fNentries->Fill(10);
 971   //fNentries->Fill(8+isSelectedPID);
 972
 973   if(fCutOnDistr && !fIsSelectedCandidate) return;
 974   //printf("\nif no cuts or cuts passed\n");
 975
 976
 977   //add distr here
 978   UInt_t pdgs[2];
 979
 980   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 981   pdgs[0]=211;
 982   pdgs[1]=321;
 983   Double_t minvD0 = part->InvMassD0();
 984   pdgs[1]=211;
 985   pdgs[0]=321;
 986   Double_t minvD0bar = part->InvMassD0bar();
 987   //cout<<"inside mass cut"<<endl;
 988
 989   Double_t invmasscut=0.03;
 990
 991   TString fillthispi="",fillthisK="",fillthis="", fillthispt="";
 992
 993   Int_t ptbin=cuts->PtBin(part->Pt());
 994   Double_t pt = part->Pt();
 995
 996   Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
 997   dz1[0]=-99; dz2[0]=-99;
 998   Double_t d0[2];
 999   Double_t decl[2]={-99,-99};
1000   Bool_t recalcvtx=kFALSE;
1001
1002
1003
1004   if(fCuts->GetIsPrimaryWithoutDaughters()){
1005     recalcvtx=kTRUE;
1006     //recalculate vertex
1007     AliAODVertex *vtxProp=0x0;
1008     vtxProp=GetPrimaryVtxSkipped(aod);
1009     if(vtxProp) {
1010       part->SetOwnPrimaryVtx(vtxProp);
1011       //Bool_t unsetvtx=kTRUE;
1012       //Calculate d0 for daughter tracks with Vtx Skipped
1013       AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1014       AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1015       esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
1016       esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
1017       delete vtxProp; vtxProp=NULL;
1018       delete esdtrack1;
1019       delete esdtrack2;
1020     }
1021
1022     d0[0]=dz1[0];
1023     d0[1]=dz2[0];
1024
1025     decl[0]=part->DecayLength2();
1026     decl[1]=part->NormalizedDecayLength2();
1027     part->UnsetOwnPrimaryVtx();
1028
1029   }
1030
1031   Double_t cosThetaStarD0 = 99;
1032   Double_t cosThetaStarD0bar = 99;
1033   Double_t cosPointingAngle = 99;
1034   Double_t normalizedDecayLength2 = -1, normalizedDecayLengthxy=-1;
1035   Double_t decayLength2 = -1, decayLengthxy=-1;
1036   Double_t ptProng[2]={-99,-99};
1037   Double_t d0Prong[2]={-99,-99};
1038
1039
1040   //disable the PID
1041   if(!fUsePid4Distr) isSelectedPID=0;
1042   if((lab>=0 && fReadMC) || (!fReadMC && isSelectedPID)){ //signal (from MC or PID)
1043
1044     //check pdg of the prongs
1045     AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1046     AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
1047
1048     if(!prong0 || !prong1) {
1049       return;
1050     }
1051     Int_t labprong[2];
1052     if(fReadMC){
1053       labprong[0]=prong0->GetLabel();
1054       labprong[1]=prong1->GetLabel();
1055     }
1056     AliAODMCParticle *mcprong=0;
1057     Int_t pdgProng[2]={0,0};
1058
1059     for (Int_t iprong=0;iprong<2;iprong++){
1060       if(fReadMC && labprong[iprong]>=0) {
1061         mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1062         pdgProng[iprong]=mcprong->GetPdgCode();
1063       }
1064     }
1065
1066     if(fSys==0){
1067       //no mass cut ditributions: ptbis
1068
1069       fillthispi="hptpiSnoMcut_";
1070       fillthispi+=ptbin;
1071
1072       fillthisK="hptKSnoMcut_";
1073       fillthisK+=ptbin;
1074
1075       if ((TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321)
1076           || (isSelectedPID==1 || isSelectedPID==3)){
1077         ((TH1F*)listout->FindObject(fillthispi))->Fill(prong0->Pt());
1078         ((TH1F*)listout->FindObject(fillthisK))->Fill(prong1->Pt());
1079       }
1080
1081       if ((TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211)
1082           || (isSelectedPID==2 || isSelectedPID==3)){
1083         ((TH1F*)listout->FindObject(fillthisK))->Fill(prong0->Pt());
1084         ((TH1F*)listout->FindObject(fillthispi))->Fill(prong1->Pt());
1085       }
1086     }
1087
1088     //no mass cut ditributions: mass
1089     fillthis="hMassS_";
1090     fillthis+=ptbin;
1091     fillthispt="histSgnPt";
1092
1093     if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)
1094         || (!fReadMC && (isSelectedPID==1 || isSelectedPID==3))){//D0
1095       ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1096       if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1097     }
1098     else { //D0bar
1099       if(fReadMC || (!fReadMC && isSelectedPID > 1)){
1100         ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1101         if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1102       }
1103     }
1104
1105     //apply cut on invariant mass on the pair
1106     if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1107
1108       cosThetaStarD0 = part->CosThetaStarD0();
1109       cosThetaStarD0bar = part->CosThetaStarD0bar();
1110       cosPointingAngle = part->CosPointingAngle();
1111       normalizedDecayLength2 = part->NormalizedDecayLength2();
1112       decayLength2 = part->DecayLength2();
1113       decayLengthxy = part->DecayLengthXY();
1114       normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1115
1116       ptProng[0]=prong0->Pt(); ptProng[1]=prong1->Pt();
1117       d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1118
1119       if(fArray==1) cout<<"LS signal: ERROR"<<endl;
1120       for (Int_t iprong=0; iprong<2; iprong++){
1121         AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1122         if (fReadMC) labprong[iprong]=prong->GetLabel();
1123
1124         //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
1125         Int_t pdgprong=0;
1126         if(fReadMC && labprong[iprong]>=0) {
1127           mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1128           pdgprong=mcprong->GetPdgCode();
1129         }
1130
1131         Bool_t isPionHere[2]={(isSelectedPID==1 || isSelectedPID==3) ? kTRUE : kFALSE,(isSelectedPID==2 || isSelectedPID==3) ? kTRUE : kFALSE};
1132
1133         if(TMath::Abs(pdgprong)==211 || isPionHere[iprong]) {
1134           //cout<<"pi"<<endl;
1135
1136           if(fSys==0){
1137             fillthispi="hptpiS_";
1138             fillthispi+=ptbin;
1139             ((TH1F*)listout->FindObject(fillthispi))->Fill(ptProng[iprong]);
1140           }
1141           fillthispi="hd0piS_";
1142           fillthispi+=ptbin;
1143           ((TH1F*)listout->FindObject(fillthispi))->Fill(d0Prong[iprong]);
1144           if(recalcvtx) {
1145
1146             fillthispi="hd0vpiS_";
1147             fillthispi+=ptbin;
1148             ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
1149           }
1150
1151         }
1152
1153         if(TMath::Abs(pdgprong)==321 || !isPionHere[iprong]) {
1154           //cout<<"kappa"<<endl;
1155           if(fSys==0){
1156             fillthisK="hptKS_";
1157             fillthisK+=ptbin;
1158             ((TH1F*)listout->FindObject(fillthisK))->Fill(ptProng[iprong]);
1159           }
1160           fillthisK="hd0KS_";
1161           fillthisK+=ptbin;
1162           ((TH1F*)listout->FindObject(fillthisK))->Fill(d0Prong[iprong]);
1163           if (recalcvtx){
1164             fillthisK="hd0vKS_";
1165             fillthisK+=ptbin;
1166             ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
1167           }
1168         }
1169
1170         if(fSys==0){
1171           fillthis="hcosthpointd0S_";
1172           fillthis+=ptbin;
1173           ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[iprong]);
1174         }
1175       } //end loop on prongs
1176
1177       fillthis="hdcaS_";
1178       fillthis+=ptbin;
1179       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1180
1181       fillthis="hcosthetapointS_";
1182       fillthis+=ptbin;
1183       ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1184
1185       fillthis="hcosthetapointxyS_";
1186       fillthis+=ptbin;
1187       ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1188
1189
1190       fillthis="hd0d0S_";
1191       fillthis+=ptbin;
1192       ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1193
1194       fillthis="hdeclS_";
1195       fillthis+=ptbin;
1196       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
1197
1198       fillthis="hnormdeclS_";
1199       fillthis+=ptbin;
1200       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1201
1202       fillthis="hdeclxyS_";
1203       fillthis+=ptbin;
1204       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1205
1206       fillthis="hnormdeclxyS_";
1207       fillthis+=ptbin;
1208       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1209
1210       fillthis="hdeclxyd0d0S_";
1211       fillthis+=ptbin;
1212       ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1213
1214       fillthis="hnormdeclxyd0d0S_";
1215       fillthis+=ptbin;
1216       ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1217
1218       if(recalcvtx) {
1219         fillthis="hdeclvS_";
1220         fillthis+=ptbin;
1221         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1222
1223         fillthis="hnormdeclvS_";
1224         fillthis+=ptbin;
1225         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1226
1227         fillthis="hd0d0vS_";
1228         fillthis+=ptbin;
1229         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1230       }
1231
1232       if(fSys==0){
1233         fillthis="hcosthetastarS_";
1234         fillthis+=ptbin;
1235         if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)) ((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1236         else {
1237           if (fReadMC || isSelectedPID>1)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1238           if(isSelectedPID==1 || isSelectedPID==3)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1239         }
1240
1241         fillthis="hcosthpointd0d0S_";
1242         fillthis+=ptbin;
1243         ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1244       }
1245
1246     } //end mass cut
1247
1248   } else{ //Background or LS
1249     //if(!fReadMC){
1250     //cout<<"is background"<<endl;
1251
1252     //no mass cut distributions: mass, ptbis
1253     fillthis="hMassB_";
1254     fillthis+=ptbin;
1255     fillthispt="histBkgPt";
1256
1257     if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))) {
1258       ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1259       if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1260     }
1261     if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
1262       ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1263       if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1264     }
1265     if(fSys==0){
1266       fillthis="hptB1prongnoMcut_";
1267       fillthis+=ptbin;
1268
1269       ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1270
1271       fillthis="hptB2prongsnoMcut_";
1272       fillthis+=ptbin;
1273       ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1274       ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1275     }
1276     //apply cut on invariant mass on the pair
1277     if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1278       if(fSys==0){
1279         ptProng[0]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt(); ptProng[1]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt();
1280         cosThetaStarD0 = part->CosThetaStarD0();
1281         cosThetaStarD0bar = part->CosThetaStarD0bar();
1282       }
1283
1284       cosPointingAngle = part->CosPointingAngle();
1285       normalizedDecayLength2 = part->NormalizedDecayLength2();
1286       decayLength2 = part->DecayLength2();
1287       decayLengthxy = part->DecayLengthXY();
1288       normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1289       d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1290
1291
1292       AliAODTrack *prongg=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1293       if(!prongg) {
1294         if(fDebug>2) cout<<"No daughter found";
1295         return;
1296       }
1297       else{
1298         if(fArray==1){
1299           if(prongg->Charge()==1) {
1300             //fTotPosPairs[ptbin]++;
1301             ((TH1F*)fOutputMass->FindObject("hpospair"))->Fill(ptbin);
1302           } else {
1303             //fTotNegPairs[ptbin]++;
1304             ((TH1F*)fOutputMass->FindObject("hnegpair"))->Fill(ptbin);
1305           }
1306         }
1307       }
1308
1309       fillthis="hd0B_";
1310       fillthis+=ptbin;
1311       ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1312       ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1313
1314       if(fReadMC){
1315         Int_t pdgMother[2]={0,0};
1316         Double_t factor[2]={1,1};
1317
1318         for(Int_t iprong=0;iprong<2;iprong++){
1319           AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1320           lab=prong->GetLabel();
1321           if(lab>=0){
1322             AliAODMCParticle* mcprong=(AliAODMCParticle*)arrMC->At(lab);
1323             if(mcprong){
1324               Int_t labmom=mcprong->GetMother();
1325               if(labmom>=0){
1326                 AliAODMCParticle* mcmother=(AliAODMCParticle*)arrMC->At(labmom);
1327                 if(mcmother) pdgMother[iprong]=mcmother->GetPdgCode();
1328               }
1329             }
1330           }
1331
1332           if(fSys==0){
1333
1334             fillthis="hd0moresB_";
1335             fillthis+=ptbin;
1336
1337             if(TMath::Abs(pdgMother[iprong])==310 || TMath::Abs(pdgMother[iprong])==130 || TMath::Abs(pdgMother[iprong])==321){ //K^0_S, K^0_L, K^+-
1338               if(ptProng[iprong]<=1)factor[iprong]=1./.7;
1339               else factor[iprong]=1./.6;
1340               fNentries->Fill(11);
1341             }
1342
1343             if(TMath::Abs(pdgMother[iprong])==3122) { //Lambda
1344               factor[iprong]=1./0.25;
1345               fNentries->Fill(12);
1346             }
1347             fillthis="hd0moresB_";
1348             fillthis+=ptbin;
1349
1350             ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[iprong],factor[iprong]);
1351
1352             if(recalcvtx){
1353               fillthis="hd0vmoresB_";
1354               fillthis+=ptbin;
1355               ((TH1F*)listout->FindObject(fillthis))->Fill(d0[iprong],factor[iprong]);
1356             }
1357           }
1358         } //loop on prongs
1359
1360         if(fSys==0){
1361           fillthis="hd0d0moresB_";
1362           fillthis+=ptbin;
1363           ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0(),factor[0]*factor[1]);
1364
1365           fillthis="hcosthetapointmoresB_";
1366           fillthis+=ptbin;
1367           ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,factor[0]*factor[1]);
1368
1369           if(recalcvtx){
1370             fillthis="hd0d0vmoresB_";
1371             fillthis+=ptbin;
1372             ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1],factor[0]*factor[1]);
1373           }
1374         }
1375       } //readMC
1376
1377       if(fSys==0){
1378         //normalise pt distr to half afterwards
1379         fillthis="hptB_";
1380         fillthis+=ptbin;
1381         ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[0]);
1382         ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[1]);
1383
1384         fillthis="hcosthetastarB_";
1385         fillthis+=ptbin;
1386         if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3)))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1387         if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1388
1389
1390         fillthis="hd0p0B_";
1391         fillthis+=ptbin;
1392         ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1393         fillthis="hd0p1B_";
1394         fillthis+=ptbin;
1395         ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1396
1397         fillthis="hcosthpointd0d0B_";
1398         fillthis+=ptbin;
1399         ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1400
1401         fillthis="hcosthpointd0B_";
1402         fillthis+=ptbin;
1403         ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[0]);
1404         ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[1]);
1405
1406
1407         if(recalcvtx){
1408
1409           fillthis="hd0vp0B_";
1410           fillthis+=ptbin;
1411           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1412           fillthis="hd0vp1B_";
1413           fillthis+=ptbin;
1414           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1415
1416           fillthis="hd0vB_";
1417           fillthis+=ptbin;
1418           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1419           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1420
1421         }
1422
1423       }
1424
1425       fillthis="hdcaB_";
1426       fillthis+=ptbin;
1427       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1428
1429       fillthis="hd0d0B_";
1430       fillthis+=ptbin;
1431       ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]*d0Prong[1]);
1432
1433       if(recalcvtx){
1434         fillthis="hd0d0vB_";
1435         fillthis+=ptbin;
1436         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1437       }
1438
1439       fillthis="hcosthetapointB_";
1440       fillthis+=ptbin;
1441       ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1442
1443       fillthis="hcosthetapointxyB_";
1444       fillthis+=ptbin;
1445       ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1446
1447       fillthis="hdeclB_";
1448       fillthis+=ptbin;
1449       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
1450
1451       fillthis="hnormdeclB_";
1452       fillthis+=ptbin;
1453       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1454
1455       fillthis="hdeclxyB_";
1456       fillthis+=ptbin;
1457       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1458
1459       fillthis="hnormdeclxyB_";
1460       fillthis+=ptbin;
1461       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1462
1463       fillthis="hdeclxyd0d0B_";
1464       fillthis+=ptbin;
1465       ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1466
1467       fillthis="hnormdeclxyd0d0B_";
1468       fillthis+=ptbin;
1469       ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1470
1471       if(recalcvtx) {
1472
1473         fillthis="hdeclvB_";
1474         fillthis+=ptbin;
1475         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1476
1477         fillthis="hnormdeclvB_";
1478         fillthis+=ptbin;
1479         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1480
1481
1482       }
1483     }//mass cut
1484   }//else (background)
1485
1486   return;
1487 }
1488
1489 //____________________________________________________________________________
1490 void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi* cuts, TList *listout){
1491   //
1492   // function used in UserExec to fill mass histograms:
1493   //
1494
1495
1496   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1497
1498   //cout<<"is selected = "<<fIsSelectedCandidate<<endl;
1499
1500   // Fill candidate variable Tree (track selection, no candidate selection)
1501   if( fWriteVariableTree && !part->HasBadDaughters()
1502       && fCuts->AreDaughtersSelected(part) && fCuts->IsSelectedPID(part) ){
1503     fCandidateVariables[0] = part->InvMassD0();
1504     fCandidateVariables[1] = part->InvMassD0bar();
1505     fCandidateVariables[2] = part->Pt();
1506     fCandidateVariables[3] = part->GetDCA();
1507     Double_t ctsD0=0. ,ctsD0bar=0.; part->CosThetaStarD0(ctsD0,ctsD0bar);
1508     fCandidateVariables[4] = ctsD0;
1509     fCandidateVariables[5] = ctsD0bar;
1510     fCandidateVariables[6] = part->Pt2Prong(0);
1511     fCandidateVariables[7] = part->Pt2Prong(1);
1512     fCandidateVariables[8] = part->Getd0Prong(0);
1513     fCandidateVariables[9] = part->Getd0Prong(1);
1514     fCandidateVariables[10] = part->Prodd0d0();
1515     fCandidateVariables[11] = part->CosPointingAngle();
1516     fCandidateVariables[12] = part->CosPointingAngleXY();
1517     fCandidateVariables[13] = part->NormalizedDecayLengthXY();
1518     fCandidateVariables[14] = fCuts->IsSelectedSpecialCuts(part);
1519     fVariablesTree->Fill();
1520   }
1521
1522   //cout<<"check cuts = "<<endl;
1523   //cuts->PrintAll();
1524   if (!fIsSelectedCandidate){
1525     //cout<<"Not Selected"<<endl;
1526     //cout<<"Rejected because "<<cuts->GetWhy()<<endl;
1527     return;
1528   }
1529
1530
1531   if(fDebug>2)  cout<<"Candidate selected"<<endl;
1532
1533   Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
1534   //printf("SELECTED\n");
1535   Int_t ptbin=cuts->PtBin(part->Pt());
1536   Double_t pt = part->Pt();
1537
1538   Double_t impparXY=part->ImpParXY()*10000.;
1539   Double_t trueImpParXY=0.;
1540   Double_t arrayForSparse[3]={invmassD0,pt,impparXY};
1541   Double_t arrayForSparseTrue[3]={invmassD0,pt,trueImpParXY};
1542
1543
1544   // AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1545   // if(!prong) {
1546   //   AliDebug(2,"No daughter found");
1547   //   return;
1548   // }
1549   // else{
1550     // if(prong->Charge()==1) {
1551     //   ((TH1F*)fDistr->FindObject("hpospair"))->Fill(fCuts->GetNPtBins()+ptbin);
1552     //   //fTotPosPairs[ptbin]++;
1553     // } else {
1554     //   ((TH1F*)fDistr->FindObject("hnegpair"))->Fill(fCuts->GetNPtBins()+ptbin);
1555     //   //fTotNegPairs[ptbin]++;
1556     // }
1557   //  }
1558
1559   // for(Int_t it=0;it<2;it++){
1560
1561   //    //request on spd points to be addes
1562   //   if(/*nSPD==2 && */part->Pt() > 5. && (TMath::Abs(invmassD0-mPDG)<0.01 || TMath::Abs(invmassD0bar-mPDG)<0.01)){
1563   //     FILE *f=fopen("4display.txt","a");
1564   //     fprintf(f,"pt: %f \n Rapidity: %f \t Period Number: %x \t Run Number: %d \t BunchCrossNumb: %d \t OrbitNumber: %d\n",part->Pt(),part->Y(421),aod->GetPeriodNumber(),aod->GetRunNumber(),aod->GetBunchCrossNumber(),aod->GetOrbitNumber());
1565   //     fclose(f);
1566   //     //printf("PrimVtx NContributors: %d \n Prongs Rel Angle: %f \n \n",ncont,relangle);
1567   //   }
1568   // }
1569
1570   TString fillthis="", fillthispt="", fillthismasspt="";
1571   Int_t pdgDgD0toKpi[2]={321,211};
1572   Int_t labD0=-1;
1573   Bool_t isPrimary=kTRUE;
1574   if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
1575
1576   //count candidates selected by cuts
1577   fNentries->Fill(1);
1578   //count true D0 selected by cuts
1579   if (fReadMC && labD0>=0) fNentries->Fill(2);
1580
1581   if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
1582
1583     arrayForSparse[0]=invmassD0; arrayForSparse[2]=impparXY;
1584
1585     if(fReadMC){
1586       if(labD0>=0) {
1587         if(fArray==1) cout<<"LS signal ERROR"<<endl;
1588
1589         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1590         Int_t pdgD0 = partD0->GetPdgCode();
1591         //      cout<<"pdg = "<<pdgD0<<endl;
1592
1593         if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
1594         if(!isPrimary)
1595           trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
1596         arrayForSparseTrue[0]=invmassD0; arrayForSparseTrue[2]=trueImpParXY;
1597
1598         if (pdgD0==421){ //D0
1599           //      cout<<"Fill S with D0"<<endl;
1600           fillthis="histSgn_";
1601           fillthis+=ptbin;
1602           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1603
1604           if(fFillPtHist){
1605             fillthismasspt="histSgnPt";
1606             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt);
1607           }
1608           if(fFillImpParHist){
1609             if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
1610             else {
1611               fHistMassPtImpParTC[2]->Fill(arrayForSparse);
1612               fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
1613             }
1614           }
1615
1616           if(fSys==0){
1617             if(TMath::Abs(invmassD0 - mPDG) < 0.027 && fFillVarHists){
1618               fillthis="histSgn27_";
1619               fillthis+=ptbin;
1620               ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1621             }
1622           }
1623         } else{ //it was a D0bar
1624           fillthis="histRfl_";
1625           fillthis+=ptbin;
1626           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1627
1628           if(fFillPtHist){
1629             fillthismasspt="histRflPt";
1630             //      cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
1631             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt);
1632           }
1633
1634         }
1635       } else {//background
1636         fillthis="histBkg_";
1637         fillthis+=ptbin;
1638         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1639
1640         if(fFillPtHist){
1641           fillthismasspt="histBkgPt";
1642           //      cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
1643           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt);
1644         }
1645         if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse);
1646
1647       }
1648
1649     }else{
1650       fillthis="histMass_";
1651       fillthis+=ptbin;
1652       //      cout<<"Filling "<<fillthis<<endl;
1653
1654       //      printf("Fill mass with D0");
1655       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1656
1657       if(fFillPtHist){
1658         fillthismasspt="histMassPt";
1659         //      cout<<"Filling "<<fillthismasspt<<endl;
1660         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt);
1661       }
1662       if(fFillImpParHist) {
1663         //      cout << "Filling fHistMassPtImpParTC[0]"<<endl;
1664         fHistMassPtImpParTC[0]->Fill(arrayForSparse);
1665       }
1666
1667     }
1668
1669   }
1670   if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
1671
1672     arrayForSparse[0]=invmassD0bar; arrayForSparse[2]=impparXY;
1673
1674     if(fReadMC){
1675       if(labD0>=0) {
1676         if(fArray==1) cout<<"LS signal ERROR"<<endl;
1677         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1678         Int_t pdgD0 = partD0->GetPdgCode();
1679         //      cout<<" pdg = "<<pdgD0<<endl;
1680
1681         if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
1682         if(!isPrimary)
1683           trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
1684         arrayForSparseTrue[0]=invmassD0bar; arrayForSparseTrue[2]=trueImpParXY;
1685
1686         if (pdgD0==-421){ //D0bar
1687           fillthis="histSgn_";
1688           fillthis+=ptbin;
1689           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1690           // if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1691           //   fillthis="histSgn27_";
1692           //   fillthis+=ptbin;
1693           //   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1694           // }
1695
1696           if(fFillPtHist){
1697             fillthismasspt="histSgnPt";
1698             //      cout<<" Filling "<< fillthismasspt << endl;
1699             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt);
1700           }
1701           if(fFillImpParHist){
1702             //      cout << " Filling impact parameter thnsparse"<<endl;
1703             if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
1704             else {
1705               fHistMassPtImpParTC[2]->Fill(arrayForSparse);
1706               fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
1707             }
1708           }
1709
1710         } else{
1711           fillthis="histRfl_";
1712           fillthis+=ptbin;
1713           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1714           if(fFillPtHist){
1715             fillthismasspt="histRflPt";
1716             //      cout << " Filling "<<fillthismasspt<<endl;
1717             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt);
1718           }
1719         }
1720       } else {//background or LS
1721         fillthis="histBkg_";
1722         fillthis+=ptbin;
1723         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1724
1725         if(fFillPtHist){
1726           fillthismasspt="histBkgPt";
1727           //      cout<<" Filling "<< fillthismasspt << endl;
1728           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt);
1729         }
1730         if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse);
1731
1732       }
1733     }else{
1734       fillthis="histMass_";
1735       fillthis+=ptbin;
1736       //      printf("Fill mass with D0bar");
1737       ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
1738
1739       if(fFillPtHist){
1740         fillthismasspt="histMassPt";
1741         //      cout<<" Filling "<< fillthismasspt << endl;
1742         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt);
1743       }
1744       if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse);
1745
1746     }
1747   }
1748
1749   return;
1750 }
1751
1752 //__________________________________________________________________________
1753 AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev){
1754   //Calculate the primary vertex w/o the daughter tracks of the candidate
1755
1756   Int_t skipped[2];
1757   Int_t nTrksToSkip=2;
1758   AliAODTrack *dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1759   if(!dgTrack){
1760     AliDebug(2,"no daughter found!");
1761     return 0x0;
1762   }
1763   skipped[0]=dgTrack->GetID();
1764   dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(1);
1765   if(!dgTrack){
1766     AliDebug(2,"no daughter found!");
1767     return 0x0;
1768   }
1769   skipped[1]=dgTrack->GetID();
1770
1771   AliESDVertex *vertexESD=0x0;
1772   AliAODVertex *vertexAOD=0x0;
1773   AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
1774
1775   //
1776   vertexer->SetSkipTracks(nTrksToSkip,skipped);
1777   vertexer->SetMinClusters(4);
1778   vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
1779   if(!vertexESD) return vertexAOD;
1780   if(vertexESD->GetNContributors()<=0) {
1781     AliDebug(2,"vertexing failed");
1782     delete vertexESD; vertexESD=NULL;
1783     return vertexAOD;
1784   }
1785
1786   delete vertexer; vertexer=NULL;
1787
1788
1789   // convert to AliAODVertex
1790   Double_t pos[3],cov[6],chi2perNDF;
1791   vertexESD->GetXYZ(pos); // position
1792   vertexESD->GetCovMatrix(cov); //covariance matrix
1793   chi2perNDF = vertexESD->GetChi2toNDF();
1794   delete vertexESD; vertexESD=NULL;
1795
1796   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1797   return vertexAOD;
1798
1799 }
1800
1801
1802 //________________________________________________________________________
1803 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1804 {
1805   // Terminate analysis
1806   //
1807   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1808
1809
1810   fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1811   if (!fOutputMass) {
1812     printf("ERROR: fOutputMass not available\n");
1813     return;
1814   }
1815   fOutputMassPt = dynamic_cast<TList*> (GetOutputData(6));
1816   if ((fFillPtHist || fFillImpParHist) && !fOutputMassPt) {
1817     printf("ERROR: fOutputMass not available\n");
1818     return;
1819   }
1820
1821   if(fFillVarHists){
1822     fDistr = dynamic_cast<TList*> (GetOutputData(2));
1823     if (!fDistr) {
1824       printf("ERROR: fDistr not available\n");
1825       return;
1826     }
1827   }
1828
1829   fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
1830
1831   if(!fNentries){
1832     printf("ERROR: fNEntries not available\n");
1833     return;
1834   }
1835   fCuts = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(4));
1836   if(!fCuts){
1837     printf("ERROR: fCuts not available\n");
1838     return;
1839   }
1840   fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(5));
1841   if (!fCounter) {
1842     printf("ERROR: fCounter not available\n");
1843     return;
1844   }
1845   Int_t nptbins=fCuts->GetNPtBins();
1846   for(Int_t ipt=0;ipt<nptbins;ipt++){
1847
1848     if(fArray==1 && fFillVarHists){
1849       fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(nptbins+ipt+1,nptbins+ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(nptbins+ipt+1,nptbins+ipt+2)); //after cuts
1850
1851
1852       if(fLsNormalization>1e-6) {
1853
1854         TString massName="histMass_";
1855         massName+=ipt;
1856         ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
1857
1858       }
1859
1860
1861       fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(ipt+1,ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(ipt+1,ipt+2));
1862       //fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1863
1864       if(fLsNormalization>1e-6) {
1865
1866         TString nameDistr="hptB_";
1867         nameDistr+=ipt;
1868         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1869         nameDistr="hdcaB_";
1870         nameDistr+=ipt;
1871         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1872         nameDistr="hcosthetastarB_";
1873         nameDistr+=ipt;
1874         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1875         nameDistr="hd0B_";
1876         nameDistr+=ipt;
1877         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1878         nameDistr="hd0d0B_";
1879         nameDistr+=ipt;
1880         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1881         nameDistr="hcosthetapointB_";
1882         nameDistr+=ipt;
1883         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1884         if(fSys==0){
1885           nameDistr="hcosthpointd0d0B_";
1886           nameDistr+=ipt;
1887           ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
1888         }
1889       }
1890     }
1891   }
1892   TString cvname,cstname;
1893
1894   if (fArray==0){
1895     cvname="D0invmass";
1896     cstname="cstat0";
1897   } else {
1898     cvname="LSinvmass";
1899     cstname="cstat1";
1900   }
1901
1902   TCanvas *cMass=new TCanvas(cvname,cvname);
1903   cMass->cd();
1904   ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
1905
1906   TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
1907   cStat->cd();
1908   cStat->SetGridy();
1909   fNentries->Draw("htext0");
1910
1911   // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
1912   // ccheck->cd();
1913
1914   return;
1915 }
1916
1917
1918 //________________________________________________________________________
1919 void AliAnalysisTaskSED0Mass::CreateImpactParameterHistos(){
1920   // Histos for impact paramter study
1921
1922   Int_t nmassbins=200;
1923   Double_t fLowmasslimit=1.5648, fUpmasslimit=2.1648;
1924   Int_t fNImpParBins=400;
1925   Double_t fLowerImpPar=-2000., fHigherImpPar=2000.;
1926   Int_t nbins[3]={nmassbins,200,fNImpParBins};
1927   Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
1928   Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
1929
1930
1931   fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
1932                                         "Mass vs. pt vs.imppar - All",
1933                                         3,nbins,xmin,xmax);
1934   fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
1935                                         "Mass vs. pt vs.imppar - promptD",
1936                                         3,nbins,xmin,xmax);
1937   fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
1938                                         "Mass vs. pt vs.imppar - DfromB",
1939                                         3,nbins,xmin,xmax);
1940   fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1941                                         "Mass vs. pt vs.true imppar -DfromB",
1942                                         3,nbins,xmin,xmax);
1943   fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
1944                                         "Mass vs. pt vs.imppar - backgr.",
1945                                         3,nbins,xmin,xmax);
1946
1947   for(Int_t i=0; i<5;i++){
1948     fOutputMassPt->Add(fHistMassPtImpParTC[i]);
1949   }
1950 }
1951
1952 //_________________________________________________________________________________________________
1953 Float_t AliAnalysisTaskSED0Mass::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD0) const {
1954   // true impact parameter calculation
1955
1956   printf(" AliAnalysisTaskSED0MassV1::GetTrueImpactParameter() \n");
1957
1958   Double_t vtxTrue[3];
1959   mcHeader->GetVertex(vtxTrue);
1960   Double_t origD[3];
1961   partD0->XvYvZv(origD);
1962   Short_t charge=partD0->Charge();
1963   Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
1964   for(Int_t iDau=0; iDau<2; iDau++){
1965     pXdauTrue[iDau]=0.;
1966     pYdauTrue[iDau]=0.;
1967     pZdauTrue[iDau]=0.;
1968   }
1969
1970   //  Int_t nDau=partD0->GetNDaughters();
1971   Int_t labelFirstDau = partD0->GetDaughter(0);
1972
1973   for(Int_t iDau=0; iDau<2; iDau++){
1974     Int_t ind = labelFirstDau+iDau;
1975     AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1976     if(!part) continue;
1977     Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
1978     if(!part){
1979       AliError("Daughter particle not found in MC array");
1980       return 99999.;
1981     }
1982     if(pdgCode==211 || pdgCode==321){
1983       pXdauTrue[iDau]=part->Px();
1984       pYdauTrue[iDau]=part->Py();
1985       pZdauTrue[iDau]=part->Pz();
1986     }
1987   }
1988
1989   Double_t d0dummy[2]={0.,0.};
1990   AliAODRecoDecayHF aodDzeroMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1991   return aodDzeroMC.ImpParXY();
1992
1993 }
1994
1995 //_________________________________________________________________________________________________
1996 Int_t AliAnalysisTaskSED0Mass::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1997   //
1998   // checking whether the mother of the particles come from a charm or a bottom quark
1999   //
2000   printf(" AliAnalysisTaskSED0MassV1::CheckOrigin() \n");
2001
2002   Int_t pdgGranma = 0;
2003   Int_t mother = 0;
2004   mother = mcPartCandidate->GetMother();
2005   Int_t istep = 0;
2006   Int_t abspdgGranma =0;
2007   Bool_t isFromB=kFALSE;
2008   Bool_t isQuarkFound=kFALSE;
2009   while (mother >0 ){
2010     istep++;
2011     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
2012     if (mcGranma){
2013       pdgGranma = mcGranma->GetPdgCode();
2014       abspdgGranma = TMath::Abs(pdgGranma);
2015       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
2016         isFromB=kTRUE;
2017       }
2018       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
2019       mother = mcGranma->GetMother();
2020     }else{
2021       AliError("Failed casting the mother particle!");
2022       break;
2023     }
2024   }
2025
2026   if(isFromB) return 5;
2027   else return 4;
2028 }