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 ClassImp(AliAnalysisTaskSEDvsMultiplicity)
  51
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
  55 AliAnalysisTaskSE(),
  56   fOutput(0),
  57   fListCuts(0),
  58   fOutputCounters(0),
  59   fListProfiles(0),
  60   fHistNEvents(0),
  61   fHistNtrEta16vsNtrEta1(0),
  62   fHistNtrCorrEta1vsNtrRawEta1(0),
  63   fHistNtrVsZvtx(0),
  64   fHistNtrCorrVsZvtx(0),
  65   fHistNtrVsNchMC(0),
  66   fHistNtrCorrVsNchMC(0),
  67   fHistNtrVsNchMCPrimary(0),
  68   fHistNtrCorrVsNchMCPrimary(0),
  69   fHistNtrVsNchMCPhysicalPrimary(0),
  70   fHistNtrCorrVsNchMCPhysicalPrimary(0),
  71   fHistGenPrimaryParticlesInelGt0(0),
  72   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
  73   fHistNtrUnCorrEvSel(0),
  74   fHistNtrCorrEvSel(0),
  75   fHistNtrCorrEvWithCand(0),
  76   fHistNtrCorrEvWithD(0),
  77   fPtVsMassVsMult(0),
  78   fPtVsMassVsMultNoPid(0),
  79   fPtVsMassVsMultUncorr(0),
  80   fPtVsMassVsMultPart(0),
  81   fPtVsMassVsMultAntiPart(0),
  82   fUpmasslimit(1.965),
  83   fLowmasslimit(1.765),
  84   fNMassBins(200),
  85   fRDCutsAnalysis(0),
  86   fCounter(0),
  87   fCounterU(0),
  88   fDoImpPar(kFALSE),
  89   fNImpParBins(400),
  90   fLowerImpPar(-2000.),
  91   fHigherImpPar(2000.),
  92   fReadMC(kFALSE),
  93   fMCOption(0),
  94   fUseBit(kTRUE),
  95   fRefMult(9.26),
  96   fPdgMeson(411)
  97 {
  98    // Default constructor
  99   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
 100   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 101 }
 102
 103 //________________________________________________________________________
 104 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts):
 105   AliAnalysisTaskSE(name),
 106   fOutput(0),
 107   fListCuts(0),
 108   fOutputCounters(0),
 109   fListProfiles(0),
 110   fHistNEvents(0),
 111   fHistNtrEta16vsNtrEta1(0),
 112   fHistNtrCorrEta1vsNtrRawEta1(0),
 113   fHistNtrVsZvtx(0),
 114   fHistNtrCorrVsZvtx(0),
 115   fHistNtrVsNchMC(0),
 116   fHistNtrCorrVsNchMC(0),
 117   fHistNtrVsNchMCPrimary(0),
 118   fHistNtrCorrVsNchMCPrimary(0),
 119   fHistNtrVsNchMCPhysicalPrimary(0),
 120   fHistNtrCorrVsNchMCPhysicalPrimary(0),
 121   fHistGenPrimaryParticlesInelGt0(0),
 122   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
 123   fHistNtrUnCorrEvSel(0),
 124   fHistNtrCorrEvSel(0),
 125   fHistNtrCorrEvWithCand(0),
 126   fHistNtrCorrEvWithD(0),
 127   fPtVsMassVsMult(0),
 128   fPtVsMassVsMultNoPid(0),
 129   fPtVsMassVsMultUncorr(0),
 130   fPtVsMassVsMultPart(0),
 131   fPtVsMassVsMultAntiPart(0),
 132   fUpmasslimit(1.965),
 133   fLowmasslimit(1.765),
 134   fNMassBins(200),
 135   fRDCutsAnalysis(cuts),
 136   fCounter(0),
 137   fCounterU(0),
 138   fDoImpPar(kFALSE),
 139   fNImpParBins(400),
 140   fLowerImpPar(-2000.),
 141   fHigherImpPar(2000.),
 142   fReadMC(kFALSE),
 143   fMCOption(0),
 144   fUseBit(kTRUE),
 145   fRefMult(9.26),
 146   fPdgMeson(pdgMeson)
 147 {
 148   //
 149   // Standard constructor
 150   //
 151   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
 152   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 153   if(fPdgMeson==413){
 154     fNMassBins=200; // FIXME
 155     SetMassLimits(0.,0.2); // FIXME
 156   }else{
 157     fNMassBins=200;
 158     SetMassLimits(fPdgMeson,0.1);
 159   }
 160   // Default constructor
 161    // Output slot #1 writes into a TList container
 162   DefineOutput(1,TList::Class());  //My private output
 163   // Output slot #2 writes cut to private output
 164   DefineOutput(2,TList::Class());
 165   // Output slot #3 writes cut to private output
 166   DefineOutput(3,TList::Class());
 167   // Output slot #4 writes cut to private output
 168   DefineOutput(4,TList::Class());
 169 }
 170 //________________________________________________________________________
 171 AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
 172 {
 173   //
 174   // Destructor
 175   //
 176   delete fOutput;
 177   delete fHistNEvents;
 178   delete fListCuts;
 179   delete fListProfiles;
 180   delete fRDCutsAnalysis;
 181   delete fCounter;
 182   delete fCounterU;
 183   for(Int_t i=0; i<4; i++) delete fMultEstimatorAvg[i];
 184   for(Int_t i=0; i<5; i++){
 185     delete fHistMassPtImpPar[i];
 186   }
 187 }
 188
 189 //_________________________________________________________________
 190 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 191   // set invariant mass limits
 192   if(uplimit>lowlimit){
 193     fLowmasslimit = lowlimit;
 194     fUpmasslimit = uplimit;
 195   }else{
 196     AliError("Wrong mass limits: upper value should be larger than lower one");
 197   }
 198 }
 199 //_________________________________________________________________
 200 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
 201   // set invariant mass limits
 202   Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
 203   SetMassLimits(mass-range,mass+range);
 204 }
 205 //________________________________________________________________________
 206 void AliAnalysisTaskSEDvsMultiplicity::Init(){
 207   //
 208   // Initialization
 209   //
 210   printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
 211
 212   fListCuts=new TList();
 213
 214   if(fPdgMeson==411){
 215      AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
 216     copycut->SetName("AnalysisCutsDplus");
 217     fListCuts->Add(copycut);
 218   }else if(fPdgMeson==421){
 219     AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
 220     copycut->SetName("AnalysisCutsDzero");
 221     fListCuts->Add(copycut);
 222   }else if(fPdgMeson==413){
 223      AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
 224     copycut->SetName("AnalysisCutsDStar");
 225     fListCuts->Add(copycut);
 226   }
 227   PostData(2,fListCuts);
 228
 229   fListProfiles = new TList();
 230   fListProfiles->SetOwner();
 231   TString period[4]={"LHC10b","LHC10c","LHC10d","LHC10e"};
 232   for(Int_t i=0; i<4; i++){
 233     if(fMultEstimatorAvg[i]){
 234       TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
 235       hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
 236       fListProfiles->Add(hprof);
 237     }
 238   }
 239   PostData(4,fListProfiles);
 240
 241   return;
 242 }
 243
 244 //________________________________________________________________________
 245 void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
 246 {
 247   // Create the output container
 248   //
 249   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
 250
 251   // Several histograms are more conveniently managed in a TList
 252   fOutput = new TList();
 253   fOutput->SetOwner();
 254   fOutput->SetName("OutputHistos");
 255
 256   fHistNtrUnCorrEvSel = new TH1F("hNtrUnCorrEvSel","Uncorrected tracklets multiplicity for selected events; Tracklets ; Entries",200,-0.5,199.5);
 257   fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Corrected tracklets multiplicity for selected events; Tracklets ; Entries",200,-0.5,199.5);
 258   fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",200,-0.5,199.5);// Total multiplicity
 259   fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",200,-0.5,199.5); //
 260   fHistNtrEta16vsNtrEta1 = new TH2F("hNtrEta16vsNtrEta1","Uncorrected Eta1.6 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<1.6",200,-0.5,199.5,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram
 261   fHistNtrCorrEta1vsNtrRawEta1 = new TH2F("hNtrCorrEta1vsNtrRawEta1","Corrected Eta1 vs Eta1.0; Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",200,-0.5,199.5,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram
 262   fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,-0.5,199.5); //
 263   fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,-0.5,199.5); //
 264
 265   fHistNtrVsNchMC = new TH2F("hNtrVsNchMC","Ntracklet vs NchMC; Nch;N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
 266   fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC","Ntracklet vs Nch; Nch;N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
 267
 268   fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch (Primary);N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
 269   fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch(Primary) ;N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
 270
 271   fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
 272   fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
 273
 274   fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",200,-0.5,199.5);
 275
 276   fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary = new TH3F("fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary", "MC: Nch (Physical Primary) vs Nch (Primary) vs Nch (Generated); Nch (Generated); Nch (Primary); Nch (Physical Primary)",200,-0.5,199.5,200,-0.5,199.5,200,-0.5,199.5);
 277
 278   fHistNtrUnCorrEvSel->Sumw2();
 279   fHistNtrCorrEvSel->Sumw2();
 280   fHistNtrCorrEvWithCand->Sumw2();
 281   fHistNtrCorrEvWithD->Sumw2();
 282   fHistGenPrimaryParticlesInelGt0->Sumw2();
 283   fOutput->Add(fHistNtrUnCorrEvSel);
 284   fOutput->Add(fHistNtrCorrEvSel);
 285   fOutput->Add(fHistNtrCorrEvWithCand);
 286   fOutput->Add(fHistNtrCorrEvWithD);
 287   fOutput->Add(fHistNtrEta16vsNtrEta1);
 288   fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1);
 289   fOutput->Add(fHistNtrVsZvtx);
 290   fOutput->Add(fHistNtrCorrVsZvtx);
 291
 292   fOutput->Add(fHistNtrVsNchMC);
 293   fOutput->Add(fHistNtrCorrVsNchMC);
 294   fOutput->Add(fHistNtrVsNchMCPrimary);
 295   fOutput->Add(fHistNtrCorrVsNchMCPrimary);
 296   fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
 297   fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
 298   fOutput->Add(fHistGenPrimaryParticlesInelGt0);
 299   fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
 300
 301
 302   fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
 303   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
 304   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
 305   fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
 306   fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
 307   fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
 308   fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
 309   fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
 310   fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
 311   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
 312   fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
 313   fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
 314   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 315   fHistNEvents->Sumw2();
 316   fHistNEvents->SetMinimum(0);
 317   fOutput->Add(fHistNEvents);
 318
 319   fPtVsMassVsMult=new TH3F("hPtVsMassvsMult", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 320
 321   fPtVsMassVsMultNoPid=new TH3F("hPtVsMassvsMultNoPid", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 322
 323   fPtVsMassVsMultUncorr=new TH3F("hPtVsMassvsMultUncorr", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 324
 325   fPtVsMassVsMultPart=new TH3F("hPtVsMassvsMultPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 326
 327   fPtVsMassVsMultAntiPart=new TH3F("hPtVsMassvsMultAntiPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 328
 329   fOutput->Add(fPtVsMassVsMult);
 330   fOutput->Add(fPtVsMassVsMultUncorr);
 331   fOutput->Add(fPtVsMassVsMultNoPid);
 332   fOutput->Add(fPtVsMassVsMultPart);
 333   fOutput->Add(fPtVsMassVsMultAntiPart);
 334
 335   if(fDoImpPar) CreateImpactParameterHistos();
 336
 337   fCounter = new AliNormalizationCounter("NormCounterCorrMult");
 338   fCounter->SetStudyMultiplicity(kTRUE,1.);
 339   fCounter->Init();
 340
 341   fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
 342   fCounterU->SetStudyMultiplicity(kTRUE,1.);
 343   fCounterU->Init();
 344
 345   fOutputCounters = new TList();
 346   fOutputCounters->SetOwner();
 347   fOutputCounters->SetName("OutputCounters");
 348   fOutputCounters->Add(fCounter);
 349   fOutputCounters->Add(fCounterU);
 350
 351   PostData(1,fOutput);
 352   PostData(2,fListCuts);
 353   PostData(3,fOutputCounters);
 354   PostData(4,fListProfiles);
 355
 356
 357
 358   return;
 359 }
 360
 361
 362
 363 //________________________________________________________________________
 364 void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
 365 {
 366   // Execute analysis for current event:
 367   // heavy flavor candidates association to MC truth
 368
 369   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 370
 371   //  AliAODTracklets* tracklets = aod->GetTracklets();
 372   //Int_t ntracklets = tracklets->GetNumberOfTracklets();
 373
 374
 375   TClonesArray *arrayCand = 0;
 376   TString arrayName="";
 377   UInt_t pdgDau[3];
 378   Int_t nDau=0;
 379   Int_t selbit=0;
 380   if(fPdgMeson==411){
 381     arrayName="Charm3Prong";
 382     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
 383     nDau=3;
 384     selbit=AliRDHFCuts::kDplusCuts;
 385   }else if(fPdgMeson==421){
 386     arrayName="D0toKpi";
 387     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
 388     nDau=2;
 389     selbit=AliRDHFCuts::kD0toKpiCuts;
 390   }else if(fPdgMeson==413){
 391     arrayName="Dstar";
 392     pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=211;
 393     nDau=3;
 394     selbit=AliRDHFCuts::kDstarCuts;
 395   }
 396
 397   if(!aod && AODEvent() && IsStandardAOD()) {
 398     // In case there is an AOD handler writing a standard AOD, use the AOD
 399     // event in memory rather than the input (ESD) event.
 400     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 401      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 402      // have to taken from the AOD event hold by the AliAODExtension
 403     AliAODHandler* aodHandler = (AliAODHandler*)
 404       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 405     if(aodHandler->GetExtensions()) {
 406       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 407       AliAODEvent *aodFromExt = ext->GetAOD();
 408       arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
 409     }
 410   } else if(aod) {
 411     arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
 412   }
 413
 414   if(!aod || !arrayCand) {
 415     printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
 416     return;
 417   }
 418
 419
 420
 421   // fix for temporary bug in ESDfilter
 422   // the AODs with null vertex pointer didn't pass the PhysSel
 423   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 424
 425   Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
 426   Int_t countTr=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
 427
 428   fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1);
 429   fHistNEvents->Fill(0); // count event
 430
 431   Double_t countTreta1corr=countTreta1;
 432   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 433   if(vtx1){
 434     if(vtx1->GetNContributors()>0){
 435       fHistNEvents->Fill(1);
 436       TProfile* estimatorAvg = GetEstimatorHistogram(aod);
 437       if(estimatorAvg){
 438         countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
 439       }
 440     }
 441   }
 442
 443   fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countTreta1corr);
 444
 445   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 446
 447   if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
 448   if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
 449   if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
 450   if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
 451
 452
 453   if(!isEvSel)return;
 454   fHistNtrEta16vsNtrEta1->Fill(countTreta1,countTr);
 455   fHistNtrCorrEta1vsNtrRawEta1->Fill(countTreta1,countTreta1corr);
 456   if(vtx1){
 457     fHistNtrVsZvtx->Fill(vtx1->GetZ(),countTreta1);
 458     fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countTreta1corr);
 459   }
 460
 461   TClonesArray *arrayMC=0;
 462   AliAODMCHeader *mcHeader=0;
 463
 464   // load MC particles
 465   if(fReadMC){
 466
 467     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 468     if(!arrayMC) {
 469       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
 470       return;
 471     }
 472     // load MC header
 473     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 474     if(!mcHeader) {
 475       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
 476       return;
 477      }
 478
 479
 480     Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
 481     Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
 482     Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
 483     if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
 484       fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary);
 485     }
 486     fHistNtrVsNchMC->Fill(nChargedMC,countTreta1);
 487     fHistNtrCorrVsNchMC->Fill(nChargedMC,countTreta1corr);
 488
 489     fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1);
 490     fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1corr);
 491
 492     fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1);
 493     fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1corr);
 494
 495     fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary);
 496   }
 497
 498   Int_t nCand = arrayCand->GetEntriesFast();
 499   Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
 500   Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 501   Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
 502   Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 503
 504   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 505     AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
 506     fHistNEvents->Fill(7);
 507     if(fUseBit && !d->HasSelectionBit(selbit)){
 508       fHistNEvents->Fill(8);
 509       continue;
 510     }
 511
 512     Double_t ptCand = d->Pt();
 513     Double_t rapid=d->Y(fPdgMeson);
 514     Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 515     if(!isFidAcc) continue;
 516     Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
 517     Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
 518     if(passTopolCuts==0) continue;
 519     nSelectedNoPID++;
 520     fHistNEvents->Fill(9);
 521     if(passAllCuts){
 522       nSelectedPID++;
 523       fHistNEvents->Fill(10);
 524     }
 525     Bool_t isPrimary=kTRUE;
 526     Int_t labD=-1;
 527     Double_t trueImpParXY=9999.;
 528     Double_t impparXY=d->ImpParXY()*10000.;
 529     Double_t dlen=0.1; //FIXME
 530     Double_t mass[2];
 531     if(fPdgMeson==411){
 532       mass[0]=d->InvMass(nDau,pdgDau);
 533       mass[1]=-1.;
 534       if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
 535     }else if(fPdgMeson==421){
 536       UInt_t pdgdaughtersD0[2]={211,321};//pi,K
 537       UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
 538       mass[0]=d->InvMass(2,pdgdaughtersD0);
 539       mass[1]=d->InvMass(2,pdgdaughtersD0bar);
 540       if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
 541     }else if(fPdgMeson==413){
 542       // FIXME
 543       AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
 544       mass[0]=temp->DeltaInvMass();
 545       mass[1]=-1.;
 546       if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.001) nSelectedInMassPeak++; //1 MeV for now... FIXME
 547     }
 548     for(Int_t iHyp=0; iHyp<2; iHyp++){
 549       if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
 550       Double_t invMass=mass[iHyp];
 551       Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
 552
 553       if(fReadMC){
 554
 555         labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
 556         Bool_t fillHisto=fDoImpPar;
 557         if(labD>=0){
 558           AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
 559           Int_t code=partD->GetPdgCode();
 560           if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
 561           if(code<0 && iHyp==0) fillHisto=kFALSE;
 562           if(code>0 && iHyp==1) fillHisto=kFALSE;
 563           if(!isPrimary){
 564             if(fPdgMeson==411){
 565               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
 566             }else if(fPdgMeson==421){
 567               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
 568             }else if(fPdgMeson==413){
 569               trueImpParXY=0.; /// FIXME
 570             }
 571             Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
 572             if(fillHisto && passAllCuts){
 573               fHistMassPtImpPar[2]->Fill(arrayForSparse);
 574               fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
 575             }
 576           }else{
 577             if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
 578           }
 579         }else{
 580           if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
 581         }
 582         if(fPdgMeson==421){
 583           if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
 584           if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
 585         }
 586       }
 587
 588       if(fPdgMeson==421){
 589         if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
 590         if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
 591       }
 592
 593       fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
 594       if(passAllCuts){
 595         fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);
 596         fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
 597         // Add separation between part antipart
 598         if(fPdgMeson==411){
 599           if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 600           else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 601         }else if(fPdgMeson==421){
 602           if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 603           if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 604         }else if(fPdgMeson==413){
 605           // FIXME ADD Dstar!!!!!!!!
 606         }
 607
 608         if(fDoImpPar){
 609           fHistMassPtImpPar[0]->Fill(arrayForSparse);
 610         }
 611
 612       }
 613
 614     }
 615   }
 616   fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
 617   fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
 618  fHistNtrUnCorrEvSel->Fill(countTreta1);
 619   fHistNtrCorrEvSel->Fill(countTreta1corr);
 620   if(nSelectedPID>0) fHistNtrCorrEvWithCand->Fill(countTreta1corr);
 621   if(nSelectedInMassPeak>0) fHistNtrCorrEvWithD->Fill(countTreta1corr);
 622
 623   PostData(1,fOutput);
 624   PostData(2,fListCuts);
 625   PostData(3,fOutputCounters);
 626
 627   return;
 628 }
 629 //________________________________________________________________________
 630 void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
 631   // Histos for impact paramter study
 632   // mass . pt , impact parameter , decay length , multiplicity
 633
 634   Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
 635   Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
 636   Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
 637
 638   fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
 639                                         "Mass vs. pt vs.imppar - All",
 640                                         5,nbins,xmin,xmax);
 641   fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
 642                                         "Mass vs. pt vs.imppar - promptD",
 643                                         5,nbins,xmin,xmax);
 644   fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
 645                                         "Mass vs. pt vs.imppar - DfromB",
 646                                         5,nbins,xmin,xmax);
 647   fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
 648                                         "Mass vs. pt vs.true imppar -DfromB",
 649                                         5,nbins,xmin,xmax);
 650   fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
 651                                         "Mass vs. pt vs.imppar - backgr.",
 652                                         5,nbins,xmin,xmax);
 653   for(Int_t i=0; i<5;i++){
 654     fOutput->Add(fHistMassPtImpPar[i]);
 655   }
 656 }
 657
 658 //________________________________________________________________________
 659 void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
 660 {
 661   // Terminate analysis
 662   //
 663   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
 664
 665   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 666   if (!fOutput) {
 667     printf("ERROR: fOutput not available\n");
 668     return;
 669   }
 670
 671   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 672   if(!fHistNEvents){
 673     printf("ERROR: fHistNEvents not available\n");
 674     return;
 675   }
 676   printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
 677
 678   return;
 679 }
 680 //_________________________________________________________________________________________________
 681 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
 682   //
 683   // checking whether the mother of the particles come from a charm or a bottom quark
 684   //
 685
 686   Int_t pdgGranma = 0;
 687   Int_t mother = 0;
 688   mother = mcPartCandidate->GetMother();
 689   Int_t istep = 0;
 690   Int_t abspdgGranma =0;
 691   Bool_t isFromB=kFALSE;
 692   Bool_t isQuarkFound=kFALSE;
 693   while (mother >0 ){
 694     istep++;
 695     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 696     if (mcGranma){
 697       pdgGranma = mcGranma->GetPdgCode();
 698       abspdgGranma = TMath::Abs(pdgGranma);
 699       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 700         isFromB=kTRUE;
 701       }
 702       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 703       mother = mcGranma->GetMother();
 704     }else{
 705       AliError("Failed casting the mother particle!");
 706       break;
 707     }
 708   }
 709
 710   if(isFromB) return 5;
 711   else return 4;
 712 }
 713
 714
 715
 716 //____________________________________________________________________________
 717 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
 718   // Get Estimator Histogram from period event->GetRunNumber();
 719   //
 720   // If you select SPD tracklets in |eta|<1 you should use type == 1
 721   //
 722
 723   Int_t runNo  = event->GetRunNumber();
 724   Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
 725   if(runNo>114930 && runNo<117223) period = 0;
 726   if(runNo>119158 && runNo<120830) period = 1;
 727   if(runNo>122373 && runNo<126438) period = 2;
 728   if(runNo>127711 && runNo<130841) period = 3;
 729   if(period<0 || period>3) return 0;
 730
 731   return fMultEstimatorAvg[period];
 732 }