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