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