PWGHF/vertexingHF/AliAnalysisTaskSEDStarSpectra.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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  * appeuear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //
  19 //
  20 //                  Base class for DStar Analysis
  21 //
  22 //
  23 //  The D* spectra study is done in pt bins:
  24 //  [0,0.5] [0.5,1] [1,2] [2,3] [3,4] [4,5] [5,6] [6,7] [7,8],
  25 //  [8,10],[10,12], [12,16], [16,20] and [20,24]
  26 //
  27 //  Cuts arew centralized in AliRDHFCutsDStartoKpipi
  28 //  Side Band and like sign background are implemented in the macro
  29 //
  30 //-----------------------------------------------------------------------
  31 //
  32 //                         Author A.Grelli
  33 //              ERC-QGP Utrecht University - a.grelli@uu.nl,
  34 //                         Author Y.Wang
  35 //        University of Heidelberg - yifei@physi.uni-heidelberg.de
  36 //                         Author C.Ivan
  37 //             ERC-QGP Utrecht University - c.ivan@uu.nl,
  38 //
  39 //-----------------------------------------------------------------------
  40
  41 #include <TSystem.h>
  42 #include <TParticle.h>
  43 #include <TH1I.h>
  44 #include "TROOT.h"
  45 #include <TDatabasePDG.h>
  46 #include <AliAnalysisDataSlot.h>
  47 #include <AliAnalysisDataContainer.h>
  48 #include "AliRDHFCutsDStartoKpipi.h"
  49 #include "AliStack.h"
  50 #include "AliMCEvent.h"
  51 #include "AliAnalysisManager.h"
  52 #include "AliAODMCHeader.h"
  53 #include "AliAODHandler.h"
  54 #include "AliLog.h"
  55 #include "AliAODVertex.h"
  56 #include "AliAODRecoDecay.h"
  57 #include "AliAODRecoDecayHF.h"
  58 #include "AliAODRecoCascadeHF.h"
  59 #include "AliAODRecoDecayHF2Prong.h"
  60 #include "AliAnalysisVertexingHF.h"
  61 #include "AliESDtrack.h"
  62 #include "AliAODMCParticle.h"
  63 #include "AliNormalizationCounter.h"
  64 #include "AliAODEvent.h"
  65 #include "AliAnalysisTaskSEDStarSpectra.h"
  66
  67 ClassImp(AliAnalysisTaskSEDStarSpectra)
  68
  69 //__________________________________________________________________________
  70 AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra():
  71   AliAnalysisTaskSE(),
  72   fEvents(0),
  73   fAnalysis(0),
  74   fD0Window(0),
  75   fPeakWindow(0),
  76   fUseMCInfo(kFALSE),
  77   fDoSearch(kFALSE),
  78   fOutput(0),
  79   fOutputAll(0),
  80   fOutputPID(0),
  81   fNSigma(3),
  82   fCuts(0),
  83   fCEvents(0),
  84   fTrueDiff2(0),
  85   fDeltaMassD1(0),
  86   fCounter(0),
  87   fDoImpParDstar(kFALSE),
  88   fNImpParBins(400),
  89   fLowerImpPar(-2000.),
  90   fHigherImpPar(2000.)
  91 {
  92   //
  93   // Default ctor
  94   //
  95   for(Int_t i=0;i<5;i++) fHistMassPtImpParTCDs[i]=0;
  96 }
  97 //___________________________________________________________________________
  98 AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
  99   AliAnalysisTaskSE(name),
 100   fEvents(0),
 101   fAnalysis(0),
 102   fD0Window(0),
 103   fPeakWindow(0),
 104   fUseMCInfo(kFALSE),
 105   fDoSearch(kFALSE),
 106   fOutput(0),
 107   fOutputAll(0),
 108   fOutputPID(0),
 109   fNSigma(3),
 110   fCuts(0),
 111   fCEvents(0),
 112   fTrueDiff2(0),
 113   fDeltaMassD1(0),
 114   fCounter(0),
 115   fDoImpParDstar(kFALSE),
 116   fNImpParBins(400),
 117   fLowerImpPar(-2000.),
 118   fHigherImpPar(2000.)
 119 {
 120   //
 121   // Constructor. Initialization of Inputs and Outputs
 122   //
 123   Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");
 124
 125   fCuts=cuts;
 126
 127   DefineOutput(1,TList::Class());  //conters
 128   DefineOutput(2,TList::Class());  //All Entries output
 129   DefineOutput(3,TList::Class());  //3sigma PID output
 130   DefineOutput(4,AliRDHFCutsDStartoKpipi::Class());   //My private output
 131   DefineOutput(5,AliNormalizationCounter::Class());   // normalization
 132 }
 133
 134 //___________________________________________________________________________
 135 AliAnalysisTaskSEDStarSpectra::~AliAnalysisTaskSEDStarSpectra() {
 136   //
 137   // destructor
 138   //
 139   Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");
 140
 141   delete fOutput;
 142   delete fOutputAll;
 143   delete fOutputPID;
 144   delete fCuts;
 145   delete fCEvents;
 146   delete fDeltaMassD1;
 147   for(Int_t i=0; i<5; i++){
 148     delete fHistMassPtImpParTCDs[i];
 149   }
 150 }
 151 //_________________________________________________
 152 void AliAnalysisTaskSEDStarSpectra::Init(){
 153   //
 154   // Initialization
 155   //
 156
 157   if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");
 158    AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
 159   // Post the data
 160   PostData(4,copyfCuts);
 161
 162   return;
 163 }
 164
 165 //_________________________________________________
 166 void AliAnalysisTaskSEDStarSpectra::UserExec(Option_t *)
 167 {
 168   // user exec
 169   if (!fInputEvent) {
 170     Error("UserExec","NO EVENT FOUND!");
 171     return;
 172   }
 173
 174   fEvents++;
 175
 176   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 177   TClonesArray *arrayDStartoD0pi=0;
 178
 179   fCEvents->Fill(1);
 180
 181   if(!aodEvent && AODEvent() && IsStandardAOD()) {
 182     // In case there is an AOD handler writing a standard AOD, use the AOD
 183     // event in memory rather than the input (ESD) event.
 184     aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 185     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 186     // have to taken from the AOD event hold by the AliAODExtension
 187     AliAODHandler* aodHandler = (AliAODHandler*)
 188       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 189     if(aodHandler->GetExtensions()) {
 190       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 191       AliAODEvent *aodFromExt = ext->GetAOD();
 192       arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 193     }
 194   } else {
 195     arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 196   }
 197
 198   // fix for temporary bug in ESDfilter
 199   // the AODs with null vertex pointer didn't pass the PhysSel
 200   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 201   fCEvents->Fill(2);
 202
 203   fCounter->StoreEvent(aodEvent,fCuts,fUseMCInfo);
 204
 205   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
 206   TString trigclass=aodEvent->GetFiredTriggerClasses();
 207   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fCEvents->Fill(5);
 208
 209   if(!fCuts->IsEventSelected(aodEvent)) {
 210     if(fCuts->GetWhyRejection()==6) // rejected for Z vertex
 211       fCEvents->Fill(6);
 212     return;
 213   }
 214
 215   Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
 216   if(!isEvSel) return;
 217
 218   // Load the event
 219   //  AliInfo(Form("Event %d",fEvents));
 220   //if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
 221
 222   // counters for efficiencies
 223   Int_t icountReco = 0;
 224
 225   //D* and D0 prongs needed to MatchToMC method
 226   Int_t pdgDgDStartoD0pi[2]={421,211};
 227   Int_t pdgDgD0toKpi[2]={321,211};
 228
 229   // AOD primary vertex
 230   AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 231   if(!vtx1) return;
 232   if(vtx1->GetNContributors()<1) return;
 233
 234   if (!arrayDStartoD0pi){
 235     AliInfo("Could not find array of HF vertices, skipping the event");
 236     return;
 237   }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 238
 239   Int_t nSelectedAna =0;
 240   Int_t nSelectedProd =0;
 241
 242   // loop over the tracks to search for candidates soft pion
 243   for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
 244
 245     // D* candidates and D0 from D*
 246     AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
 247     if(!dstarD0pi->GetSecondaryVtx()) continue;
 248     AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 249     if (!theD0particle) continue;
 250
 251     Int_t isDStar = 0;
 252     TClonesArray *mcArray = 0;
 253     AliAODMCHeader *mcHeader=0;
 254
 255     Bool_t isPrimary=kTRUE;
 256     Float_t truePt=0.;
 257     Float_t xDecay=0.;
 258     Float_t yDecay=0.;
 259     Float_t zDecay=0.;
 260     Float_t pdgCode=-2;
 261     Float_t trueImpParXY=0.;
 262
 263     // mc analysis
 264     if(fUseMCInfo){
 265     //MC array need for maching
 266       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 267       if (!mcArray) {
 268         AliError("Could not find Monte-Carlo in AOD");
 269         return;
 270       }
 271       // load MC header
 272       mcHeader =  (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 273       if(!mcHeader) {
 274         printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
 275         return;
 276       }
 277       // find associated MC particle for D* ->D0toKpi
 278       Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
 279       if(mcLabel>=0){
 280
 281         AliAODMCParticle *partDSt = (AliAODMCParticle*)mcArray->At(mcLabel);
 282         Int_t checkOrigin = CheckOrigin(mcArray,partDSt);
 283         if(checkOrigin==5) isPrimary=kFALSE;
 284         AliAODMCParticle *dg0 = (AliAODMCParticle*)mcArray->At(partDSt->GetDaughter(0));
 285         AliAODMCParticle *dg01 = (AliAODMCParticle*)mcArray->At(dg0->GetDaughter(0));
 286         truePt=dg0->Pt();
 287         xDecay=dg01->Xv();
 288         yDecay=dg01->Yv();
 289         zDecay=dg01->Zv();
 290         pdgCode=TMath::Abs(partDSt->GetPdgCode());
 291         if(!isPrimary){
 292           trueImpParXY=GetTrueImpactParameterD0(mcHeader,mcArray,dg0)*1000.;
 293         }
 294         isDStar = 1;
 295       }else{
 296         pdgCode=-1;
 297       }
 298     }
 299
 300     Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
 301
 302     // quality selction on tracks and region of interest
 303     Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
 304     if(!isTkSelected) continue;
 305
 306     if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
 307
 308
 309     //histos for impact par studies - D0!!!
 310     Double_t ptCand = dstarD0pi->Get2Prong()->Pt();
 311     Double_t invMass=dstarD0pi->InvMassD0();
 312     Double_t impparXY=dstarD0pi->Get2Prong()->ImpParXY()*10000.;
 313
 314     Double_t arrayForSparse[3]={invMass,ptCand,impparXY};
 315     Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};
 316
 317   // set the D0 search window bin by bin - useful to calculate side band bkg
 318     if (ptbin==0){
 319       if(fAnalysis==1){
 320         fD0Window=0.035;
 321         fPeakWindow=0.03;
 322       }else{
 323         fD0Window=0.020;
 324         fPeakWindow=0.0018;
 325       }
 326     }
 327     if (ptbin==1){
 328       if(fAnalysis==1){
 329         fD0Window=0.035;
 330         fPeakWindow=0.03;
 331       }else{
 332         fD0Window=0.020;
 333         fPeakWindow=0.0018;
 334       }
 335     }
 336     if (ptbin==2){
 337       if(fAnalysis==1){
 338         fD0Window=0.035;
 339         fPeakWindow=0.03;
 340       }else{
 341         fD0Window=0.020;
 342         fPeakWindow=0.0018;
 343       }
 344     }
 345     if (ptbin==3){
 346       if(fAnalysis==1){
 347         fD0Window=0.035;
 348         fPeakWindow=0.03;
 349       }else{
 350         fD0Window=0.022;
 351         fPeakWindow=0.0016;
 352       }
 353     }
 354     if (ptbin==4){
 355       if(fAnalysis==1){
 356         fD0Window=0.035;
 357         fPeakWindow=0.03;
 358       }else{
 359         fD0Window=0.026;
 360         fPeakWindow=0.0014;
 361       }
 362     }
 363     if (ptbin==5){
 364       if(fAnalysis==1){
 365         fD0Window=0.045;
 366         fPeakWindow=0.03;
 367       }else{
 368         fD0Window=0.026;
 369         fPeakWindow=0.0014;
 370       }
 371     }
 372    if (ptbin==6){
 373       if(fAnalysis==1){
 374         fD0Window=0.045;
 375         fPeakWindow=0.03;
 376       }else{
 377         fD0Window=0.026;
 378         fPeakWindow=0.006;
 379       }
 380     }
 381     if (ptbin==7){
 382       if(fAnalysis==1){
 383         fD0Window=0.055;
 384         fPeakWindow=0.03;
 385       }else{
 386         fD0Window=0.026;
 387         fPeakWindow=0.006;
 388       }
 389     }
 390     if (ptbin>7){
 391       if(fAnalysis==1){
 392         fD0Window=0.074;
 393         fPeakWindow=0.03;
 394       }else{
 395         fD0Window=0.026;
 396         fPeakWindow=0.006;
 397       }
 398     }
 399
 400     nSelectedProd++;
 401     nSelectedAna++;
 402
 403     // check that we are close to signal in the DeltaM - here to save time for PbPb
 404     Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 405     Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
 406     Double_t invmassDelta = dstarD0pi->DeltaInvMass();
 407
 408     if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>fPeakWindow) continue;
 409
 410
 411     Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
 412
 413     // after cuts
 414     if(fDoImpParDstar && isSelected){
 415       fHistMassPtImpParTCDs[0]->Fill(arrayForSparse);
 416       if(isPrimary) fHistMassPtImpParTCDs[1]->Fill(arrayForSparse);
 417       else{
 418         fHistMassPtImpParTCDs[2]->Fill(arrayForSparse);
 419         fHistMassPtImpParTCDs[3]->Fill(arrayForSparseTrue);
 420       }
 421     }
 422
 423     // fill PID
 424     FillSpectrum(dstarD0pi,isDStar,fCuts,isSelected,fOutputPID);
 425     SideBandBackground(dstarD0pi,fCuts,isSelected, fOutputPID);
 426     //WrongSignForDStar(dstarD0pi,fCuts,fOutputPID);
 427
 428     //swich off the PID selection
 429     fCuts->SetUsePID(kFALSE);
 430     Int_t isSelectedNoPID=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
 431     fCuts->SetUsePID(kTRUE);
 432
 433     FillSpectrum(dstarD0pi,isDStar,fCuts,isSelectedNoPID,fOutputPID);
 434     SideBandBackground(dstarD0pi,fCuts,isSelectedNoPID, fOutputPID);
 435
 436     // rare D search ------
 437     if(fDoSearch){
 438       TLorentzVector lorentzTrack1(0,0,0,0); // lorentz 4 vector
 439       TLorentzVector lorentzTrack2(0,0,0,0); // lorentz 4 vector
 440
 441       for (Int_t i=0; i<aodEvent->GetNTracks(); i++){
 442
 443         AliAODTrack* aodTrack = aodEvent->GetTrack(i);
 444
 445         if(dstarD0pi->Charge() == aodTrack->Charge()) continue;
 446         if((!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
 447         if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>0.02) continue;
 448
 449         //build the D1 mass
 450         Double_t mass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 451
 452         lorentzTrack1.SetPxPyPzE( dstarD0pi->Px(),dstarD0pi->Py(), dstarD0pi->Pz(), dstarD0pi->E(413) );
 453         lorentzTrack2.SetPxPyPzE( aodTrack->Px(),aodTrack->Py(), aodTrack->Pz(),aodTrack->E(mass) );
 454
 455         //D1 mass
 456         Double_t d1mass = ((lorentzTrack1+lorentzTrack2).M());
 457         //mass difference - at 0.4117 and 0.4566
 458         fDeltaMassD1->Fill(d1mass-dstarD0pi->InvMassDstarKpipi());
 459       }
 460     }
 461
 462     if(isDStar == 1) {
 463       fTrueDiff2->Fill(dstarD0pi->Pt(),dstarD0pi->DeltaInvMass());
 464     }
 465
 466   }
 467
 468   fCounter->StoreCandidates(aodEvent,nSelectedProd,kTRUE);
 469   fCounter->StoreCandidates(aodEvent,nSelectedAna,kFALSE);
 470
 471   AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
 472
 473   PostData(1,fOutput);
 474   PostData(2,fOutputAll);
 475   PostData(3,fOutputPID);
 476   PostData(5,fCounter);
 477
 478 }
 479 //________________________________________ terminate ___________________________
 480 void AliAnalysisTaskSEDStarSpectra::Terminate(Option_t*)
 481 {
 482   // The Terminate() function is the last function to be called during
 483   // a query. It always runs on the client, it can be used to present
 484   // the results graphically or save the results to file.
 485
 486   //Info("Terminate","");
 487   AliAnalysisTaskSE::Terminate();
 488
 489   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 490   if (!fOutput) {
 491     printf("ERROR: fOutput not available\n");
 492     return;
 493   }
 494
 495   fCEvents        = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
 496   fDeltaMassD1     = dynamic_cast<TH1F*>(fOutput->FindObject("fDeltaMassD1"));
 497   fTrueDiff2      = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
 498
 499   fOutputAll = dynamic_cast<TList*> (GetOutputData(1));
 500   if (!fOutputAll) {
 501     printf("ERROR: fOutputAll not available\n");
 502     return;
 503   }
 504   fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
 505   if (!fOutputPID) {
 506     printf("ERROR: fOutputPID not available\n");
 507     return;
 508   }
 509
 510
 511   return;
 512 }
 513 //___________________________________________________________________________
 514 void AliAnalysisTaskSEDStarSpectra::UserCreateOutputObjects() {
 515  // output
 516   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 517
 518   //slot #1
 519   //OpenFile(1);
 520   fOutput = new TList();
 521   fOutput->SetOwner();
 522   fOutput->SetName("chist0");
 523
 524   fOutputAll = new TList();
 525   fOutputAll->SetOwner();
 526   fOutputAll->SetName("listAll");
 527
 528   fOutputPID = new TList();
 529   fOutputPID->SetOwner();
 530   fOutputPID->SetName("listPID");
 531
 532   // define histograms
 533   DefineHistograms();
 534
 535  //Counter for Normalization
 536  fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
 537  fCounter->Init();
 538
 539  if(fDoImpParDstar) CreateImpactParameterHistos();
 540
 541   PostData(1,fOutput);
 542   PostData(2,fOutputAll);
 543   PostData(3,fOutputPID);
 544
 545   return;
 546 }
 547 //___________________________________ hiostograms _______________________________________
 548 void  AliAnalysisTaskSEDStarSpectra::DefineHistograms(){
 549   // Create histograms
 550
 551   fCEvents = new TH1F("fCEvents","conter",11,0,11);
 552   fCEvents->SetStats(kTRUE);
 553   fCEvents->GetXaxis()->SetTitle("1");
 554   fCEvents->GetYaxis()->SetTitle("counts");
 555   fOutput->Add(fCEvents);
 556
 557   fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
 558   fOutput->Add(fTrueDiff2);
 559
 560   fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);
 561   fOutput->Add(fDeltaMassD1);
 562
 563   const Int_t nhist=14;
 564   TString nameMass=" ", nameSgn=" ", nameBkg=" ";
 565
 566   for(Int_t i=-2;i<nhist;i++){
 567     nameMass="histDeltaMass_";
 568     nameMass+=i+1;
 569     nameSgn="histDeltaSgn_";
 570     nameSgn+=i+1;
 571     nameBkg="histDeltaBkg_";
 572     nameBkg+=i+1;
 573
 574     if (i==-2) {
 575       nameMass="histDeltaMass";
 576       nameSgn="histDeltaSgn";
 577       nameBkg="histDeltaBkg";
 578     }
 579
 580     TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
 581     TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
 582     TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
 583
 584     nameMass="histD0Mass_";
 585     nameMass+=i+1;
 586     nameSgn="histD0Sgn_";
 587     nameSgn+=i+1;
 588     nameBkg="histD0Bkg_";
 589     nameBkg+=i+1;
 590
 591     if (i==-2) {
 592       nameMass="histD0Mass";
 593       nameSgn="histD0Sgn";
 594       nameBkg="histD0Bkg";
 595     }
 596
 597     TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
 598     TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
 599     TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
 600
 601     nameMass="histDstarMass_";
 602     nameMass+=i+1;
 603     nameSgn="histDstarSgn_";
 604     nameSgn+=i+1;
 605     nameBkg="histDstarBkg_";
 606     nameBkg+=i+1;
 607
 608     if (i==-2) {
 609       nameMass="histDstarMass";
 610       nameSgn="histDstarSgn";
 611       nameBkg="histDstarBkg";
 612     }
 613
 614     TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
 615     TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
 616     TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
 617
 618     nameMass="histSideBandMass_";
 619     nameMass+=i+1;
 620     if (i==-2) {
 621       nameMass="histSideBandMass";
 622     }
 623
 624     TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
 625
 626     nameMass="histWrongSignMass_";
 627     nameMass+=i+1;
 628     if (i==-2) {
 629       nameMass="histWrongSignMass";
 630     }
 631
 632     TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
 633
 634
 635     spectrumMass->Sumw2();
 636     spectrumSgn->Sumw2();
 637     spectrumBkg->Sumw2();
 638
 639     spectrumMass->SetLineColor(6);
 640     spectrumSgn->SetLineColor(2);
 641     spectrumBkg->SetLineColor(4);
 642
 643     spectrumMass->SetMarkerStyle(20);
 644     spectrumSgn->SetMarkerStyle(20);
 645     spectrumBkg->SetMarkerStyle(20);
 646     spectrumMass->SetMarkerSize(0.6);
 647     spectrumSgn->SetMarkerSize(0.6);
 648     spectrumBkg->SetMarkerSize(0.6);
 649     spectrumMass->SetMarkerColor(6);
 650     spectrumSgn->SetMarkerColor(2);
 651     spectrumBkg->SetMarkerColor(4);
 652
 653     spectrumD0Mass->Sumw2();
 654     spectrumD0Sgn->Sumw2();
 655     spectrumD0Bkg->Sumw2();
 656
 657     spectrumD0Mass->SetLineColor(6);
 658     spectrumD0Sgn->SetLineColor(2);
 659     spectrumD0Bkg->SetLineColor(4);
 660
 661     spectrumD0Mass->SetMarkerStyle(20);
 662     spectrumD0Sgn->SetMarkerStyle(20);
 663     spectrumD0Bkg->SetMarkerStyle(20);
 664     spectrumD0Mass->SetMarkerSize(0.6);
 665     spectrumD0Sgn->SetMarkerSize(0.6);
 666     spectrumD0Bkg->SetMarkerSize(0.6);
 667     spectrumD0Mass->SetMarkerColor(6);
 668     spectrumD0Sgn->SetMarkerColor(2);
 669     spectrumD0Bkg->SetMarkerColor(4);
 670
 671     spectrumDstarMass->Sumw2();
 672     spectrumDstarSgn->Sumw2();
 673     spectrumDstarBkg->Sumw2();
 674
 675     spectrumDstarMass->SetLineColor(6);
 676     spectrumDstarSgn->SetLineColor(2);
 677     spectrumDstarBkg->SetLineColor(4);
 678
 679     spectrumDstarMass->SetMarkerStyle(20);
 680     spectrumDstarSgn->SetMarkerStyle(20);
 681     spectrumDstarBkg->SetMarkerStyle(20);
 682     spectrumDstarMass->SetMarkerSize(0.6);
 683     spectrumDstarSgn->SetMarkerSize(0.6);
 684     spectrumDstarBkg->SetMarkerSize(0.6);
 685     spectrumDstarMass->SetMarkerColor(6);
 686     spectrumDstarSgn->SetMarkerColor(2);
 687     spectrumDstarBkg->SetMarkerColor(4);
 688
 689     spectrumSideBandMass->Sumw2();
 690     spectrumSideBandMass->SetLineColor(4);
 691     spectrumSideBandMass->SetMarkerStyle(20);
 692     spectrumSideBandMass->SetMarkerSize(0.6);
 693     spectrumSideBandMass->SetMarkerColor(4);
 694
 695     spectrumWrongSignMass->Sumw2();
 696     spectrumWrongSignMass->SetLineColor(4);
 697     spectrumWrongSignMass->SetMarkerStyle(20);
 698     spectrumWrongSignMass->SetMarkerSize(0.6);
 699     spectrumWrongSignMass->SetMarkerColor(4);
 700
 701     TH1F* allMass = (TH1F*)spectrumMass->Clone();
 702     TH1F* allSgn  = (TH1F*)spectrumSgn->Clone();
 703     TH1F* allBkg  = (TH1F*)spectrumBkg->Clone();
 704
 705     TH1F* pidMass = (TH1F*)spectrumMass->Clone();
 706     TH1F* pidSgn  = (TH1F*)spectrumSgn->Clone();
 707     TH1F* pidBkg  = (TH1F*)spectrumBkg->Clone();
 708
 709     fOutputAll->Add(allMass);
 710     fOutputAll->Add(allSgn);
 711     fOutputAll->Add(allBkg);
 712
 713     fOutputPID->Add(pidMass);
 714     fOutputPID->Add(pidSgn);
 715     fOutputPID->Add(pidBkg);
 716
 717     TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
 718     TH1F* allD0Sgn  = (TH1F*)spectrumD0Sgn->Clone();
 719     TH1F* allD0Bkg  = (TH1F*)spectrumD0Bkg->Clone();
 720
 721     TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();
 722     TH1F* pidD0Sgn  = (TH1F*)spectrumD0Sgn->Clone();
 723     TH1F* pidD0Bkg  = (TH1F*)spectrumD0Bkg->Clone();
 724
 725     fOutputAll->Add(allD0Mass);
 726     fOutputAll->Add(allD0Sgn);
 727     fOutputAll->Add(allD0Bkg);
 728
 729     fOutputPID->Add(pidD0Mass);
 730     fOutputPID->Add(pidD0Sgn);
 731     fOutputPID->Add(pidD0Bkg);
 732
 733     TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
 734     TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
 735     TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
 736
 737     TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();
 738     TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
 739     TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
 740
 741     fOutputAll->Add(allDstarMass);
 742     fOutputAll->Add(allDstarSgn);
 743     fOutputAll->Add(allDstarBkg);
 744
 745     fOutputPID->Add(pidDstarMass);
 746     fOutputPID->Add(pidDstarSgn);
 747     fOutputPID->Add(pidDstarBkg);
 748
 749     TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
 750     TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
 751
 752     fOutputAll->Add(allSideBandMass);
 753     fOutputPID->Add(pidSideBandMass);
 754
 755     TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
 756     TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
 757
 758     fOutputAll->Add(allWrongSignMass);
 759     fOutputPID->Add(pidWrongSignMass);
 760
 761   }
 762
 763   // pt spectra
 764   nameMass="ptMass";
 765   nameSgn="ptSgn";
 766   nameBkg="ptBkg";
 767
 768   TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",200,0,10);
 769   TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
 770   TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
 771
 772   ptspectrumMass->Sumw2();
 773   ptspectrumSgn->Sumw2();
 774   ptspectrumBkg->Sumw2();
 775
 776   ptspectrumMass->SetLineColor(6);
 777   ptspectrumSgn->SetLineColor(2);
 778   ptspectrumBkg->SetLineColor(4);
 779
 780   ptspectrumMass->SetMarkerStyle(20);
 781   ptspectrumSgn->SetMarkerStyle(20);
 782   ptspectrumBkg->SetMarkerStyle(20);
 783   ptspectrumMass->SetMarkerSize(0.6);
 784   ptspectrumSgn->SetMarkerSize(0.6);
 785   ptspectrumBkg->SetMarkerSize(0.6);
 786   ptspectrumMass->SetMarkerColor(6);
 787   ptspectrumSgn->SetMarkerColor(2);
 788   ptspectrumBkg->SetMarkerColor(4);
 789
 790   TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
 791   TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
 792   TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
 793
 794   TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();
 795   TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();
 796   TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();
 797
 798   fOutputAll->Add(ptallMass);
 799   fOutputAll->Add(ptallSgn);
 800   fOutputAll->Add(ptallBkg);
 801
 802   fOutputPID->Add(ptpidMass);
 803   fOutputPID->Add(ptpidSgn);
 804   fOutputPID->Add(ptpidBkg);
 805
 806   // eta spectra
 807   nameMass="etaMass";
 808   nameSgn="etaSgn";
 809   nameBkg="etaBkg";
 810
 811   TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
 812   TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
 813   TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
 814
 815   etaspectrumMass->Sumw2();
 816   etaspectrumSgn->Sumw2();
 817   etaspectrumBkg->Sumw2();
 818
 819   etaspectrumMass->SetLineColor(6);
 820   etaspectrumSgn->SetLineColor(2);
 821   etaspectrumBkg->SetLineColor(4);
 822
 823   etaspectrumMass->SetMarkerStyle(20);
 824   etaspectrumSgn->SetMarkerStyle(20);
 825   etaspectrumBkg->SetMarkerStyle(20);
 826   etaspectrumMass->SetMarkerSize(0.6);
 827   etaspectrumSgn->SetMarkerSize(0.6);
 828   etaspectrumBkg->SetMarkerSize(0.6);
 829   etaspectrumMass->SetMarkerColor(6);
 830   etaspectrumSgn->SetMarkerColor(2);
 831   etaspectrumBkg->SetMarkerColor(4);
 832
 833   TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
 834   TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
 835   TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
 836
 837   TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();
 838   TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();
 839   TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();
 840
 841   fOutputAll->Add(etaallMass);
 842   fOutputAll->Add(etaallSgn);
 843   fOutputAll->Add(etaallBkg);
 844
 845   fOutputPID->Add(etapidMass);
 846   fOutputPID->Add(etapidSgn);
 847   fOutputPID->Add(etapidBkg);
 848
 849   return;
 850 }
 851 //________________________________________________________________________
 852 void AliAnalysisTaskSEDStarSpectra::FillSpectrum(AliAODRecoCascadeHF *part, Int_t isDStar, AliRDHFCutsDStartoKpipi *cuts,Int_t isSel, TList *listout){
 853   //
 854   // Fill histos for D* spectrum
 855   //
 856
 857   if(!isSel) return;
 858
 859   // D0 window
 860   Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 861   Double_t invmassD0   = part->InvMassD0();
 862   if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return;
 863
 864
 865   Int_t ptbin=cuts->PtBin(part->Pt());
 866   Double_t pt = part->Pt();
 867   Double_t eta = part->Eta();
 868
 869   Double_t invmassDelta = part->DeltaInvMass();
 870   Double_t invmassDstar = part->InvMassDstarKpipi();
 871
 872   TString fillthis="";
 873   Bool_t massInRange=kFALSE;
 874
 875   Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
 876
 877   // delta M(Kpipi)-M(Kpi)
 878   if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<fPeakWindow) massInRange=kTRUE;
 879
 880   if(fUseMCInfo) {
 881     if(isDStar==1) {
 882       fillthis="histD0Sgn_";
 883       fillthis+=ptbin;
 884       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 885       fillthis="histD0Sgn";
 886       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 887       fillthis="histDstarSgn_";
 888       fillthis+=ptbin;
 889       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 890       fillthis="histDstarSgn";
 891       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 892       fillthis="histDeltaSgn_";
 893       fillthis+=ptbin;
 894       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 895       fillthis="histDeltaSgn";
 896       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 897     if (massInRange) {
 898         fillthis="ptSgn";
 899         ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
 900         fillthis="etaSgn";
 901         ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
 902       }
 903     }
 904     else {//background
 905       fillthis="histD0Bkg_";
 906       fillthis+=ptbin;
 907       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 908       fillthis="histD0Bkg";
 909       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 910       fillthis="histDstarBkg_";
 911       fillthis+=ptbin;
 912       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 913       fillthis="histDstarBkg";
 914       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 915       fillthis="histDeltaBkg_";
 916       fillthis+=ptbin;
 917       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 918       fillthis="histDeltaBkg";
 919       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 920      if (massInRange) {
 921         fillthis="ptBkg";
 922         ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
 923         fillthis="etaBkg";
 924         ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
 925       }
 926     }
 927   }
 928   //no MC info, just cut selection
 929   fillthis="histD0Mass_";
 930   fillthis+=ptbin;
 931   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 932   fillthis="histD0Mass";
 933   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 934   fillthis="histDstarMass_";
 935   fillthis+=ptbin;
 936   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 937   fillthis="histDstarMass";
 938   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 939   fillthis="histDeltaMass_";
 940   fillthis+=ptbin;
 941   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 942   fillthis="histDeltaMass";
 943   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 944
 945   if (massInRange) {
 946     fillthis="ptMass";
 947     ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
 948     fillthis="etaMass";
 949     ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
 950   }
 951
 952   return;
 953 }
 954 //______________________________ side band background for D*___________________________________
 955 void AliAnalysisTaskSEDStarSpectra::SideBandBackground(AliAODRecoCascadeHF *part,  AliRDHFCutsDStartoKpipi *cuts, Int_t isSel, TList *listout){
 956
 957   //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
 958   // (expected detector resolution) on the left and right frm the D0 mass. Each band
 959   //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
 960
 961   if(!isSel) return;
 962
 963   Int_t ptbin=cuts->PtBin(part->Pt());
 964
 965   Bool_t massInRange=kFALSE;
 966
 967   // select the side bands intervall
 968   Double_t invmassD0    = part->InvMassD0();
 969   if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){
 970
 971     // for pt and eta
 972     Double_t invmassDelta = part->DeltaInvMass();
 973     if (TMath::Abs(invmassDelta-0.14557)<fPeakWindow) massInRange=kTRUE;
 974
 975     TString fillthis="";
 976     fillthis="histSideBandMass_";
 977     fillthis+=ptbin;
 978     ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 979     fillthis="histSideBandMass";
 980     ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 981
 982   }
 983 }
 984 //________________________________________________________________________________________________________________
 985 void AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(AliAODRecoCascadeHF *part,  AliRDHFCutsDStartoKpipi *cuts, TList *listout){
 986   //
 987   // assign the wrong charge to the soft pion to create background
 988   //
 989   Int_t ptbin=cuts->PtBin(part->Pt());
 990
 991   Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 992   Double_t invmassD0   = part->InvMassD0();
 993   if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return;
 994
 995   AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();
 996
 997   Int_t okD0WrongSign,okD0barWrongSign;
 998   Double_t wrongMassD0=0.;
 999
1000   Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
1001    if (!isSelected){
1002     return;
1003   }
1004
1005   okD0WrongSign =  1;
1006   okD0barWrongSign = 1;
1007
1008   //if is D*+ than assume D0bar
1009   if(part->Charge()>0 && (isSelected ==1)) {
1010     okD0WrongSign = 0;
1011   }
1012   if(part->Charge()<0 && (isSelected ==2)){
1013     okD0barWrongSign = 0;
1014   }
1015
1016   // assign the wrong mass in case the cuts return both D0 and D0bar
1017   if(part->Charge()>0 && (isSelected ==3)) {
1018     okD0WrongSign = 0;
1019   } else if(part->Charge()<0 && (isSelected ==3)){
1020     okD0barWrongSign = 0;
1021   }
1022
1023   //wrong D0 inv mass
1024   if(okD0WrongSign!=0){
1025     wrongMassD0 = theD0particle->InvMassD0();
1026   }else if(okD0WrongSign==0){
1027     wrongMassD0 = theD0particle->InvMassD0bar();
1028   }
1029
1030   if(TMath::Abs(wrongMassD0-1.865)<fD0Window){
1031
1032     // wrong D* inv mass
1033     Double_t e[3];
1034     if (part->Charge()>0){
1035       e[0]=theD0particle->EProng(0,321);
1036       e[1]=theD0particle->EProng(1,211);
1037     }else{
1038       e[0]=theD0particle->EProng(0,211);
1039       e[1]=theD0particle->EProng(1,321);
1040     }
1041     e[2]=part->EProng(0,211);
1042
1043     Double_t esum = e[0]+e[1]+e[2];
1044     Double_t pds = part->P();
1045
1046     Double_t   wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);
1047
1048     TString fillthis="";
1049     fillthis="histWrongSignMass_";
1050     fillthis+=ptbin;
1051     ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1052     fillthis="histWrongSignMass";
1053     ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1054
1055   }
1056 }
1057 //-------------------------------------------------------------------------------
1058 Int_t AliAnalysisTaskSEDStarSpectra::CheckOrigin(TClonesArray* arrayMC, const AliAODMCParticle *mcPartCandidate) const {
1059   //
1060   // checking whether the mother of the particles come from a charm or a bottom quark
1061   //
1062
1063   Int_t pdgGranma = 0;
1064   Int_t mother = 0;
1065   mother = mcPartCandidate->GetMother();
1066   Int_t istep = 0;
1067   Int_t abspdgGranma =0;
1068   Bool_t isFromB=kFALSE;
1069   Bool_t isQuarkFound=kFALSE;
1070   while (mother >0 ){
1071     istep++;
1072     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1073     if (mcGranma){
1074       pdgGranma = mcGranma->GetPdgCode();
1075       abspdgGranma = TMath::Abs(pdgGranma);
1076       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1077         isFromB=kTRUE;
1078       }
1079       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1080       mother = mcGranma->GetMother();
1081     }else{
1082       AliError("Failed casting the mother particle!");
1083       break;
1084     }
1085   }
1086
1087   if(isFromB) return 5;
1088   else return 4;
1089 }
1090 //-------------------------------------------------------------------------------------
1091 Float_t AliAnalysisTaskSEDStarSpectra::GetTrueImpactParameterD0(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
1092   // true impact parameter calculation
1093
1094   Double_t vtxTrue[3];
1095   mcHeader->GetVertex(vtxTrue);
1096   Double_t origD[3];
1097   partDp->XvYvZv(origD);
1098   Short_t charge=partDp->Charge();
1099   Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1100   Int_t labelFirstDau = partDp->GetDaughter(0);
1101
1102   Int_t nDau=partDp->GetNDaughters();
1103
1104   Int_t theDau=0;
1105   if(nDau==2){
1106     for(Int_t iDau=0; iDau<2; iDau++){
1107       Int_t ind = labelFirstDau+iDau;
1108       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1109       if(!part){
1110         AliError("Daughter particle not found in MC array");
1111         return 99999.;
1112       }
1113       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1114       if(pdgCode==211 || pdgCode==321){
1115         pXdauTrue[theDau]=part->Px();
1116         pYdauTrue[theDau]=part->Py();
1117         pZdauTrue[theDau]=part->Pz();
1118         ++theDau;
1119       }
1120     }
1121   }
1122   if(theDau!=2){
1123     AliError("Wrong number of decay prongs");
1124     return 99999.;
1125   }
1126
1127   Double_t d0dummy[3]={0.,0.,0.};
1128   AliAODRecoDecayHF aodD0MC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1129   return aodD0MC.ImpParXY();
1130
1131 }
1132 //______________________________________________________-
1133 void AliAnalysisTaskSEDStarSpectra::CreateImpactParameterHistos(){
1134   // Histos for impact paramter study
1135
1136   Int_t nbins[3]={400,200,fNImpParBins};
1137   Double_t xmin[3]={1.75,0.,fLowerImpPar};
1138   Double_t xmax[3]={1.98,20.,fHigherImpPar};
1139
1140   fHistMassPtImpParTCDs[0]=new THnSparseF("hMassPtImpParAll",
1141                                         "Mass vs. pt vs.imppar - All",
1142                                         3,nbins,xmin,xmax);
1143   fHistMassPtImpParTCDs[1]=new THnSparseF("hMassPtImpParPrompt",
1144                                         "Mass vs. pt vs.imppar - promptD",
1145                                         3,nbins,xmin,xmax);
1146   fHistMassPtImpParTCDs[2]=new THnSparseF("hMassPtImpParBfeed",
1147                                         "Mass vs. pt vs.imppar - DfromB",
1148                                         3,nbins,xmin,xmax);
1149   fHistMassPtImpParTCDs[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1150                                         "Mass vs. pt vs.true imppar -DfromB",
1151                                         3,nbins,xmin,xmax);
1152   fHistMassPtImpParTCDs[4]=new THnSparseF("hMassPtImpParBkg",
1153                                         "Mass vs. pt vs.imppar - backgr.",
1154                                         3,nbins,xmin,xmax);
1155
1156   for(Int_t i=0; i<5;i++){
1157     fOutput->Add(fHistMassPtImpParTCDs[i]);
1158   }
1159 }
1160