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