PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2012, 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 (2Prongs)
  21 // and hadrons correlations
  22 //
  23 // Authors: Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
  24 // Fabio Colamaria, fabio.colamaria@ba.infn.it (correlations)
  25 /////////////////////////////////////////////////////////////
  26
  27 #include <Riostream.h>
  28 #include <TClonesArray.h>
  29 #include <TCanvas.h>
  30 #include <TNtuple.h>
  31 #include <TTree.h>
  32 #include <TList.h>
  33 #include <TH1F.h>
  34 #include <TH2F.h>
  35 #include <THnSparse.h>
  36 #include <TDatabasePDG.h>
  37
  38 #include <AliAnalysisDataSlot.h>
  39 #include <AliAnalysisDataContainer.h>
  40 #include "AliAnalysisManager.h"
  41 #include "AliESDtrack.h"
  42 #include "AliVertexerTracks.h"
  43 #include "AliAODHandler.h"
  44 #include "AliAODEvent.h"
  45 #include "AliAODVertex.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODMCHeader.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliAODRecoDecayHF2Prong.h"
  50 #include "AliAODRecoCascadeHF.h"
  51 #include "AliAnalysisVertexingHF.h"
  52 #include "AliAnalysisTaskSE.h"
  53 #include "AliAnalysisTaskSED0Correlations.h"
  54 #include "AliNormalizationCounter.h"
  55 #include "AliVertexingHFUtils.h"
  56
  57 using std::cout;
  58 using std::endl;
  59
  60 ClassImp(AliAnalysisTaskSED0Correlations)
  61
  62
  63 //________________________________________________________________________
  64 AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations():
  65 AliAnalysisTaskSE(),
  66   fNPtBinsCorr(0),
  67   fBinLimsCorr(),
  68   fPtThreshLow(),
  69   fPtThreshUp(),
  70   fEvents(0),
  71   fAlreadyFilled(kFALSE),
  72   fOutputMass(0),
  73   fOutputCorr(0),
  74   fOutputStudy(0),
  75   fNentries(0),
  76   fCutsD0(0),
  77   fCutsTracks(0),
  78   fCorrelatorTr(0),
  79   fCorrelatorKc(0),
  80   fCorrelatorK0(0),
  81   fReadMC(0),
  82   fRecoTr(kTRUE),
  83   fRecoD0(kTRUE),
  84   fSelEvType(kFALSE),
  85   fMixing(kFALSE),
  86   fCounter(0),
  87   fNPtBins(1),
  88   fFillOnlyD0D0bar(0),
  89   fIsSelectedCandidate(0),
  90   fSys(0),
  91   fEtaForCorrel(0),
  92   fIsRejectSDDClusters(0),
  93   fFillGlobal(kFALSE),
  94   fMultEv(0.),
  95   fSoftPiCut(kTRUE),
  96   fMEAxisThresh(kFALSE),
  97   fKaonCorr(kFALSE)
  98 {
  99   // Default constructor
 100
 101 }
 102
 103 //________________________________________________________________________
 104 AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const char *name,AliRDHFCutsD0toKpi* cutsD0, AliHFAssociatedTrackCuts* cutsTrk):
 105   AliAnalysisTaskSE(name),
 106   fNPtBinsCorr(0),
 107   fBinLimsCorr(),
 108   fPtThreshLow(),
 109   fPtThreshUp(),
 110   fEvents(0),
 111   fAlreadyFilled(kFALSE),
 112   fOutputMass(0),
 113   fOutputCorr(0),
 114   fOutputStudy(0),
 115   fNentries(0),
 116   fCutsD0(0),
 117   fCutsTracks(cutsTrk),
 118   fCorrelatorTr(0),
 119   fCorrelatorKc(0),
 120   fCorrelatorK0(0),
 121   fReadMC(0),
 122   fRecoTr(kTRUE),
 123   fRecoD0(kTRUE),
 124   fSelEvType(kFALSE),
 125   fMixing(kFALSE),
 126   fCounter(0),
 127   fNPtBins(1),
 128   fFillOnlyD0D0bar(0),
 129   fIsSelectedCandidate(0),
 130   fSys(0),
 131   fEtaForCorrel(0),
 132   fIsRejectSDDClusters(0),
 133   fFillGlobal(kFALSE),
 134   fMultEv(0.),
 135   fSoftPiCut(kTRUE),
 136   fMEAxisThresh(kFALSE),
 137   fKaonCorr(kFALSE)
 138 {
 139   // Default constructor
 140
 141   fNPtBins=cutsD0->GetNPtBins();
 142
 143   fCutsD0=cutsD0;
 144
 145   // Output slot #1 writes into a TList container (mass with cuts)
 146   DefineOutput(1,TList::Class());  //My private output
 147   // Output slot #2 writes into a TH1F container (number of events)
 148   DefineOutput(2,TH1F::Class());  //My private output
 149   // Output slot #3 writes into a AliRDHFD0toKpi container (cuts)
 150   DefineOutput(3,AliRDHFCutsD0toKpi::Class());  //My private output
 151   // Output slot #4 writes Normalization Counter
 152   DefineOutput(4,AliNormalizationCounter::Class());
 153   // Output slot #5 writes into a TList container (correl output)
 154   DefineOutput(5,TList::Class());  //My private output
 155   // Output slot #6 writes into a TList container (correl advanced)
 156   DefineOutput(6,TList::Class());  //My private output
 157   // Output slot #7 writes into a AliHFAssociatedTrackCuts container (cuts)
 158   DefineOutput(7,AliHFAssociatedTrackCuts::Class());  //My private output
 159 }
 160
 161 //________________________________________________________________________
 162 AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const AliAnalysisTaskSED0Correlations &source):
 163   AliAnalysisTaskSE(source),
 164   fNPtBinsCorr(source.fNPtBinsCorr),
 165   fBinLimsCorr(source.fBinLimsCorr),
 166   fPtThreshLow(source.fPtThreshLow),
 167   fPtThreshUp(source.fPtThreshUp),
 168   fEvents(source.fEvents),
 169   fAlreadyFilled(source.fAlreadyFilled),
 170   fOutputMass(source.fOutputMass),
 171   fOutputCorr(source.fOutputCorr),
 172   fOutputStudy(source.fOutputStudy),
 173   fNentries(source.fNentries),
 174   fCutsD0(source.fCutsD0),
 175   fCutsTracks(source.fCutsTracks),
 176   fCorrelatorTr(source.fCorrelatorTr),
 177   fCorrelatorKc(source.fCorrelatorKc),
 178   fCorrelatorK0(source.fCorrelatorK0),
 179   fReadMC(source.fReadMC),
 180   fRecoTr(source.fRecoTr),
 181   fRecoD0(source.fRecoD0),
 182   fSelEvType(source.fSelEvType),
 183   fMixing(source.fMixing),
 184   fCounter(source.fCounter),
 185   fNPtBins(source.fNPtBins),
 186   fFillOnlyD0D0bar(source.fFillOnlyD0D0bar),
 187   fIsSelectedCandidate(source.fIsSelectedCandidate),
 188   fSys(source.fSys),
 189   fEtaForCorrel(source.fEtaForCorrel),
 190   fIsRejectSDDClusters(source.fIsRejectSDDClusters),
 191   fFillGlobal(source.fFillGlobal),
 192   fMultEv(source.fMultEv),
 193   fSoftPiCut(source.fSoftPiCut),
 194   fMEAxisThresh(source.fMEAxisThresh),
 195   fKaonCorr(source.fKaonCorr)
 196 {
 197   // Copy constructor
 198 }
 199
 200 //________________________________________________________________________
 201 AliAnalysisTaskSED0Correlations::~AliAnalysisTaskSED0Correlations()
 202 {
 203   if (fOutputMass) {
 204     delete fOutputMass;
 205     fOutputMass = 0;
 206   }
 207   if (fOutputCorr) {
 208     delete fOutputCorr;
 209     fOutputCorr = 0;
 210   }
 211   if (fOutputStudy) {
 212     delete fOutputStudy;
 213     fOutputStudy = 0;
 214   }
 215   if (fCutsD0) {
 216     delete fCutsD0;
 217     fCutsD0 = 0;
 218   }
 219   if (fNentries){
 220     delete fNentries;
 221     fNentries = 0;
 222   }
 223   if (fCorrelatorTr) {
 224     delete fCorrelatorTr;
 225     fCorrelatorTr = 0;
 226   }
 227   if (fCorrelatorKc) {
 228     delete fCorrelatorKc;
 229     fCorrelatorKc = 0;
 230   }
 231   if (fCorrelatorK0) {
 232     delete fCorrelatorK0;
 233     fCorrelatorK0 = 0;
 234   }
 235   if (fCounter){
 236     delete fCounter;
 237     fCounter=0;
 238   }
 239 }
 240
 241 //______________________________________________________________________________
 242 AliAnalysisTaskSED0Correlations& AliAnalysisTaskSED0Correlations::operator=(const AliAnalysisTaskSED0Correlations& orig)
 243 {
 244 // Assignment
 245   if (&orig == this) return *this; //if address is the same (same object), returns itself
 246
 247   AliAnalysisTaskSE::operator=(orig); //Uses the AliAnalysisTaskSE operator to assign the inherited part of the class
 248   fNPtBinsCorr = orig.fNPtBinsCorr;
 249   fBinLimsCorr = orig.fBinLimsCorr;
 250   fPtThreshLow = orig.fPtThreshLow;
 251   fPtThreshUp = orig.fPtThreshUp;
 252   fEvents = orig.fEvents;
 253   fAlreadyFilled = orig.fAlreadyFilled;
 254   fOutputMass = orig.fOutputMass;
 255   fOutputCorr = orig.fOutputCorr;
 256   fOutputStudy = orig.fOutputStudy;
 257   fNentries = orig.fNentries;
 258   fCutsD0 = orig.fCutsD0;
 259   fCutsTracks = orig.fCutsTracks;
 260   fCorrelatorTr = orig.fCorrelatorTr;
 261   fCorrelatorKc = orig.fCorrelatorKc;
 262   fCorrelatorK0 = orig.fCorrelatorK0;
 263   fReadMC = orig.fReadMC;
 264   fRecoTr = orig.fRecoTr;
 265   fRecoD0 = orig.fRecoD0;
 266   fSelEvType = orig.fSelEvType;
 267   fMixing = orig.fMixing;
 268   fCounter = orig.fCounter;
 269   fNPtBins = orig.fNPtBins;
 270   fFillOnlyD0D0bar = orig.fFillOnlyD0D0bar;
 271   fIsSelectedCandidate = orig.fIsSelectedCandidate;
 272   fSys = orig.fSys;
 273   fEtaForCorrel = orig.fEtaForCorrel;
 274   fIsRejectSDDClusters = orig.fIsRejectSDDClusters;
 275   fFillGlobal = orig.fFillGlobal;
 276   fMultEv = orig.fMultEv;
 277   fSoftPiCut = orig.fSoftPiCut;
 278   fMEAxisThresh = orig.fMEAxisThresh;
 279   fKaonCorr = orig.fKaonCorr;
 280
 281   return *this; //returns pointer of the class
 282 }
 283
 284 //________________________________________________________________________
 285 void AliAnalysisTaskSED0Correlations::Init()
 286 {
 287   // Initialization
 288
 289   if(fDebug > 1) printf("AnalysisTaskSED0Correlations::Init() \n");
 290
 291   //Copy of cuts objects
 292   AliRDHFCutsD0toKpi* copyfCutsD0 = new AliRDHFCutsD0toKpi(*fCutsD0);
 293   const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
 294   copyfCutsD0->SetName(nameoutput);
 295
 296   //needer to clear completely the objects inside with Clear() method
 297   // Post the data
 298   PostData(3,copyfCutsD0);
 299   PostData(7,fCutsTracks);
 300
 301   return;
 302 }
 303
 304 //________________________________________________________________________
 305 void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
 306 {
 307
 308   // Create the output container
 309   //
 310   if(fDebug > 1) printf("AnalysisTaskSED0Correlations::UserCreateOutputObjects() \n");
 311
 312   //HFCorrelator creation and definition
 313   fCorrelatorTr = new AliHFCorrelator("CorrelatorTr",fCutsTracks,fSys,fCutsD0);//fSys=0 use multiplicity, =1 use centrality
 314   fCorrelatorKc = new AliHFCorrelator("CorrelatorKc",fCutsTracks,fSys,fCutsD0);
 315   fCorrelatorK0 = new AliHFCorrelator("CorrelatorK0",fCutsTracks,fSys,fCutsD0);
 316   fCorrelatorTr->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);// set the Delta Phi Interval you want (in this case -0.5Pi to 1.5 Pi)
 317   fCorrelatorKc->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);
 318   fCorrelatorK0->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);
 319   fCorrelatorTr->SetEventMixing(fMixing);// sets the analysis on a single event (kFALSE) or mixed events (kTRUE)
 320   fCorrelatorKc->SetEventMixing(fMixing);
 321   fCorrelatorK0->SetEventMixing(fMixing);
 322   fCorrelatorTr->SetAssociatedParticleType(1);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
 323   fCorrelatorKc->SetAssociatedParticleType(2);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
 324   fCorrelatorK0->SetAssociatedParticleType(3);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
 325   fCorrelatorTr->SetApplyDisplacementCut(2); //0: don't calculate d0; 1: return d0; 2: return d0/d0err
 326   fCorrelatorKc->SetApplyDisplacementCut(2);
 327   fCorrelatorK0->SetApplyDisplacementCut(0);
 328   fCorrelatorTr->SetUseMC(fReadMC);// sets Montecarlo flag
 329   fCorrelatorKc->SetUseMC(fReadMC);
 330   fCorrelatorK0->SetUseMC(fReadMC);
 331   fCorrelatorTr->SetUseReco(fRecoTr);// sets (if MC analysis) wheter to analyze Reco or Kinem tracks
 332   fCorrelatorKc->SetUseReco(fRecoTr);
 333   fCorrelatorK0->SetUseReco(fRecoTr);
 334   fCorrelatorKc->SetPIDmode(2); //switch for K+/- PID option
 335   Bool_t pooldefTr = fCorrelatorTr->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
 336   Bool_t pooldefKc = fCorrelatorKc->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
 337   Bool_t pooldefK0 = fCorrelatorK0->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
 338   if(!pooldefTr) AliInfo("Warning:: Event pool not defined properly");
 339   if(!pooldefKc) AliInfo("Warning:: Event pool not defined properly");
 340   if(!pooldefK0) AliInfo("Warning:: Event pool not defined properly");
 341
 342   // Several histograms are more conveniently managed in a TList
 343   fOutputMass = new TList();
 344   fOutputMass->SetOwner();
 345   fOutputMass->SetName("listMass");
 346
 347   fOutputCorr = new TList();
 348   fOutputCorr->SetOwner();
 349   fOutputCorr->SetName("correlationslist");
 350
 351   fOutputStudy = new TList();
 352   fOutputStudy->SetOwner();
 353   fOutputStudy->SetName("MCstudyplots");
 354
 355   TString nameMass=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassWg=" ",nameSgnWg=" ", nameBkgWg=" ", nameRflWg=" ";
 356
 357 //for origin c case (or for data)
 358   for(Int_t i=0;i<fCutsD0->GetNPtBins();i++){
 359
 360     nameMass="histMass_"; if(fReadMC) nameMass+="c_";
 361     nameMass+=i;
 362     nameMassWg="histMass_WeigD0Eff_"; if(fReadMC) nameMassWg+="c_";
 363     nameMassWg+=i;
 364     nameSgn="histSgn_"; if(fReadMC) nameSgn+="c_";
 365     nameSgn+=i;
 366     nameSgnWg="histSgn_WeigD0Eff_"; if(fReadMC) nameSgnWg+="c_";
 367     nameSgnWg+=i;
 368     nameBkg="histBkg_"; if(fReadMC) nameBkg+="c_";
 369     nameBkg+=i;
 370     nameBkgWg="histBkg_WeigD0Eff_"; if(fReadMC) nameBkgWg+="c_";
 371     nameBkgWg+=i;
 372     nameRfl="histRfl_"; if(fReadMC) nameRfl+="c_";
 373     nameRfl+=i;
 374     nameRflWg="histRfl_WeigD0Eff_"; if(fReadMC) nameRflWg+="c_";
 375     nameRflWg+=i;
 376
 377     //histograms of invariant mass distributions
 378
 379     //MC signal
 380     if(fReadMC){
 381       TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass c - MC; M [GeV]; Entries",120,1.5648,2.1648);
 382       TH1F* tmpStWg = new TH1F(nameSgnWg.Data(), "D^{0} invariant mass c - MC; M [GeV] - weight 1/D0eff; Entries",120,1.5648,2.1648);
 383       tmpSt->Sumw2();
 384       tmpStWg->Sumw2();
 385
 386       //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 387       TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass c - MC; M [GeV]; Entries",120,1.5648,2.1648);
 388       TH1F* tmpRtWg = new TH1F(nameRflWg.Data(), "Reflected signal invariant mass c - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
 389       TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass c - MC; M [GeV]; Entries",120,1.5648,2.1648);
 390       TH1F* tmpBtWg = new TH1F(nameBkgWg.Data(), "Background invariant mass c - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
 391       tmpBt->Sumw2();
 392       tmpBtWg->Sumw2();
 393       tmpRt->Sumw2();
 394       tmpRtWg->Sumw2();
 395       fOutputMass->Add(tmpSt);
 396       fOutputMass->Add(tmpStWg);
 397       fOutputMass->Add(tmpRt);
 398       fOutputMass->Add(tmpRtWg);
 399       fOutputMass->Add(tmpBt);
 400       fOutputMass->Add(tmpBtWg);
 401     }
 402
 403     //mass
 404     TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass c; M [GeV]; Entries",120,1.5648,2.1648);
 405     tmpMt->Sumw2();
 406     fOutputMass->Add(tmpMt);
 407     //mass weighted by 1/D0eff
 408     TH1F* tmpMtwg = new TH1F(nameMassWg.Data(),"D^{0} invariant mass c - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
 409     tmpMtwg->Sumw2();
 410     fOutputMass->Add(tmpMtwg);
 411   }
 412
 413 //for origin b case (no Bkg and Mass histos, here for weights you should use c+b efficiencies, while on data (on MC they're useless))
 414   for(Int_t i=0;i<fCutsD0->GetNPtBins();i++){
 415
 416     nameSgn="histSgn_b_";
 417     nameSgn+=i;
 418     nameSgnWg="histSgn_WeigD0Eff_b_";
 419     nameSgnWg+=i;
 420     nameRfl="histRfl_b_";
 421     nameRfl+=i;
 422     nameRflWg="histRfl_WeigD0Eff_b_";
 423     nameRflWg+=i;
 424
 425     //histograms of invariant mass distributions
 426
 427     //MC signal
 428     if(fReadMC){
 429       TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass b - MC; M [GeV]; Entries",120,1.5648,2.1648);
 430       TH1F* tmpStWg = new TH1F(nameSgnWg.Data(), "D^{0} invariant mass b - MC; M [GeV] - weight 1/D0eff; Entries",120,1.5648,2.1648);
 431       tmpSt->Sumw2();
 432       tmpStWg->Sumw2();
 433
 434       //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 435       TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass b - MC; M [GeV]; Entries",120,1.5648,2.1648);
 436       TH1F* tmpRtWg = new TH1F(nameRflWg.Data(), "Reflected signal invariant mass b - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
 437       tmpRt->Sumw2();
 438       tmpRtWg->Sumw2();
 439       fOutputMass->Add(tmpSt);
 440       fOutputMass->Add(tmpStWg);
 441       fOutputMass->Add(tmpRt);
 442       fOutputMass->Add(tmpRtWg);
 443     }
 444   }
 445
 446   const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
 447
 448   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", 20,-0.5,19.5);
 449
 450   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 451   fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
 452   fReadMC ? fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected") : fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
 453   fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
 454   fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
 455   fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
 456   if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
 457   if(fReadMC && fSys==0){
 458     fNentries->GetXaxis()->SetBinLabel(12,"K");
 459     fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
 460   }
 461   fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
 462   fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
 463   if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
 464   if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
 465   fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
 466   fNentries->GetXaxis()->SetBinLabel(19,"nEventsSelected");
 467   if(fReadMC) fNentries->GetXaxis()->SetBinLabel(20,"nEvsWithProdMech");
 468   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
 469
 470   fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(4)->GetContainer()->GetName()));
 471   fCounter->Init();
 472
 473   CreateCorrelationsObjs(); //creates histos for correlations analysis
 474
 475   // Post the data
 476   PostData(1,fOutputMass);
 477   PostData(2,fNentries);
 478   PostData(4,fCounter);
 479   PostData(5,fOutputCorr);
 480   PostData(6,fOutputStudy);
 481
 482   return;
 483 }
 484
 485 //________________________________________________________________________
 486 void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
 487 {
 488   // Execute analysis for current event:
 489   // heavy flavor candidates association to MC truth
 490   //cout<<"I'm in UserExec"<<endl;
 491
 492
 493   //cuts order
 494   //     printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 495   //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 496   //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 497   //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 498   //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 499   //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 500   //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 501   //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 502   //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 503
 504   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 505   fEvents++;
 506
 507   TString bname="D0toKpi";
 508
 509   TClonesArray *inputArray=0;
 510
 511   fMultEv = 0.; //reset event multiplicity
 512
 513   if(!aod && AODEvent() && IsStandardAOD()) {
 514     // In case there is an AOD handler writing a standard AOD, use the AOD
 515     // event in memory rather than the input (ESD) event.
 516     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 517     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 518     // have to taken from the AOD event hold by the AliAODExtension
 519     AliAODHandler* aodHandler = (AliAODHandler*)
 520       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 521
 522     if(aodHandler->GetExtensions()) {
 523       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 524       AliAODEvent* aodFromExt = ext->GetAOD();
 525       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
 526     }
 527   } else if(aod) {
 528     inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
 529   }
 530
 531   if(!inputArray || !aod) {
 532     printf("AliAnalysisTaskSED0Correlations::UserExec: input branch not found!\n");
 533     return;
 534   }
 535
 536   // fix for temporary bug in ESDfilter
 537   // the AODs with null vertex pointer didn't pass the PhysSel
 538   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 539
 540   TClonesArray *mcArray = 0;
 541   AliAODMCHeader *mcHeader = 0;
 542
 543   if(fReadMC) {
 544     // load MC particles
 545     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 546     if(!mcArray) {
 547       printf("AliAnalysisTaskSED0Correlations::UserExec: MC particles branch not found!\n");
 548       return;
 549     }
 550
 551     // load MC header
 552     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 553     if(!mcHeader) {
 554       printf("AliAnalysisTaskSED0Correlations::UserExec: MC header branch not found!\n");
 555       return;
 556     }
 557   }
 558
 559   //histogram filled with 1 for every AOD
 560   fNentries->Fill(0);
 561   fCounter->StoreEvent(aod,fCutsD0,fReadMC);
 562
 563   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
 564   TString trigclass=aod->GetFiredTriggerClasses();
 565   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
 566
 567   if(!fCutsD0->IsEventSelected(aod)) {
 568     if(fCutsD0->GetWhyRejection()==1) // rejected for pileup
 569       fNentries->Fill(13);
 570     if(fSys==1 && (fCutsD0->GetWhyRejection()==2 || fCutsD0->GetWhyRejection()==3)) fNentries->Fill(15);
 571     if(fCutsD0->GetWhyRejection()==7) fNentries->Fill(17);
 572     return;
 573   }
 574
 575   fNentries->Fill(18); //event selected after selection
 576
 577   //Setting PIDResponse for associated tracks
 578   fCorrelatorTr->SetPidAssociated();
 579   fCorrelatorKc->SetPidAssociated();
 580   fCorrelatorK0->SetPidAssociated();
 581
 582   //Selection on production type (MC)
 583   if(fReadMC && fSelEvType){
 584
 585     Bool_t isMCeventgood = kFALSE;
 586
 587     Int_t eventType = mcHeader->GetEventType();
 588     Int_t NMCevents = fCutsTracks->GetNofMCEventType();
 589
 590     for(Int_t k=0; k<NMCevents; k++){
 591       Int_t * MCEventType = fCutsTracks->GetMCEventType();
 592
 593       if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
 594       ((TH1D*)fOutputStudy->FindObject("EventTypeMC"))->Fill(eventType);
 595     }
 596
 597     if(NMCevents && !isMCeventgood){
 598       if(fDebug>2)std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
 599       return;
 600     }
 601     fNentries->Fill(19); //event with particular production type
 602
 603   } //end of selection
 604
 605
 606   // Check the Nb of SDD clusters
 607   if (fIsRejectSDDClusters) {
 608     Bool_t skipEvent = kFALSE;
 609     Int_t ntracks = 0;
 610     if (aod) ntracks = aod->GetNumberOfTracks();
 611     for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
 612       //    ... get the track
 613       AliAODTrack * track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrack));
 614       if(!track){
 615         AliWarning("Error in casting to AOD track. Not a standard AOD?");
 616         continue;
 617       }
 618       if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
 619         skipEvent=kTRUE;
 620         fNentries->Fill(16);
 621         break;
 622       }
 623     }
 624     if (skipEvent) return;
 625   }
 626
 627   //HFCorrelators initialization (for this event)
 628   fCorrelatorTr->SetAODEvent(aod); // set the AOD event from which you are processing
 629   fCorrelatorKc->SetAODEvent(aod);
 630   fCorrelatorK0->SetAODEvent(aod);
 631   Bool_t correlatorONTr = fCorrelatorTr->Initialize(); // initialize the pool for event mixing
 632   Bool_t correlatorONKc = fCorrelatorKc->Initialize();
 633   Bool_t correlatorONK0 = fCorrelatorK0->Initialize();
 634   if(!correlatorONTr) {AliInfo("AliHFCorrelator (tracks) didn't initialize the pool correctly or processed a bad event"); return;}
 635   if(!correlatorONKc) {AliInfo("AliHFCorrelator (charged K) didn't initialize the pool correctly or processed a bad event"); return;}
 636   if(!correlatorONK0) {AliInfo("AliHFCorrelator (K0) didn't initialize the pool correctly or processed a bad event"); return;}
 637
 638   if(fReadMC) {
 639     fCorrelatorTr->SetMCArray(mcArray); // set the TClonesArray *fmcArray for analysis on monte carlo
 640     fCorrelatorKc->SetMCArray(mcArray);
 641     fCorrelatorK0->SetMCArray(mcArray);
 642   }
 643
 644   // AOD primary vertex
 645   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 646
 647   //vtx1->Print();
 648   TString primTitle = vtx1->GetTitle();
 649   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
 650     fNentries->Fill(3);
 651   }
 652
 653   //Reset flag for tracks distributions fill
 654   fAlreadyFilled=kFALSE;
 655
 656   //***** Loop over D0 candidates *****
 657   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 658   if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 659
 660   if(fFillGlobal) { //loop on V0 and tracks for each event, to fill Pt distr. and InvMass distr.
 661
 662     TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
 663     Int_t pdgCodes[2] = {211,211};
 664     Int_t idArrayV0[v0array->GetEntriesFast()][2];
 665     for(int iV0=0; iV0<v0array->GetEntriesFast(); iV0++) {
 666       for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
 667       AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
 668       if(SelectV0(v0,vtx1,2,idArrayV0)) { //option 2 = for mass inv plots only
 669         if(fReadMC && fRecoTr && (v0->MatchToMC(310,mcArray,2,pdgCodes)<0)) continue; //310 = K0s, 311 = K0 generico!!
 670         ((TH2F*)fOutputStudy->FindObject("hK0MassInv"))->Fill(v0->MassK0Short(),v0->Pt()); //invariant mass plot
 671         ((TH1F*)fOutputStudy->FindObject("hist_Pt_K0_AllEv"))->Fill(v0->Pt()); //pT distribution (in all events), K0 case
 672       }
 673     }
 674
 675     for(Int_t itrack=0; itrack<aod->GetNumberOfTracks(); itrack++) { // loop on tacks
 676       AliAODTrack * track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrack));
 677       if(!track){
 678         AliWarning("Error in casting to AOD track. Not a standard AOD?");
 679         continue;
 680       }
 681       //rejection of tracks
 682       if(track->GetID() < 0) continue; //discard negative ID tracks
 683       if(track->Pt() < fPtThreshLow.at(0) || track->Pt() > fPtThreshUp.at(0)) continue; //discard tracks outside pt range for hadrons/K
 684       if(!fCutsTracks->IsHadronSelected(track) || !fCutsTracks->CheckHadronKinematic(track->Pt(),0.1)) continue; //0.1 = dummy (d0 value, no cut on it for me)
 685       //pT distribution (in all events), charg and hadr cases
 686       ((TH1F*)fOutputStudy->FindObject("hist_Pt_Charg_AllEv"))->Fill(track->Pt());
 687       if(fCutsTracks->CheckKaonCompatibility(track,kFALSE,0,2)) ((TH1F*)fOutputStudy->FindObject("hist_Pt_Kcharg_AllEv"))->Fill(track->Pt());
 688     }
 689
 690   } //end of loops for global plot fill
 691
 692   Int_t nSelectedloose=0,nSelectedtight=0;
 693
 694   //Fill Event Multiplicity (needed only in Reco)
 695   fMultEv = (Double_t)(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.));
 696
 697   //RecoD0 case ************************************************
 698   if(fRecoD0) {
 699
 700     for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 701       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
 702
 703       if(d->Pt()<2.) continue; //to save time and merging memory...
 704
 705       if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
 706         fNentries->Fill(2);
 707         continue; //skip the D0 from Dstar
 708       }
 709
 710       if (fCutsD0->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
 711         nSelectedloose++;
 712         nSelectedtight++;
 713         if(fSys==0){
 714           if(fCutsD0->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
 715         }
 716         Int_t ptbin=fCutsD0->PtBin(d->Pt());
 717         if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 718
 719         fIsSelectedCandidate=fCutsD0->IsSelected(d,AliRDHFCuts::kAll,aod); //D0 selected
 720         if(!fIsSelectedCandidate) continue;
 721
 722         //D0 infos
 723         Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(d->Phi());
 724                  phiD0 = fCorrelatorKc->SetCorrectPhiRange(d->Phi());  //bad usage, but returns a Double_t...
 725                  phiD0 = fCorrelatorK0->SetCorrectPhiRange(d->Phi());
 726         fCorrelatorTr->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta()); // sets the parameters of the trigger particles that are needed
 727         fCorrelatorKc->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
 728         fCorrelatorK0->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
 729         fCorrelatorTr->SetD0Properties(d,fIsSelectedCandidate); //sets special properties for D0
 730         fCorrelatorKc->SetD0Properties(d,fIsSelectedCandidate);
 731         fCorrelatorK0->SetD0Properties(d,fIsSelectedCandidate);
 732
 733         if(!fReadMC) {
 734           if (TMath::Abs(d->Eta())<fEtaForCorrel) {
 735             CalculateCorrelations(d); //correlations on real data
 736             if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv);
 737           }
 738         } else { //correlations on MC -> association of selected D0 to MCinfo with MCtruth
 739           if (TMath::Abs(d->Eta())<fEtaForCorrel) {
 740             Int_t pdgDgD0toKpi[2]={321,211};
 741             Int_t labD0 = d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not
 742             if (labD0>-1) {
 743               CalculateCorrelations(d,labD0,mcArray);
 744               if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv); //Fill multiplicity histo
 745             }
 746           }
 747         }
 748
 749         FillMassHists(d,mcArray,fCutsD0,fOutputMass);
 750       }
 751     }
 752   }
 753   //End RecoD0 case ************************************************
 754
 755   //MCKineD0 case ************************************************
 756   if(fReadMC && !fRecoD0) {
 757
 758     for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { //Loop over all the tracks of MCArray
 759       AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 760       if (!mcPart) {
 761         AliWarning("Particle not found in tree, skipping");
 762         continue;
 763       }
 764
 765       if(TMath::Abs(mcPart->GetPdgCode()) == 421){  // THIS IS A D0
 766         if (fCutsD0->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y()) ) {
 767           nSelectedloose++;
 768           nSelectedtight++;
 769
 770           //Removal of cases in which D0 decay is not in Kpi!
 771           if(mcPart->GetNDaughters()!=2) continue;
 772           AliAODMCParticle* mcDau1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
 773           AliAODMCParticle* mcDau2 = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)));
 774           if(!mcDau1 || !mcDau2) continue;
 775           Int_t pdg1 = TMath::Abs(mcDau1->GetPdgCode());
 776           Int_t pdg2 = TMath::Abs(mcDau2->GetPdgCode());
 777           if(!((pdg1 == 211 && pdg2 == 321) || (pdg2 == 211 && pdg1 == 321))) continue;
 778           if(TMath::Abs(mcDau1->Eta())>0.8||TMath::Abs(mcDau2->Eta())>0.8) continue;
 779             //Check momentum conservation (to exclude 4-prong decays with tracks outside y=1.5)
 780             Double_t p1[3]  = {mcDau1->Px(),mcDau1->Py(),mcDau1->Pz()};
 781             Double_t p2[3]  = {mcDau2->Px(),mcDau2->Py(),mcDau2->Pz()};
 782             Double_t pD0[3] = {mcPart->Px(),mcPart->Py(),mcPart->Pz()};
 783             if(TMath::Abs( (p1[0]+p2[0]-pD0[0])*(p1[0]+p2[0]-pD0[0]) + (p1[1]+p2[1]-pD0[1])*(p1[1]+p2[1]-pD0[1]) + (p1[2]+p2[2]-pD0[2])*(p1[2]+p2[2]-pD0[2]) )>0.1) continue;
 784
 785           if(fSys==0) fNentries->Fill(6);
 786           Int_t ptbin=fCutsD0->PtBin(mcPart->Pt());
 787           if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 788
 789           //D0 infos
 790           Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(mcPart->Phi());
 791                    phiD0 = fCorrelatorKc->SetCorrectPhiRange(mcPart->Phi());  //bad usage, but returns a Double_t...
 792                    phiD0 = fCorrelatorK0->SetCorrectPhiRange(mcPart->Phi());
 793           fCorrelatorTr->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta()); // sets the parameters of the trigger particles that are needed
 794           fCorrelatorKc->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
 795           fCorrelatorK0->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
 796           //fCorrelatorTr->SetD0Properties(mcPart,fIsSelectedCandidate); //needed for D* soft pions rejection, useless in MCKine
 797           //fCorrelatorKc->SetD0Properties(mcPart,fIsSelectedCandidate);
 798           //fCorrelatorK0->SetD0Properties(mcPart,fIsSelectedCandidate);
 799
 800           if (TMath::Abs(mcPart->Eta())<fEtaForCorrel) {
 801
 802             //Removal of D0 from D* feeddown! This solves also the problem of soft pions, now excluded
 803          /*   Int_t mother = mcPart->GetMother();
 804             AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
 805             if(!mcMoth) continue;
 806             if(TMath::Abs(mcMoth->GetPdgCode())==413) continue;*/
 807
 808             if (mcPart->GetPdgCode()==421) fIsSelectedCandidate = 1;
 809             else fIsSelectedCandidate = 2;
 810
 811             TString fillthis="histSgn_";
 812             if(CheckD0Origin(mcArray,mcPart)==4) fillthis+="c_";
 813             else if(CheckD0Origin(mcArray,mcPart)==5) fillthis+="b_";
 814             else continue;
 815             fillthis+=ptbin;
 816             ((TH1F*)(fOutputMass->FindObject(fillthis)))->Fill(1.864);
 817
 818             CalculateCorrelationsMCKine(mcPart,mcArray);
 819             if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv); //Fill multiplicity histo
 820           }
 821         }
 822       }
 823     }
 824   }
 825   //End MCKineD0 case ************************************************
 826
 827   if(fMixing /* && fAlreadyFilled*/) { // update the pool for Event Mixing, if: enabled,  event is ok, at least a SelD0 found! (fAlreadyFilled's role!)
 828     Bool_t updatedTr = fCorrelatorTr->PoolUpdate();
 829     Bool_t updatedKc = fCorrelatorKc->PoolUpdate();
 830     Bool_t updatedK0 = fCorrelatorK0->PoolUpdate();
 831     if(!updatedTr || !updatedKc || !updatedK0) AliInfo("Pool was not updated");
 832   }
 833
 834   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
 835   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
 836
 837   // Post the data
 838   PostData(1,fOutputMass);
 839   PostData(2,fNentries);
 840   PostData(4,fCounter);
 841   PostData(5,fOutputCorr);
 842   PostData(6,fOutputStudy);
 843
 844   return;
 845 }
 846
 847 //____________________________________________________________________________
 848 void AliAnalysisTaskSED0Correlations::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
 849   //
 850   // function used in UserExec to fill mass histograms:
 851   //
 852   if (!fIsSelectedCandidate) return;
 853
 854   if(fDebug>2)  cout<<"Candidate selected"<<endl;
 855
 856   Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
 857   Int_t ptbin = cuts->PtBin(part->Pt());
 858
 859   TString fillthis="";
 860   Int_t pdgDgD0toKpi[2]={321,211};
 861   Int_t labD0=-1;
 862   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)
 863
 864   //count candidates selected by cuts
 865   fNentries->Fill(1);
 866   //count true D0 selected by cuts
 867   if (fReadMC && labD0>=0) fNentries->Fill(2);
 868
 869   if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
 870
 871     if(fReadMC){
 872       if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==4) {
 873         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 874         Int_t pdgD0 = partD0->GetPdgCode();
 875         if (pdgD0==421){ //D0
 876           fillthis="histSgn_c_";
 877           fillthis+=ptbin;
 878           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 879           fillthis="histSgn_WeigD0Eff_c_";
 880           fillthis+=ptbin;
 881           Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 882           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 883         } else{ //it was a D0bar
 884           fillthis="histRfl_c_";
 885           fillthis+=ptbin;
 886           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 887           fillthis="histRfl_WeigD0Eff_c_";
 888           fillthis+=ptbin;
 889           Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 890           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 891         }
 892       } else if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==5) {
 893           AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 894           Int_t pdgD0 = partD0->GetPdgCode();
 895           if (pdgD0==421){ //D0
 896             fillthis="histSgn_b_";
 897             fillthis+=ptbin;
 898             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 899             fillthis="histSgn_WeigD0Eff_b_";
 900             fillthis+=ptbin;
 901             Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
 902             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 903           } else{ //it was a D0bar
 904             fillthis="histRfl_b_";
 905             fillthis+=ptbin;
 906             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 907             fillthis="histRfl_WeigD0Eff_b_";
 908             fillthis+=ptbin;
 909             Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
 910             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 911           }
 912       } else {//background
 913         fillthis="histBkg_c_";
 914         fillthis+=ptbin;
 915         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 916         fillthis="histBkg_WeigD0Eff_c_";
 917         fillthis+=ptbin;
 918         Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 919         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 920       }
 921     }else{
 922       fillthis="histMass_";
 923       fillthis+=ptbin;
 924       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 925       fillthis="histMass_WeigD0Eff_";
 926       fillthis+=ptbin;
 927       Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 928       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 929     }
 930
 931   }
 932   if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
 933
 934     if(fReadMC){
 935       if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==4) {
 936         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 937         Int_t pdgD0 = partD0->GetPdgCode();
 938         if (pdgD0==-421){ //D0
 939           fillthis="histSgn_c_";
 940           fillthis+=ptbin;
 941           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 942           fillthis="histSgn_WeigD0Eff_c_";
 943           fillthis+=ptbin;
 944           Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 945           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 946         } else{ //it was a D0bar
 947           fillthis="histRfl_c_";
 948           fillthis+=ptbin;
 949           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 950           fillthis="histRfl_WeigD0Eff_c_";
 951           fillthis+=ptbin;
 952           Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 953           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 954         }
 955       } else if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==5) {
 956           AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 957           Int_t pdgD0 = partD0->GetPdgCode();
 958           if (pdgD0==-421){ //D0
 959             fillthis="histSgn_b_";
 960             fillthis+=ptbin;
 961             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 962             fillthis="histSgn_WeigD0Eff_b_";
 963             fillthis+=ptbin;
 964             Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
 965             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 966           } else{ //it was a D0bar
 967             fillthis="histRfl_b_";
 968             fillthis+=ptbin;
 969             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 970             fillthis="histRfl_WeigD0Eff_b_";
 971             fillthis+=ptbin;
 972             Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
 973             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 974           }
 975       } else {//background
 976         fillthis="histBkg_c_";
 977         fillthis+=ptbin;
 978         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 979         fillthis="histBkg_WeigD0Eff_c_";
 980         fillthis+=ptbin;
 981         Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 982         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 983       }
 984     }else{
 985       fillthis="histMass_";
 986       fillthis+=ptbin;
 987       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 988       fillthis="histMass_WeigD0Eff_";
 989       fillthis+=ptbin;
 990       Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 991       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 992     }
 993
 994   }
 995
 996   return;
 997 }
 998
 999 //________________________________________________________________________
1000 void AliAnalysisTaskSED0Correlations::Terminate(Option_t */*option*/)
1001 {
1002   // Terminate analysis
1003   //
1004   if(fDebug > 1) printf("AnalysisTaskSED0Correlations: Terminate() \n");
1005
1006   fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1007   if (!fOutputMass) {
1008     printf("ERROR: fOutputMass not available\n");
1009     return;
1010   }
1011
1012   fNentries = dynamic_cast<TH1F*>(GetOutputData(2));
1013
1014   if(!fNentries){
1015     printf("ERROR: fNEntries not available\n");
1016     return;
1017   }
1018
1019   fCutsD0 = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(3));
1020   if(!fCutsD0){
1021     printf("ERROR: fCuts not available\n");
1022     return;
1023   }
1024
1025   fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(4));
1026   if (!fCounter) {
1027     printf("ERROR: fCounter not available\n");
1028     return;
1029   }
1030   fOutputCorr = dynamic_cast<TList*> (GetOutputData(5));
1031   if (!fOutputCorr) {
1032     printf("ERROR: fOutputCorr not available\n");
1033     return;
1034   }
1035   fOutputStudy = dynamic_cast<TList*> (GetOutputData(6));
1036   if (!fOutputStudy) {
1037     printf("ERROR: fOutputStudy not available\n");
1038     return;
1039   }
1040   fCutsTracks = dynamic_cast<AliHFAssociatedTrackCuts*>(GetOutputData(7));
1041   if(!fCutsTracks){
1042     printf("ERROR: fCutsTracks not available\n");
1043     return;
1044   }
1045
1046   return;
1047 }
1048
1049 //_________________________________________________________________________________________________
1050 Int_t AliAnalysisTaskSED0Correlations::CheckD0Origin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1051   //
1052   // checking whether the mother of the particles come from a charm or a bottom quark
1053   //
1054   printf("AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
1055
1056   Int_t pdgGranma = 0;
1057   Int_t mother = 0;
1058   mother = mcPartCandidate->GetMother();
1059   Int_t abspdgGranma =0;
1060   Bool_t isFromB=kFALSE;
1061   Bool_t isQuarkFound=kFALSE;
1062
1063   while (mother > 0){
1064     AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1065     if (mcMoth){
1066       pdgGranma = mcMoth->GetPdgCode();
1067       abspdgGranma = TMath::Abs(pdgGranma);
1068       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1069         isFromB=kTRUE;
1070       }
1071       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1072       mother = mcMoth->GetMother();
1073     }else{
1074       AliError("Failed casting the mother particle!");
1075       break;
1076     }
1077   }
1078
1079   if(isQuarkFound) {
1080     if(isFromB) return 5;
1081     else return 4;
1082   }
1083   else return 1;
1084 }
1085
1086 //________________________________________________________________________
1087 void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
1088 //
1089
1090   TString namePlot = "";
1091
1092   //These for limits in THnSparse (one per bin, same limits).
1093   //Vars: DeltaPhi, InvMass, PtTrack, Displacement, DeltaEta
1094   Int_t nBinsPhi[5] = {32,150,6,3,16};
1095   Double_t binMinPhi[5] = {-TMath::Pi()/2.,1.6,0.,0.,-1.6};  //is the minimum for all the bins
1096   Double_t binMaxPhi[5] = {3.*TMath::Pi()/2.,2.2,3.0,3.,1.6};  //is the maximum for all the bins
1097
1098   //Vars: DeltaPhi, InvMass, DeltaEta
1099   Int_t nBinsMix[4] = {32,150,16,6};
1100   Double_t binMinMix[4] = {-TMath::Pi()/2.,1.6,-1.6,0.};  //is the minimum for all the bins
1101   Double_t binMaxMix[4] = {3.*TMath::Pi()/2.,2.2,1.6,3.};  //is the maximum for all the bins
1102
1103   for(Int_t i=0;i<fNPtBinsCorr;i++){
1104
1105     if(!fMixing) {
1106       //THnSparse plots: correlations for various invariant mass (MC and data)
1107       namePlot="hPhi_K0_Bin";
1108       namePlot+=i;
1109
1110       THnSparseF *hPhiK = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1111       hPhiK->Sumw2();
1112       fOutputCorr->Add(hPhiK);
1113
1114       namePlot="hPhi_Kcharg_Bin";
1115       namePlot+=i;
1116
1117       THnSparseF *hPhiH = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1118       hPhiH->Sumw2();
1119       fOutputCorr->Add(hPhiH);
1120
1121       namePlot="hPhi_Charg_Bin";
1122       namePlot+=i;
1123
1124       THnSparseF *hPhiC = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1125       hPhiC->Sumw2();
1126       fOutputCorr->Add(hPhiC);
1127
1128       //histos for c/b origin for D0 (MC only)
1129       if (fReadMC) {
1130
1131         //generic origin for tracks
1132         namePlot="hPhi_K0_From_c_Bin";
1133         namePlot+=i;
1134
1135         THnSparseF *hPhiK_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1136         hPhiK_c->Sumw2();
1137         fOutputCorr->Add(hPhiK_c);
1138
1139         namePlot="hPhi_Kcharg_From_c_Bin";
1140         namePlot+=i;
1141
1142         THnSparseF *hPhiH_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1143         hPhiH_c->Sumw2();
1144         fOutputCorr->Add(hPhiH_c);
1145
1146         namePlot="hPhi_Charg_From_c_Bin";
1147         namePlot+=i;
1148
1149         THnSparseF *hPhiC_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1150         hPhiC_c->Sumw2();
1151         fOutputCorr->Add(hPhiC_c);
1152
1153         namePlot="hPhi_K0_From_b_Bin";
1154         namePlot+=i;
1155
1156         THnSparseF *hPhiK_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1157         hPhiK_b->Sumw2();
1158         fOutputCorr->Add(hPhiK_b);
1159
1160         namePlot="hPhi_Kcharg_From_b_Bin";
1161         namePlot+=i;
1162
1163         THnSparseF *hPhiH_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1164         hPhiH_b->Sumw2();
1165         fOutputCorr->Add(hPhiH_b);
1166
1167         namePlot="hPhi_Charg_From_b_Bin";
1168         namePlot+=i;
1169
1170         THnSparseF *hPhiC_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1171         hPhiC_b->Sumw2();
1172         fOutputCorr->Add(hPhiC_b);
1173
1174         //HF-only tracks (c for c->D0, b for b->D0)
1175         namePlot="hPhi_K0_HF_From_c_Bin";
1176         namePlot+=i;
1177
1178         THnSparseF *hPhiK_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1179         hPhiK_HF_c->Sumw2();
1180         fOutputCorr->Add(hPhiK_HF_c);
1181
1182         namePlot="hPhi_Kcharg_HF_From_c_Bin";
1183         namePlot+=i;
1184
1185         THnSparseF *hPhiH_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1186         hPhiH_HF_c->Sumw2();
1187         fOutputCorr->Add(hPhiH_HF_c);
1188
1189         namePlot="hPhi_Charg_HF_From_c_Bin";
1190         namePlot+=i;
1191
1192         THnSparseF *hPhiC_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1193         hPhiC_HF_c->Sumw2();
1194         fOutputCorr->Add(hPhiC_HF_c);
1195
1196         namePlot="hPhi_K0_HF_From_b_Bin";
1197         namePlot+=i;
1198
1199         THnSparseF *hPhiK_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1200         hPhiK_HF_b->Sumw2();
1201         fOutputCorr->Add(hPhiK_HF_b);
1202
1203         namePlot="hPhi_Kcharg_HF_From_b_Bin";
1204         namePlot+=i;
1205
1206         THnSparseF *hPhiH_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1207         hPhiH_HF_b->Sumw2();
1208         fOutputCorr->Add(hPhiH_HF_b);
1209
1210         namePlot="hPhi_Charg_HF_From_b_Bin";
1211         namePlot+=i;
1212
1213         THnSparseF *hPhiC_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1214         hPhiC_HF_b->Sumw2();
1215         fOutputCorr->Add(hPhiC_HF_b);
1216
1217         namePlot="hPhi_K0_NonHF_Bin";
1218         namePlot+=i;
1219
1220         THnSparseF *hPhiK_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1221         hPhiK_Non->Sumw2();
1222         fOutputCorr->Add(hPhiK_Non);
1223
1224         namePlot="hPhi_Kcharg_NonHF_Bin";
1225         namePlot+=i;
1226
1227         THnSparseF *hPhiH_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1228         hPhiH_Non->Sumw2();
1229         fOutputCorr->Add(hPhiH_Non);
1230
1231         namePlot="hPhi_Charg_NonHF_Bin";
1232         namePlot+=i;
1233
1234         THnSparseF *hPhiC_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1235         hPhiC_Non->Sumw2();
1236         fOutputCorr->Add(hPhiC_Non);
1237       }
1238
1239       //leading hadron correlations
1240       namePlot="hPhi_Lead_Bin";
1241       namePlot+=i;
1242
1243       THnSparseF *hCorrLead = new THnSparseF(namePlot.Data(), "Leading particle correlations; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1244       hCorrLead->Sumw2();
1245       fOutputCorr->Add(hCorrLead);
1246
1247       if (fReadMC) {
1248         namePlot="hPhi_Lead_From_c_Bin";
1249         namePlot+=i;
1250
1251         THnSparseF *hCorrLead_c = new THnSparseF(namePlot.Data(), "Leading particle correlations - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1252         hCorrLead_c->Sumw2();
1253         fOutputCorr->Add(hCorrLead_c);
1254
1255         namePlot="hPhi_Lead_From_b_Bin";
1256         namePlot+=i;
1257
1258         THnSparseF *hCorrLead_b = new THnSparseF(namePlot.Data(), "Leading particle correlations - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1259         hCorrLead_b->Sumw2();
1260         fOutputCorr->Add(hCorrLead_b);
1261
1262         namePlot="hPhi_Lead_HF_From_c_Bin";
1263         namePlot+=i;
1264
1265         THnSparseF *hCorrLead_HF_c = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1266         hCorrLead_HF_c->Sumw2();
1267         fOutputCorr->Add(hCorrLead_HF_c);
1268
1269         namePlot="hPhi_Lead_HF_From_b_Bin";
1270         namePlot+=i;
1271
1272         THnSparseF *hCorrLead_HF_b = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1273         hCorrLead_HF_b->Sumw2();
1274         fOutputCorr->Add(hCorrLead_HF_b);
1275
1276         namePlot="hPhi_Lead_NonHF_Bin";
1277         namePlot+=i;
1278
1279         THnSparseF *hCorrLead_Non = new THnSparseF(namePlot.Data(), "Leading particle correlations - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1280         hCorrLead_Non->Sumw2();
1281         fOutputCorr->Add(hCorrLead_Non);
1282       }
1283
1284       //pT weighted correlations
1285       namePlot="hPhi_Weig_Bin";
1286       namePlot+=i;
1287
1288       THnSparseF *hCorrWeig = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted); #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1289       fOutputCorr->Add(hCorrWeig);
1290
1291       if (fReadMC) {
1292         namePlot="hPhi_Weig_From_c_Bin";
1293         namePlot+=i;
1294
1295         THnSparseF *hCorrWeig_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1296         fOutputCorr->Add(hCorrWeig_c);
1297
1298         namePlot="hPhi_Weig_From_b_Bin";
1299         namePlot+=i;
1300
1301         THnSparseF *hCorrWeig_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1302         fOutputCorr->Add(hCorrWeig_b);
1303
1304         namePlot="hPhi_Weig_HF_From_c_Bin";
1305         namePlot+=i;
1306
1307         THnSparseF *hCorrWeig_HF_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1308         fOutputCorr->Add(hCorrWeig_HF_c);
1309
1310         namePlot="hPhi_Weig_HF_From_b_Bin";
1311         namePlot+=i;
1312
1313         THnSparseF *hCorrWeig_HF_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1314         fOutputCorr->Add(hCorrWeig_HF_b);
1315
1316         namePlot="hPhi_Weig_NonHF_Bin";
1317         namePlot+=i;
1318
1319         THnSparseF *hCorrWeig_Non = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1320         fOutputCorr->Add(hCorrWeig_Non);
1321       }
1322
1323     //pT distribution histos
1324     namePlot = "hist_Pt_Charg_Bin"; namePlot+=i;
1325     TH1F *hPtC = new TH1F(namePlot.Data(), "Charged track pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1326     hPtC->SetMinimum(0);
1327     fOutputStudy->Add(hPtC);
1328
1329     namePlot = "hist_Pt_Kcharg_Bin"; namePlot+=i;
1330     TH1F *hPtH = new TH1F(namePlot.Data(), "Hadrons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1331     hPtH->SetMinimum(0);
1332     fOutputStudy->Add(hPtH);
1333
1334     namePlot = "hist_Pt_K0_Bin"; namePlot+=i;
1335     TH1F *hPtK = new TH1F(namePlot.Data(), "Kaons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1336     hPtK->SetMinimum(0);
1337     fOutputStudy->Add(hPtK);
1338
1339     //D* feeddown pions rejection histos
1340     namePlot = "hDstarPions_Bin"; namePlot+=i;
1341     TH2F *hDstarPions = new TH2F(namePlot.Data(), "Tracks rejected for D* inv.mass cut; # Tracks",2,0.,2.,150,1.6,2.2);
1342     hDstarPions->GetXaxis()->SetBinLabel(1,"Not rejected");
1343     hDstarPions->GetXaxis()->SetBinLabel(2,"Rejected");
1344     hDstarPions->SetMinimum(0);
1345     fOutputStudy->Add(hDstarPions);
1346
1347     //Events multiplicity
1348     Double_t yAxisMult[13] = {0, 4, 8, 12, 16, 20, 28, 36, 44, 100};
1349     namePlot = "hMultiplEvt_Bin"; namePlot+=i;
1350     TH1F *hMultEv = new TH1F(namePlot.Data(), "Event multiplicity",9,yAxisMult);
1351     hMultEv->SetMinimum(0);
1352     fOutputStudy->Add(hMultEv);
1353
1354     }
1355
1356     if(fMixing) {
1357       //THnSparse plots for event mixing!
1358       namePlot="hPhi_K0_Bin";
1359       namePlot+=i;namePlot+="_EvMix";
1360
1361       THnSparseF *hPhiK_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1362       hPhiK_EvMix->Sumw2();
1363       fOutputCorr->Add(hPhiK_EvMix);
1364
1365       namePlot="hPhi_Kcharg_Bin";
1366       namePlot+=i;namePlot+="_EvMix";
1367
1368       THnSparseF *hPhiH_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1369       hPhiH_EvMix->Sumw2();
1370       fOutputCorr->Add(hPhiH_EvMix);
1371
1372       namePlot="hPhi_Charg_Bin";
1373       namePlot+=i;namePlot+="_EvMix";
1374
1375       THnSparseF *hPhiC_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1376       hPhiC_EvMix->Sumw2();
1377       fOutputCorr->Add(hPhiC_EvMix);
1378     }
1379   }
1380
1381   //out of bin loop
1382   if(!fMixing) {
1383     TH1F *hCountC = new TH1F("hist_Count_Charg", "Charged track counter; # Tracks",100,0.,100.);
1384     hCountC->SetMinimum(0);
1385     fOutputStudy->Add(hCountC);
1386
1387     TH1F *hCountH = new TH1F("hist_Count_Kcharg", "Hadrons counter; # Tracks",20,0.,20.);
1388     hCountH->SetMinimum(0);
1389     fOutputStudy->Add(hCountH);
1390
1391     TH1F *hCountK = new TH1F("hist_Count_K0", "Kaons counter; # Tracks",20,0.,20.);
1392     hCountK->SetMinimum(0);
1393     fOutputStudy->Add(hCountK);
1394   }
1395
1396   if (fReadMC) {
1397     TH1D *hEventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
1398     fOutputStudy->Add(hEventTypeMC);
1399   }
1400
1401   if (fFillGlobal) { //all-events plots
1402     //pt distributions
1403     TH1F *hPtCAll = new TH1F("hist_Pt_Charg_AllEv", "Charged track pT (All); p_{T} (GeV/c)",240,0.,12.);
1404     hPtCAll->SetMinimum(0);
1405     fOutputStudy->Add(hPtCAll);
1406
1407     TH1F *hPtHAll = new TH1F("hist_Pt_Kcharg_AllEv", "Hadrons pT (All); p_{T} (GeV/c)",240,0.,12.);
1408     hPtHAll->SetMinimum(0);
1409     fOutputStudy->Add(hPtHAll);
1410
1411     TH1F *hPtKAll = new TH1F("hist_Pt_K0_AllEv", "Kaons pT (All); p_{T} (GeV/c)",240,0.,12.);
1412     hPtKAll->SetMinimum(0);
1413     fOutputStudy->Add(hPtKAll);
1414
1415     //K0 Invariant Mass plots
1416     TH2F *hK0MassInv = new TH2F("hK0MassInv", "K0 invariant mass; Invariant mass (MeV/c^{2}); pT (GeV/c)",200,0.4,0.6,100,0.,10.);
1417     hK0MassInv->SetMinimum(0);
1418     fOutputStudy->Add(hK0MassInv);
1419   }
1420
1421   if(!fMixing) {
1422     //phi distributions
1423     TH1F *hPhiDistCAll = new TH1F("hist_PhiDistr_Charg", "Charged track phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1424     hPhiDistCAll->SetMinimum(0);
1425     fOutputStudy->Add(hPhiDistCAll);
1426
1427     TH1F *hPhiDistHAll = new TH1F("hist_PhiDistr_Kcharg", "Hadrons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1428     hPhiDistHAll->SetMinimum(0);
1429     fOutputStudy->Add(hPhiDistHAll);
1430
1431     TH1F *hPhiDistKAll = new TH1F("hist_PhiDistr_K0", "Kaons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1432     hPhiDistKAll->SetMinimum(0);
1433     fOutputStudy->Add(hPhiDistKAll);
1434
1435     TH1F *hPhiDistDAll = new TH1F("hist_PhiDistr_D0", "D^{0} phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1436     hPhiDistDAll->SetMinimum(0);
1437     fOutputStudy->Add(hPhiDistDAll);
1438   }
1439
1440   //for MC analysis only
1441   for(Int_t i=0;i<fNPtBinsCorr;i++) {
1442
1443     if (fReadMC && !fMixing) {
1444
1445       //displacement histos
1446       namePlot="histDispl_K0_Bin"; namePlot+=i;
1447       TH1F *hDisplK = new TH1F(namePlot.Data(), "Kaons Displacement; DCA",150,0.,0.15);
1448       hDisplK->SetMinimum(0);
1449       fOutputStudy->Add(hDisplK);
1450
1451       namePlot="histDispl_K0_HF_Bin";  namePlot+=i;
1452       TH1F *hDisplK_HF = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only); DCA",150,0.,0.15);
1453       hDisplK_HF->SetMinimum(0);
1454       fOutputStudy->Add(hDisplK_HF);
1455
1456       namePlot="histDispl_Kcharg_Bin"; namePlot+=i;
1457       TH1F *hDisplHadr = new TH1F(namePlot.Data(), "Hadrons Displacement; DCA",150,0.,0.15);
1458       hDisplHadr->SetMinimum(0);
1459       fOutputStudy->Add(hDisplHadr);
1460
1461       namePlot="histDispl_Kcharg_HF_Bin";  namePlot+=i;
1462       TH1F *hDisplHadr_HF = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only); DCA",150,0.,0.15);
1463       hDisplHadr_HF->SetMinimum(0);
1464       fOutputStudy->Add(hDisplHadr_HF);
1465
1466       namePlot="histDispl_Charg_Bin"; namePlot+=i;
1467       TH1F *hDisplCharg = new TH1F(namePlot.Data(), "Charged tracks Displacement; DCA",150,0.,0.15);
1468       hDisplCharg->SetMinimum(0);
1469       fOutputStudy->Add(hDisplCharg);
1470
1471       namePlot="histDispl_Charg_HF_Bin";  namePlot+=i;
1472       TH1F *hDisplCharg_HF = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only); DCA",150,0.,0.15);
1473       hDisplCharg_HF->SetMinimum(0);
1474       fOutputStudy->Add(hDisplCharg_HF);
1475
1476       namePlot="histDispl_K0_From_c_Bin"; namePlot+=i;
1477       TH1F *hDisplK_c = new TH1F(namePlot.Data(), "Kaons Displacement - c origin; DCA",150,0.,0.15);
1478       hDisplK_c->SetMinimum(0);
1479       fOutputStudy->Add(hDisplK_c);
1480
1481       namePlot="histDispl_K0_HF_From_c_Bin";  namePlot+=i;
1482       TH1F *hDisplK_HF_c = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1483       hDisplK_HF_c->SetMinimum(0);
1484       fOutputStudy->Add(hDisplK_HF_c);
1485
1486       namePlot="histDispl_Kcharg_From_c_Bin"; namePlot+=i;
1487       TH1F *hDisplHadr_c = new TH1F(namePlot.Data(), "Hadrons Displacement - c origin; DCA",150,0.,0.15);
1488       hDisplHadr_c->SetMinimum(0);
1489       fOutputStudy->Add(hDisplHadr_c);
1490
1491       namePlot="histDispl_Kcharg_HF_From_c_Bin";  namePlot+=i;
1492       TH1F *hDisplHadr_HF_c = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1493       hDisplHadr_HF_c->SetMinimum(0);
1494       fOutputStudy->Add(hDisplHadr_HF_c);
1495
1496       namePlot="histDispl_Charg_From_c_Bin"; namePlot+=i;
1497       TH1F *hDisplCharg_c = new TH1F(namePlot.Data(), "Charged tracks Displacement - c origin; DCA",150,0.,0.15);
1498       hDisplCharg_c->Sumw2();
1499       hDisplCharg_c->SetMinimum(0);
1500       fOutputStudy->Add(hDisplCharg_c);
1501
1502       namePlot="histDispl_Charg_HF_From_c_Bin";  namePlot+=i;
1503       TH1F *hDisplCharg_HF_c = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1504       hDisplCharg_HF_c->SetMinimum(0);
1505       fOutputStudy->Add(hDisplCharg_HF_c);
1506
1507       namePlot="histDispl_K0_From_b_Bin"; namePlot+=i;
1508       TH1F *hDisplK_b = new TH1F(namePlot.Data(), "Kaons Displacement - b origin; DCA",150,0.,0.15);
1509       hDisplK_b->SetMinimum(0);
1510       fOutputStudy->Add(hDisplK_b);
1511
1512       namePlot="histDispl_K0_HF_From_b_Bin";  namePlot+=i;
1513       TH1F *hDisplK_HF_b = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1514       hDisplK_HF_b->SetMinimum(0);
1515       fOutputStudy->Add(hDisplK_HF_b);
1516
1517       namePlot="histDispl_Kcharg_From_b_Bin"; namePlot+=i;
1518       TH1F *hDisplHadr_b = new TH1F(namePlot.Data(), "Hadrons Displacement - b origin; DCA",150,0.,0.15);
1519       hDisplHadr_b->SetMinimum(0);
1520       fOutputStudy->Add(hDisplHadr_b);
1521
1522       namePlot="histDispl_Kcharg_HF_From_b_Bin";  namePlot+=i;
1523       TH1F *hDisplHadr_HF_b = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1524       hDisplHadr_HF_b->SetMinimum(0);
1525       fOutputStudy->Add(hDisplHadr_HF_b);
1526
1527       namePlot="histDispl_Charg_From_b_Bin"; namePlot+=i;
1528       TH1F *hDisplCharg_b = new TH1F(namePlot.Data(), "Charged tracks Displacement - b origin; DCA",150,0.,0.15);
1529       hDisplCharg_b->SetMinimum(0);
1530       fOutputStudy->Add(hDisplCharg_b);
1531
1532       namePlot="histDispl_Charg_HF_From_b_Bin";  namePlot+=i;
1533       TH1F *hDisplCharg_HF_b = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1534       hDisplCharg_HF_b->SetMinimum(0);
1535       fOutputStudy->Add(hDisplCharg_HF_b);
1536
1537       //origin of tracks histos
1538       namePlot="histOrig_Charg_Bin";  namePlot+=i;
1539       TH1F *hOrigin_Charm = new TH1F(namePlot.Data(), "Origin of charged tracks",9,0.,9.);
1540       hOrigin_Charm->SetMinimum(0);
1541       hOrigin_Charm->GetXaxis()->SetBinLabel(1,"Not HF");
1542       hOrigin_Charm->GetXaxis()->SetBinLabel(2,"D->#");
1543       hOrigin_Charm->GetXaxis()->SetBinLabel(3,"D->X->#");
1544       hOrigin_Charm->GetXaxis()->SetBinLabel(4,"c hadr.");
1545       hOrigin_Charm->GetXaxis()->SetBinLabel(5,"B->#");
1546       hOrigin_Charm->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1547       hOrigin_Charm->GetXaxis()->SetBinLabel(7,"B->D->#");
1548       hOrigin_Charm->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1549       hOrigin_Charm->GetXaxis()->SetBinLabel(9,"b hadr.");
1550       fOutputStudy->Add(hOrigin_Charm);
1551
1552       namePlot="histOrig_Kcharg_Bin";  namePlot+=i;
1553       TH1F *hOrigin_Kcharg = new TH1F(namePlot.Data(), "Origin of hadrons",9,0.,9.);
1554       hOrigin_Kcharg->SetMinimum(0);
1555       hOrigin_Kcharg->GetXaxis()->SetBinLabel(1,"Not HF");
1556       hOrigin_Kcharg->GetXaxis()->SetBinLabel(2,"D->#");
1557       hOrigin_Kcharg->GetXaxis()->SetBinLabel(3,"D->X->#");
1558       hOrigin_Kcharg->GetXaxis()->SetBinLabel(4,"c hadr.");
1559       hOrigin_Kcharg->GetXaxis()->SetBinLabel(5,"B->#");
1560       hOrigin_Kcharg->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1561       hOrigin_Kcharg->GetXaxis()->SetBinLabel(7,"B->D->#");
1562       hOrigin_Kcharg->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1563       hOrigin_Kcharg->GetXaxis()->SetBinLabel(9,"b hadr.");
1564       fOutputStudy->Add(hOrigin_Kcharg);
1565
1566       namePlot="histOrig_K0_Bin";  namePlot+=i;
1567       TH1F *hOrigin_K = new TH1F(namePlot.Data(), "Origin of kaons",9,0.,9.);
1568       hOrigin_K->SetMinimum(0);
1569       hOrigin_K->GetXaxis()->SetBinLabel(1,"Not HF");
1570       hOrigin_K->GetXaxis()->SetBinLabel(2,"D->#");
1571       hOrigin_K->GetXaxis()->SetBinLabel(3,"D->X->#");
1572       hOrigin_K->GetXaxis()->SetBinLabel(4,"c hadr.");
1573       hOrigin_K->GetXaxis()->SetBinLabel(5,"B->#");
1574       hOrigin_K->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1575       hOrigin_K->GetXaxis()->SetBinLabel(7,"B->D->#");
1576       hOrigin_K->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1577       hOrigin_K->GetXaxis()->SetBinLabel(9,"b hadr.");
1578       fOutputStudy->Add(hOrigin_K);
1579     }
1580
1581     if (fReadMC) {
1582       //origin of D0 histos
1583       namePlot="histOrig_D0_Bin";  namePlot+=i;
1584       TH1F *hOrigin_D0 = new TH1F(namePlot.Data(), "Origin of D0",2,0.,2.);
1585       hOrigin_D0->SetMinimum(0);
1586       hOrigin_D0->GetXaxis()->SetBinLabel(1,"From c");
1587       hOrigin_D0->GetXaxis()->SetBinLabel(2,"From b");
1588       fOutputStudy->Add(hOrigin_D0);
1589
1590       //primary tracks (Kine & Reco)
1591       namePlot="hPhysPrim_Bin";  namePlot+=i;
1592       TH1F *hPhysPrim = new TH1F(namePlot.Data(), "Origin of hadrons",2,0.,2.);
1593       hPhysPrim->SetMinimum(0);
1594       hPhysPrim->GetXaxis()->SetBinLabel(1,"OK");
1595       hPhysPrim->GetXaxis()->SetBinLabel(2,"NO");
1596       fOutputStudy->Add(hPhysPrim);
1597     }
1598   }
1599 }
1600
1601 //________________________________________________________________________
1602 void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Prong* d, Int_t labD0, TClonesArray* mcArray) {
1603 //
1604 // Method for correlations D0-hadrons study
1605 //
1606   Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
1607   Double_t mD0, mD0bar;
1608   Int_t origD0 = 0, PDGD0 = 0, ptbin = 0;
1609   d->InvMassD0(mD0,mD0bar);
1610   Double_t mInv[2] = {mD0, mD0bar};
1611   ptbin = PtBinCorr(d->Pt());
1612
1613   if(ptbin < 0) return;
1614
1615   //Fill of D0 phi distribution
1616   if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi());
1617
1618   //Origin of D0
1619   TString orig="";
1620   if(fReadMC) {
1621     origD0=CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0));
1622     PDGD0 = ((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode();
1623     switch (CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0))) {
1624       case 4:
1625         orig = "_From_c";
1626         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1627         break;
1628       case 5:
1629         orig = "_From_b";
1630         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1631         break;
1632       default:
1633         return;
1634     }
1635   }
1636
1637   Double_t highPt = 0; Double_t lead[4] = {0,0,0,1};  //infos for leading particle (pt,deltaphi)
1638
1639   //loop over the tracks in the pool
1640   Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1641   Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1642   Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1643
1644   Int_t NofEventsinPool = 1;
1645   if(fMixing) {
1646     NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
1647     if(!execPoolTr) {
1648       AliInfo("Mixed event analysis: track pool is not ready");
1649       NofEventsinPool = 0;
1650     }
1651   }
1652
1653   //Charged tracks
1654   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1655     Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
1656     if(!analyzetracksTr) {
1657       AliInfo("AliHFCorrelator::Cannot process the track array");
1658       continue;
1659     }
1660
1661     for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
1662
1663       Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
1664       if(!runcorrelation) continue;
1665
1666       AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
1667
1668       if(!fMixing) {
1669         if(fSoftPiCut && !track->CheckSoftPi()) { //removal of soft pions
1670           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1671           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1672           continue;
1673         } else { //not a soft pion
1674           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1675           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1676         }
1677         Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1678         if(track->GetID() == idDaughs[0] || track->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1679       }
1680       if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1681
1682       if(fReadMC) {
1683         AliAODMCParticle* trkKine = (AliAODMCParticle*)mcArray->At(track->GetLabel());
1684         if (!trkKine) continue;
1685         if (!trkKine->IsPhysicalPrimary()) {
1686           ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(1.);
1687           continue; //reject the Reco track if correspondent Kine track is not primary
1688         } else ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(0.);
1689       }
1690
1691       Double_t effTr = track->GetWeight(); //extract track efficiency
1692       Double_t effD0 = 1.;
1693       if(fReadMC) {
1694         if(origD0==4) effD0 = fCutsTracks->GetTrigWeight(d->Pt(),fMultEv);
1695         if(origD0==5) effD0 = fCutsTracks->GetTrigWeightB(d->Pt(),fMultEv);
1696       } else effD0 = fCutsTracks->GetTrigWeight(d->Pt(),fMultEv);
1697       Double_t eff = effTr*effD0;
1698
1699       FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack,1./eff); //fills for charged tracks
1700
1701       if(!fMixing) N_Charg++;
1702
1703       //retrieving leading info...
1704       if(track->Pt() > highPt) {
1705         if(fReadMC && track->GetLabel()<1) continue;
1706         if(fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) continue;
1707         lead[0] = fCorrelatorTr->GetDeltaPhi();
1708         lead[1] = fCorrelatorTr->GetDeltaEta();
1709         lead[2] = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
1710         if(fReadMC) {
1711           if(origD0==4) lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeight(d->Pt(),fMultEv)); //weight is 1./efficiency
1712           if(origD0==5) lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeightB(d->Pt(),fMultEv)); //weight is 1./efficiency
1713         } else lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeight(d->Pt(),fMultEv));
1714         highPt = track->Pt();
1715       }
1716
1717     } // end of tracks loop
1718   } //end of event loop for fCorrelatorTr
1719
1720  if(fKaonCorr) { //loops for Kcharg and K0
1721
1722   if(fMixing) {
1723     NofEventsinPool = fCorrelatorKc->GetNofEventsInPool();
1724     if(!execPoolKc) {
1725       AliInfo("Mixed event analysis: K+/- pool is not ready");
1726       NofEventsinPool = 0;
1727     }
1728   }
1729
1730   //Charged Kaons loop
1731   for (Int_t jMix = 0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1732     Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
1733     if(!analyzetracksKc) {
1734       AliInfo("AliHFCorrelator::Cannot process the K+/- array");
1735       continue;
1736     }
1737
1738     for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
1739
1740       Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
1741       if(!runcorrelation) continue;
1742
1743       AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
1744
1745       if(!fMixing) {
1746         if(fSoftPiCut && !kCharg->CheckSoftPi()) { //removal of soft pions
1747           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1748           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1749           continue;
1750         } else {
1751           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1752           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1753         }
1754         Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1755         if(kCharg->GetID() == idDaughs[0] || kCharg->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1756       }
1757       if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1758
1759       FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
1760
1761       if(!fMixing) N_KCharg++;
1762
1763     } // end of charged kaons loop
1764   } //end of event loop for fCorrelatorKc
1765
1766   if(fMixing) {
1767     NofEventsinPool = fCorrelatorK0->GetNofEventsInPool();
1768     if(!execPoolK0) {
1769       AliInfo("Mixed event analysis: K0 pool is not ready");
1770       NofEventsinPool = 0;
1771     }
1772   }
1773
1774   //K0 loop
1775   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1776     Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
1777     if(!analyzetracksK0) {
1778       AliInfo("AliHFCorrelator::Cannot process the K0 array");
1779       continue;
1780     }
1781
1782     for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
1783
1784       Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
1785       if(!runcorrelation) continue;
1786
1787       AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
1788
1789       if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1790
1791       FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
1792
1793       if(!fMixing) N_Kaons++;
1794
1795     } // end of charged kaons loop
1796   } //end of event loop for fCorrelatorK0
1797
1798  } //end of 'if(fKaonCorr)'
1799
1800   Double_t fillSpLeadD0[4] = {lead[0],mD0,lead[1],0.4}; //dummy value for threshold of leading!
1801   Double_t fillSpLeadD0bar[4] = {lead[0],mD0bar,lead[1],0.4};
1802
1803   //leading track correlations fill
1804   if(!fMixing) {
1805     if(fReadMC) {
1806       if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0
1807         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]); //c and b D0
1808         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]); //c or b D0
1809         if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);
1810         if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);
1811         if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]);  //non HF
1812       }
1813       if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421 && fIsSelectedCandidate>1) { //D0bar
1814         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1815         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]); //c or b D0
1816         if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1817         if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1818         if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);  //non HF
1819       }
1820     } else {
1821         if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]);
1822         if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1823     }
1824
1825     //Fill of count histograms
1826     if (!fAlreadyFilled) {
1827       ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
1828       ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
1829       ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
1830     }
1831   }
1832
1833   fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
1834
1835 }
1836
1837 //________________________________________________________________________
1838 void AliAnalysisTaskSED0Correlations::CalculateCorrelationsMCKine(AliAODMCParticle* d, TClonesArray* mcArray) {
1839 //
1840 // Method for correlations D0-hadrons study
1841 //
1842   Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
1843   Double_t mD0 = 1.864, mD0bar = 1.864;
1844   Double_t mInv[2] = {mD0, mD0bar};
1845   Int_t origD0 = 0, PDGD0 = 0;
1846   Int_t ptbin = PtBinCorr(d->Pt());
1847
1848   if(ptbin < 0) return;
1849
1850   //Fill of D0 phi distribution
1851   if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi());
1852
1853   //Origin of D0
1854   TString orig="";
1855   origD0=CheckD0Origin(mcArray,d);
1856   PDGD0 = d->GetPdgCode();
1857   switch (CheckD0Origin(mcArray,d)) {
1858     case 4:
1859       orig = "_From_c";
1860       ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1861       break;
1862     case 5:
1863       orig = "_From_b";
1864       ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1865       break;
1866     default:
1867       return;
1868   }
1869
1870   Double_t highPt = 0; Double_t lead[3] = {0,0,0};  //infos for leading particle (pt,deltaphi)
1871
1872   //loop over the tracks in the pool
1873   Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1874   Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1875   Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1876
1877   Int_t NofEventsinPool = 1;
1878   if(fMixing) {
1879     NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
1880     if(!execPoolTr) {
1881       AliInfo("Mixed event analysis: track pool is not ready");
1882       NofEventsinPool = 0;
1883     }
1884   }
1885
1886   //Charged tracks
1887   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1888
1889     Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
1890     if(!analyzetracksTr) {
1891       AliInfo("AliHFCorrelator::Cannot process the track array");
1892       continue;
1893     }
1894
1895     for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
1896
1897       Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
1898       if(!runcorrelation) continue;
1899
1900       AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
1901       if(track->GetLabel()<0) continue;
1902       if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1903       if(track->Pt() < 0.3 || TMath::Abs(track->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1904       if(!fMixing) N_Charg++;
1905
1906       AliAODMCParticle *trkMC = (AliAODMCParticle*)mcArray->At(track->GetLabel());
1907       if(!trkMC) continue;
1908
1909       if (!trkMC->IsPhysicalPrimary()) {  //reject material budget, or other fake tracks
1910         ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(1.);
1911         continue;
1912       } else ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(0.);
1913
1914       if (IsDDaughter(d,trkMC,mcArray)) continue;
1915       if (fSoftPiCut && IsSoftPion_MCKine(d,trkMC,mcArray)) continue; //remove soft pions (if requestes, e.g. for templates)
1916
1917       FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack); //fills for charged tracks
1918
1919       //retrieving leading info...
1920       if(track->Pt() > highPt) {
1921         lead[0] = fCorrelatorTr->GetDeltaPhi();
1922         lead[1] = fCorrelatorTr->GetDeltaEta();
1923         lead[2] = fReadMC ? CheckTrackOrigin(mcArray,trkMC) : 0;
1924         highPt = track->Pt();
1925       }
1926
1927     } // end of tracks loop
1928   } //end of event loop for fCorrelatorTr
1929
1930  if(fKaonCorr) { //loops for Kcharg and K0
1931
1932   if(fMixing) {
1933     NofEventsinPool = fCorrelatorKc->GetNofEventsInPool();
1934     if(!execPoolKc) {
1935       AliInfo("Mixed event analysis: K+/- pool is not ready");
1936       NofEventsinPool = 0;
1937     }
1938   }
1939
1940   //Charged Kaons loop
1941   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1942     Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
1943     if(!analyzetracksKc) {
1944       AliInfo("AliHFCorrelator::Cannot process the K+/- array");
1945       continue;
1946     }
1947
1948     for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
1949
1950       Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
1951       if(!runcorrelation) continue;
1952
1953       AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
1954       if(kCharg->GetLabel()<1) continue;
1955       if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1956       if(TMath::Abs(kCharg->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1957       if(!fMixing) N_KCharg++;
1958
1959       AliAODMCParticle *kChargMC = (AliAODMCParticle*)mcArray->At(kCharg->GetLabel());
1960       if(!kChargMC) continue;
1961
1962       if (IsDDaughter(d,kChargMC,mcArray)) continue;
1963       FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
1964
1965     } // end of charged kaons loop
1966   } //end of event loop for fCorrelatorKc
1967
1968   if(fMixing) {
1969     NofEventsinPool = fCorrelatorK0->GetNofEventsInPool();
1970     if(!execPoolK0) {
1971       AliInfo("Mixed event analysis: K0 pool is not ready");
1972       NofEventsinPool = 0;
1973     }
1974   }
1975
1976   //K0 loop
1977   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1978     Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
1979     if(!analyzetracksK0) {
1980       AliInfo("AliHFCorrelator::Cannot process the K0 array");
1981       continue;
1982     }
1983
1984     for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
1985
1986       Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
1987       if(!runcorrelation) continue;
1988
1989       AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
1990       if(k0->GetLabel()<1) continue;
1991       if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1992       if(TMath::Abs(k0->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1993
1994       AliAODMCParticle *k0MC = (AliAODMCParticle*)mcArray->At(k0->GetLabel());
1995       if(!k0MC) continue;
1996
1997       if (IsDDaughter(d,k0MC,mcArray)) continue;
1998       FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
1999
2000       if(!fMixing) N_Kaons++;
2001
2002     } // end of charged kaons loop
2003   } //end of event loop for fCorrelatorK0
2004
2005  } //end of 'if(fKaonCorr)'
2006
2007   Double_t fillSpLeadMC[4] = {lead[0],mD0,lead[1],0.4}; //mD0 = mD0bar = 1.864
2008
2009   //leading track correlations fill
2010   if(!fMixing) {
2011     if(d->GetPdgCode()==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0
2012       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC); //c and b D0
2013       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
2014       if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
2015       if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
2016       if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC);  //non HF
2017     }
2018     if(d->GetPdgCode()==-421 && fIsSelectedCandidate>1) { //D0bar
2019       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC);
2020       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
2021       if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
2022       if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
2023       if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC);  //non HF
2024     }
2025
2026     //Fill of count histograms
2027     if (!fAlreadyFilled) {
2028       ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
2029       ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
2030       ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
2031     }
2032
2033     fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
2034   }
2035
2036 }
2037
2038 //________________________________________________________________________
2039 void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, Double_t mInv[], Int_t origD0, Int_t PdgD0, AliReducedParticle* track, Int_t ptbin, Int_t type, Double_t wg) {
2040   //
2041   //fills the THnSparse for correlations, calculating the variables
2042   //
2043
2044   //Initialization of variables
2045   Double_t mD0, mD0bar, deltaphi = 0., deltaeta = 0.;
2046   mD0 = mInv[0];
2047   mD0bar = mInv[1];
2048
2049   if (fReadMC && track->GetLabel()<1) return;
2050   if (fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) return;
2051   Double_t ptTrack = track->Pt();
2052   Double_t d0Track = type!=kK0 ? track->GetImpPar() : 0.;
2053   Double_t phiTr = track->Phi();
2054   Double_t origTr = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
2055
2056   TString part = "", orig = "";
2057
2058   switch (type) {
2059     case(kTrack): {
2060       part = "Charg";
2061       deltaphi = fCorrelatorTr->GetDeltaPhi();
2062       deltaeta = fCorrelatorTr->GetDeltaEta();
2063       break;
2064     }
2065     case(kKCharg): {
2066       part = "Kcharg";
2067       deltaphi = fCorrelatorKc->GetDeltaPhi();
2068       deltaeta = fCorrelatorKc->GetDeltaEta();
2069       break;
2070     }
2071     case(kK0): {
2072       part = "K0";
2073       deltaphi = fCorrelatorK0->GetDeltaPhi();
2074       deltaeta = fCorrelatorK0->GetDeltaEta();
2075       break;
2076     }
2077   }
2078
2079   if(fMixing == kSE) {
2080
2081     //Fixes limits; needed to include overflow into THnSparse projections!
2082     Double_t pTorig = track->Pt();
2083     Double_t d0orig = track->GetImpPar();
2084     Double_t ptLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
2085     Double_t displLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
2086     Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
2087     if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
2088     if(d0Track > displLim_Sparse) d0Track = (displLim_Sparse-0.001);
2089     if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
2090     if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
2091
2092     //variables for filling histos
2093     Double_t fillSpPhiD0[5] = {deltaphi,mD0,ptTrack,d0Track,deltaeta};
2094     Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,ptTrack,d0Track,deltaeta};
2095     Double_t fillSpWeigD0[4] = {deltaphi,mD0,deltaeta,ptTrack};
2096     Double_t fillSpWeigD0bar[4] = {deltaphi,mD0bar,deltaeta,ptTrack};
2097
2098     if(fReadMC == 0) {
2099       //sparse fill for data (tracks, K+-, K0) + weighted
2100       if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
2101         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2102         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2103       }
2104       if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
2105         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2106         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2107       }
2108       if(!fAlreadyFilled) {
2109         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
2110         ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
2111       }
2112     }
2113
2114     if(fReadMC) {
2115
2116       if(origD0==4) {orig = "_From_c";} else {orig = "_From_b";}
2117
2118       //sparse fill for data (tracks, K+-, K0) + weighted
2119       if(PdgD0==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0 (from MCTruth)
2120          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2121          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2122          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2123          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2124          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2125          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2126          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2127          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2128          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2129          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2130       }
2131       if(PdgD0==-421 && fIsSelectedCandidate>1) { //D0bar
2132          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2133          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2134          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2135          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2136          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2137          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2138          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2139          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2140          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2141          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2142       }
2143       if(!fAlreadyFilled) {
2144         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_Bin%d",part.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
2145         if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF_Bin%d",part.Data(),ptbin)))->Fill(d0orig);
2146         if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig);
2147         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
2148         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
2149         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_%s_Bin%d",part.Data(),ptbin)))->Fill(origTr);
2150         ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
2151       }
2152     }//end MC case
2153
2154   } //end of SE fill
2155
2156   if(fMixing == kME) {
2157
2158     //Fixes limits; needed to include overflow into THnSparse projections!
2159     Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(2)->GetXmax();
2160     if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
2161     if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
2162     Double_t ptLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(3)->GetXmax(); //all plots have same axes...
2163     if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
2164
2165     //variables for filling histos
2166     Double_t fillSpPhiD0[4] = {deltaphi,mD0,deltaeta,0.4}; //dummy for ME threshold! unless explicitly set by flag...
2167     Double_t fillSpPhiD0bar[4] = {deltaphi,mD0bar,deltaeta,0.4};
2168     if(fMEAxisThresh) {
2169       fillSpPhiD0[3] = ptTrack;
2170       fillSpPhiD0bar[3] = ptTrack;
2171     }
2172
2173     if(fReadMC == 0) {
2174       //sparse fill for data (tracks, K+-, K0)
2175       if(fIsSelectedCandidate == 1||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2176       if(fIsSelectedCandidate == 2||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2177     }
2178     if(fReadMC == 1) {
2179       //sparse fill for data (tracks, K+-, K0)
2180       if(PdgD0==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3))  ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2181       if(PdgD0==-421 && fIsSelectedCandidate>1) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2182     }//end MC case
2183
2184   } //end of ME fill
2185
2186   return;
2187 }
2188
2189 //_________________________________________________________________________________________________
2190 Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
2191   //
2192   // checks on particle (#) origin:
2193   // 0) Not HF
2194   // 1) D->#
2195   // 2) D->X->#
2196   // 3) c hadronization
2197   // 4) B->#
2198   // 5) B->X-># (X!=D)
2199   // 6) B->D->#
2200   // 7) B->D->X->#
2201   // 8) b hadronization
2202   //
2203   if(fDebug>2) printf("AliAnalysisTaskSED0Correlations::CheckTrkOrigin() \n");
2204
2205   Int_t pdgGranma = 0;
2206   Int_t mother = 0;
2207   mother = mcPartCandidate->GetMother();
2208   Int_t istep = 0;
2209   Int_t abspdgGranma =0;
2210   Bool_t isFromB=kFALSE;
2211   Bool_t isDdaugh=kFALSE;
2212   Bool_t isDchaindaugh=kFALSE;
2213   Bool_t isBdaugh=kFALSE;
2214   Bool_t isBchaindaugh=kFALSE;
2215   Bool_t isQuarkFound=kFALSE;
2216
2217   if (mother<0) return -1;
2218   while (mother >= 0){
2219     istep++;
2220     AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
2221     if (mcMoth){
2222       pdgGranma = mcMoth->GetPdgCode();
2223       abspdgGranma = TMath::Abs(pdgGranma);
2224       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
2225         isBchaindaugh=kTRUE;
2226         if(istep==1) isBdaugh=kTRUE;
2227       }
2228       if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)){
2229         isDchaindaugh=kTRUE;
2230         if(istep==1) isDdaugh=kTRUE;
2231       }
2232       if(abspdgGranma==4 || abspdgGranma==5) {isQuarkFound=kTRUE; if(abspdgGranma==5) isFromB = kTRUE;}
2233       mother = mcMoth->GetMother();
2234     }else{
2235       AliError("Failed casting the mother particle!");
2236       return -1;
2237     }
2238   }
2239
2240   //decides what to return based on the flag status
2241   if(isQuarkFound) {
2242     if(!isFromB) {  //charm
2243       if(isDdaugh) return 1; //charm immediate
2244       else if(isDchaindaugh) return 2; //charm chain
2245       else return 3; //charm hadronization
2246     }
2247     else { //beauty
2248       if(isBdaugh) return 4; //b immediate
2249       else if(isBchaindaugh) { //b chain
2250         if(isDchaindaugh) {
2251           if(isDdaugh) return 6; //d immediate
2252           return 7; //d chain
2253           }
2254         else return 5; //b, not d
2255       }
2256       else return 8; //b hadronization
2257     }
2258   }
2259   else if(!isDdaugh && !isDchaindaugh && !isBdaugh && !isBchaindaugh) return 0; //no HF decay
2260      //in this case, it's !isQuarkFound, but not in 100% cases it's a non HF particle!
2261      //rarely you can find a D/B meson which comes from a -1! It isn't a Non-HF, in that case! And I'll return -1...
2262
2263   return -1; //some problem spotted
2264 }
2265 //________________________________________________________________________
2266 Bool_t AliAnalysisTaskSED0Correlations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* mcArray) const {
2267
2268   //Daughter removal in MCKine case
2269   Bool_t isDaughter = kFALSE;
2270   Int_t labelD0 = d->GetLabel();
2271
2272   Int_t mother = track->GetMother();
2273
2274   //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
2275   while (mother > 0){
2276     AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother)); //it's the mother of the track!
2277     if (mcMoth){
2278       if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
2279       mother = mcMoth->GetMother(); //goes back by one
2280     } else{
2281       AliError("Failed casting the mother particle!");
2282       break;
2283     }
2284   }
2285
2286   return isDaughter;
2287 }
2288
2289 //________________________________________________________________________
2290 Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
2291   //
2292   //give the pt bin where the pt lies.
2293   //
2294   Int_t ptbin=-1;
2295   if(pt<fBinLimsCorr.at(0)) return ptbin; //out of bounds
2296
2297   Int_t i = 0;
2298   while(pt>fBinLimsCorr.at(i)) {ptbin=i; i++;}
2299
2300   return ptbin;
2301 }
2302
2303 //---------------------------------------------------------------------------
2304 Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t opt, Int_t idArrayV0[][2]) const
2305 {
2306   //
2307   // Selection for K0 hypotheses
2308   // options: 1 = selects mass invariant about 3 sigma inside the peak + threshold of 0.3 GeV
2309   //          2 = no previous selections
2310
2311   if(!fCutsTracks->IsKZeroSelected(v0,vtx)) return kFALSE;
2312
2313   AliAODTrack *v0Daug1 = (AliAODTrack*)v0->GetDaughter(0);
2314   AliAODTrack *v0Daug2 = (AliAODTrack*)v0->GetDaughter(1);
2315
2316   if(opt==1) { //additional cuts for correlations (V0 has to be closer than 3 sigma from K0 mass)
2317     if(TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.004) return kFALSE;
2318   }
2319
2320   //This part removes double counting for swapped tracks!
2321   Int_t i = 0;  //while loop (until the last-written entry pair of ID!
2322   while(idArrayV0[i][0]!=-2 && idArrayV0[i][1]!=-2) {
2323     if((v0Daug1->GetID()==idArrayV0[i][0] && v0Daug2->GetID()==idArrayV0[i][1])||
2324        (v0Daug1->GetID()==idArrayV0[i][1] && v0Daug2->GetID()==idArrayV0[i][0])) return kFALSE;
2325     i++;
2326   }
2327   idArrayV0[i][0]=v0Daug1->GetID();
2328   idArrayV0[i][1]=v0Daug2->GetID();
2329
2330   return kTRUE;
2331 }
2332
2333 //---------------------------------------------------------------------------
2334 Bool_t AliAnalysisTaskSED0Correlations::IsSoftPion_MCKine(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* arrayMC) const
2335 {
2336   //
2337   // Removes soft pions in Kine
2338
2339   //Daughter removal in MCKine case
2340   Bool_t isSoftPi = kFALSE;
2341   Int_t labelD0 = d->GetLabel();
2342
2343   Int_t mother = track->GetMother();
2344   if(mother<0) return isSoftPi; //safety check
2345
2346   AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother)); //it's the mother of the track!
2347   if(!mcMoth){
2348     return isSoftPi;
2349   }
2350   if(TMath::Abs(mcMoth->GetPdgCode())==413 && mcMoth->GetNDaughters()==2) { //mother is D* with 2 daughs
2351     Int_t labdau1 = mcMoth->GetDaughter(0);
2352     Int_t labdau2 = mcMoth->GetDaughter(1);
2353     AliAODMCParticle* dau1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labdau1));
2354     AliAODMCParticle* dau2 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labdau2));
2355     if(!dau1 || !dau2) return isSoftPi; //safety check
2356     if(dau1->GetLabel()==labelD0 || dau2->GetLabel()==labelD0) { //one of the daughs is the D0 trigger
2357       if((TMath::Abs(dau1->GetPdgCode())==421 && TMath::Abs(dau2->GetPdgCode())==211)||(TMath::Abs(dau1->GetPdgCode())==211 && TMath::Abs(dau2->GetPdgCode())==421)) {
2358         isSoftPi = kTRUE; //ok, soft pion was found
2359         return isSoftPi;
2360       }
2361     }
2362   }
2363
2364   return isSoftPi;
2365 }
2366
2367 //________________________________________________________________________
2368 void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
2369
2370   cout << "--------------------------\n";
2371   cout << "PtBins = " << fNPtBinsCorr << "\n";
2372   cout << "PtBin limits--------------\n";
2373   for (int i=0; i<fNPtBinsCorr; i++) {
2374     cout << "Bin "<<i+1<<" = "<<fBinLimsCorr.at(i)<<" to "<<fBinLimsCorr.at(i+1)<<"\n";
2375   }
2376   cout << "\n--------------------------\n";
2377   cout << "PtBin tresh. tracks low---\n";
2378   for (int i=0; i<fNPtBinsCorr; i++) {
2379     cout << fPtThreshLow.at(i) << "; ";
2380   }
2381   cout << "PtBin tresh. tracks up----\n";
2382   for (int i=0; i<fNPtBinsCorr; i++) {
2383     cout << fPtThreshUp.at(i) << "; ";
2384   }
2385   cout << "\n--------------------------\n";
2386   cout << "D0 Eta cut for Correl = "<<fEtaForCorrel<<"\n";
2387   cout << "--------------------------\n";
2388   cout << "MC Truth = "<<fReadMC<<" - MC Reco: Trk = "<<fRecoTr<<", D0 = "<<fRecoD0<<"\n";
2389   cout << "--------------------------\n";
2390   cout << "Sel of Event tpye (PP,GS,FE,...)= "<<fSelEvType<<"\n";
2391   cout << "--------------------------\n";
2392   cout << "Ev Mixing = "<<fMixing<<"\n";
2393   cout << "--------------------------\n";
2394   cout << "ME thresh axis = "<<fMEAxisThresh<<"\n";
2395   cout << "--------------------------\n";
2396   cout << "Soft Pi Cut = "<<fSoftPiCut<<"\n";
2397   cout << "--------------------------\n";
2398 }
2399