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 //
  17 //                  Base class for DStar Analysis
  18 //
  19 //
  20 //  The D* spectra study is done in pt bins:
  21 //  [0.7,1] [1,2] [2,3] [3,5] [5,8] [8,12],[12,18]
  22 //
  23 //  Optimized cuts used and TPC PID is on request (flag in che .C)
  24 //  Cuts option of analysis: 0 Heidelberg ; 1 Utrecht
  25 //  Side Band and like sign background are implemented in the macro
  26 //
  27 //-----------------------------------------------------------------------
  28 //
  29 //                         Author A.Grelli
  30 //              ERC-QGP Utrecht University - a.grelli@uu.nl,
  31 //                         Author Y.Wang
  32 //        University of Heidelberg - yifei@physi.uni-heidelberg.de
  33 //                         Author C.Ivan
  34 //             ERC-QGP Utrecht University - c.ivan@uu.nl,
  35 //
  36 //-----------------------------------------------------------------------
  37
  38 #include <TSystem.h>
  39 #include <TParticle.h>
  40 #include <TH1I.h>
  41 #include "TROOT.h"
  42
  43
  44 #include <AliAnalysisDataSlot.h>
  45 #include <AliAnalysisDataContainer.h>
  46 #include "AliRDHFCutsDStartoKpipi.h"
  47 #include "AliPID.h"
  48 #include "AliTPCPIDResponse.h"
  49 //#include "AliAODPidHF.h"
  50 #include "AliStack.h"
  51 #include "AliMCEvent.h"
  52 #include "AliAnalysisManager.h"
  53 #include "AliAODMCHeader.h"
  54 #include "AliAODHandler.h"
  55 #include "AliLog.h"
  56 #include "AliAODVertex.h"
  57 //#include "AliAODJet.h"
  58 #include "AliAODRecoDecay.h"
  59 #include "AliAODRecoDecayHF.h"
  60 #include "AliAODRecoCascadeHF.h"
  61 #include "AliAODRecoDecayHF2Prong.h"
  62 #include "AliAnalysisVertexingHF.h"
  63 #include "AliESDtrack.h"
  64 //#include "AliVertexerTracks.h"
  65 #include "AliAODMCParticle.h"
  66 #include "AliAnalysisTaskSE.h"
  67 #include "AliAnalysisTaskSEDStarSpectra.h"
  68
  69 ClassImp(AliAnalysisTaskSEDStarSpectra)
  70
  71 //__________________________________________________________________________
  72 AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra():
  73   AliAnalysisTaskSE(),
  74   fEvents(0),
  75   fAnalysis(0),
  76   fD0Window(0),
  77   fPeakWindow(0),
  78   fUseMCInfo(kFALSE),
  79   fOutput(0),
  80   fOutputSpectrum(0),
  81   fOutputAll(0),
  82   fOutputPID3(0),
  83   fOutputPID2(0),
  84   fOutputPID1(0),
  85   fNSigma(3),
  86   fPID(kTRUE),
  87   fCuts(0),
  88   fCEvents(0),
  89   fTrueDiff2(0)
  90 {
  91   //
  92   // Default ctor
  93   //
  94 }
  95 //___________________________________________________________________________
  96 AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
  97   AliAnalysisTaskSE(name),
  98   fEvents(0),
  99   fAnalysis(0),
 100   fD0Window(0),
 101   fPeakWindow(0),
 102   fUseMCInfo(kFALSE),
 103   fOutput(0),
 104   fOutputSpectrum(0),
 105   fOutputAll(0),
 106   fOutputPID3(0),
 107   fOutputPID2(0),
 108   fOutputPID1(0),
 109   fNSigma(3),
 110   fPID(kTRUE),
 111   fCuts(0),
 112   fCEvents(0),
 113   fTrueDiff2(0)
 114 {
 115   //
 116   // Constructor. Initialization of Inputs and Outputs
 117   //
 118   Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");
 119
 120   fCuts=cuts;
 121
 122   DefineOutput(1,TList::Class());  //conters
 123   DefineOutput(2,TList::Class());  //Spectrum output
 124   DefineOutput(3,TList::Class());  //All Entries output
 125   DefineOutput(4,TList::Class());  //3sigma PID output
 126   DefineOutput(5,TList::Class());  //2sigma PID output
 127   DefineOutput(6,TList::Class());  //1sigma PID output
 128   DefineOutput(7,AliRDHFCutsDStartoKpipi::Class());  //My private output
 129
 130 }
 131
 132 //___________________________________________________________________________
 133 AliAnalysisTaskSEDStarSpectra::~AliAnalysisTaskSEDStarSpectra() {
 134   //
 135   // destructor
 136   //
 137   Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");
 138
 139   if (fOutput) {
 140     delete fOutput;
 141     fOutput = 0;
 142   }
 143   if (fOutputSpectrum) {
 144     delete fOutputSpectrum;
 145     fOutputSpectrum = 0;
 146   }
 147   if (fOutputAll) {
 148     delete fOutputAll;
 149     fOutputAll = 0;
 150   }
 151   if (fOutputPID3) {
 152     delete fOutputPID3;
 153     fOutputPID3 = 0;
 154   }
 155   if (fOutputPID2) {
 156     delete fOutputPID2;
 157     fOutputPID2 = 0;
 158   }
 159   if (fOutputPID1) {
 160     delete fOutputPID1;
 161     fOutputPID1 = 0;
 162   }
 163   if (fCuts) {
 164     delete fCuts;
 165     fCuts = 0;
 166   }
 167   if(fCEvents){
 168     delete fCEvents;
 169     fCEvents =0;
 170   }
 171 }
 172 //_________________________________________________
 173 void AliAnalysisTaskSEDStarSpectra::Init(){
 174   //
 175   // Initialization
 176   //
 177
 178   if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");
 179    AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
 180   // Post the data
 181   PostData(7,copyfCuts);
 182
 183   return;
 184 }
 185
 186 //_________________________________________________
 187 void AliAnalysisTaskSEDStarSpectra::UserExec(Option_t *)
 188 {
 189   // user exec
 190   if (!fInputEvent) {
 191     Error("UserExec","NO EVENT FOUND!");
 192     return;
 193   }
 194
 195   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 196   TClonesArray *arrayDStartoD0pi=0;
 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
 216   // fix for temporary bug in ESDfilter
 217   // the AODs with null vertex pointer didn't pass the PhysSel
 218   if(!aodEvent->GetPrimaryVertex()) return;
 219
 220
 221   fCEvents->Fill(1);
 222   // Load the event
 223   fEvents++;
 224   AliInfo(Form("Event %d",fEvents));
 225   if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
 226
 227   // counters for efficiencies
 228   Int_t icountReco = 0;
 229
 230   //D* and D0 prongs needed to MatchToMC method
 231   Int_t pdgDgDStartoD0pi[2]={421,211};
 232   Int_t pdgDgD0toKpi[2]={321,211};
 233
 234   // AOD primary vertex
 235   AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 236   if(!vtx1) return;
 237   if(vtx1->GetNContributors()<1) return;
 238
 239   if (!arrayDStartoD0pi){
 240     AliInfo("Could not find array of HF vertices, skipping the event");
 241     return;
 242   }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 243
 244   // loop over the tracks to search for candidates soft pion
 245
 246   for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
 247
 248     // D* candidates and D0 from D*
 249     AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
 250     AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 251     if (!theD0particle) continue;
 252
 253     Int_t isDStar = 0;
 254     // mc analysis
 255     if(fUseMCInfo){
 256     //MC array need for maching
 257       TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 258       if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 259       // find associated MC particle for D* ->D0toKpi
 260       Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
 261       if(mcLabel>=0) isDStar = 1;
 262     }
 263
 264     // soft pion candidate
 265     AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
 266
 267     Double_t pt = dstarD0pi->Pt();
 268     Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
 269     if(!isTkSelected) continue;
 270
 271     // cut in acceptance for the soft pion and for the D0 daughters  - TO BE REMOVED ONCE WILL BE IN THE CUTTING CLASS
 272     Bool_t okTracks = SingleTrackSelections(theD0particle, track2);
 273     if (!okTracks) continue;
 274
 275     Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
 276
 277     // set the D0 search window bin by bin
 278     if (ptbin==0){
 279       if(fAnalysis==1){
 280         fD0Window=0.015;
 281         fPeakWindow=0.0018;
 282       }else{
 283         fD0Window=0.020;
 284         fPeakWindow=0.0018;
 285       }
 286     }
 287     if (ptbin==1){
 288       if(fAnalysis==1){
 289         fD0Window=0.015;
 290         fPeakWindow=0.0018;
 291       }else{
 292         fD0Window=0.020;
 293         fPeakWindow=0.0018;
 294       }
 295     }
 296     if (ptbin==2){
 297       if(fAnalysis==1){
 298         fD0Window=0.018;
 299         fPeakWindow=0.0018;
 300       }else{
 301         fD0Window=0.020;
 302         fPeakWindow=0.0018;
 303       }
 304     }
 305     if (ptbin==3){
 306       if(fAnalysis==1){
 307         fD0Window=0.036;
 308         fPeakWindow=0.0018;
 309       }else{
 310         fD0Window=0.022;
 311         fPeakWindow=0.0016;
 312       }
 313     }
 314     if (ptbin==4){
 315       if(fAnalysis==1){
 316         fD0Window=0.036;
 317         fPeakWindow=0.0016;
 318       }else{
 319         fD0Window=0.026;
 320         fPeakWindow=0.0014;
 321       }
 322     }
 323     if (ptbin>=5){
 324       if(fAnalysis==1){
 325         fD0Window=0.062;
 326         fPeakWindow=0.0014;
 327       }else{
 328         fD0Window=0.026;
 329         fPeakWindow=0.0014;
 330       }
 331     }
 332
 333     FillSpectrum(dstarD0pi,isDStar,1,3,fCuts,fOutputPID3);
 334     FillSpectrum(dstarD0pi,isDStar,1,2,fCuts,fOutputPID2);
 335     FillSpectrum(dstarD0pi,isDStar,1,1,fCuts,fOutputPID1);
 336     FillSpectrum(dstarD0pi,isDStar,fPID,fNSigma,fCuts,fOutputSpectrum);
 337     //FillSpectrum(dstarD0pi,isDStar,0,0,fCuts,fOutputAll);
 338
 339     SideBandBackground(dstarD0pi,1,3,fCuts,fOutputPID3);
 340     SideBandBackground(dstarD0pi,1,2,fCuts,fOutputPID2);
 341     SideBandBackground(dstarD0pi,1,1,fCuts,fOutputPID1);
 342     SideBandBackground(dstarD0pi,fPID,fNSigma,fCuts,fOutputSpectrum);
 343     //SideBandBackground(dstarD0pi,0,0,fCuts,fOutputAll);
 344
 345     WrongSignForDStar(dstarD0pi,1,3,fCuts,fOutputPID3);
 346     WrongSignForDStar(dstarD0pi,1,2,fCuts,fOutputPID2);
 347     WrongSignForDStar(dstarD0pi,1,1,fCuts,fOutputPID1);
 348     WrongSignForDStar(dstarD0pi,fPID,fNSigma,fCuts,fOutputSpectrum);
 349     //WrongSignForDStar(dstarD0pi,0,0,fCuts,fOutputAll);
 350
 351     if(isDStar == 1) {
 352       fTrueDiff2->Fill(pt,dstarD0pi->DeltaInvMass());
 353     }
 354
 355   }
 356
 357   AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
 358
 359   PostData(1,fOutput);
 360   PostData(2,fOutputSpectrum);
 361   PostData(3,fOutputAll);
 362   PostData(4,fOutputPID3);
 363   PostData(5,fOutputPID2);
 364   PostData(6,fOutputPID1);
 365
 366
 367 }
 368 //________________________________________ terminate ___________________________
 369 void AliAnalysisTaskSEDStarSpectra::Terminate(Option_t*)
 370 {
 371   // The Terminate() function is the last function to be called during
 372   // a query. It always runs on the client, it can be used to present
 373   // the results graphically or save the results to file.
 374
 375   Info("Terminate","");
 376   AliAnalysisTaskSE::Terminate();
 377
 378   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 379   if (!fOutput) {
 380     printf("ERROR: fOutput not available\n");
 381     return;
 382   }
 383
 384   fCEvents        = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
 385   fTrueDiff2      = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
 386
 387   fOutputSpectrum = dynamic_cast<TList*> (GetOutputData(2));
 388   if (!fOutputSpectrum) {
 389     printf("ERROR: fOutputSpectrum not available\n");
 390     return;
 391   }
 392   fOutputAll = dynamic_cast<TList*> (GetOutputData(3));
 393   if (!fOutputAll) {
 394     printf("ERROR: fOutputAll not available\n");
 395     return;
 396   }
 397   fOutputPID3 = dynamic_cast<TList*> (GetOutputData(4));
 398   if (!fOutputPID3) {
 399     printf("ERROR: fOutputPID3 not available\n");
 400     return;
 401   }
 402   fOutputPID2 = dynamic_cast<TList*> (GetOutputData(5));
 403   if (!fOutputPID2) {
 404     printf("ERROR: fOutputPID2 not available\n");
 405     return;
 406   }
 407   fOutputPID1 = dynamic_cast<TList*> (GetOutputData(6));
 408   if (!fOutputPID1) {
 409     printf("ERROR: fOutputPID1 not available\n");
 410     return;
 411   }
 412
 413   return;
 414 }
 415 //___________________________________________________________________________
 416 void AliAnalysisTaskSEDStarSpectra::UserCreateOutputObjects() {
 417  // output
 418   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 419
 420   //slot #1
 421   //OpenFile(1);
 422   fOutput = new TList();
 423   fOutput->SetOwner();
 424   fOutput->SetName("chist0");
 425
 426
 427   fOutputSpectrum = new TList();
 428   fOutputSpectrum->SetOwner();
 429   fOutputSpectrum->SetName("listSpectrum");
 430
 431   fOutputPID3 = new TList();
 432   fOutputPID3->SetOwner();
 433   fOutputPID3->SetName("listPID3");
 434
 435   fOutputPID2 = new TList();
 436   fOutputPID2->SetOwner();
 437   fOutputPID2->SetName("listPID2");
 438
 439   fOutputPID1 = new TList();
 440   fOutputPID1->SetOwner();
 441   fOutputPID1->SetName("listPID1");
 442
 443   fOutputAll = new TList();
 444   fOutputAll->SetOwner();
 445   fOutputAll->SetName("listAll");
 446
 447   // define histograms
 448   DefineHistograms();
 449
 450   PostData(1,fOutput);
 451   PostData(2,fOutputSpectrum);
 452   PostData(3,fOutputAll);
 453   PostData(4,fOutputPID3);
 454   PostData(5,fOutputPID2);
 455   PostData(6,fOutputPID1);
 456
 457   return;
 458 }
 459 //___________________________________ hiostograms _______________________________________
 460 void  AliAnalysisTaskSEDStarSpectra::DefineHistograms(){
 461
 462   fCEvents = new TH1F("fCEvents","conter",10,0,10);
 463   fCEvents->SetStats(kTRUE);
 464   fCEvents->GetXaxis()->SetTitle("1");
 465   fCEvents->GetYaxis()->SetTitle("counts");
 466   fOutput->Add(fCEvents);
 467
 468   fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
 469   fOutput->Add(fTrueDiff2);
 470
 471   const Int_t nhist=6;
 472   TString nameMass=" ", nameSgn=" ", nameBkg=" ";
 473
 474   for(Int_t i=-2;i<nhist;i++){
 475     nameMass="histDeltaMass_";
 476     nameMass+=i+1;
 477     nameSgn="histDeltaSgn_";
 478     nameSgn+=i+1;
 479     nameBkg="histDeltaBkg_";
 480     nameBkg+=i+1;
 481
 482     if (i==-2) {
 483       nameMass="histDeltaMass";
 484       nameSgn="histDeltaSgn";
 485       nameBkg="histDeltaBkg";
 486     }
 487
 488     TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
 489     TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
 490     TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
 491
 492     nameMass="histD0Mass_";
 493     nameMass+=i+1;
 494     nameSgn="histD0Sgn_";
 495     nameSgn+=i+1;
 496     nameBkg="histD0Bkg_";
 497     nameBkg+=i+1;
 498
 499     if (i==-2) {
 500       nameMass="histD0Mass";
 501       nameSgn="histD0Sgn";
 502       nameBkg="histD0Bkg";
 503     }
 504
 505     TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
 506     TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
 507     TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
 508
 509     nameMass="histDstarMass_";
 510     nameMass+=i+1;
 511     nameSgn="histDstarSgn_";
 512     nameSgn+=i+1;
 513     nameBkg="histDstarBkg_";
 514     nameBkg+=i+1;
 515
 516     if (i==-2) {
 517       nameMass="histDstarMass";
 518       nameSgn="histDstarSgn";
 519       nameBkg="histDstarBkg";
 520     }
 521
 522     TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
 523     TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
 524     TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
 525
 526     nameMass="histSideBandMass_";
 527     nameMass+=i+1;
 528     if (i==-2) {
 529       nameMass="histSideBandMass";
 530     }
 531
 532     TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
 533
 534     nameMass="histWrongSignMass_";
 535     nameMass+=i+1;
 536     if (i==-2) {
 537       nameMass="histWrongSignMass";
 538     }
 539
 540     TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
 541
 542
 543     spectrumMass->Sumw2();
 544     spectrumSgn->Sumw2();
 545     spectrumBkg->Sumw2();
 546
 547     spectrumMass->SetLineColor(6);
 548     spectrumSgn->SetLineColor(2);
 549     spectrumBkg->SetLineColor(4);
 550
 551     spectrumMass->SetMarkerStyle(20);
 552     spectrumSgn->SetMarkerStyle(20);
 553     spectrumBkg->SetMarkerStyle(20);
 554     spectrumMass->SetMarkerSize(0.6);
 555     spectrumSgn->SetMarkerSize(0.6);
 556     spectrumBkg->SetMarkerSize(0.6);
 557     spectrumMass->SetMarkerColor(6);
 558     spectrumSgn->SetMarkerColor(2);
 559     spectrumBkg->SetMarkerColor(4);
 560
 561     spectrumD0Mass->Sumw2();
 562     spectrumD0Sgn->Sumw2();
 563     spectrumD0Bkg->Sumw2();
 564
 565     spectrumD0Mass->SetLineColor(6);
 566     spectrumD0Sgn->SetLineColor(2);
 567     spectrumD0Bkg->SetLineColor(4);
 568
 569     spectrumD0Mass->SetMarkerStyle(20);
 570     spectrumD0Sgn->SetMarkerStyle(20);
 571     spectrumD0Bkg->SetMarkerStyle(20);
 572     spectrumD0Mass->SetMarkerSize(0.6);
 573     spectrumD0Sgn->SetMarkerSize(0.6);
 574     spectrumD0Bkg->SetMarkerSize(0.6);
 575     spectrumD0Mass->SetMarkerColor(6);
 576     spectrumD0Sgn->SetMarkerColor(2);
 577     spectrumD0Bkg->SetMarkerColor(4);
 578
 579     spectrumDstarMass->Sumw2();
 580     spectrumDstarSgn->Sumw2();
 581     spectrumDstarBkg->Sumw2();
 582
 583     spectrumDstarMass->SetLineColor(6);
 584     spectrumDstarSgn->SetLineColor(2);
 585     spectrumDstarBkg->SetLineColor(4);
 586
 587     spectrumDstarMass->SetMarkerStyle(20);
 588     spectrumDstarSgn->SetMarkerStyle(20);
 589     spectrumDstarBkg->SetMarkerStyle(20);
 590     spectrumDstarMass->SetMarkerSize(0.6);
 591     spectrumDstarSgn->SetMarkerSize(0.6);
 592     spectrumDstarBkg->SetMarkerSize(0.6);
 593     spectrumDstarMass->SetMarkerColor(6);
 594     spectrumDstarSgn->SetMarkerColor(2);
 595     spectrumDstarBkg->SetMarkerColor(4);
 596
 597     spectrumSideBandMass->Sumw2();
 598     spectrumSideBandMass->SetLineColor(4);
 599     spectrumSideBandMass->SetMarkerStyle(20);
 600     spectrumSideBandMass->SetMarkerSize(0.6);
 601     spectrumSideBandMass->SetMarkerColor(4);
 602
 603     spectrumWrongSignMass->Sumw2();
 604     spectrumWrongSignMass->SetLineColor(4);
 605     spectrumWrongSignMass->SetMarkerStyle(20);
 606     spectrumWrongSignMass->SetMarkerSize(0.6);
 607     spectrumWrongSignMass->SetMarkerColor(4);
 608
 609     TH1F* allMass = (TH1F*)spectrumMass->Clone();
 610     TH1F* allSgn  = (TH1F*)spectrumSgn->Clone();
 611     TH1F* allBkg  = (TH1F*)spectrumBkg->Clone();
 612
 613     TH1F* pid3Mass = (TH1F*)spectrumMass->Clone();
 614     TH1F* pid3Sgn  = (TH1F*)spectrumSgn->Clone();
 615     TH1F* pid3Bkg  = (TH1F*)spectrumBkg->Clone();
 616
 617     TH1F* pid2Mass = (TH1F*)spectrumMass->Clone();
 618     TH1F* pid2Sgn  = (TH1F*)spectrumSgn->Clone();
 619     TH1F* pid2Bkg  = (TH1F*)spectrumBkg->Clone();
 620
 621     TH1F* pid1Mass = (TH1F*)spectrumMass->Clone();
 622     TH1F* pid1Sgn  = (TH1F*)spectrumSgn->Clone();
 623     TH1F* pid1Bkg  = (TH1F*)spectrumBkg->Clone();
 624
 625     fOutputSpectrum->Add(spectrumMass);
 626     fOutputSpectrum->Add(spectrumSgn);
 627     fOutputSpectrum->Add(spectrumBkg);
 628
 629     fOutputAll->Add(allMass);
 630     fOutputAll->Add(allSgn);
 631     fOutputAll->Add(allBkg);
 632
 633     fOutputPID3->Add(pid3Mass);
 634     fOutputPID3->Add(pid3Sgn);
 635     fOutputPID3->Add(pid3Bkg);
 636
 637     fOutputPID2->Add(pid2Mass);
 638     fOutputPID2->Add(pid2Sgn);
 639     fOutputPID2->Add(pid2Bkg);
 640
 641     fOutputPID1->Add(pid1Mass);
 642     fOutputPID1->Add(pid1Sgn);
 643     fOutputPID1->Add(pid1Bkg);
 644
 645     TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
 646     TH1F* allD0Sgn  = (TH1F*)spectrumD0Sgn->Clone();
 647     TH1F* allD0Bkg  = (TH1F*)spectrumD0Bkg->Clone();
 648
 649     TH1F* pid3D0Mass = (TH1F*)spectrumD0Mass->Clone();
 650     TH1F* pid3D0Sgn  = (TH1F*)spectrumD0Sgn->Clone();
 651     TH1F* pid3D0Bkg  = (TH1F*)spectrumD0Bkg->Clone();
 652
 653     TH1F* pid2D0Mass = (TH1F*)spectrumD0Mass->Clone();
 654     TH1F* pid2D0Sgn  = (TH1F*)spectrumD0Sgn->Clone();
 655     TH1F* pid2D0Bkg  = (TH1F*)spectrumD0Bkg->Clone();
 656
 657     TH1F* pid1D0Mass = (TH1F*)spectrumD0Mass->Clone();
 658     TH1F* pid1D0Sgn  = (TH1F*)spectrumD0Sgn->Clone();
 659     TH1F* pid1D0Bkg  = (TH1F*)spectrumD0Bkg->Clone();
 660
 661     fOutputSpectrum->Add(spectrumD0Mass);
 662     fOutputSpectrum->Add(spectrumD0Sgn);
 663     fOutputSpectrum->Add(spectrumD0Bkg);
 664
 665     fOutputAll->Add(allD0Mass);
 666     fOutputAll->Add(allD0Sgn);
 667     fOutputAll->Add(allD0Bkg);
 668
 669     fOutputPID3->Add(pid3D0Mass);
 670     fOutputPID3->Add(pid3D0Sgn);
 671     fOutputPID3->Add(pid3D0Bkg);
 672
 673     fOutputPID2->Add(pid2D0Mass);
 674     fOutputPID2->Add(pid2D0Sgn);
 675     fOutputPID2->Add(pid2D0Bkg);
 676
 677     fOutputPID1->Add(pid1D0Mass);
 678     fOutputPID1->Add(pid1D0Sgn);
 679     fOutputPID1->Add(pid1D0Bkg);
 680
 681     TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
 682     TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
 683     TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
 684
 685     TH1F* pid3DstarMass = (TH1F*)spectrumDstarMass->Clone();
 686     TH1F* pid3DstarSgn = (TH1F*)spectrumDstarSgn->Clone();
 687     TH1F* pid3DstarBkg = (TH1F*)spectrumDstarBkg->Clone();
 688
 689     TH1F* pid2DstarMass = (TH1F*)spectrumDstarMass->Clone();
 690     TH1F* pid2DstarSgn = (TH1F*)spectrumDstarSgn->Clone();
 691     TH1F* pid2DstarBkg = (TH1F*)spectrumDstarBkg->Clone();
 692
 693     TH1F* pid1DstarMass = (TH1F*)spectrumDstarMass->Clone();
 694     TH1F* pid1DstarSgn = (TH1F*)spectrumDstarSgn->Clone();
 695     TH1F* pid1DstarBkg = (TH1F*)spectrumDstarBkg->Clone();
 696
 697     fOutputSpectrum->Add(spectrumDstarMass);
 698     fOutputSpectrum->Add(spectrumDstarSgn);
 699     fOutputSpectrum->Add(spectrumDstarBkg);
 700
 701     fOutputAll->Add(allDstarMass);
 702     fOutputAll->Add(allDstarSgn);
 703     fOutputAll->Add(allDstarBkg);
 704
 705     fOutputPID3->Add(pid3DstarMass);
 706     fOutputPID3->Add(pid3DstarSgn);
 707     fOutputPID3->Add(pid3DstarBkg);
 708
 709     fOutputPID2->Add(pid2DstarMass);
 710     fOutputPID2->Add(pid2DstarSgn);
 711     fOutputPID2->Add(pid2DstarBkg);
 712
 713     fOutputPID1->Add(pid1DstarMass);
 714     fOutputPID1->Add(pid1DstarSgn);
 715     fOutputPID1->Add(pid1DstarBkg);
 716
 717     TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
 718     TH1F* pid3SideBandMass = (TH1F*)spectrumSideBandMass->Clone();
 719     TH1F* pid2SideBandMass = (TH1F*)spectrumSideBandMass->Clone();
 720     TH1F* pid1SideBandMass = (TH1F*)spectrumSideBandMass->Clone();
 721
 722     fOutputSpectrum->Add(spectrumSideBandMass);
 723     fOutputAll->Add(allSideBandMass);
 724     fOutputPID3->Add(pid3SideBandMass);
 725     fOutputPID2->Add(pid2SideBandMass);
 726     fOutputPID1->Add(pid1SideBandMass);
 727
 728     TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
 729     TH1F* pid3WrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
 730     TH1F* pid2WrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
 731     TH1F* pid1WrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
 732
 733     fOutputSpectrum->Add(spectrumWrongSignMass);
 734     fOutputAll->Add(allWrongSignMass);
 735     fOutputPID3->Add(pid3WrongSignMass);
 736     fOutputPID2->Add(pid2WrongSignMass);
 737     fOutputPID1->Add(pid1WrongSignMass);
 738
 739   }
 740
 741   // pt spectra
 742   nameMass="ptMass";
 743   nameSgn="ptSgn";
 744   nameBkg="ptBkg";
 745
 746   TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",200,0,10);
 747   TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
 748   TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
 749
 750   ptspectrumMass->Sumw2();
 751   ptspectrumSgn->Sumw2();
 752   ptspectrumBkg->Sumw2();
 753
 754   ptspectrumMass->SetLineColor(6);
 755   ptspectrumSgn->SetLineColor(2);
 756   ptspectrumBkg->SetLineColor(4);
 757
 758   ptspectrumMass->SetMarkerStyle(20);
 759   ptspectrumSgn->SetMarkerStyle(20);
 760   ptspectrumBkg->SetMarkerStyle(20);
 761   ptspectrumMass->SetMarkerSize(0.6);
 762   ptspectrumSgn->SetMarkerSize(0.6);
 763   ptspectrumBkg->SetMarkerSize(0.6);
 764   ptspectrumMass->SetMarkerColor(6);
 765   ptspectrumSgn->SetMarkerColor(2);
 766   ptspectrumBkg->SetMarkerColor(4);
 767
 768   TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
 769   TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
 770   TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
 771
 772   TH1F* ptpid3Mass = (TH1F*)ptspectrumMass->Clone();
 773   TH1F* ptpid3Sgn = (TH1F*)ptspectrumSgn->Clone();
 774   TH1F* ptpid3Bkg = (TH1F*)ptspectrumBkg->Clone();
 775
 776   TH1F* ptpid2Mass = (TH1F*)ptspectrumMass->Clone();
 777   TH1F* ptpid2Sgn = (TH1F*)ptspectrumSgn->Clone();
 778   TH1F* ptpid2Bkg = (TH1F*)ptspectrumBkg->Clone();
 779
 780   TH1F* ptpid1Mass = (TH1F*)ptspectrumMass->Clone();
 781   TH1F* ptpid1Sgn = (TH1F*)ptspectrumSgn->Clone();
 782   TH1F* ptpid1Bkg = (TH1F*)ptspectrumBkg->Clone();
 783
 784   fOutputSpectrum->Add(ptspectrumMass);
 785   fOutputSpectrum->Add(ptspectrumSgn);
 786   fOutputSpectrum->Add(ptspectrumBkg);
 787
 788   fOutputAll->Add(ptallMass);
 789   fOutputAll->Add(ptallSgn);
 790   fOutputAll->Add(ptallBkg);
 791
 792   fOutputPID3->Add(ptpid3Mass);
 793   fOutputPID3->Add(ptpid3Sgn);
 794   fOutputPID3->Add(ptpid3Bkg);
 795
 796   fOutputPID2->Add(ptpid2Mass);
 797   fOutputPID2->Add(ptpid2Sgn);
 798   fOutputPID2->Add(ptpid2Bkg);
 799
 800   fOutputPID1->Add(ptpid1Mass);
 801   fOutputPID1->Add(ptpid1Sgn);
 802   fOutputPID1->Add(ptpid1Bkg);
 803
 804   // eta spectra
 805   nameMass="etaMass";
 806   nameSgn="etaSgn";
 807   nameBkg="etaBkg";
 808
 809   TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
 810   TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
 811   TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
 812
 813   etaspectrumMass->Sumw2();
 814   etaspectrumSgn->Sumw2();
 815   etaspectrumBkg->Sumw2();
 816
 817   etaspectrumMass->SetLineColor(6);
 818   etaspectrumSgn->SetLineColor(2);
 819   etaspectrumBkg->SetLineColor(4);
 820
 821   etaspectrumMass->SetMarkerStyle(20);
 822   etaspectrumSgn->SetMarkerStyle(20);
 823   etaspectrumBkg->SetMarkerStyle(20);
 824   etaspectrumMass->SetMarkerSize(0.6);
 825   etaspectrumSgn->SetMarkerSize(0.6);
 826   etaspectrumBkg->SetMarkerSize(0.6);
 827   etaspectrumMass->SetMarkerColor(6);
 828   etaspectrumSgn->SetMarkerColor(2);
 829   etaspectrumBkg->SetMarkerColor(4);
 830
 831   TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
 832   TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
 833   TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
 834
 835   TH1F* etapid3Mass = (TH1F*)etaspectrumMass->Clone();
 836   TH1F* etapid3Sgn = (TH1F*)etaspectrumSgn->Clone();
 837   TH1F* etapid3Bkg = (TH1F*)etaspectrumBkg->Clone();
 838
 839   TH1F* etapid2Mass = (TH1F*)etaspectrumMass->Clone();
 840   TH1F* etapid2Sgn = (TH1F*)etaspectrumSgn->Clone();
 841   TH1F* etapid2Bkg = (TH1F*)etaspectrumBkg->Clone();
 842
 843   TH1F* etapid1Mass = (TH1F*)etaspectrumMass->Clone();
 844   TH1F* etapid1Sgn = (TH1F*)etaspectrumSgn->Clone();
 845   TH1F* etapid1Bkg = (TH1F*)etaspectrumBkg->Clone();
 846
 847   fOutputSpectrum->Add(etaspectrumMass);
 848   fOutputSpectrum->Add(etaspectrumSgn);
 849   fOutputSpectrum->Add(etaspectrumBkg);
 850
 851   fOutputAll->Add(etaallMass);
 852   fOutputAll->Add(etaallSgn);
 853   fOutputAll->Add(etaallBkg);
 854
 855   fOutputPID3->Add(etapid3Mass);
 856   fOutputPID3->Add(etapid3Sgn);
 857   fOutputPID3->Add(etapid3Bkg);
 858
 859   fOutputPID2->Add(etapid2Mass);
 860   fOutputPID2->Add(etapid2Sgn);
 861   fOutputPID2->Add(etapid2Bkg);
 862
 863   fOutputPID1->Add(etapid1Mass);
 864   fOutputPID1->Add(etapid1Sgn);
 865   fOutputPID1->Add(etapid1Bkg);
 866
 867   return;
 868 }
 869 //________________________________________________________________________
 870 void AliAnalysisTaskSEDStarSpectra::FillSpectrum(AliAODRecoCascadeHF *part, Int_t isDStar, Bool_t PIDon, Int_t nSigma, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
 871   //
 872   // Fill histos for D* spectrum
 873   //
 874
 875   Int_t ptbin=cuts->PtBin(part->Pt());
 876
 877   Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
 878   if (!isSelected){
 879     return;
 880   }
 881
 882   Double_t invmassD0   = part->InvMassD0();
 883   if (TMath::Abs(invmassD0-1.865)>fD0Window) return;
 884
 885   Double_t pt = part->Pt();
 886   Double_t eta = part->Eta();
 887
 888   AliAODTrack *softPi = (AliAODTrack*)part->GetBachelor();
 889
 890   //PID of D0 daughters
 891   AliAODTrack *pos = (AliAODTrack*)part->Get2Prong()->GetDaughter(0);
 892   AliAODTrack *neg = (AliAODTrack*)part->Get2Prong()->GetDaughter(1);
 893
 894   Bool_t isPID = kTRUE;
 895
 896   //  Double_t cutsN = -1;
 897
 898   //if(part->Charge()<0){
 899   //  cutsN = (pos->Pt()-neg->Pt())/(part->Get2Prong()->Pt());
 900   // }
 901   //if(part->Charge()>0){
 902   //  cutsN = (neg->Pt()-pos->Pt())/(part->Get2Prong()->Pt());
 903   // }
 904
 905   //  if(ptbin==1 || ptbin==2){
 906   //  if(cutsN>0.7) return;
 907   // }
 908   //if(ptbin==3){
 909   //  if(cutsN>0.8) return;
 910   // }
 911
 912
 913   if(fAnalysis==1 && ptbin==1){
 914     if(part->Get2Prong()->DecayLength()>0.1 || part->Get2Prong()->DecayLength()<0.015)return; //(0.05)
 915    }
 916   if(fAnalysis==1 && ptbin==2){
 917     if(part->Get2Prong()->DecayLength()>0.1 || part->Get2Prong()->DecayLength()<0.015)return;
 918   }
 919
 920   if (PIDon) {
 921     if(fDebug > 1) printf("AnalysisTaskSEDStar::TPCPIDon \n");
 922     if(fDebug > 1) printf("AnalysisTaskSEDStar::NSigmaTPC: %d\n", nSigma);
 923
 924     if (part->Charge()>0){
 925       if(!SelectPID(pos, AliPID::kPion, nSigma)) return;//pion+
 926       if(!SelectPID(neg, AliPID::kKaon, nSigma)) return;//kaon-
 927     }else{
 928       if(!SelectPID(pos, AliPID::kKaon, nSigma)) return;//kaon+
 929       if(!SelectPID(neg, AliPID::kPion, nSigma)) return;//pion-
 930     }
 931     if (ptbin>1){
 932       isPID =SelectTOFPID(part->Get2Prong(), softPi);
 933       if(!isPID) return;
 934     }
 935   }
 936
 937   Double_t invmassDelta = part->DeltaInvMass();
 938   Double_t invmassDstar = part->InvMassDstarKpipi();
 939
 940   TString fillthis="";
 941   Bool_t massInRange=kFALSE;
 942   if (TMath::Abs(invmassDelta-0.14557)<fPeakWindow) massInRange=kTRUE;
 943
 944
 945   if(fUseMCInfo) {
 946     if(isDStar==1) {
 947       fillthis="histD0Sgn_";
 948       fillthis+=ptbin;
 949       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 950       fillthis="histD0Sgn";
 951       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 952       fillthis="histDstarSgn_";
 953       fillthis+=ptbin;
 954       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 955       fillthis="histDstarSgn";
 956       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 957       fillthis="histDeltaSgn_";
 958       fillthis+=ptbin;
 959       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 960       fillthis="histDeltaSgn";
 961       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 962       //if (massInRange) {
 963       if(ptbin<=1){
 964         fillthis="ptSgn";
 965         ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
 966         fillthis="etaSgn";
 967         ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
 968       }
 969     }
 970     else {//background
 971       fillthis="histD0Bkg_";
 972       fillthis+=ptbin;
 973       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 974       fillthis="histD0Bkg";
 975       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 976       fillthis="histDstarBkg_";
 977       fillthis+=ptbin;
 978       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 979       fillthis="histDstarBkg";
 980       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
 981       fillthis="histDeltaBkg_";
 982       fillthis+=ptbin;
 983       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 984       fillthis="histDeltaBkg";
 985       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
 986       //if (massInRange) {
 987       if(ptbin<=1){
 988         fillthis="ptBkg";
 989         ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
 990         fillthis="etaBkg";
 991         ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
 992       }
 993     }
 994   }
 995   //no MC info, just cut selection
 996   fillthis="histD0Mass_";
 997   fillthis+=ptbin;
 998   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 999   fillthis="histD0Mass";
1000   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1001   fillthis="histDstarMass_";
1002   fillthis+=ptbin;
1003   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
1004   fillthis="histDstarMass";
1005   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
1006   fillthis="histDeltaMass_";
1007   fillthis+=ptbin;
1008   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
1009   fillthis="histDeltaMass";
1010   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
1011
1012   if (massInRange) {
1013     fillthis="ptMass";
1014     ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
1015     fillthis="etaMass";
1016     ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
1017   }
1018
1019   return;
1020 }
1021 //______________________________ side band background for D*___________________________________
1022 void AliAnalysisTaskSEDStarSpectra::SideBandBackground(AliAODRecoCascadeHF *part, Bool_t PIDon, Int_t nSigma,  AliRDHFCutsDStartoKpipi *cuts, TList *listout){
1023
1024   //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1025   // (expected detector resolution) on the left and right frm the D0 mass. Each band
1026   //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
1027
1028   Int_t ptbin=cuts->PtBin(part->Pt());
1029
1030   Bool_t massInRange=kFALSE;
1031
1032   Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
1033   if (!isSelected){
1034     return;
1035   }
1036   //  Double_t pt = part->Pt();
1037
1038   AliAODTrack *softPi = (AliAODTrack*)part->GetBachelor();
1039
1040   // select the side bands intervall
1041   Double_t invmassD0    = part->InvMassD0();
1042   if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){
1043
1044     // for pt and eta
1045     Double_t invmassDelta = part->DeltaInvMass();
1046     if (TMath::Abs(invmassDelta-0.14557)<fPeakWindow) massInRange=kTRUE;
1047
1048     //PID of D0 daughters
1049     AliAODTrack *pos = (AliAODTrack*)part->Get2Prong()->GetDaughter(0);
1050     AliAODTrack *neg = (AliAODTrack*)part->Get2Prong()->GetDaughter(1);
1051
1052     //Double_t cutsN = -1;
1053
1054
1055     /* if(part->Charge()<0){
1056       cutsN = (pos->Pt()-neg->Pt())/(part->Get2Prong()->Pt());
1057     }
1058     if(part->Charge()>0){
1059       cutsN = (neg->Pt()-pos->Pt())/(part->Get2Prong()->Pt());
1060     }
1061
1062     if(ptbin==1 || ptbin==2){
1063       if(cutsN>0.5) return;
1064     }
1065     if(ptbin==3){
1066       if(cutsN>0.7) return;
1067     }
1068     */
1069     if(fAnalysis==1 && ptbin==1){
1070       if(part->Get2Prong()->DecayLength()>0.08 || part->Get2Prong()->DecayLength()<0.022)return; //(0.05)
1071     }
1072     if(fAnalysis==1 && ptbin==2){
1073       if(part->Get2Prong()->DecayLength()>0.08 || part->Get2Prong()->DecayLength()<0.017)return;
1074     }
1075
1076     Bool_t isPID = kTRUE;
1077
1078     if (PIDon) {
1079       if(fDebug > 1) printf("AnalysisTaskSEDStar::TPCPIDon \n");
1080       if(fDebug > 1) printf("AnalysisTaskSEDStar::NSigmaTPC: %d\n", nSigma);
1081
1082       if (part->Charge()>0){
1083         if(!SelectPID(pos, AliPID::kPion, nSigma)) return;//pion+
1084         if(!SelectPID(neg, AliPID::kKaon, nSigma)) return;//kaon-
1085       }else{
1086         if(!SelectPID(pos, AliPID::kKaon, nSigma)) return;//kaon+
1087         if(!SelectPID(neg, AliPID::kPion, nSigma)) return;//pion-
1088       }
1089       if (ptbin>2){
1090         isPID =SelectTOFPID(part->Get2Prong(), softPi);
1091         if(!isPID) return;
1092       }
1093     }
1094
1095     TString fillthis="";
1096     fillthis="histSideBandMass_";
1097     fillthis+=ptbin;
1098     ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
1099     fillthis="histSideBandMass";
1100     ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
1101
1102   }
1103 }
1104 //________________________________________________________________________________________________________________
1105 void AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(AliAODRecoCascadeHF *part, Bool_t PIDon, Int_t nSigma,  AliRDHFCutsDStartoKpipi *cuts, TList *listout){
1106   //
1107   // assign the wrong charge to the soft pion to create background
1108   //
1109   Int_t ptbin=cuts->PtBin(part->Pt());
1110
1111   AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();
1112   AliAODTrack *theSoftPi = (AliAODTrack*)part->GetBachelor();
1113
1114   Int_t okD0WrongSign,okD0barWrongSign;
1115   Double_t wrongMassD0=0.;
1116
1117   Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
1118    if (!isSelected){
1119     return;
1120   }
1121
1122   okD0WrongSign =  1;
1123   okD0barWrongSign = 1;
1124
1125   //if is D*+ than assume D0bar
1126   if(part->Charge()>0 && (isSelected ==1)) {
1127     okD0WrongSign = 0;
1128   }
1129   if(part->Charge()<0 && (isSelected ==2)){
1130     okD0barWrongSign = 0;
1131   }
1132
1133   // assign the wrong mass in case the cuts return both D0 and D0bar
1134   if(part->Charge()>0 && (isSelected ==3)) {
1135     okD0WrongSign = 0;
1136   } else if(part->Charge()<0 && (isSelected ==3)){
1137     okD0barWrongSign = 0;
1138   }
1139
1140   //wrong D0 inv mass
1141   if(okD0WrongSign!=0){
1142     wrongMassD0 = theD0particle->InvMassD0();
1143   }else if(okD0WrongSign==0){
1144     wrongMassD0 = theD0particle->InvMassD0bar();
1145   }
1146
1147   if(TMath::Abs(wrongMassD0-1.865)<fD0Window){
1148
1149     //PID of D0 daughters
1150     AliAODTrack *pos = (AliAODTrack*)part->Get2Prong()->GetDaughter(0);
1151     AliAODTrack *neg = (AliAODTrack*)part->Get2Prong()->GetDaughter(1);
1152
1153     Bool_t isPID = kTRUE;
1154
1155
1156     if(fAnalysis==1 && ptbin==1){
1157       if(part->Get2Prong()->DecayLength()>0.08 || part->Get2Prong()->DecayLength()<0.022)return; //(0.05)
1158     }
1159     if(fAnalysis==1 && ptbin==2){
1160       if(part->Get2Prong()->DecayLength()>0.08 || part->Get2Prong()->DecayLength()<0.017)return;
1161     }
1162
1163     if (PIDon) {
1164       if(fDebug > 1) printf("AnalysisTaskSEDStar::TPCPIDon \n");
1165       if(fDebug > 1) printf("AnalysisTaskSEDStar::NSigmaTPC: %d\n", nSigma);
1166
1167       if (part->Charge()>0){
1168         if(!SelectPID(pos, AliPID::kPion, nSigma)) return;//pion+
1169         if(!SelectPID(neg, AliPID::kKaon, nSigma)) return;//kaon-
1170       }else{
1171         if(!SelectPID(pos, AliPID::kKaon, nSigma)) return;//kaon+
1172         if(!SelectPID(neg, AliPID::kPion, nSigma)) return;//pion-
1173       }
1174       if (ptbin>2){
1175         isPID =SelectTOFPID(theD0particle, theSoftPi);
1176         if(!isPID) return;
1177       }
1178     }
1179
1180     // wrong D* inv mass
1181     Double_t e[3];
1182     if (part->Charge()>0){
1183       e[0]=theD0particle->EProng(0,321);
1184       e[1]=theD0particle->EProng(1,211);
1185     }else{
1186       e[0]=theD0particle->EProng(0,211);
1187       e[1]=theD0particle->EProng(1,321);
1188     }
1189     e[2]=part->EProng(0,211);
1190
1191     Double_t esum = e[0]+e[1]+e[2];
1192     Double_t pds = part->P();
1193
1194     Double_t   wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);
1195
1196     TString fillthis="";
1197     fillthis="histWrongSignMass_";
1198     fillthis+=ptbin;
1199     ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1200     fillthis="histWrongSignMass";
1201     ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1202
1203   }
1204 }
1205
1206 //_____________________________________________SINGLE TRACK PRE-SELECTION___________________________________________
1207 Bool_t AliAnalysisTaskSEDStarSpectra::SingleTrackSelections(const AliAODRecoDecayHF2Prong* theD0particle, const AliAODTrack *track2){
1208
1209   // Preselection on D0 daughters and the soft pion
1210
1211   // cut in acceptance for the soft pion and for the D0 daughters
1212   Bool_t acceptanceProng0 = (TMath::Abs(theD0particle->EtaProng(0))<= 0.9);
1213   Bool_t acceptanceProng1 = (TMath::Abs(theD0particle->EtaProng(1))<= 0.9);
1214   // soft pion acceptance ... is it fine 0.9?????
1215   Bool_t acceptanceProng2 = (TMath::Abs(track2->Eta())<= 1.0);
1216
1217   if (!(acceptanceProng0 && acceptanceProng1 && acceptanceProng2)) return kFALSE;
1218   AliDebug(2,"D* reco daughters in acceptance");
1219
1220   return kTRUE;
1221 }
1222
1223 //_____________________________ pid _______________________________________-
1224 Bool_t AliAnalysisTaskSEDStarSpectra::SelectPID(const AliAODTrack *track,  AliPID::EParticleType type, Double_t nsig){//type(0-4): {e,mu,pi,K,p}
1225   //
1226   // Method to extract the PID for the pion/kaon. The particle type for PID can be set by user
1227   // At the moment only TPC PID.
1228   //
1229
1230   //TPC
1231
1232   Bool_t isParticle=kTRUE;
1233   if ((track->GetStatus()&AliESDtrack::kTPCpid )==0) return isParticle;
1234   AliAODPid *pid = track->GetDetPid();
1235   static AliTPCPIDResponse theTPCpid;
1236   Double_t nsigma = theTPCpid.GetNumberOfSigmas(track->P(),pid->GetTPCsignal(),track->GetTPCClusterMap().CountBits(), type);
1237   if (TMath::Abs(nsigma)>nsig) isParticle=kFALSE;
1238
1239   return isParticle;
1240 }
1241 //-------------------------------------------------
1242 Bool_t AliAnalysisTaskSEDStarSpectra::SelectTOFPID(const AliAODRecoDecayHF2Prong* d, const AliAODTrack *tracksoft){
1243
1244
1245   // ######### SPECIAL PID CUTS #################################
1246   Int_t isKaon[2]={0,0};
1247   Bool_t isD0D0barPID[2]={kTRUE,kTRUE};
1248   for(Int_t daught=0;daught<2;daught++){
1249
1250
1251     AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
1252     // AliESDtrack *esdtrack=new
1253     AliESDtrack((AliVTrack*)d->GetDaughter(daught));
1254     if(!(aodtrack->GetStatus()&AliESDtrack::kTOFpid)){
1255       isKaon[daught]=0;
1256       //delete esdtrack;
1257       return kTRUE;
1258     }
1259     if(!(aodtrack->GetStatus()&AliESDtrack::kTOFout)){
1260       isKaon[daught]=0;
1261       //        delete esdtrack;
1262       return kTRUE;
1263     }
1264     if(!(aodtrack->GetStatus()&AliESDtrack::kTIME)){
1265       isKaon[daught]=0;
1266       //delete esdtrack;
1267       return kTRUE;
1268     }
1269     if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
1270       isKaon[daught]=0;
1271       //        delete esdtrack;
1272       return kTRUE;
1273     }
1274     if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
1275       isKaon[daught]=0;
1276       //        delete esdtrack;
1277       return kTRUE;
1278     }
1279
1280     AliAODPid *pid=aodtrack->GetDetPid();
1281     if(!pid) {
1282       isKaon[daught]=0;
1283       //delete esdtrack;
1284       return kTRUE;
1285     }
1286     Double_t tofSig=pid->GetTOFsignal();
1287
1288     Double_t times[5];
1289     //      esdtrack->GetIntegratedTimes(times);
1290     pid->GetIntegratedTimes(times);
1291     //fHistCheck->Fill(esdtrack->P(),esdtrack->GetTOFsignal()-times[3]); //
1292     //3 is the kaon
1293
1294     //  printf("Test momentum VS time %f, %f \n",aodtrack->P(),tofSig-times[3]);
1295     if(TMath::Abs(tofSig-times[3])>3.*160.){
1296       isKaon[daught]=2;
1297       if(aodtrack->Charge()==-1){
1298         isD0D0barPID[0]=kFALSE;
1299       }
1300       else isD0D0barPID[1]=kFALSE;
1301     }
1302     else {
1303       isKaon[daught]=1;
1304
1305       if(aodtrack->P()<1.5){
1306         if(aodtrack->Charge()==-1){
1307           isD0D0barPID[1]=kFALSE;
1308         }
1309         else isD0D0barPID[0]=kFALSE;
1310
1311       }
1312       //delete esdtrack;
1313     }
1314   }
1315
1316   Double_t psCharge = tracksoft->Charge();
1317
1318   if(psCharge>0 &&  !isD0D0barPID[0]) return kFALSE;
1319   if(psCharge<0 &&  !isD0D0barPID[1]) return kFALSE;
1320
1321   return kTRUE;
1322 }