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