PWG3/vertexingHF/AliAnalysisTaskSEDStarSpectra.cxx

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