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 // Jeremy Wilkinson, jwilkinson@physi.uni-heidelberg.de (weighted Bayesian
  27 ////////////////////////////////////////////////////////////
  28
  29 #include <Riostream.h>
  30 #include <TClonesArray.h>
  31 #include <TCanvas.h>
  32 #include <TNtuple.h>
  33 #include <TTree.h>
  34 #include <TList.h>
  35 #include <TH1F.h>
  36 #include <TH2F.h>
  37 #include <TDatabasePDG.h>
  38
  39 #include <AliAnalysisDataSlot.h>
  40 #include <AliAnalysisDataContainer.h>
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDtrack.h"
  43 #include "AliVertexerTracks.h"
  44 #include "AliAODHandler.h"
  45 #include "AliAODEvent.h"
  46 #include "AliAODVertex.h"
  47 #include "AliAODTrack.h"
  48 #include "AliAODMCHeader.h"
  49 #include "AliAODMCParticle.h"
  50 #include "AliAODRecoDecayHF2Prong.h"
  51 #include "AliAODRecoCascadeHF.h"
  52 #include "AliAnalysisVertexingHF.h"
  53 #include "AliAnalysisTaskSE.h"
  54 #include "AliAnalysisTaskSED0Mass.h"
  55 #include "AliNormalizationCounter.h"
  56
  57 using std::cout;
  58 using std::endl;
  59
  60 ClassImp(AliAnalysisTaskSED0Mass)
  61
  62
  63 //________________________________________________________________________
  64 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
  65 AliAnalysisTaskSE(),
  66   fOutputMass(0),
  67   fOutputMassPt(0),
  68   fOutputMassY(0),
  69   fDistr(0),
  70   fNentries(0),
  71   fCuts(0),
  72   fArray(0),
  73   fReadMC(0),
  74   fCutOnDistr(0),
  75   fUsePid4Distr(0),
  76   fCounter(0),
  77   fNPtBins(1),
  78   fLsNormalization(1.),
  79   fFillOnlyD0D0bar(0),
  80   fDaughterTracks(),
  81   fIsSelectedCandidate(0),
  82   fFillVarHists(kTRUE),
  83   fSys(0),
  84   fIsRejectSDDClusters(0),
  85   fFillPtHist(kTRUE),
  86   fFillYHist(kFALSE),
  87   fFillImpParHist(kFALSE),
  88   fUseSelectionBit(kTRUE),
  89   fWriteVariableTree(kFALSE),
  90   fVariablesTree(0),
  91   fCandidateVariables(),
  92   fPIDCheck(kFALSE),
  93   fDrawDetSignal(kFALSE),
  94   fDetSignal(0)
  95 {
  96   // Default constructor
  97   for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
  98
  99 }
 100
 101 //________________________________________________________________________
 102 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
 103   AliAnalysisTaskSE(name),
 104   fOutputMass(0),
 105   fOutputMassPt(0),
 106   fOutputMassY(0),
 107   fDistr(0),
 108   fNentries(0),
 109   fCuts(0),
 110   fArray(0),
 111   fReadMC(0),
 112   fCutOnDistr(0),
 113   fUsePid4Distr(0),
 114   fCounter(0),
 115   fNPtBins(1),
 116   fLsNormalization(1.),
 117   fFillOnlyD0D0bar(0),
 118   fDaughterTracks(),
 119   fIsSelectedCandidate(0),
 120   fFillVarHists(kTRUE),
 121   fSys(0),
 122   fIsRejectSDDClusters(0),
 123   fFillPtHist(kTRUE),
 124   fFillYHist(kFALSE),
 125   fFillImpParHist(kFALSE),
 126   fUseSelectionBit(kTRUE),
 127   fWriteVariableTree(kFALSE),
 128   fVariablesTree(0),
 129   fCandidateVariables(),
 130   fPIDCheck(kFALSE),
 131   fDrawDetSignal(kFALSE),
 132   fDetSignal(0)
 133 {
 134   // Default constructor
 135
 136   fNPtBins=cuts->GetNPtBins();
 137
 138   fCuts=cuts;
 139   for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
 140
 141   // Output slot #1 writes into a TList container (mass with cuts)
 142   DefineOutput(1,TList::Class());  //My private output
 143   // Output slot #2 writes into a TList container (distributions)
 144   DefineOutput(2,TList::Class());  //My private output
 145   // Output slot #3 writes into a TH1F container (number of events)
 146   DefineOutput(3,TH1F::Class());  //My private output
 147   // Output slot #4 writes into a TList container (cuts)
 148   DefineOutput(4,AliRDHFCutsD0toKpi::Class());  //My private output
 149   // Output slot #5 writes Normalization Counter
 150   DefineOutput(5,AliNormalizationCounter::Class());
 151   // Output slot #6 stores the mass vs pt and impact parameter distributions
 152   DefineOutput(6,TList::Class());  //My private output
 153   // Output slot #7 keeps a tree of the candidate variables after track selection
 154   DefineOutput(7,TTree::Class());  //My private outpu
 155   // Output slot #8 writes into a TList container (Detector signals)
 156   DefineOutput(8, TList::Class()); //My private output
 157   // Output slot #9 stores the mass vs rapidity (y) distributions
 158   DefineOutput(9, TList::Class()); //My private output
 159 }
 160
 161 //________________________________________________________________________
 162 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
 163 {
 164   if (fOutputMass) {
 165     delete fOutputMass;
 166     fOutputMass = 0;
 167   }
 168   if (fOutputMassPt) {
 169     delete fOutputMassPt;
 170     fOutputMassPt = 0;
 171   }
 172   if (fOutputMassY) {
 173     delete fOutputMassY;
 174     fOutputMassY = 0;
 175   }
 176   if (fDistr) {
 177     delete fDistr;
 178     fDistr = 0;
 179   }
 180   if (fCuts) {
 181     delete fCuts;
 182     fCuts = 0;
 183   }
 184   for(Int_t i=0; i<5; i++){
 185     if(fHistMassPtImpParTC[i]) delete fHistMassPtImpParTC[i];
 186   }
 187   if (fNentries){
 188     delete fNentries;
 189     fNentries = 0;
 190   }
 191   if(fCounter){
 192     delete fCounter;
 193     fCounter=0;
 194   }
 195   if(fVariablesTree){
 196     delete fVariablesTree;
 197     fVariablesTree = 0;
 198   }
 199   if (fDetSignal) {
 200     delete fDetSignal;
 201     fDetSignal = 0;
 202   }
 203
 204 }
 205
 206 //________________________________________________________________________
 207 void AliAnalysisTaskSED0Mass::Init()
 208 {
 209   // Initialization
 210
 211   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
 212
 213
 214   AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
 215   const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 216   copyfCuts->SetName(nameoutput);
 217   // Post the data
 218   PostData(4,copyfCuts);
 219
 220
 221   return;
 222 }
 223
 224 //________________________________________________________________________
 225 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
 226 {
 227
 228   // Create the output container
 229   //
 230   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
 231
 232   // Several histograms are more conveniently managed in a TList
 233   fOutputMass = new TList();
 234   fOutputMass->SetOwner();
 235   fOutputMass->SetName("listMass");
 236
 237   fOutputMassPt = new TList();
 238   fOutputMassPt->SetOwner();
 239   fOutputMassPt->SetName("listMassPt");
 240
 241   fOutputMassY = new TList();
 242   fOutputMassY->SetOwner();
 243   fOutputMassY->SetName("listMassY");
 244
 245   fDistr = new TList();
 246   fDistr->SetOwner();
 247   fDistr->SetName("distributionslist");
 248
 249   fDetSignal = new TList();
 250   fDetSignal->SetOwner();
 251   fDetSignal->SetName("detectorsignals");
 252
 253   TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
 254   TString nameMassPt="", nameSgnPt="", nameBkgPt="", nameRflPt="";
 255   TString nameMassY="", nameSgnY="", nameBkgY="", nameRflY="";
 256   Int_t nbins2dPt=60; Float_t binInPt=0., binFinPt=30.;
 257   Int_t nbins2dY=60; Float_t binInY=-1.5, binFinY=1.5;
 258
 259   for(Int_t i=0;i<fCuts->GetNPtBins();i++){
 260
 261     nameMass="histMass_";
 262     nameMass+=i;
 263     nameSgn27="histSgn27_";
 264     nameSgn27+=i;
 265     nameSgn="histSgn_";
 266     nameSgn+=i;
 267     nameBkg="histBkg_";
 268     nameBkg+=i;
 269     nameRfl="histRfl_";
 270     nameRfl+=i;
 271     nameMassNocutsS="hMassS_";
 272     nameMassNocutsS+=i;
 273     nameMassNocutsB="hMassB_";
 274     nameMassNocutsB+=i;
 275
 276     //histograms of cut variable distributions
 277     if(fFillVarHists){
 278       if(fReadMC){
 279
 280         namedistr="hNclsD0vsptS_";
 281         namedistr+=i;
 282         TH2F *hNclsD0vsptS = new TH2F(namedistr.Data(),"N cls distrubution [S];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
 283         namedistr="hNclsD0barvsptS_";
 284         namedistr+=i;
 285         TH2F *hNclsD0barvsptS = new TH2F(namedistr.Data(),"N cls distrubution [S];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
 286
 287         namedistr="hNITSpointsD0vsptS_";
 288         namedistr+=i;
 289         TH2F *hNITSpointsD0vsptS = new TH2F(namedistr.Data(),"N ITS points distrubution [S];p_{T} [GeV/c];N points",200,0.,20.,7,0.,7.);
 290
 291         namedistr="hNSPDpointsD0S_";
 292         namedistr+=i;
 293         TH1I *hNSPDpointsD0S = new TH1I(namedistr.Data(),"N SPD points distrubution [S]; ;N tracks",4,0,4);
 294         hNSPDpointsD0S->GetXaxis()->SetBinLabel(1, "no SPD");
 295         hNSPDpointsD0S->GetXaxis()->SetBinLabel(2, "kOnlyFirst");
 296         hNSPDpointsD0S->GetXaxis()->SetBinLabel(3, "kOnlySecond");
 297         hNSPDpointsD0S->GetXaxis()->SetBinLabel(4, "kBoth");
 298
 299         namedistr="hptD0S_";
 300         namedistr+=i;
 301         TH1F *hptD0S = new TH1F(namedistr.Data(), "p_{T} distribution [S];p_{T} [GeV/c]",200,0.,20.);
 302         namedistr="hptD0barS_";
 303         namedistr+=i;
 304         TH1F *hptD0barS = new TH1F(namedistr.Data(), "p_{T} distribution [S];p_{T} [GeV/c]",200,0.,20.);
 305
 306         namedistr="hphiD0S_";
 307         namedistr+=i;
 308         TH1F *hphiD0S = new TH1F(namedistr.Data(), "#phi distribution [S];#phi [rad]",100,0.,2*TMath::Pi());
 309         namedistr="hphiD0barS_";
 310         namedistr+=i;
 311         TH1F *hphiD0barS = new TH1F(namedistr.Data(), "#phi distribution [S];#phi [rad]",100,0.,2*TMath::Pi());
 312
 313
 314         namedistr="hetaphiD0candidateS_";
 315         namedistr+=i;
 316         TH2F *hetaphiD0candidateS = new TH2F(namedistr.Data(), "D^{0} candidates #eta #phi distribution [S];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
 317         namedistr="hetaphiD0barcandidateS_";
 318         namedistr+=i;
 319         TH2F *hetaphiD0barcandidateS = new TH2F(namedistr.Data(), "anti-D^{0} candidates #eta #phi distribution [S];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
 320
 321         namedistr="hetaphiD0candidatesignalregionS_";
 322         namedistr+=i;
 323         TH2F *hetaphiD0candidatesignalregionS = new TH2F(namedistr.Data(), "D^{0} candidates #eta #phi distribution [S] [mass cut];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
 324         namedistr="hetaphiD0barcandidatesignalregionS_";
 325         namedistr+=i;
 326         TH2F *hetaphiD0barcandidatesignalregionS = new TH2F(namedistr.Data(), "anti-D^{0} candidates #eta #phi distribution [S] [mass cut];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
 327
 328         //  dca
 329         namedistr="hdcaS_";
 330         namedistr+=i;
 331         TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 332         // impact parameter
 333         namedistr="hd0piS_";
 334         namedistr+=i;
 335         TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
 336
 337         namedistr="hd0KS_";
 338         namedistr+=i;
 339         TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
 340         namedistr="hd0d0S_";
 341         namedistr+=i;
 342         TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 343
 344         //decay lenght
 345         namedistr="hdeclS_";
 346         namedistr+=i;
 347         TH1F *hdeclengthS = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm]",200,0,0.015);
 348
 349         namedistr="hnormdeclS_";
 350         namedistr+=i;
 351         TH1F *hnormdeclengthS = new TH1F(namedistr.Data(), "Normalized Decay Length^{2} distribution;(Decay Length/Err)^{2} ",200,0,12.);
 352
 353         namedistr="hdeclxyS_";
 354         namedistr+=i;
 355         TH1F* hdeclxyS=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
 356         namedistr="hnormdeclxyS_";
 357         namedistr+=i;
 358         TH1F* hnormdeclxyS=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
 359
 360         namedistr="hdeclxyd0d0S_";
 361         namedistr+=i;
 362         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);
 363
 364         namedistr="hnormdeclxyd0d0S_";
 365         namedistr+=i;
 366         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);
 367
 368         //  costhetapoint
 369         namedistr="hcosthetapointS_";
 370         namedistr+=i;
 371         TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 372
 373         namedistr="hcosthetapointxyS_";
 374         namedistr+=i;
 375         TH1F *hcosthetapointxyS = new TH1F(namedistr.Data(), "cos#theta_{Point} XYdistribution;cos#theta_{Point}",300,0.,1.);
 376
 377         TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
 378         tmpMS->Sumw2();
 379
 380         fDistr->Add(hNclsD0vsptS);
 381         fDistr->Add(hNclsD0barvsptS);
 382         fDistr->Add(hNITSpointsD0vsptS);
 383         fDistr->Add(hNSPDpointsD0S);
 384         fDistr->Add(hptD0S);
 385         fDistr->Add(hphiD0S);
 386         fDistr->Add(hptD0barS);
 387         fDistr->Add(hphiD0barS);
 388         fDistr->Add(hetaphiD0candidateS);
 389         fDistr->Add(hetaphiD0candidatesignalregionS);
 390         fDistr->Add(hetaphiD0barcandidateS);
 391         fDistr->Add(hetaphiD0barcandidatesignalregionS);
 392
 393         fDistr->Add(hdcaS);
 394
 395         fDistr->Add(hd0piS);
 396         fDistr->Add(hd0KS);
 397
 398         fDistr->Add(hd0d0S);
 399
 400         fDistr->Add(hcosthetapointS);
 401
 402         fDistr->Add(hcosthetapointxyS);
 403
 404         fDistr->Add(hdeclengthS);
 405
 406         fDistr->Add(hnormdeclengthS);
 407
 408         fDistr->Add(hdeclxyS);
 409
 410         fDistr->Add(hnormdeclxyS);
 411
 412         fDistr->Add(hdeclxyd0d0S);
 413         fDistr->Add(hnormdeclxyd0d0S);
 414
 415         fDistr->Add(tmpMS);
 416       }
 417
 418
 419       //Ncls, phi, pt distrubutions
 420
 421       namedistr="hNclsD0vsptB_";
 422       namedistr+=i;
 423       TH2F *hNclsD0vsptB = new TH2F(namedistr.Data(),"N cls distrubution [B];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
 424       namedistr="hNclsD0barvsptB_";
 425       namedistr+=i;
 426       TH2F *hNclsD0barvsptB = new TH2F(namedistr.Data(),"N cls distrubution [B];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
 427
 428       namedistr="hNITSpointsD0vsptB_";
 429       namedistr+=i;
 430       TH2F *hNITSpointsD0vsptB = new TH2F(namedistr.Data(),"N ITS points distrubution [B];p_{T} [GeV/c];N points",200,0.,20.,7,0.,7.);
 431
 432       namedistr="hNSPDpointsD0B_";
 433       namedistr+=i;
 434       TH1I *hNSPDpointsD0B = new TH1I(namedistr.Data(),"N SPD points distrubution [B]; ;N tracks",4,0,4);
 435       hNSPDpointsD0B->GetXaxis()->SetBinLabel(1, "no SPD");
 436       hNSPDpointsD0B->GetXaxis()->SetBinLabel(2, "kOnlyFirst");
 437       hNSPDpointsD0B->GetXaxis()->SetBinLabel(3, "kOnlySecond");
 438       hNSPDpointsD0B->GetXaxis()->SetBinLabel(4, "kBoth");
 439
 440       namedistr="hptD0B_";
 441       namedistr+=i;
 442       TH1F *hptD0B = new TH1F(namedistr.Data(), "p_{T} distribution [B];p_{T} [GeV/c]",200,0.,20.);
 443       namedistr="hptD0barB_";
 444       namedistr+=i;
 445       TH1F *hptD0barB = new TH1F(namedistr.Data(), "p_{T} distribution [B];p_{T} [GeV/c]",200,0.,20.);
 446
 447       namedistr="hphiD0B_";
 448       namedistr+=i;
 449       TH1F *hphiD0B = new TH1F(namedistr.Data(), "#phi distribution [B];#phi [rad]",100,0.,2*TMath::Pi());
 450       namedistr="hphiD0barB_";
 451       namedistr+=i;
 452       TH1F *hphiD0barB = new TH1F(namedistr.Data(), "#phi distribution [B];#phi [rad]",100,0.,2*TMath::Pi());
 453
 454       namedistr="hetaphiD0candidateB_";
 455       namedistr+=i;
 456       TH2F *hetaphiD0candidateB = new TH2F(namedistr.Data(), "D^{0} candidates #eta #phi distribution [B];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
 457       namedistr="hetaphiD0barcandidateB_";
 458       namedistr+=i;
 459       TH2F *hetaphiD0barcandidateB = new TH2F(namedistr.Data(), "anti-D^{0} candidates #eta #phi distribution [B];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
 460
 461       namedistr="hetaphiD0candidatesignalregionB_";
 462       namedistr+=i;
 463       TH2F *hetaphiD0candidatesignalregionB = new TH2F(namedistr.Data(), "D^{0} candidates #eta #phi distribution [B] [mass cut];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
 464       namedistr="hetaphiD0barcandidatesignalregionB_";
 465       namedistr+=i;
 466       TH2F *hetaphiD0barcandidatesignalregionB = new TH2F(namedistr.Data(), "anti-D^{0} candidates #eta #phi distribution [B] [mass cut];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
 467
 468       //  dca
 469       namedistr="hdcaB_";
 470       namedistr+=i;
 471       TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 472
 473       // impact parameter
 474       namedistr="hd0B_";
 475       namedistr+=i;
 476       TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
 477
 478       namedistr="hd0d0B_";
 479       namedistr+=i;
 480       TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 481
 482       //decay lenght
 483       namedistr="hdeclB_";
 484       namedistr+=i;
 485       TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm^{2}]",200,0,0.015);
 486
 487       namedistr="hnormdeclB_";
 488       namedistr+=i;
 489       TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;(Decay Length/Err)^{2} ",200,0,12.);
 490
 491       namedistr="hdeclxyB_";
 492       namedistr+=i;
 493       TH1F* hdeclxyB=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
 494       namedistr="hnormdeclxyB_";
 495       namedistr+=i;
 496       TH1F* hnormdeclxyB=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
 497
 498       namedistr="hdeclxyd0d0B_";
 499       namedistr+=i;
 500       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);
 501
 502       namedistr="hnormdeclxyd0d0B_";
 503       namedistr+=i;
 504       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);
 505
 506       //  costhetapoint
 507       namedistr="hcosthetapointB_";
 508       namedistr+=i;
 509       TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 510
 511       namedistr="hcosthetapointxyB_";
 512       namedistr+=i;
 513       TH1F *hcosthetapointxyB = new TH1F(namedistr.Data(), "cos#theta_{Point} XY distribution;cos#theta_{Point} XY",300,0.,1.);
 514
 515       TH1F* tmpMB = new TH1F(nameMassNocutsB.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
 516       tmpMB->Sumw2();
 517
 518
 519       fDistr->Add(hNclsD0vsptB);
 520       fDistr->Add(hNclsD0barvsptB);
 521       fDistr->Add(hNITSpointsD0vsptB);
 522       fDistr->Add(hNSPDpointsD0B);
 523       fDistr->Add(hptD0B);
 524       fDistr->Add(hphiD0B);
 525       fDistr->Add(hptD0barB);
 526       fDistr->Add(hphiD0barB);
 527       fDistr->Add(hetaphiD0candidateB);
 528       fDistr->Add(hetaphiD0candidatesignalregionB);
 529       fDistr->Add(hetaphiD0barcandidateB);
 530       fDistr->Add(hetaphiD0barcandidatesignalregionB);
 531
 532       fDistr->Add(hdcaB);
 533
 534       fDistr->Add(hd0B);
 535
 536       fDistr->Add(hd0d0B);
 537
 538       fDistr->Add(hcosthetapointB);
 539
 540       fDistr->Add(hcosthetapointxyB);
 541
 542       fDistr->Add(hdeclengthB);
 543
 544       fDistr->Add(hnormdeclengthB);
 545
 546       fDistr->Add(hdeclxyB);
 547
 548       fDistr->Add(hnormdeclxyB);
 549
 550       fDistr->Add(hdeclxyd0d0B);
 551       fDistr->Add(hnormdeclxyd0d0B);
 552
 553       fDistr->Add(tmpMB);
 554
 555       //histograms filled only when the secondary vertex is recalculated w/o the daughter tracks (as requested in the cut object)
 556
 557       if(fCuts->GetIsPrimaryWithoutDaughters()){
 558         if(fReadMC){
 559           namedistr="hd0vpiS_";
 560           namedistr+=i;
 561           TH1F *hd0vpiS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions)(vtx w/o these tracks);d0(#pi) [cm]",200,-0.1,0.1);
 562
 563           namedistr="hd0vKS_";
 564           namedistr+=i;
 565           TH1F *hd0vKS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons) (vtx w/o these tracks);d0(K) [cm]",200,-0.1,0.1);
 566
 567           namedistr="hd0d0vS_";
 568           namedistr+=i;
 569           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);
 570
 571           namedistr="hdeclvS_";
 572           namedistr+=i;
 573           TH1F *hdeclengthvS = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
 574
 575           namedistr="hnormdeclvS_";
 576           namedistr+=i;
 577           TH1F *hnormdeclengthvS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
 578
 579           fDistr->Add(hd0vpiS);
 580           fDistr->Add(hd0vKS);
 581           fDistr->Add(hd0d0vS);
 582           fDistr->Add(hdeclengthvS);
 583           fDistr->Add(hnormdeclengthvS);
 584
 585         }
 586
 587         namedistr="hd0vmoresB_";
 588         namedistr+=i;
 589         TH1F *hd0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
 590
 591         namedistr="hd0d0vmoresB_";
 592         namedistr+=i;
 593         TH1F *hd0d0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
 594
 595
 596         namedistr="hd0vB_";
 597         namedistr+=i;
 598         TH1F *hd0vB = new TH1F(namedistr.Data(), "Impact parameter distribution (vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
 599
 600         namedistr="hd0vp0B_";
 601         namedistr+=i;
 602         TH1F *hd0vp0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong + ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
 603
 604         namedistr="hd0vp1B_";
 605         namedistr+=i;
 606         TH1F *hd0vp1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong - ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
 607
 608         namedistr="hd0d0vB_";
 609         namedistr+=i;
 610         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);
 611
 612         namedistr="hdeclvB_";
 613         namedistr+=i;
 614         TH1F *hdeclengthvB = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
 615
 616         namedistr="hnormdeclvB_";
 617         namedistr+=i;
 618         TH1F *hnormdeclengthvB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
 619
 620         fDistr->Add(hd0vB);
 621         fDistr->Add(hd0vp0B);
 622         fDistr->Add(hd0vp1B);
 623         fDistr->Add(hd0vmoresB);
 624
 625         fDistr->Add(hd0d0vB);
 626         fDistr->Add(hd0d0vmoresB);
 627
 628         fDistr->Add(hdeclengthvB);
 629
 630         fDistr->Add(hnormdeclengthvB);
 631       }
 632
 633
 634     }
 635     //histograms of invariant mass distributions
 636
 637
 638     //MC signal
 639     if(fReadMC){
 640       TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",600,1.5648,2.1648);
 641
 642       TH1F *tmpSl=(TH1F*)tmpSt->Clone();
 643       tmpSt->Sumw2();
 644       tmpSl->Sumw2();
 645
 646       //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 647       TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",600,1.5648,2.1648);
 648       //TH1F *tmpRl=(TH1F*)tmpRt->Clone();
 649       TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",600,1.5648,2.1648);
 650       //TH1F *tmpBl=(TH1F*)tmpBt->Clone();
 651       tmpBt->Sumw2();
 652       //tmpBl->Sumw2();
 653       tmpRt->Sumw2();
 654       //tmpRl->Sumw2();
 655
 656       fOutputMass->Add(tmpSt);
 657       fOutputMass->Add(tmpRt);
 658       fOutputMass->Add(tmpBt);
 659
 660     }
 661
 662     //mass
 663     TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",600,1.5648,2.1648);
 664     //TH1F *tmpMl=(TH1F*)tmpMt->Clone();
 665     tmpMt->Sumw2();
 666     //tmpMl->Sumw2();
 667     //distribution w/o cuts large range
 668     // TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
 669
 670     fOutputMass->Add(tmpMt);
 671
 672     if(fSys==0){ //histograms filled only in pp to save time in PbPb
 673       if(fFillVarHists){
 674
 675         if(fReadMC){
 676           //  pT
 677           namedistr="hptpiS_";
 678           namedistr+=i;
 679           TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
 680
 681           namedistr="hptKS_";
 682           namedistr+=i;
 683           TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
 684
 685           //  costhetastar
 686           namedistr="hcosthetastarS_";
 687           namedistr+=i;
 688           TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 689
 690           //pT no mass cut
 691
 692           namedistr="hptpiSnoMcut_";
 693           namedistr+=i;
 694           TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
 695
 696           namedistr="hptKSnoMcut_";
 697           namedistr+=i;
 698           TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
 699
 700           fDistr->Add(hptpiS);
 701           fDistr->Add(hptKS);
 702           fDistr->Add(hcosthetastarS);
 703
 704           fDistr->Add(hptpiSnoMcut);
 705           fDistr->Add(hptKSnoMcut);
 706
 707           //  costhetapoint vs d0 or d0d0
 708           namedistr="hcosthpointd0S_";
 709           namedistr+=i;
 710           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);
 711           namedistr="hcosthpointd0d0S_";
 712           namedistr+=i;
 713           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);
 714
 715           fDistr->Add(hcosthpointd0S);
 716           fDistr->Add(hcosthpointd0d0S);
 717
 718           //to compare with AliAnalysisTaskCharmFraction
 719           TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",600,1.5648,2.1648);
 720           TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
 721           tmpS27t->Sumw2();
 722           tmpS27l->Sumw2();
 723
 724           fOutputMass->Add(tmpS27t);
 725           fOutputMass->Add(tmpS27l);
 726
 727         }
 728
 729         //  pT
 730         namedistr="hptB_";
 731         namedistr+=i;
 732         TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 733
 734         //  costhetastar
 735         namedistr="hcosthetastarB_";
 736         namedistr+=i;
 737         TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 738
 739         //pT no mass cut
 740         namedistr="hptB1prongnoMcut_";
 741         namedistr+=i;
 742         TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 743
 744         namedistr="hptB2prongsnoMcut_";
 745         namedistr+=i;
 746         TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 747
 748         fDistr->Add(hptB);
 749         fDistr->Add(hcosthetastarB);
 750
 751         fDistr->Add(hptB1pnoMcut);
 752         fDistr->Add(hptB2pnoMcut);
 753
 754         //impact parameter of negative/positive track
 755         namedistr="hd0p0B_";
 756         namedistr+=i;
 757         TH1F *hd0p0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.1,0.1);
 758
 759         namedistr="hd0p1B_";
 760         namedistr+=i;
 761         TH1F *hd0p1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong -);d0 [cm]",200,-0.1,0.1);
 762
 763         //impact parameter corrected for strangeness
 764         namedistr="hd0moresB_";
 765         namedistr+=i;
 766         TH1F *hd0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
 767
 768         namedistr="hd0d0moresB_";
 769         namedistr+=i;
 770         TH1F *hd0d0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
 771
 772
 773         namedistr="hcosthetapointmoresB_";
 774         namedistr+=i;
 775         TH1F *hcosthetapointmoresB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 776
 777         //  costhetapoint vs d0 or d0d0
 778         namedistr="hcosthpointd0B_";
 779         namedistr+=i;
 780         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);
 781
 782         namedistr="hcosthpointd0d0B_";
 783         namedistr+=i;
 784         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);
 785
 786         fDistr->Add(hd0p0B);
 787         fDistr->Add(hd0p1B);
 788
 789         fDistr->Add(hd0moresB);
 790         fDistr->Add(hd0d0moresB);
 791         fDistr->Add(hcosthetapointmoresB);
 792
 793
 794         fDistr->Add(hcosthpointd0B);
 795
 796
 797         fDistr->Add(hcosthpointd0d0B);
 798       }
 799
 800     } //end pp histo
 801   }
 802
 803
 804   //for Like sign analysis
 805
 806   if(fArray==1){
 807     namedistr="hpospair";
 808     TH1F* hpospair=new TH1F(namedistr.Data(),"Number of positive pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
 809     namedistr="hnegpair";
 810     TH1F* hnegpair=new TH1F(namedistr.Data(),"Number of negative pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
 811     fOutputMass->Add(hpospair);
 812     fOutputMass->Add(hnegpair);
 813   }
 814
 815
 816   // 2D Pt distributions and impact parameter histograms
 817   if(fFillPtHist) {
 818
 819     nameMassPt="histMassPt";
 820     nameSgnPt="histSgnPt";
 821     nameBkgPt="histBkgPt";
 822     nameRflPt="histRflPt";
 823
 824     //MC signal
 825     if(fReadMC){
 826       TH2F* tmpStPt = new TH2F(nameSgnPt.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.16484,nbins2dPt,binInPt,binFinPt);
 827       TH2F *tmpSlPt=(TH2F*)tmpStPt->Clone();
 828       tmpStPt->Sumw2();
 829       tmpSlPt->Sumw2();
 830
 831       //Reflection: histo filled with D0MassV1 which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 832       TH2F* tmpRtPt = new TH2F(nameRflPt.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
 833       TH2F* tmpBtPt = new TH2F(nameBkgPt.Data(), "Background invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
 834       tmpBtPt->Sumw2();
 835       tmpRtPt->Sumw2();
 836
 837       fOutputMassPt->Add(tmpStPt);
 838       fOutputMassPt->Add(tmpRtPt);
 839       fOutputMassPt->Add(tmpBtPt);
 840
 841       //       cout<<endl<<endl<<"***************************************"<<endl;
 842       //       cout << " created and added to the list "<< nameSgnPt.Data() <<" "<< tmpStPt <<
 843       //        ", "<<nameRflPt.Data() <<" " << tmpRtPt<<", "<<nameBkgPt.Data()<< " " << tmpBtPt <<endl;
 844       //       cout<<"***************************************"<<endl<<endl;
 845     }
 846
 847     TH2F* tmpMtPt = new TH2F(nameMassPt.Data(),"D^{0} invariant mass; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
 848     tmpMtPt->Sumw2();
 849
 850     fOutputMassPt->Add(tmpMtPt);
 851   }
 852
 853   if(fFillImpParHist) CreateImpactParameterHistos();
 854
 855   // 2D Y distributions
 856
 857   if(fFillYHist) {
 858     for(Int_t i=0;i<fCuts->GetNPtBins();i++){
 859       nameMassY="histMassY_";
 860       nameMassY+=i;
 861       nameSgnY="histSgnY_";
 862       nameSgnY+=i;
 863       nameBkgY="histBkgY_";
 864       nameBkgY+=i;
 865       nameRflY="histRflY_";
 866       nameRflY+=i;
 867       //MC signal
 868       if(fReadMC){
 869         TH2F* tmpStY = new TH2F(nameSgnY.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; y",600,1.5648,2.16484,nbins2dY,binInY,binFinY);
 870         tmpStY->Sumw2();
 871         //Reflection: histo filled with D0MassV1 which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 872         TH2F* tmpRtY = new TH2F(nameRflY.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; y",600,1.5648,2.1648,nbins2dY,binInY,binFinY);
 873         TH2F* tmpBtY = new TH2F(nameBkgY.Data(), "Background invariant mass - MC; M [GeV]; Entries; y",600,1.5648,2.1648,nbins2dY,binInY,binFinY);
 874         tmpBtY->Sumw2();
 875         tmpRtY->Sumw2();
 876
 877         fOutputMassY->Add(tmpStY);
 878         fOutputMassY->Add(tmpRtY);
 879         fOutputMassY->Add(tmpBtY);
 880       }
 881       TH2F* tmpMtY = new TH2F(nameMassY.Data(),"D^{0} invariant mass; M [GeV]; Entries; y",600,1.5648,2.1648,nbins2dY,binInY,binFinY);
 882       tmpMtY->Sumw2();
 883       fOutputMassY->Add(tmpMtY);
 884     }
 885   }
 886
 887
 888   const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
 889
 890   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);
 891
 892   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 893   fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
 894   if(fReadMC) fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
 895   else fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
 896   fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
 897   fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
 898   fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
 899   if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
 900   if(fFillVarHists || fPIDCheck){
 901     fNentries->GetXaxis()->SetBinLabel(8,"PID=0");
 902     fNentries->GetXaxis()->SetBinLabel(9,"PID=1");
 903     fNentries->GetXaxis()->SetBinLabel(10,"PID=2");
 904     fNentries->GetXaxis()->SetBinLabel(11,"PID=3");
 905   }
 906   if(fReadMC && fSys==0){
 907     fNentries->GetXaxis()->SetBinLabel(12,"K");
 908     fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
 909   }
 910   fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
 911   fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
 912   if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
 913   if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
 914   fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
 915   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
 916
 917   fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
 918   fCounter->Init();
 919
 920   //
 921   // Output slot 7 : tree of the candidate variables after track selection
 922   //
 923   nameoutput = GetOutputSlot(7)->GetContainer()->GetName();
 924   fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
 925   Int_t nVar = 15;
 926   fCandidateVariables = new Double_t [nVar];
 927   TString * fCandidateVariableNames = new TString[nVar];
 928   fCandidateVariableNames[0] = "massD0";
 929   fCandidateVariableNames[1] = "massD0bar";
 930   fCandidateVariableNames[2] = "pt";
 931   fCandidateVariableNames[3] = "dca";
 932   fCandidateVariableNames[4] = "costhsD0";
 933   fCandidateVariableNames[5] = "costhsD0bar";
 934   fCandidateVariableNames[6] = "ptk";
 935   fCandidateVariableNames[7] = "ptpi";
 936   fCandidateVariableNames[8] = "d0k";
 937   fCandidateVariableNames[9] = "d0pi";
 938   fCandidateVariableNames[10] = "d0xd0";
 939   fCandidateVariableNames[11] = "costhp";
 940   fCandidateVariableNames[12] = "costhpxy";
 941   fCandidateVariableNames[13] = "lxy";
 942   fCandidateVariableNames[14] = "specialcuts";
 943   for(Int_t ivar=0; ivar<nVar; ivar++){
 944     fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/d",fCandidateVariableNames[ivar].Data()));
 945   }
 946
 947   //
 948   // Output slot 8 : List for detector response histograms
 949   //
 950   if (fDrawDetSignal) {
 951     TH2F *TOFSigBefPID = new TH2F("TOFSigBefPID", "TOF signal of daughters before PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
 952     TH2F *TOFSigAftPID = new TH2F("TOFSigAftPID", "TOF signal after PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
 953
 954     TH2F *TPCSigBefPID = new TH2F("TPCSigBefPID", "TPC dE/dx before PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
 955     TH2F *TPCSigAftPID = new TH2F("TPCSigAftPID", "TPC dE/dx after PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
 956
 957     fDetSignal->Add(TOFSigBefPID);
 958     fDetSignal->Add(TOFSigAftPID);
 959     fDetSignal->Add(TPCSigBefPID);
 960     fDetSignal->Add(TPCSigAftPID);
 961   }
 962
 963   // Post the data
 964   PostData(1,fOutputMass);
 965   PostData(2,fDistr);
 966   PostData(3,fNentries);
 967   PostData(5,fCounter);
 968   PostData(6,fOutputMassPt);
 969   PostData(7,fVariablesTree);
 970   PostData(8, fDetSignal);
 971   PostData(9,fOutputMassY);
 972
 973   return;
 974 }
 975
 976 //________________________________________________________________________
 977 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
 978 {
 979   // Execute analysis for current event:
 980   // heavy flavor candidates association to MC truth
 981   //cout<<"I'm in UserExec"<<endl;
 982
 983
 984   //cuts order
 985   //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 986   //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 987   //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 988   //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 989   //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 990   //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 991   //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 992   //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 993   //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 994
 995
 996   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 997
 998   TString bname;
 999   if(fArray==0){ //D0 candidates
1000     // load D0->Kpi candidates
1001     //cout<<"D0 candidates"<<endl;
1002     bname="D0toKpi";
1003   } else { //LikeSign candidates
1004     // load like sign candidates
1005     //cout<<"LS candidates"<<endl;
1006     bname="LikeSign2Prong";
1007   }
1008
1009   TClonesArray *inputArray=0;
1010
1011   if(!aod && AODEvent() && IsStandardAOD()) {
1012     // In case there is an AOD handler writing a standard AOD, use the AOD
1013     // event in memory rather than the input (ESD) event.
1014     aod = dynamic_cast<AliAODEvent*> (AODEvent());
1015     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
1016     // have to taken from the AOD event hold by the AliAODExtension
1017     AliAODHandler* aodHandler = (AliAODHandler*)
1018       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
1019
1020     if(aodHandler->GetExtensions()) {
1021       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
1022       AliAODEvent* aodFromExt = ext->GetAOD();
1023       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
1024     }
1025   } else if(aod) {
1026     inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
1027   }
1028
1029
1030   if(!inputArray || !aod) {
1031     printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
1032     return;
1033   }
1034
1035   // fix for temporary bug in ESDfilter
1036   // the AODs with null vertex pointer didn't pass the PhysSel
1037   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
1038
1039
1040   TClonesArray *mcArray = 0;
1041   AliAODMCHeader *mcHeader = 0;
1042
1043   if(fReadMC) {
1044     // load MC particles
1045     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1046     if(!mcArray) {
1047       printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
1048       return;
1049     }
1050
1051     // load MC header
1052     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1053     if(!mcHeader) {
1054       printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
1055       return;
1056     }
1057   }
1058
1059   //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
1060
1061   //histogram filled with 1 for every AOD
1062   fNentries->Fill(0);
1063   fCounter->StoreEvent(aod,fCuts,fReadMC);
1064   //fCounter->StoreEvent(aod,fReadMC);
1065
1066   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
1067   TString trigclass=aod->GetFiredTriggerClasses();
1068   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
1069
1070   if(!fCuts->IsEventSelected(aod)) {
1071     if(fCuts->GetWhyRejection()==1) // rejected for pileup
1072       fNentries->Fill(13);
1073     if(fSys==1 && (fCuts->GetWhyRejection()==2 || fCuts->GetWhyRejection()==3)) fNentries->Fill(15);
1074     if(fCuts->GetWhyRejection()==7) fNentries->Fill(17);
1075     return;
1076   }
1077
1078   // Check the Nb of SDD clusters
1079   if (fIsRejectSDDClusters) {
1080     Bool_t skipEvent = kFALSE;
1081     Int_t ntracks = 0;
1082     if (aod) ntracks = aod->GetNumberOfTracks();
1083     for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
1084       //    ... get the track
1085       AliAODTrack * track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrack));
1086       if(!track) AliFatal("Not a standard AOD");
1087       if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
1088         skipEvent=kTRUE;
1089         fNentries->Fill(16);
1090         break;
1091       }
1092     }
1093     if (skipEvent) return;
1094   }
1095
1096   // AOD primary vertex
1097   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
1098
1099   Bool_t isGoodVtx=kFALSE;
1100
1101   //vtx1->Print();
1102   TString primTitle = vtx1->GetTitle();
1103   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
1104     isGoodVtx=kTRUE;
1105     fNentries->Fill(3);
1106   }
1107
1108   // loop over candidates
1109   Int_t nInD0toKpi = inputArray->GetEntriesFast();
1110   if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
1111
1112   // FILE *f=fopen("4display.txt","a");
1113   // fprintf(f,"Number of D0->Kpi: %d\n",nInD0toKpi);
1114   Int_t nSelectedloose=0,nSelectedtight=0;
1115   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
1116     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
1117
1118     if(fUseSelectionBit && d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
1119         fNentries->Fill(2);
1120         continue; //skip the D0 from Dstar
1121       }
1122
1123     // Bool_t unsetvtx=kFALSE;
1124     // if(!d->GetOwnPrimaryVtx()) {
1125     //   d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
1126     //   unsetvtx=kTRUE;
1127     // }
1128
1129
1130     if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
1131       nSelectedloose++;
1132       nSelectedtight++;
1133       if(fSys==0){
1134         if(fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
1135       }
1136       Int_t ptbin=fCuts->PtBin(d->Pt());
1137       if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
1138       fIsSelectedCandidate=fCuts->IsSelected(d,AliRDHFCuts::kAll,aod); //selected
1139       if(fFillVarHists) {
1140         //if(!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate)) {
1141         fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(0),0);
1142         fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(1),1);
1143         //check daughters
1144         if(!fDaughterTracks.UncheckedAt(0) || !fDaughterTracks.UncheckedAt(1)) {
1145           AliDebug(1,"at least one daughter not found!");
1146           fNentries->Fill(5);
1147           fDaughterTracks.Clear();
1148           continue;
1149         }
1150         //}
1151         FillVarHists(aod,d,mcArray,fCuts,fDistr);
1152       }
1153
1154       if (fDrawDetSignal) {
1155         DrawDetSignal(d, fDetSignal);
1156       }
1157
1158       FillMassHists(d,mcArray,mcHeader,fCuts,fOutputMass);
1159       if (fPIDCheck) {
1160         Int_t isSelectedPIDfill = 3;
1161         if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(d); //0 rejected,1 D0,2 Dobar, 3 both
1162
1163         if (isSelectedPIDfill == 0)fNentries->Fill(7);
1164         if (isSelectedPIDfill == 1)fNentries->Fill(8);
1165         if (isSelectedPIDfill == 2)fNentries->Fill(9);
1166         if (isSelectedPIDfill == 3)fNentries->Fill(10);
1167       }
1168
1169     }
1170
1171     fDaughterTracks.Clear();
1172     //if(unsetvtx) d->UnsetOwnPrimaryVtx();
1173   } //end for prongs
1174   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1175   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1176   // Post the data
1177   PostData(1,fOutputMass);
1178   PostData(2,fDistr);
1179   PostData(3,fNentries);
1180   PostData(5,fCounter);
1181   PostData(6,fOutputMassPt);
1182   PostData(7,fVariablesTree);
1183   PostData(8, fDetSignal);
1184
1185   return;
1186 }
1187
1188 //____________________________________________________________________________
1189 void AliAnalysisTaskSED0Mass::DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal)
1190 {
1191   //
1192   // Function called in UserExec for drawing detector signal histograms:
1193   //
1194   fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(0), 0);
1195   fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(1), 1);
1196
1197   AliESDtrack *esdtrack1 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1198   AliESDtrack *esdtrack2 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1199
1200
1201   //For filling detector histograms
1202   Int_t isSelectedPIDfill = 3;
1203   if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1204
1205
1206   //fill "before PID" histos with every daughter
1207   ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1208   ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1209   ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1210   ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1211
1212   if (isSelectedPIDfill != 0)  { //fill "After PID" with everything that's not rejected
1213     ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1214     ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1215     ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1216     ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1217
1218   }
1219
1220   delete esdtrack1;
1221   delete esdtrack2;
1222
1223   return;
1224 }
1225
1226 //____________________________________________________________________________
1227 void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
1228   //
1229   // function used in UserExec to fill variable histograms:
1230   //
1231
1232
1233   Int_t pdgDgD0toKpi[2]={321,211};
1234   Int_t lab=-9999;
1235   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)
1236   //Double_t pt = d->Pt(); //mother pt
1237   Int_t isSelectedPID=3;
1238   if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1239   if (!fPIDCheck) {  //if fPIDCheck, this is already filled elsewhere
1240     if (isSelectedPID==0)fNentries->Fill(7);
1241     if (isSelectedPID==1)fNentries->Fill(8);
1242     if (isSelectedPID==2)fNentries->Fill(9);
1243     if (isSelectedPID==3)fNentries->Fill(10);
1244     //fNentries->Fill(8+isSelectedPID);
1245   }
1246
1247   if(fCutOnDistr && !fIsSelectedCandidate) return;
1248   //printf("\nif no cuts or cuts passed\n");
1249
1250
1251   //add distr here
1252   UInt_t pdgs[2];
1253
1254   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1255   pdgs[0]=211;
1256   pdgs[1]=321;
1257   Double_t minvD0 = part->InvMassD0();
1258   pdgs[1]=211;
1259   pdgs[0]=321;
1260   Double_t minvD0bar = part->InvMassD0bar();
1261   //cout<<"inside mass cut"<<endl;
1262
1263   Double_t invmasscut=0.03;
1264
1265   TString fillthispi="",fillthisK="",fillthis="", fillthispt="", fillthisetaphi="";
1266
1267   Int_t ptbin=cuts->PtBin(part->Pt());
1268   Double_t pt = part->Pt();
1269
1270   Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
1271   dz1[0]=-99; dz2[0]=-99;
1272   Double_t d0[2];
1273   Double_t decl[2]={-99,-99};
1274   Bool_t recalcvtx=kFALSE;
1275
1276
1277
1278   if(fCuts->GetIsPrimaryWithoutDaughters()){
1279     recalcvtx=kTRUE;
1280     //recalculate vertex
1281     AliAODVertex *vtxProp=0x0;
1282     vtxProp=GetPrimaryVtxSkipped(aod);
1283     if(vtxProp) {
1284       part->SetOwnPrimaryVtx(vtxProp);
1285       //Bool_t unsetvtx=kTRUE;
1286       //Calculate d0 for daughter tracks with Vtx Skipped
1287       AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1288       AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1289       esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
1290       esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
1291       delete vtxProp; vtxProp=NULL;
1292       delete esdtrack1;
1293       delete esdtrack2;
1294     }
1295
1296     d0[0]=dz1[0];
1297     d0[1]=dz2[0];
1298
1299     decl[0]=part->DecayLength2();
1300     decl[1]=part->NormalizedDecayLength2();
1301     part->UnsetOwnPrimaryVtx();
1302
1303   }
1304
1305   Double_t cosThetaStarD0 = 99;
1306   Double_t cosThetaStarD0bar = 99;
1307   Double_t cosPointingAngle = 99;
1308   Double_t normalizedDecayLength2 = -1, normalizedDecayLengthxy=-1;
1309   Double_t decayLength2 = -1, decayLengthxy=-1;
1310   Double_t ptProng[2]={-99,-99};
1311   Double_t d0Prong[2]={-99,-99};
1312   Double_t etaD = 99.;
1313   Double_t phiD = 99.;
1314
1315
1316   //disable the PID
1317   if(!fUsePid4Distr) isSelectedPID=0;
1318   if((lab>=0 && fReadMC) || (!fReadMC && isSelectedPID)){ //signal (from MC or PID)
1319
1320     //check pdg of the prongs
1321     AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1322     AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
1323
1324     if(!prong0 || !prong1) {
1325       return;
1326     }
1327     Int_t labprong[2];
1328     if(fReadMC){
1329       labprong[0]=prong0->GetLabel();
1330       labprong[1]=prong1->GetLabel();
1331     }
1332     AliAODMCParticle *mcprong=0;
1333     Int_t pdgProng[2]={0,0};
1334
1335     for (Int_t iprong=0;iprong<2;iprong++){
1336       if(fReadMC && labprong[iprong]>=0) {
1337         mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1338         pdgProng[iprong]=mcprong->GetPdgCode();
1339       }
1340     }
1341
1342     if(fSys==0){
1343       //no mass cut ditributions: ptbis
1344
1345       fillthispi="hptpiSnoMcut_";
1346       fillthispi+=ptbin;
1347
1348       fillthisK="hptKSnoMcut_";
1349       fillthisK+=ptbin;
1350
1351       if ((TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321)
1352           || (isSelectedPID==1 || isSelectedPID==3)){
1353         ((TH1F*)listout->FindObject(fillthispi))->Fill(prong0->Pt());
1354         ((TH1F*)listout->FindObject(fillthisK))->Fill(prong1->Pt());
1355       }
1356
1357       if ((TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211)
1358           || (isSelectedPID==2 || isSelectedPID==3)){
1359         ((TH1F*)listout->FindObject(fillthisK))->Fill(prong0->Pt());
1360         ((TH1F*)listout->FindObject(fillthispi))->Fill(prong1->Pt());
1361       }
1362     }
1363
1364     //no mass cut ditributions: mass
1365
1366     etaD = part->Eta();
1367     phiD = part->Phi();
1368
1369
1370     fillthis="hMassS_";
1371     fillthis+=ptbin;
1372     fillthispt="histSgnPt";
1373
1374     if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)
1375         || (!fReadMC && (isSelectedPID==1 || isSelectedPID==3))){//D0
1376       ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1377       if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1378
1379       fillthisetaphi="hetaphiD0candidateS_";
1380       fillthisetaphi+=ptbin;
1381       ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1382
1383       if(TMath::Abs(minvD0-mPDG)<0.05){
1384         fillthisetaphi="hetaphiD0candidatesignalregionS_";
1385         fillthisetaphi+=ptbin;
1386         ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1387       }
1388
1389     }
1390     else { //D0bar
1391       if(fReadMC || (!fReadMC && isSelectedPID > 1)){
1392         ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1393         if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1394
1395         fillthisetaphi="hetaphiD0barcandidateS_";
1396         fillthisetaphi+=ptbin;
1397         ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1398
1399         if(TMath::Abs(minvD0bar-mPDG)<0.05){
1400           fillthisetaphi="hetaphiD0barcandidatesignalregionS_";
1401           fillthisetaphi+=ptbin;
1402           ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1403         }
1404
1405       }
1406     }
1407
1408     //apply cut on invariant mass on the pair
1409     if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1410
1411       cosThetaStarD0 = part->CosThetaStarD0();
1412       cosThetaStarD0bar = part->CosThetaStarD0bar();
1413       cosPointingAngle = part->CosPointingAngle();
1414       normalizedDecayLength2 = part->NormalizedDecayLength2();
1415       decayLength2 = part->DecayLength2();
1416       decayLengthxy = part->DecayLengthXY();
1417       normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1418
1419       ptProng[0]=prong0->Pt(); ptProng[1]=prong1->Pt();
1420       d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1421
1422       if(fArray==1) cout<<"LS signal: ERROR"<<endl;
1423       for (Int_t iprong=0; iprong<2; iprong++){
1424         AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1425         if (fReadMC) labprong[iprong]=prong->GetLabel();
1426
1427         //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
1428         Int_t pdgprong=0;
1429         if(fReadMC && labprong[iprong]>=0) {
1430           mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1431           pdgprong=mcprong->GetPdgCode();
1432         }
1433
1434         Bool_t isPionHere[2]={(isSelectedPID==1 || isSelectedPID==3) ? kTRUE : kFALSE,(isSelectedPID==2 || isSelectedPID==3) ? kTRUE : kFALSE};
1435
1436         if(TMath::Abs(pdgprong)==211 || isPionHere[iprong]) {
1437           //cout<<"pi"<<endl;
1438
1439           if(fSys==0){
1440             fillthispi="hptpiS_";
1441             fillthispi+=ptbin;
1442             ((TH1F*)listout->FindObject(fillthispi))->Fill(ptProng[iprong]);
1443           }
1444
1445           fillthispi="hd0piS_";
1446           fillthispi+=ptbin;
1447           ((TH1F*)listout->FindObject(fillthispi))->Fill(d0Prong[iprong]);
1448           if(recalcvtx) {
1449
1450             fillthispi="hd0vpiS_";
1451             fillthispi+=ptbin;
1452             ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
1453           }
1454
1455         }
1456
1457         if(TMath::Abs(pdgprong)==321 || !isPionHere[iprong]) {
1458           //cout<<"kappa"<<endl;
1459           if(fSys==0){
1460             fillthisK="hptKS_";
1461             fillthisK+=ptbin;
1462             ((TH1F*)listout->FindObject(fillthisK))->Fill(ptProng[iprong]);
1463           }
1464
1465
1466           fillthisK="hd0KS_";
1467           fillthisK+=ptbin;
1468           ((TH1F*)listout->FindObject(fillthisK))->Fill(d0Prong[iprong]);
1469           if (recalcvtx){
1470             fillthisK="hd0vKS_";
1471             fillthisK+=ptbin;
1472             ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
1473           }
1474         }
1475
1476         if(fSys==0){
1477           fillthis="hcosthpointd0S_";
1478           fillthis+=ptbin;
1479           ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[iprong]);
1480         }
1481       } //end loop on prongs
1482
1483       fillthis="hdcaS_";
1484       fillthis+=ptbin;
1485       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1486
1487       fillthis="hcosthetapointS_";
1488       fillthis+=ptbin;
1489       ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1490
1491       fillthis="hcosthetapointxyS_";
1492       fillthis+=ptbin;
1493       ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1494
1495
1496       fillthis="hd0d0S_";
1497       fillthis+=ptbin;
1498       ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1499
1500       fillthis="hdeclS_";
1501       fillthis+=ptbin;
1502       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
1503
1504       fillthis="hnormdeclS_";
1505       fillthis+=ptbin;
1506       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1507
1508       fillthis="hdeclxyS_";
1509       fillthis+=ptbin;
1510       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1511
1512       fillthis="hnormdeclxyS_";
1513       fillthis+=ptbin;
1514       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1515
1516       fillthis="hdeclxyd0d0S_";
1517       fillthis+=ptbin;
1518       ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1519
1520       fillthis="hnormdeclxyd0d0S_";
1521       fillthis+=ptbin;
1522       ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1523
1524       if(recalcvtx) {
1525         fillthis="hdeclvS_";
1526         fillthis+=ptbin;
1527         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1528
1529         fillthis="hnormdeclvS_";
1530         fillthis+=ptbin;
1531         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1532
1533         fillthis="hd0d0vS_";
1534         fillthis+=ptbin;
1535         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1536       }
1537
1538       if(fSys==0){
1539         fillthis="hcosthetastarS_";
1540         fillthis+=ptbin;
1541         if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)) ((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1542         else {
1543           if (fReadMC || isSelectedPID>1)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1544           if(isSelectedPID==1 || isSelectedPID==3)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1545         }
1546
1547         fillthis="hcosthpointd0d0S_";
1548         fillthis+=ptbin;
1549         ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1550       }
1551
1552       if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)){
1553         for(Int_t it=0; it<2; it++){
1554           fillthis="hptD0S_";
1555           fillthis+=ptbin;
1556           ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1557           fillthis="hphiD0S_";
1558           fillthis+=ptbin;
1559           ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1560           Int_t nPointsITS = 0;
1561           for (Int_t il=0; il<6; il++){
1562             if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1563           }
1564           fillthis="hNITSpointsD0vsptS_";
1565           fillthis+=ptbin;
1566           ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(),nPointsITS);
1567           fillthis="hNSPDpointsD0S_";
1568           fillthis+=ptbin;
1569           if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
1570             ((TH1I*)listout->FindObject(fillthis))->Fill(0);
1571           }
1572           if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
1573             ((TH1I*)listout->FindObject(fillthis))->Fill(1);
1574           }
1575           if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
1576             ((TH1I*)listout->FindObject(fillthis))->Fill(2);
1577           }
1578           if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
1579             ((TH1I*)listout->FindObject(fillthis))->Fill(3);
1580           }
1581           fillthis="hNclsD0vsptS_";
1582           fillthis+=ptbin;
1583           Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1584           Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->GetTPCNcls();
1585           ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1586         }
1587       }
1588       else {
1589         if (fReadMC || isSelectedPID>1){
1590           for(Int_t it=0; it<2; it++){
1591             fillthis="hptD0barS_";
1592             fillthis+=ptbin;
1593             ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1594             fillthis="hphiD0barS_";
1595             fillthis+=ptbin;
1596             ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1597             fillthis="hNclsD0barvsptS_";
1598             fillthis+=ptbin;
1599             Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1600             Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
1601             ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1602           }
1603         }
1604         if(isSelectedPID==1 || isSelectedPID==3){
1605           for(Int_t it=0; it<2; it++){
1606             fillthis="hptD0S_";
1607             fillthis+=ptbin;
1608             ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1609             fillthis="hphiD0S_";
1610             fillthis+=ptbin;
1611             ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1612             Int_t nPointsITS = 0;
1613             for (Int_t il=0; il<6; il++){
1614               if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1615             }
1616             fillthis="hNITSpointsD0vsptS_";
1617             fillthis+=ptbin;
1618             ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(), nPointsITS);
1619             fillthis="hNSPDpointsD0S_";
1620             fillthis+=ptbin;
1621             if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
1622               ((TH1I*)listout->FindObject(fillthis))->Fill(0);
1623             }
1624             if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
1625               ((TH1I*)listout->FindObject(fillthis))->Fill(1);
1626             }
1627             if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
1628               ((TH1I*)listout->FindObject(fillthis))->Fill(2);
1629             }
1630             if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
1631               ((TH1I*)listout->FindObject(fillthis))->Fill(3);
1632             }
1633             fillthis="hNclsD0vsptS_";
1634             fillthis+=ptbin;
1635             Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1636             Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->GetTPCNcls();
1637             ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1638           }
1639         }
1640       }
1641
1642
1643       } //end mass cut
1644
1645   } else{ //Background or LS
1646     //if(!fReadMC){
1647     //cout<<"is background"<<endl;
1648
1649     etaD = part->Eta();
1650     phiD = part->Phi();
1651
1652     //no mass cut distributions: mass, ptbis
1653     fillthis="hMassB_";
1654     fillthis+=ptbin;
1655     fillthispt="histBkgPt";
1656
1657     if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))) {
1658       ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1659       if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1660
1661       fillthisetaphi="hetaphiD0candidateB_";
1662       fillthisetaphi+=ptbin;
1663       ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1664
1665       if(TMath::Abs(minvD0-mPDG)<0.05){
1666         fillthisetaphi="hetaphiD0candidatesignalregionB_";
1667         fillthisetaphi+=ptbin;
1668         ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1669       }
1670     }
1671     if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
1672       ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1673       if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1674
1675       fillthisetaphi="hetaphiD0barcandidateB_";
1676       fillthisetaphi+=ptbin;
1677       ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1678
1679       if(TMath::Abs(minvD0bar-mPDG)<0.05){
1680         fillthisetaphi="hetaphiD0barcandidatesignalregionB_";
1681         fillthisetaphi+=ptbin;
1682         ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1683       }
1684
1685     }
1686     if(fSys==0){
1687       fillthis="hptB1prongnoMcut_";
1688       fillthis+=ptbin;
1689
1690       ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1691
1692       fillthis="hptB2prongsnoMcut_";
1693       fillthis+=ptbin;
1694       ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1695       ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1696     }
1697
1698
1699       //apply cut on invariant mass on the pair
1700       if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1701       if(fSys==0){
1702         ptProng[0]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt(); ptProng[1]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt();
1703         cosThetaStarD0 = part->CosThetaStarD0();
1704         cosThetaStarD0bar = part->CosThetaStarD0bar();
1705       }
1706
1707       cosPointingAngle = part->CosPointingAngle();
1708       normalizedDecayLength2 = part->NormalizedDecayLength2();
1709       decayLength2 = part->DecayLength2();
1710       decayLengthxy = part->DecayLengthXY();
1711       normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1712       d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1713
1714
1715       AliAODTrack *prongg=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1716       if(!prongg) {
1717         if(fDebug>2) cout<<"No daughter found";
1718         return;
1719       }
1720       else{
1721         if(fArray==1){
1722           if(prongg->Charge()==1) {
1723             //fTotPosPairs[ptbin]++;
1724             ((TH1F*)fOutputMass->FindObject("hpospair"))->Fill(ptbin);
1725           } else {
1726             //fTotNegPairs[ptbin]++;
1727             ((TH1F*)fOutputMass->FindObject("hnegpair"))->Fill(ptbin);
1728           }
1729         }
1730       }
1731
1732       //fill pt and phi distrib for prongs with M cut
1733
1734       if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))){
1735           for(Int_t it=0; it<2; it++){
1736             fillthis="hptD0B_";
1737             fillthis+=ptbin;
1738             ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1739             fillthis="hphiD0B_";
1740             fillthis+=ptbin;
1741             ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1742
1743             Int_t nPointsITS = 0;
1744             for (Int_t il=0; il<6; il++){
1745               if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1746             }
1747             fillthis="hNITSpointsD0vsptB_";
1748             fillthis+=ptbin;
1749             ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(), nPointsITS);
1750             fillthis="hNSPDpointsD0B_";
1751             fillthis+=ptbin;
1752             if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
1753               ((TH1I*)listout->FindObject(fillthis))->Fill(0);
1754             }
1755             if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
1756               ((TH1I*)listout->FindObject(fillthis))->Fill(1);
1757             }
1758             if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
1759               ((TH1I*)listout->FindObject(fillthis))->Fill(2);
1760             }
1761             if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
1762               ((TH1I*)listout->FindObject(fillthis))->Fill(3);
1763             }
1764             fillthis="hNclsD0vsptB_";
1765             fillthis+=ptbin;
1766             Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1767             Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
1768             ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1769           }
1770
1771
1772         }
1773
1774         if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
1775           for(Int_t it=0; it<2; it++){
1776             fillthis="hptD0barB_";
1777             fillthis+=ptbin;
1778             ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1779             fillthis="hphiD0barB_";
1780             fillthis+=ptbin;
1781             ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1782             fillthis="hNclsD0barvsptB_";
1783             fillthis+=ptbin;
1784             Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1785             Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
1786             ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1787           }
1788         }
1789
1790       fillthis="hd0B_";
1791       fillthis+=ptbin;
1792       ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1793       ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1794
1795       if(fReadMC){
1796         Int_t pdgMother[2]={0,0};
1797         Double_t factor[2]={1,1};
1798
1799         for(Int_t iprong=0;iprong<2;iprong++){
1800           AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1801           lab=prong->GetLabel();
1802           if(lab>=0){
1803             AliAODMCParticle* mcprong=(AliAODMCParticle*)arrMC->At(lab);
1804             if(mcprong){
1805               Int_t labmom=mcprong->GetMother();
1806               if(labmom>=0){
1807                 AliAODMCParticle* mcmother=(AliAODMCParticle*)arrMC->At(labmom);
1808                 if(mcmother) pdgMother[iprong]=mcmother->GetPdgCode();
1809               }
1810             }
1811           }
1812
1813           if(fSys==0){
1814
1815             fillthis="hd0moresB_";
1816             fillthis+=ptbin;
1817
1818             if(TMath::Abs(pdgMother[iprong])==310 || TMath::Abs(pdgMother[iprong])==130 || TMath::Abs(pdgMother[iprong])==321){ //K^0_S, K^0_L, K^+-
1819               if(ptProng[iprong]<=1)factor[iprong]=1./.7;
1820               else factor[iprong]=1./.6;
1821               fNentries->Fill(11);
1822             }
1823
1824             if(TMath::Abs(pdgMother[iprong])==3122) { //Lambda
1825               factor[iprong]=1./0.25;
1826               fNentries->Fill(12);
1827             }
1828             fillthis="hd0moresB_";
1829             fillthis+=ptbin;
1830
1831             ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[iprong],factor[iprong]);
1832
1833             if(recalcvtx){
1834               fillthis="hd0vmoresB_";
1835               fillthis+=ptbin;
1836               ((TH1F*)listout->FindObject(fillthis))->Fill(d0[iprong],factor[iprong]);
1837             }
1838           }
1839         } //loop on prongs
1840
1841         if(fSys==0){
1842           fillthis="hd0d0moresB_";
1843           fillthis+=ptbin;
1844           ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0(),factor[0]*factor[1]);
1845
1846           fillthis="hcosthetapointmoresB_";
1847           fillthis+=ptbin;
1848           ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,factor[0]*factor[1]);
1849
1850           if(recalcvtx){
1851             fillthis="hd0d0vmoresB_";
1852             fillthis+=ptbin;
1853             ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1],factor[0]*factor[1]);
1854           }
1855         }
1856       } //readMC
1857
1858       if(fSys==0){
1859         //normalise pt distr to half afterwards
1860         fillthis="hptB_";
1861         fillthis+=ptbin;
1862         ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[0]);
1863         ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[1]);
1864
1865         fillthis="hcosthetastarB_";
1866         fillthis+=ptbin;
1867         if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3)))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1868         if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1869
1870
1871         fillthis="hd0p0B_";
1872         fillthis+=ptbin;
1873         ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1874         fillthis="hd0p1B_";
1875         fillthis+=ptbin;
1876         ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1877
1878         fillthis="hcosthpointd0d0B_";
1879         fillthis+=ptbin;
1880         ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1881
1882         fillthis="hcosthpointd0B_";
1883         fillthis+=ptbin;
1884         ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[0]);
1885         ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[1]);
1886
1887
1888         if(recalcvtx){
1889
1890           fillthis="hd0vp0B_";
1891           fillthis+=ptbin;
1892           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1893           fillthis="hd0vp1B_";
1894           fillthis+=ptbin;
1895           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1896
1897           fillthis="hd0vB_";
1898           fillthis+=ptbin;
1899           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1900           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1901
1902         }
1903
1904       }
1905
1906       fillthis="hdcaB_";
1907       fillthis+=ptbin;
1908       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1909
1910       fillthis="hd0d0B_";
1911       fillthis+=ptbin;
1912       ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]*d0Prong[1]);
1913
1914       if(recalcvtx){
1915         fillthis="hd0d0vB_";
1916         fillthis+=ptbin;
1917         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1918       }
1919
1920       fillthis="hcosthetapointB_";
1921       fillthis+=ptbin;
1922       ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1923
1924       fillthis="hcosthetapointxyB_";
1925       fillthis+=ptbin;
1926       ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1927
1928       fillthis="hdeclB_";
1929       fillthis+=ptbin;
1930       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
1931
1932       fillthis="hnormdeclB_";
1933       fillthis+=ptbin;
1934       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1935
1936       fillthis="hdeclxyB_";
1937       fillthis+=ptbin;
1938       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1939
1940       fillthis="hnormdeclxyB_";
1941       fillthis+=ptbin;
1942       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1943
1944       fillthis="hdeclxyd0d0B_";
1945       fillthis+=ptbin;
1946       ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1947
1948       fillthis="hnormdeclxyd0d0B_";
1949       fillthis+=ptbin;
1950       ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1951
1952
1953       if(recalcvtx) {
1954
1955         fillthis="hdeclvB_";
1956         fillthis+=ptbin;
1957         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1958
1959         fillthis="hnormdeclvB_";
1960         fillthis+=ptbin;
1961         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1962
1963
1964       }
1965     }//mass cut
1966   }//else (background)
1967
1968   return;
1969 }
1970
1971 //____________________________________________________________________________
1972 void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi* cuts, TList *listout){
1973   //
1974   // function used in UserExec to fill mass histograms:
1975   //
1976
1977
1978   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1979
1980   //cout<<"is selected = "<<fIsSelectedCandidate<<endl;
1981
1982   // Fill candidate variable Tree (track selection, no candidate selection)
1983   if( fWriteVariableTree && !part->HasBadDaughters()
1984       && fCuts->AreDaughtersSelected(part) && fCuts->IsSelectedPID(part) ){
1985     fCandidateVariables[0] = part->InvMassD0();
1986     fCandidateVariables[1] = part->InvMassD0bar();
1987     fCandidateVariables[2] = part->Pt();
1988     fCandidateVariables[3] = part->GetDCA();
1989     Double_t ctsD0=0. ,ctsD0bar=0.; part->CosThetaStarD0(ctsD0,ctsD0bar);
1990     fCandidateVariables[4] = ctsD0;
1991     fCandidateVariables[5] = ctsD0bar;
1992     fCandidateVariables[6] = part->Pt2Prong(0);
1993     fCandidateVariables[7] = part->Pt2Prong(1);
1994     fCandidateVariables[8] = part->Getd0Prong(0);
1995     fCandidateVariables[9] = part->Getd0Prong(1);
1996     fCandidateVariables[10] = part->Prodd0d0();
1997     fCandidateVariables[11] = part->CosPointingAngle();
1998     fCandidateVariables[12] = part->CosPointingAngleXY();
1999     fCandidateVariables[13] = part->NormalizedDecayLengthXY();
2000     fCandidateVariables[14] = fCuts->IsSelectedSpecialCuts(part);
2001     fVariablesTree->Fill();
2002   }
2003
2004   //cout<<"check cuts = "<<endl;
2005   //cuts->PrintAll();
2006   if (!fIsSelectedCandidate){
2007     //cout<<"Not Selected"<<endl;
2008     //cout<<"Rejected because "<<cuts->GetWhy()<<endl;
2009     return;
2010   }
2011
2012
2013   if(fDebug>2)  cout<<"Candidate selected"<<endl;
2014
2015   Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
2016   //printf("SELECTED\n");
2017   Int_t ptbin=cuts->PtBin(part->Pt());
2018   Double_t pt = part->Pt();
2019   Double_t y = part->YD0();
2020
2021   Double_t impparXY=part->ImpParXY()*10000.;
2022   Double_t trueImpParXY=0.;
2023   Double_t arrayForSparse[3]={invmassD0,pt,impparXY};
2024   Double_t arrayForSparseTrue[3]={invmassD0,pt,trueImpParXY};
2025
2026
2027   // AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2028   // if(!prong) {
2029   //   AliDebug(2,"No daughter found");
2030   //   return;
2031   // }
2032   // else{
2033     // if(prong->Charge()==1) {
2034     //   ((TH1F*)fDistr->FindObject("hpospair"))->Fill(fCuts->GetNPtBins()+ptbin);
2035     //   //fTotPosPairs[ptbin]++;
2036     // } else {
2037     //   ((TH1F*)fDistr->FindObject("hnegpair"))->Fill(fCuts->GetNPtBins()+ptbin);
2038     //   //fTotNegPairs[ptbin]++;
2039     // }
2040   //  }
2041
2042   // for(Int_t it=0;it<2;it++){
2043
2044   //    //request on spd points to be addes
2045   //   if(/*nSPD==2 && */part->Pt() > 5. && (TMath::Abs(invmassD0-mPDG)<0.01 || TMath::Abs(invmassD0bar-mPDG)<0.01)){
2046   //     FILE *f=fopen("4display.txt","a");
2047   //     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());
2048   //     fclose(f);
2049   //     //printf("PrimVtx NContributors: %d \n Prongs Rel Angle: %f \n \n",ncont,relangle);
2050   //   }
2051   // }
2052
2053   TString fillthis="", fillthispt="", fillthismasspt="", fillthismassy="";
2054   Int_t pdgDgD0toKpi[2]={321,211};
2055   Int_t labD0=-1;
2056   Bool_t isPrimary=kTRUE;
2057   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)
2058
2059   //Define weights for Bayesian (if appropriate)
2060
2061   Double_t weigD0=1.;
2062   Double_t weigD0bar=1.;
2063   if (fCuts->GetCombPID() && (fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)) {
2064     weigD0=fCuts->GetWeightsNegative()[AliPID::kKaon] * fCuts->GetWeightsPositive()[AliPID::kPion];
2065     weigD0bar=fCuts->GetWeightsPositive()[AliPID::kKaon] * fCuts->GetWeightsNegative()[AliPID::kPion];
2066     if (weigD0 > 1.0 || weigD0 < 0.) {weigD0 = 0.;}
2067     if (weigD0bar > 1.0 || weigD0bar < 0.) {weigD0bar = 0.;} //Prevents filling with weight > 1, or < 0
2068   }
2069
2070   //count candidates selected by cuts
2071   fNentries->Fill(1);
2072   //count true D0 selected by cuts
2073   if (fReadMC && labD0>=0) fNentries->Fill(2);
2074
2075   if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
2076
2077     arrayForSparse[0]=invmassD0; arrayForSparse[2]=impparXY;
2078
2079     if(fReadMC){
2080       if(labD0>=0) {
2081         if(fArray==1) cout<<"LS signal ERROR"<<endl;
2082
2083         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
2084         Int_t pdgD0 = partD0->GetPdgCode();
2085         //      cout<<"pdg = "<<pdgD0<<endl;
2086
2087         if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
2088         if(!isPrimary)
2089           trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
2090         arrayForSparseTrue[0]=invmassD0; arrayForSparseTrue[2]=trueImpParXY;
2091
2092         if (pdgD0==421){ //D0
2093           //      cout<<"Fill S with D0"<<endl;
2094           fillthis="histSgn_";
2095           fillthis+=ptbin;
2096           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2097
2098           if(fFillPtHist){
2099             fillthismasspt="histSgnPt";
2100             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2101           }
2102           if(fFillImpParHist){
2103             if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0);
2104             else {
2105               fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0);
2106               fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0);
2107             }
2108           }
2109
2110           if(fFillYHist){
2111             fillthismassy="histSgnY_";
2112             fillthismassy+=ptbin;
2113             ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2114           }
2115
2116           if(fSys==0){
2117             if(TMath::Abs(invmassD0 - mPDG) < 0.027 && fFillVarHists){
2118               fillthis="histSgn27_";
2119               fillthis+=ptbin;
2120               ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2121             }
2122           }
2123         } else{ //it was a D0bar
2124           fillthis="histRfl_";
2125           fillthis+=ptbin;
2126           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2127
2128           if(fFillPtHist){
2129             fillthismasspt="histRflPt";
2130             //      cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
2131             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2132           }
2133
2134           if(fFillYHist){
2135             fillthismassy="histRflY_";
2136             fillthismassy+=ptbin;
2137             //      cout << " Filling "<<fillthismassy<<" D0bar"<<endl;
2138             ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2139           }
2140
2141         }
2142       } else {//background
2143         fillthis="histBkg_";
2144         fillthis+=ptbin;
2145         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2146
2147         if(fFillPtHist){
2148           fillthismasspt="histBkgPt";
2149           //      cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
2150           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2151         }
2152         if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0);
2153
2154         if(fFillYHist){
2155           fillthismassy="histBkgY_";
2156           fillthismassy+=ptbin;
2157           //      cout << " Filling "<<fillthismassy<<" D0bar"<<endl;
2158           ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2159         }
2160
2161       }
2162
2163     }else{
2164       fillthis="histMass_";
2165       fillthis+=ptbin;
2166       //      cout<<"Filling "<<fillthis<<endl;
2167
2168       //      printf("Fill mass with D0");
2169        ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2170
2171
2172       if(fFillPtHist){
2173         fillthismasspt="histMassPt";
2174         //      cout<<"Filling "<<fillthismasspt<<endl;
2175         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2176       }
2177       if(fFillImpParHist) {
2178         //      cout << "Filling fHistMassPtImpParTC[0]"<<endl;
2179         fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0);
2180       }
2181
2182       if(fFillYHist){
2183         fillthismassy="histMassY_";
2184         fillthismassy+=ptbin;
2185         //      cout<<"Filling "<<fillthismassy<<endl;
2186         ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2187       }
2188
2189     }
2190
2191   }
2192   if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
2193
2194     arrayForSparse[0]=invmassD0bar; arrayForSparse[2]=impparXY;
2195
2196     if(fReadMC){
2197       if(labD0>=0) {
2198         if(fArray==1) cout<<"LS signal ERROR"<<endl;
2199         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
2200         Int_t pdgD0 = partD0->GetPdgCode();
2201         //      cout<<" pdg = "<<pdgD0<<endl;
2202
2203         if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
2204         if(!isPrimary)
2205           trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
2206         arrayForSparseTrue[0]=invmassD0bar; arrayForSparseTrue[2]=trueImpParXY;
2207
2208         if (pdgD0==-421){ //D0bar
2209           fillthis="histSgn_";
2210           fillthis+=ptbin;
2211           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
2212           // if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
2213           //   fillthis="histSgn27_";
2214           //   fillthis+=ptbin;
2215           //   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
2216           // }
2217
2218           if(fFillPtHist){
2219             fillthismasspt="histSgnPt";
2220             //      cout<<" Filling "<< fillthismasspt << endl;
2221             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2222           }
2223           if(fFillImpParHist){
2224             //      cout << " Filling impact parameter thnsparse"<<endl;
2225             if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0bar);
2226             else {
2227               fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0bar);
2228               fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0bar);
2229             }
2230           }
2231
2232           if(fFillYHist){
2233             fillthismassy="histSgnY_";
2234             fillthismassy+=ptbin;
2235             //      cout<<" Filling "<< fillthismassy << endl;
2236             ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2237           }
2238
2239         } else{
2240           fillthis="histRfl_";
2241           fillthis+=ptbin;
2242           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
2243           if(fFillPtHist){
2244             fillthismasspt="histRflPt";
2245             //      cout << " Filling "<<fillthismasspt<<endl;
2246             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2247           }
2248           if(fFillYHist){
2249             fillthismassy="histRflY_";
2250             fillthismassy+=ptbin;
2251             //      cout << " Filling "<<fillthismassy<<endl;
2252             ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2253           }
2254         }
2255       } else {//background or LS
2256         fillthis="histBkg_";
2257         fillthis+=ptbin;
2258         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
2259
2260         if(fFillPtHist){
2261           fillthismasspt="histBkgPt";
2262           //      cout<<" Filling "<< fillthismasspt << endl;
2263           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2264         }
2265         if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0bar);
2266         if(fFillYHist){
2267           fillthismassy="histBkgY_";
2268           fillthismassy+=ptbin;
2269           //      cout<<" Filling "<< fillthismassy << endl;
2270           ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2271         }
2272       }
2273     }else{
2274       fillthis="histMass_";
2275       fillthis+=ptbin;
2276       //      printf("Fill mass with D0bar");
2277
2278         ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar,weigD0bar);
2279
2280
2281       if(fFillPtHist){
2282         fillthismasspt="histMassPt";
2283         //      cout<<" Filling "<< fillthismasspt << endl;
2284         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2285       }
2286       if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0bar);
2287       if(fFillYHist){
2288         fillthismassy="histMassY_";
2289         fillthismassy+=ptbin;
2290         //      cout<<" Filling "<< fillthismassy << endl;
2291         ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2292       }
2293     }
2294   }
2295
2296   return;
2297 }
2298
2299 //__________________________________________________________________________
2300 AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev){
2301   //Calculate the primary vertex w/o the daughter tracks of the candidate
2302
2303   Int_t skipped[2];
2304   Int_t nTrksToSkip=2;
2305   AliAODTrack *dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2306   if(!dgTrack){
2307     AliDebug(2,"no daughter found!");
2308     return 0x0;
2309   }
2310   skipped[0]=dgTrack->GetID();
2311   dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(1);
2312   if(!dgTrack){
2313     AliDebug(2,"no daughter found!");
2314     return 0x0;
2315   }
2316   skipped[1]=dgTrack->GetID();
2317
2318   AliESDVertex *vertexESD=0x0;
2319   AliAODVertex *vertexAOD=0x0;
2320   AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
2321
2322   //
2323   vertexer->SetSkipTracks(nTrksToSkip,skipped);
2324   vertexer->SetMinClusters(4);
2325   vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
2326   if(!vertexESD) return vertexAOD;
2327   if(vertexESD->GetNContributors()<=0) {
2328     AliDebug(2,"vertexing failed");
2329     delete vertexESD; vertexESD=NULL;
2330     return vertexAOD;
2331   }
2332
2333   delete vertexer; vertexer=NULL;
2334
2335
2336   // convert to AliAODVertex
2337   Double_t pos[3],cov[6],chi2perNDF;
2338   vertexESD->GetXYZ(pos); // position
2339   vertexESD->GetCovMatrix(cov); //covariance matrix
2340   chi2perNDF = vertexESD->GetChi2toNDF();
2341   delete vertexESD; vertexESD=NULL;
2342
2343   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2344   return vertexAOD;
2345
2346 }
2347
2348
2349 //________________________________________________________________________
2350 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
2351 {
2352   // Terminate analysis
2353   //
2354   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
2355
2356
2357   fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
2358   if (!fOutputMass) {
2359     printf("ERROR: fOutputMass not available\n");
2360     return;
2361   }
2362   fOutputMassPt = dynamic_cast<TList*> (GetOutputData(6));
2363   if ((fFillPtHist || fFillImpParHist) && !fOutputMassPt) {
2364     printf("ERROR: fOutputMass not available\n");
2365     return;
2366   }
2367
2368   if(fFillVarHists){
2369     fDistr = dynamic_cast<TList*> (GetOutputData(2));
2370     if (!fDistr) {
2371       printf("ERROR: fDistr not available\n");
2372       return;
2373     }
2374   }
2375
2376   fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
2377
2378   if(!fNentries){
2379     printf("ERROR: fNEntries not available\n");
2380     return;
2381   }
2382   fCuts = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(4));
2383   if(!fCuts){
2384     printf("ERROR: fCuts not available\n");
2385     return;
2386   }
2387   fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(5));
2388   if (!fCounter) {
2389     printf("ERROR: fCounter not available\n");
2390     return;
2391   }
2392   if (fDrawDetSignal) {
2393     fDetSignal = dynamic_cast<TList*>(GetOutputData(8));
2394     if (!fDetSignal) {
2395       printf("ERROR: fDetSignal not available\n");
2396       return;
2397     }
2398   }
2399   if(fFillYHist){
2400     fOutputMassY = dynamic_cast<TList*> (GetOutputData(9));
2401     if (fFillYHist && !fOutputMassY) {
2402       printf("ERROR: fOutputMassY not available\n");
2403       return;
2404     }
2405   }
2406
2407   Int_t nptbins=fCuts->GetNPtBins();
2408   for(Int_t ipt=0;ipt<nptbins;ipt++){
2409
2410     if(fArray==1 && fFillVarHists){
2411       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
2412
2413
2414       if(fLsNormalization>1e-6) {
2415
2416         TString massName="histMass_";
2417         massName+=ipt;
2418         ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
2419
2420       }
2421
2422
2423       fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(ipt+1,ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(ipt+1,ipt+2));
2424       //fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
2425
2426       if(fLsNormalization>1e-6) {
2427
2428         TString nameDistr="hdcaB_";
2429         nameDistr+=ipt;
2430         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2431         nameDistr="hd0B_";
2432         nameDistr+=ipt;
2433         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2434         nameDistr="hd0d0B_";
2435         nameDistr+=ipt;
2436         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2437         nameDistr="hcosthetapointB_";
2438         nameDistr+=ipt;
2439         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2440         if(fSys==0){
2441           nameDistr="hptB_";
2442           nameDistr+=ipt;
2443           ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2444           nameDistr="hcosthetastarB_";
2445           nameDistr+=ipt;
2446           ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2447           nameDistr="hcosthpointd0d0B_";
2448           nameDistr+=ipt;
2449           ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
2450         }
2451       }
2452     }
2453   }
2454   TString cvname,cstname;
2455
2456   if (fArray==0){
2457     cvname="D0invmass";
2458     cstname="cstat0";
2459   } else {
2460     cvname="LSinvmass";
2461     cstname="cstat1";
2462   }
2463
2464   TCanvas *cMass=new TCanvas(cvname,cvname);
2465   cMass->cd();
2466   ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
2467
2468   TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
2469   cStat->cd();
2470   cStat->SetGridy();
2471   fNentries->Draw("htext0");
2472
2473   // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
2474   // ccheck->cd();
2475
2476   return;
2477 }
2478
2479
2480 //________________________________________________________________________
2481 void AliAnalysisTaskSED0Mass::CreateImpactParameterHistos(){
2482   // Histos for impact paramter study
2483
2484   Int_t nmassbins=200;
2485   Double_t fLowmasslimit=1.5648, fUpmasslimit=2.1648;
2486   Int_t fNImpParBins=400;
2487   Double_t fLowerImpPar=-2000., fHigherImpPar=2000.;
2488   Int_t nbins[3]={nmassbins,200,fNImpParBins};
2489   Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
2490   Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
2491
2492
2493   fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
2494                                         "Mass vs. pt vs.imppar - All",
2495                                         3,nbins,xmin,xmax);
2496   fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
2497                                         "Mass vs. pt vs.imppar - promptD",
2498                                         3,nbins,xmin,xmax);
2499   fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
2500                                         "Mass vs. pt vs.imppar - DfromB",
2501                                         3,nbins,xmin,xmax);
2502   fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
2503                                         "Mass vs. pt vs.true imppar -DfromB",
2504                                         3,nbins,xmin,xmax);
2505   fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
2506                                         "Mass vs. pt vs.imppar - backgr.",
2507                                         3,nbins,xmin,xmax);
2508
2509   for(Int_t i=0; i<5;i++){
2510     fOutputMassPt->Add(fHistMassPtImpParTC[i]);
2511   }
2512 }
2513
2514 //_________________________________________________________________________________________________
2515 Float_t AliAnalysisTaskSED0Mass::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD0) const {
2516   // true impact parameter calculation
2517
2518   printf(" AliAnalysisTaskSED0MassV1::GetTrueImpactParameter() \n");
2519
2520   Double_t vtxTrue[3];
2521   mcHeader->GetVertex(vtxTrue);
2522   Double_t origD[3];
2523   partD0->XvYvZv(origD);
2524   Short_t charge=partD0->Charge();
2525   Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
2526   for(Int_t iDau=0; iDau<2; iDau++){
2527     pXdauTrue[iDau]=0.;
2528     pYdauTrue[iDau]=0.;
2529     pZdauTrue[iDau]=0.;
2530   }
2531
2532   //  Int_t nDau=partD0->GetNDaughters();
2533   Int_t labelFirstDau = partD0->GetDaughter(0);
2534
2535   for(Int_t iDau=0; iDau<2; iDau++){
2536     Int_t ind = labelFirstDau+iDau;
2537     AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
2538     if(!part) continue;
2539     Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
2540     if(!part){
2541       AliError("Daughter particle not found in MC array");
2542       return 99999.;
2543     }
2544     if(pdgCode==211 || pdgCode==321){
2545       pXdauTrue[iDau]=part->Px();
2546       pYdauTrue[iDau]=part->Py();
2547       pZdauTrue[iDau]=part->Pz();
2548     }
2549   }
2550
2551   Double_t d0dummy[2]={0.,0.};
2552   AliAODRecoDecayHF aodDzeroMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
2553   return aodDzeroMC.ImpParXY();
2554
2555 }
2556
2557 //_________________________________________________________________________________________________
2558 Int_t AliAnalysisTaskSED0Mass::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
2559   //
2560   // checking whether the mother of the particles come from a charm or a bottom quark
2561   //
2562   printf(" AliAnalysisTaskSED0MassV1::CheckOrigin() \n");
2563
2564   Int_t pdgGranma = 0;
2565   Int_t mother = 0;
2566   mother = mcPartCandidate->GetMother();
2567   Int_t istep = 0;
2568   Int_t abspdgGranma =0;
2569   Bool_t isFromB=kFALSE;
2570   Bool_t isQuarkFound=kFALSE;
2571   while (mother >0 ){
2572     istep++;
2573     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
2574     if (mcGranma){
2575       pdgGranma = mcGranma->GetPdgCode();
2576       abspdgGranma = TMath::Abs(pdgGranma);
2577       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
2578         isFromB=kTRUE;
2579       }
2580       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
2581       mother = mcGranma->GetMother();
2582     }else{
2583       AliError("Failed casting the mother particle!");
2584       break;
2585     }
2586   }
2587
2588   if(isFromB) return 5;
2589   else return 4;
2590 }