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