PWGHF/vertexingHF/AliAnalysisTaskSEDStarSpectra.cxx

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