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