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