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