PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //*************************************************************************
  19 // Class AliAnalysisTaskSEDvsMultiplicity
  20 // AliAnalysisTaskSE for the D meson vs. multiplcity analysis
  21 // Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TClonesArray.h>
  25 #include <TCanvas.h>
  26 #include <TList.h>
  27 #include <TString.h>
  28 #include <TDatabasePDG.h>
  29 #include <TH1F.h>
  30 #include <TH2F.h>
  31 #include <TH3F.h>
  32 #include <THnSparse.h>
  33 #include <TProfile.h>
  34 #include "AliAnalysisManager.h"
  35 #include "AliRDHFCuts.h"
  36 #include "AliRDHFCutsDplustoKpipi.h"
  37 #include "AliRDHFCutsDStartoKpipi.h"
  38 #include "AliRDHFCutsD0toKpi.h"
  39 #include "AliAODHandler.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODVertex.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODRecoDecayHF.h"
  44 #include "AliAODRecoCascadeHF.h"
  45 #include "AliAnalysisVertexingHF.h"
  46 #include "AliAnalysisTaskSE.h"
  47 #include "AliAnalysisTaskSEDvsMultiplicity.h"
  48 #include "AliNormalizationCounter.h"
  49 #include "AliVertexingHFUtils.h"
  50 #include "AliAODVZERO.h"
  51 ClassImp(AliAnalysisTaskSEDvsMultiplicity)
  52
  53
  54 //________________________________________________________________________
  55 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
  56 AliAnalysisTaskSE(),
  57   fOutput(0),
  58   fListCuts(0),
  59   fOutputCounters(0),
  60   fListProfiles(0),
  61   fHistNEvents(0),
  62   fHistNtrEta16vsNtrEta1(0),
  63   fHistNtrCorrEta1vsNtrRawEta1(0),
  64   fHistNtrVsZvtx(0),
  65   fHistNtrCorrVsZvtx(0),
  66   fHistNtrVsNchMC(0),
  67   fHistNtrCorrVsNchMC(0),
  68   fHistNtrVsNchMCPrimary(0),
  69   fHistNtrCorrVsNchMCPrimary(0),
  70   fHistNtrVsNchMCPhysicalPrimary(0),
  71   fHistNtrCorrVsNchMCPhysicalPrimary(0),
  72   fHistGenPrimaryParticlesInelGt0(0),
  73   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
  74   fHistNtrUnCorrEvSel(0),
  75   fHistNtrCorrEvSel(0),
  76   fHistNtrCorrEvWithCand(0),
  77   fHistNtrCorrEvWithD(0),
  78   fPtVsMassVsMult(0),
  79   fPtVsMassVsMultNoPid(0),
  80   fPtVsMassVsMultUncorr(0),
  81   fPtVsMassVsMultPart(0),
  82   fPtVsMassVsMultAntiPart(0),
  83   fUpmasslimit(1.965),
  84   fLowmasslimit(1.765),
  85   fNMassBins(200),
  86   fRDCutsAnalysis(0),
  87   fCounter(0),
  88   fCounterU(0),
  89   fDoImpPar(kFALSE),
  90   fNImpParBins(400),
  91   fLowerImpPar(-2000.),
  92   fHigherImpPar(2000.),
  93   fReadMC(kFALSE),
  94   fMCOption(0),
  95   fUseBit(kTRUE),
  96   fSubtractTrackletsFromDau(kFALSE),
  97   fUseNchWeight(kFALSE),
  98   fHistoMCNch(0),
  99   fHistoMeasNch(0),
 100   fRefMult(9.26),
 101   fPdgMeson(411),
 102   fMultiplicityEstimator(kNtrk10)
 103 {
 104    // Default constructor
 105   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
 106   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 107 }
 108
 109 //________________________________________________________________________
 110 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts):
 111   AliAnalysisTaskSE(name),
 112   fOutput(0),
 113   fListCuts(0),
 114   fOutputCounters(0),
 115   fListProfiles(0),
 116   fHistNEvents(0),
 117   fHistNtrEta16vsNtrEta1(0),
 118   fHistNtrCorrEta1vsNtrRawEta1(0),
 119   fHistNtrVsZvtx(0),
 120   fHistNtrCorrVsZvtx(0),
 121   fHistNtrVsNchMC(0),
 122   fHistNtrCorrVsNchMC(0),
 123   fHistNtrVsNchMCPrimary(0),
 124   fHistNtrCorrVsNchMCPrimary(0),
 125   fHistNtrVsNchMCPhysicalPrimary(0),
 126   fHistNtrCorrVsNchMCPhysicalPrimary(0),
 127   fHistGenPrimaryParticlesInelGt0(0),
 128   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
 129   fHistNtrUnCorrEvSel(0),
 130   fHistNtrCorrEvSel(0),
 131   fHistNtrCorrEvWithCand(0),
 132   fHistNtrCorrEvWithD(0),
 133   fPtVsMassVsMult(0),
 134   fPtVsMassVsMultNoPid(0),
 135   fPtVsMassVsMultUncorr(0),
 136   fPtVsMassVsMultPart(0),
 137   fPtVsMassVsMultAntiPart(0),
 138   fUpmasslimit(1.965),
 139   fLowmasslimit(1.765),
 140   fNMassBins(200),
 141   fRDCutsAnalysis(cuts),
 142   fCounter(0),
 143   fCounterU(0),
 144   fDoImpPar(kFALSE),
 145   fNImpParBins(400),
 146   fLowerImpPar(-2000.),
 147   fHigherImpPar(2000.),
 148   fReadMC(kFALSE),
 149   fMCOption(0),
 150   fUseBit(kTRUE),
 151   fSubtractTrackletsFromDau(kFALSE),
 152   fUseNchWeight(kFALSE),
 153   fHistoMCNch(0),
 154   fHistoMeasNch(0),
 155   fRefMult(9.26),
 156   fPdgMeson(pdgMeson),
 157   fMultiplicityEstimator(kNtrk10)
 158 {
 159   //
 160   // Standard constructor
 161   //
 162   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
 163   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 164   if(fPdgMeson==413){
 165     fNMassBins=200; // FIXME
 166     SetMassLimits(0.,0.2); // FIXME
 167   }else{
 168     fNMassBins=200;
 169     SetMassLimits(fPdgMeson,0.1);
 170   }
 171   // Default constructor
 172    // Otput slot #1 writes into a TList container
 173   DefineOutput(1,TList::Class());  //My private output
 174   // Output slot #2 writes cut to private output
 175   DefineOutput(2,TList::Class());
 176   // Output slot #3 writes cut to private output
 177   DefineOutput(3,TList::Class());
 178   // Output slot #4 writes cut to private output
 179   DefineOutput(4,TList::Class());
 180 }
 181 //________________________________________________________________________
 182 AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
 183 {
 184   //
 185   // Destructor
 186   //
 187   delete fOutput;
 188   delete fHistNEvents;
 189   delete fListCuts;
 190   delete fListProfiles;
 191   delete fRDCutsAnalysis;
 192   delete fCounter;
 193   delete fCounterU;
 194   for(Int_t i=0; i<4; i++) delete fMultEstimatorAvg[i];
 195   for(Int_t i=0; i<5; i++){
 196     delete fHistMassPtImpPar[i];
 197   }
 198   if(fHistoMCNch) delete fHistoMCNch;
 199   if(fHistoMeasNch) delete fHistoMeasNch;
 200 }
 201
 202 //_________________________________________________________________
 203 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 204   // set invariant mass limits
 205   if(uplimit>lowlimit){
 206     fLowmasslimit = lowlimit;
 207     fUpmasslimit = uplimit;
 208   }else{
 209     AliError("Wrong mass limits: upper value should be larger than lower one");
 210   }
 211 }
 212 //_________________________________________________________________
 213 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
 214   // set invariant mass limits
 215   Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
 216   SetMassLimits(mass-range,mass+range);
 217 }
 218 //________________________________________________________________________
 219 void AliAnalysisTaskSEDvsMultiplicity::Init(){
 220   //
 221   // Initialization
 222   //
 223   printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
 224
 225   if(fUseNchWeight && !fReadMC){ AliFatal("Nch weights can only be used in MC mode"); return; }
 226   if(fUseNchWeight && !fHistoMCNch){ AliFatal("Nch weights can only be used without histogram"); return; }
 227
 228   fListCuts=new TList();
 229
 230   if(fPdgMeson==411){
 231      AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
 232     copycut->SetName("AnalysisCutsDplus");
 233     fListCuts->Add(copycut);
 234   }else if(fPdgMeson==421){
 235     AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
 236     copycut->SetName("AnalysisCutsDzero");
 237     fListCuts->Add(copycut);
 238   }else if(fPdgMeson==413){
 239      AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
 240     copycut->SetName("AnalysisCutsDStar");
 241     fListCuts->Add(copycut);
 242   }
 243   PostData(2,fListCuts);
 244
 245   fListProfiles = new TList();
 246   fListProfiles->SetOwner();
 247   TString period[4]={"LHC10b","LHC10c","LHC10d","LHC10e"};
 248   for(Int_t i=0; i<4; i++){
 249     if(fMultEstimatorAvg[i]){
 250       TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
 251       hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
 252       fListProfiles->Add(hprof);
 253     }
 254   }
 255   PostData(4,fListProfiles);
 256
 257   return;
 258 }
 259
 260 //________________________________________________________________________
 261 void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
 262 {
 263   // Create the output container
 264   //
 265   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
 266
 267   // Several histograms are more conveniently managed in a TList
 268   fOutput = new TList();
 269   fOutput->SetOwner();
 270   fOutput->SetName("OutputHistos");
 271
 272   Int_t nMultBins = 200;
 273   Float_t firstMultBin = -0.5;
 274   Float_t lastMultBin = 199.5;
 275   if(fMultiplicityEstimator==kVZERO) {
 276     nMultBins = 600;
 277     lastMultBin = 599.5;
 278   }
 279
 280   fHistNtrUnCorrEvSel = new TH1F("hNtrUnCorrEvSel","Uncorrected tracklets multiplicity for selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
 281   fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Corrected tracklets multiplicity for selected events; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);
 282   fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin);// Total multiplicity
 283   fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",nMultBins,firstMultBin,lastMultBin); //
 284   fHistNtrEta16vsNtrEta1 = new TH2F("hNtrEta16vsNtrEta1","Uncorrected Eta1.6 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<1.6",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 1.6 vs eta 1.0 histogram
 285   fHistNtrCorrEta1vsNtrRawEta1 = new TH2F("hNtrCorrEta1vsNtrRawEta1","Corrected Eta1 vs Eta1.0; Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //eta 1.6 vs eta 1.0 histogram
 286   fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,nMultBins,firstMultBin,lastMultBin); //
 287   fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,nMultBins,firstMultBin,lastMultBin); //
 288
 289   fHistNtrVsNchMC = new TH2F("hNtrVsNchMC","Ntracklet vs NchMC; Nch;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 290   fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC","Ntracklet vs Nch; Nch;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 291
 292   fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch (Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 293   fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch(Primary) ;N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 294
 295   fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 296   fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin); //
 297
 298   fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",nMultBins,firstMultBin,lastMultBin);
 299
 300   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary = new TH3F("fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary", "MC: Nch (Physical Primary) vs Nch (Primary) vs Nch (Generated); Nch (Generated); Nch (Primary); Nch (Physical Primary)",nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin,nMultBins,firstMultBin,lastMultBin);
 301
 302   fHistNtrUnCorrEvSel->Sumw2();
 303   fHistNtrCorrEvSel->Sumw2();
 304   fHistNtrCorrEvWithCand->Sumw2();
 305   fHistNtrCorrEvWithD->Sumw2();
 306   fHistGenPrimaryParticlesInelGt0->Sumw2();
 307   fOutput->Add(fHistNtrUnCorrEvSel);
 308   fOutput->Add(fHistNtrCorrEvSel);
 309   fOutput->Add(fHistNtrCorrEvWithCand);
 310   fOutput->Add(fHistNtrCorrEvWithD);
 311   fOutput->Add(fHistNtrEta16vsNtrEta1);
 312   fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1);
 313   fOutput->Add(fHistNtrVsZvtx);
 314   fOutput->Add(fHistNtrCorrVsZvtx);
 315
 316   fOutput->Add(fHistNtrVsNchMC);
 317   fOutput->Add(fHistNtrCorrVsNchMC);
 318   fOutput->Add(fHistNtrVsNchMCPrimary);
 319   fOutput->Add(fHistNtrCorrVsNchMCPrimary);
 320   fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
 321   fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
 322   fOutput->Add(fHistGenPrimaryParticlesInelGt0);
 323   fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
 324
 325
 326   fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
 327   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
 328   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
 329   fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
 330   fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
 331   fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
 332   fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
 333   fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
 334   fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
 335   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
 336   fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
 337   fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
 338   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 339   fHistNEvents->Sumw2();
 340   fHistNEvents->SetMinimum(0);
 341   fOutput->Add(fHistNEvents);
 342
 343   fPtVsMassVsMult=new TH3F("hPtVsMassvsMult", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 344
 345   fPtVsMassVsMultNoPid=new TH3F("hPtVsMassvsMultNoPid", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 346
 347   fPtVsMassVsMultUncorr=new TH3F("hPtVsMassvsMultUncorr", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 348
 349   fPtVsMassVsMultPart=new TH3F("hPtVsMassvsMultPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 350
 351   fPtVsMassVsMultAntiPart=new TH3F("hPtVsMassvsMultAntiPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",nMultBins,firstMultBin,lastMultBin,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 352
 353   fOutput->Add(fPtVsMassVsMult);
 354   fOutput->Add(fPtVsMassVsMultUncorr);
 355   fOutput->Add(fPtVsMassVsMultNoPid);
 356   fOutput->Add(fPtVsMassVsMultPart);
 357   fOutput->Add(fPtVsMassVsMultAntiPart);
 358
 359   if(fDoImpPar) CreateImpactParameterHistos();
 360
 361   fCounter = new AliNormalizationCounter("NormCounterCorrMult");
 362   fCounter->SetStudyMultiplicity(kTRUE,1.);
 363   fCounter->Init();
 364
 365   fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
 366   fCounterU->SetStudyMultiplicity(kTRUE,1.);
 367   fCounterU->Init();
 368
 369   fOutputCounters = new TList();
 370   fOutputCounters->SetOwner();
 371   fOutputCounters->SetName("OutputCounters");
 372   fOutputCounters->Add(fCounter);
 373   fOutputCounters->Add(fCounterU);
 374
 375   PostData(1,fOutput);
 376   PostData(2,fListCuts);
 377   PostData(3,fOutputCounters);
 378   PostData(4,fListProfiles);
 379
 380   if(fUseNchWeight) CreateMeasuredNchHisto();
 381
 382   return;
 383 }
 384
 385
 386
 387 //________________________________________________________________________
 388 void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
 389 {
 390   // Execute analysis for current event:
 391   // heavy flavor candidates association to MC truth
 392
 393   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 394
 395   //  AliAODTracklets* tracklets = aod->GetTracklets();
 396   //Int_t ntracklets = tracklets->GetNumberOfTracklets();
 397
 398
 399   TClonesArray *arrayCand = 0;
 400   TString arrayName="";
 401   UInt_t pdgDau[3];
 402   Int_t nDau=0;
 403   Int_t selbit=0;
 404   if(fPdgMeson==411){
 405     arrayName="Charm3Prong";
 406     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
 407     nDau=3;
 408     selbit=AliRDHFCuts::kDplusCuts;
 409   }else if(fPdgMeson==421){
 410     arrayName="D0toKpi";
 411     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
 412     nDau=2;
 413     selbit=AliRDHFCuts::kD0toKpiCuts;
 414   }else if(fPdgMeson==413){
 415     arrayName="Dstar";
 416     pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=211;
 417     nDau=3;
 418     selbit=AliRDHFCuts::kDstarCuts;
 419   }
 420
 421   if(!aod && AODEvent() && IsStandardAOD()) {
 422     // In case there is an AOD handler writing a standard AOD, use the AOD
 423     // event in memory rather than the input (ESD) event.
 424     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 425      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 426      // have to taken from the AOD event hold by the AliAODExtension
 427     AliAODHandler* aodHandler = (AliAODHandler*)
 428       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 429     if(aodHandler->GetExtensions()) {
 430       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 431       AliAODEvent *aodFromExt = ext->GetAOD();
 432       arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
 433     }
 434   } else if(aod) {
 435     arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
 436   }
 437
 438   if(!aod || !arrayCand) {
 439     printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
 440     return;
 441   }
 442
 443
 444
 445   // fix for temporary bug in ESDfilter
 446   // the AODs with null vertex pointer didn't pass the PhysSel
 447   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 448
 449   Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
 450   Int_t countTr=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
 451
 452   Int_t vzeroMult=0;
 453   AliAODVZERO *vzeroAOD = (AliAODVZERO*)aod->GetVZEROData();
 454   if(vzeroAOD) vzeroMult = vzeroAOD->GetMTotV0A() +  vzeroAOD->GetMTotV0C();
 455
 456   Int_t countMult = countTreta1;
 457   if(fMultiplicityEstimator==kNtrk10to16) { countMult = countTr - countTreta1; }
 458   if(fMultiplicityEstimator==kVZERO) { countMult = vzeroMult; }
 459
 460
 461   fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countMult);
 462   fHistNEvents->Fill(0); // count event
 463
 464   Double_t countTreta1corr=countTreta1;
 465   Double_t countCorr=countMult;
 466   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 467   //  if(vtx1){
 468   // FIX ME: No correction to the VZERO !!
 469   if(vtx1 && (fMultiplicityEstimator!=kVZERO)){
 470     if(vtx1->GetNContributors()>0){
 471       fHistNEvents->Fill(1);
 472       TProfile* estimatorAvg = GetEstimatorHistogram(aod);
 473       if(estimatorAvg){
 474         countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
 475         countCorr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countMult,vtx1->GetZ(),fRefMult);
 476       }
 477     }
 478   }
 479
 480
 481   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 482
 483   if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
 484   if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
 485   if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
 486   if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
 487
 488
 489   if(!isEvSel)return;
 490   fHistNtrEta16vsNtrEta1->Fill(countTreta1,countTr);
 491   fHistNtrCorrEta1vsNtrRawEta1->Fill(countTreta1,countTreta1corr);
 492   if(vtx1){
 493     fHistNtrVsZvtx->Fill(vtx1->GetZ(),countMult);
 494     fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countCorr);
 495   }
 496
 497   TClonesArray *arrayMC=0;
 498   AliAODMCHeader *mcHeader=0;
 499
 500   Double_t nchWeight=1.0;
 501
 502   // load MC particles
 503   if(fReadMC){
 504
 505     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 506     if(!arrayMC) {
 507       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
 508       return;
 509     }
 510     // load MC header
 511     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 512     if(!mcHeader) {
 513       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
 514       return;
 515      }
 516
 517
 518     Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
 519     Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
 520     Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
 521     if(fUseNchWeight){
 522       Double_t tmpweight = 1.0;
 523       if(nChargedMCPhysicalPrimary<=0) tmpweight = 0.0;
 524       else{
 525         Double_t pMeas=fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nChargedMCPhysicalPrimary));
 526         //      printf(" pMeas=%2.2f  and histo MCNch %s \n",pMeas,fHistoMCNch);
 527         Double_t pMC=fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nChargedMCPhysicalPrimary));
 528         tmpweight = pMeas/pMC;
 529       }
 530       nchWeight *= tmpweight;
 531       AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,nchWeight));
 532     }
 533     if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
 534       fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary,nchWeight);
 535     }
 536
 537     fHistNtrVsNchMC->Fill(nChargedMC,countMult,nchWeight);
 538     fHistNtrCorrVsNchMC->Fill(nChargedMC,countCorr,nchWeight);
 539
 540     fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countMult,nchWeight);
 541     fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countCorr,nchWeight);
 542
 543     fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countMult,nchWeight);
 544     fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countCorr,nchWeight);
 545
 546     fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary,nchWeight);
 547   }
 548
 549   Int_t nCand = arrayCand->GetEntriesFast();
 550   Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
 551   Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 552   Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
 553   Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 554   Double_t aveMult=0.;
 555   Double_t nSelCand=0.;
 556   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 557     AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
 558     fHistNEvents->Fill(7);
 559     if(fUseBit && !d->HasSelectionBit(selbit)){
 560       fHistNEvents->Fill(8);
 561       continue;
 562     }
 563
 564     Double_t ptCand = d->Pt();
 565     Double_t rapid=d->Y(fPdgMeson);
 566     Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 567     if(!isFidAcc) continue;
 568     Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
 569     Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
 570     if(passTopolCuts==0) continue;
 571     nSelectedNoPID++;
 572     fHistNEvents->Fill(9);
 573     if(passAllCuts){
 574       nSelectedPID++;
 575       fHistNEvents->Fill(10);
 576     }
 577     Double_t multForCand = countCorr;
 578     if(fSubtractTrackletsFromDau){
 579       for(Int_t iDau=0; iDau<nDau; iDau++){
 580         AliAODTrack *t = (AliAODTrack*)d->GetDaughter(iDau);
 581         if(!t) continue;
 582         if(t->HasPointOnITSLayer(0) && t->HasPointOnITSLayer(1)){
 583           if(multForCand>0) multForCand-=1;
 584         }
 585       }
 586     }
 587     Bool_t isPrimary=kTRUE;
 588     Int_t labD=-1;
 589     Double_t trueImpParXY=9999.;
 590     Double_t impparXY=d->ImpParXY()*10000.;
 591     Double_t dlen=0.1; //FIXME
 592     Double_t mass[2];
 593     if(fPdgMeson==411){
 594       mass[0]=d->InvMass(nDau,pdgDau);
 595       mass[1]=-1.;
 596       if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
 597     }else if(fPdgMeson==421){
 598       UInt_t pdgdaughtersD0[2]={211,321};//pi,K
 599       UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
 600       mass[0]=d->InvMass(2,pdgdaughtersD0);
 601       mass[1]=d->InvMass(2,pdgdaughtersD0bar);
 602       if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
 603     }else if(fPdgMeson==413){
 604       // FIXME
 605       AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
 606       mass[0]=temp->DeltaInvMass();
 607       mass[1]=-1.;
 608       if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.001) nSelectedInMassPeak++; //1 MeV for now... FIXME
 609     }
 610     for(Int_t iHyp=0; iHyp<2; iHyp++){
 611       if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
 612       Double_t invMass=mass[iHyp];
 613       Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,multForCand};
 614
 615       if(fReadMC){
 616
 617         labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
 618         Bool_t fillHisto=fDoImpPar;
 619         if(labD>=0){
 620           AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
 621           Int_t code=partD->GetPdgCode();
 622           if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
 623           if(code<0 && iHyp==0) fillHisto=kFALSE;
 624           if(code>0 && iHyp==1) fillHisto=kFALSE;
 625           if(!isPrimary){
 626             if(fPdgMeson==411){
 627               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
 628             }else if(fPdgMeson==421){
 629               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
 630             }else if(fPdgMeson==413){
 631               trueImpParXY=0.; /// FIXME
 632             }
 633             Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,multForCand};
 634             if(fillHisto && passAllCuts){
 635               fHistMassPtImpPar[2]->Fill(arrayForSparse);
 636               fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
 637             }
 638           }else{
 639             if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
 640           }
 641         }else{
 642           if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
 643         }
 644         if(fPdgMeson==421){
 645           if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
 646           if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
 647         }
 648       }
 649
 650       if(fPdgMeson==421){
 651         if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
 652         if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
 653       }
 654
 655       fPtVsMassVsMultNoPid->Fill(multForCand,invMass,ptCand);
 656
 657       if(fPdgMeson==421){
 658         if(iHyp==0 && !(passAllCuts&1)) continue; // candidate not passing as D0
 659         if(iHyp==1 && !(passAllCuts&2)) continue; // candidate not passing as D0bar
 660       }
 661       if(passAllCuts){
 662         aveMult+=multForCand;
 663         nSelCand+=1.;
 664         fPtVsMassVsMult->Fill(multForCand,invMass,ptCand,nchWeight);
 665         fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand,nchWeight);
 666         // Add separation between part antipart
 667         if(fPdgMeson==411){
 668           if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
 669           else fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
 670         }else if(fPdgMeson==421){
 671           if(passAllCuts&1) fPtVsMassVsMultPart->Fill(multForCand,invMass,ptCand,nchWeight);
 672           if(passAllCuts&2) fPtVsMassVsMultAntiPart->Fill(multForCand,invMass,ptCand,nchWeight);
 673         }else if(fPdgMeson==413){
 674           // FIXME ADD Dstar!!!!!!!!
 675         }
 676
 677         if(fDoImpPar){
 678           fHistMassPtImpPar[0]->Fill(arrayForSparse);
 679         }
 680
 681       }
 682
 683     }
 684   }
 685   if(fSubtractTrackletsFromDau && nSelCand>0){
 686     aveMult/=nSelCand;
 687     fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)(aveMult+0.5001));
 688   }else{
 689     fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countTreta1corr);
 690   }
 691
 692
 693   fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
 694   fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
 695   fHistNtrUnCorrEvSel->Fill(countTreta1,nchWeight);
 696   fHistNtrCorrEvSel->Fill(countTreta1corr,nchWeight);
 697   if(nSelectedPID>0) fHistNtrCorrEvWithCand->Fill(countTreta1corr,nchWeight);
 698   if(nSelectedInMassPeak>0) fHistNtrCorrEvWithD->Fill(countTreta1corr,nchWeight);
 699
 700   PostData(1,fOutput);
 701   PostData(2,fListCuts);
 702   PostData(3,fOutputCounters);
 703
 704   return;
 705 }
 706
 707 //________________________________________________________________________
 708 void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
 709   // Histos for impact paramter study
 710   // mass . pt , impact parameter , decay length , multiplicity
 711
 712   Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
 713   Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
 714   Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
 715
 716   fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
 717                                         "Mass vs. pt vs.imppar - All",
 718                                         5,nbins,xmin,xmax);
 719   fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
 720                                         "Mass vs. pt vs.imppar - promptD",
 721                                         5,nbins,xmin,xmax);
 722   fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
 723                                         "Mass vs. pt vs.imppar - DfromB",
 724                                         5,nbins,xmin,xmax);
 725   fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
 726                                         "Mass vs. pt vs.true imppar -DfromB",
 727                                         5,nbins,xmin,xmax);
 728   fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
 729                                         "Mass vs. pt vs.imppar - backgr.",
 730                                         5,nbins,xmin,xmax);
 731   for(Int_t i=0; i<5;i++){
 732     fOutput->Add(fHistMassPtImpPar[i]);
 733   }
 734 }
 735
 736 //________________________________________________________________________
 737 void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
 738 {
 739   // Terminate analysis
 740   //
 741   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
 742
 743   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 744   if (!fOutput) {
 745     printf("ERROR: fOutput not available\n");
 746     return;
 747   }
 748
 749   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 750   if(!fHistNEvents){
 751     printf("ERROR: fHistNEvents not available\n");
 752     return;
 753   }
 754   printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
 755
 756   return;
 757 }
 758 //_________________________________________________________________________________________________
 759 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
 760   //
 761   // checking whether the mother of the particles come from a charm or a bottom quark
 762   //
 763
 764   Int_t pdgGranma = 0;
 765   Int_t mother = 0;
 766   mother = mcPartCandidate->GetMother();
 767   Int_t istep = 0;
 768   Int_t abspdgGranma =0;
 769   Bool_t isFromB=kFALSE;
 770   Bool_t isQuarkFound=kFALSE;
 771   while (mother >0 ){
 772     istep++;
 773     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 774     if (mcGranma){
 775       pdgGranma = mcGranma->GetPdgCode();
 776       abspdgGranma = TMath::Abs(pdgGranma);
 777       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 778         isFromB=kTRUE;
 779       }
 780       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 781       mother = mcGranma->GetMother();
 782     }else{
 783       AliError("Failed casting the mother particle!");
 784       break;
 785     }
 786   }
 787
 788   if(isFromB) return 5;
 789   else return 4;
 790 }
 791
 792
 793
 794 //____________________________________________________________________________
 795 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
 796   // Get Estimator Histogram from period event->GetRunNumber();
 797   //
 798   // If you select SPD tracklets in |eta|<1 you should use type == 1
 799   //
 800
 801   Int_t runNo  = event->GetRunNumber();
 802   Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
 803   if(runNo>114930 && runNo<117223) period = 0;
 804   if(runNo>119158 && runNo<120830) period = 1;
 805   if(runNo>122373 && runNo<126438) period = 2;
 806   if(runNo>127711 && runNo<130841) period = 3;
 807   if(period<0 || period>3) return 0;
 808
 809   return fMultEstimatorAvg[period];
 810 }
 811
 812 //__________________________________________________________________________________________________
 813 void AliAnalysisTaskSEDvsMultiplicity::CreateMeasuredNchHisto(){
 814   // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
 815   Double_t nchbins[66]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
 816                         10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
 817                         20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
 818                         30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
 819                         40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
 820                         50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
 821                         60.50,62.50,64.50,66.50,68.50,70.50};
 822   Double_t pch[65]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
 823                     0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
 824                     0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
 825                     0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
 826                     0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
 827                     0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
 828                     0.000296,0.000265,0.000193,0.00016,0.000126};
 829
 830   if(fHistoMeasNch) delete fHistoMeasNch;
 831   fHistoMeasNch=new TH1F("hMeaseNch","",65,nchbins);
 832   for(Int_t i=0; i<65; i++){
 833     fHistoMeasNch->SetBinContent(i+1,pch[i]);
 834     fHistoMeasNch->SetBinError(i+1,0.);
 835   }
 836 }