PWG3/vertexingHF/AliAnalysisTaskSEDStarJets.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 in Jets Analysis
  18 //
  19 //  The D* (+ and -) is reconstructed inside jets. Two different cuts are
  20 //  implemented:
  21 //
  22 //  1) C.Ivan D* cuts adapted for correlation analysis
  23 //  2) Topological cut enforcing the correlation D0-softPion pt + relaxed
  24 //     CosThetaStar. This second should be better for correlations.
  25 //
  26 //  USAGE:
  27 //
  28 //  The analysis is performed separately for D*+ and D*-. A flag in the .C is
  29 //  taking care to decide which analysis.
  30 //
  31 //  The cut number 2 can be activaded with a flag in the .C (not active in this version)
  32 //  Cuts 1) are the default.
  33 //
  34 //  The task requires reconstructed jets in the AODs
  35 //
  36 //-----------------------------------------------------------------------
  37 //                         Author A.Grelli
  38 //              ERC-QGP Utrecht University - a.grelli@uu.nl
  39 //-----------------------------------------------------------------------
  40
  41 #include <TDatabasePDG.h>
  42 #include <TParticle.h>
  43 #include <TVector3.h>
  44
  45 #include "AliAnalysisTaskSEDStarJets.h"
  46 #include "AliStack.h"
  47 #include "AliMCEvent.h"
  48 #include "AliAODMCHeader.h"
  49 #include "AliAnalysisManager.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 "AliAODRecoDecayHF2Prong.h"
  57 #include "AliESDtrack.h"
  58 #include "AliAODMCParticle.h"
  59
  60 ClassImp(AliAnalysisTaskSEDStarJets)
  61
  62 //__________________________________________________________________________
  63 AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets() :
  64   AliAnalysisTaskSE(),
  65   fCountReco(0),
  66   fCountRecoAcc(0),
  67   fCountRecoITSClusters(0),
  68   fCountRecoPPR(0),
  69   fCountDStar(0),
  70   fEvents(0),
  71   fMinITSClusters(0),
  72   fComputeD0(kTRUE),
  73   fUseMCInfo(kTRUE),
  74   ftopologicalCut(kFALSE),
  75   fRequireNormalization(kTRUE),
  76   fLorentzTrack1(0,0,0,0),
  77   fLorentzTrack2(0,0,0,0),
  78   fLorentzTrack3(0,0,0,0),
  79   fLorentzTrack4(0,0,0,0),
  80   fOutput(0),
  81   fD0ptvsSoftPtSignal(0),
  82   fD0ptvsSoftPt(0),
  83   ftrigger(0),
  84   fPtPion(0),
  85   fInvMass(0),
  86   fRECOPtDStar(0),
  87   fDStar(0),
  88   fDiff(0),
  89   fDiffSideBand(0),
  90   fDStarMass(0),
  91   fPhi(0),
  92   fPhiBkg(0),
  93   fTrueDiff(0),
  94   fResZ(0),
  95   fResZBkg(0),
  96   fcharmpt(0),
  97   fdstarE(0),
  98   fEjet(0),
  99   fPhijet(0),
 100   fEtaJet(0),
 101   fdstarpt(0)
 102 {
 103   //
 104   // Default ctor
 105   //
 106 }
 107 //___________________________________________________________________________
 108 AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name) :
 109   AliAnalysisTaskSE(name),
 110   fCountReco(0),
 111   fCountRecoAcc(0),
 112   fCountRecoITSClusters(0),
 113   fCountRecoPPR(0),
 114   fCountDStar(0),
 115   fEvents(0),
 116   fMinITSClusters(0),
 117   fComputeD0(kTRUE),
 118   fUseMCInfo(kTRUE),
 119   ftopologicalCut(kFALSE),
 120   fRequireNormalization(kTRUE),
 121   fLorentzTrack1(0,0,0,0),
 122   fLorentzTrack2(0,0,0,0),
 123   fLorentzTrack3(0,0,0,0),
 124   fLorentzTrack4(0,0,0,0),
 125   fOutput(0),
 126   fD0ptvsSoftPtSignal(0),
 127   fD0ptvsSoftPt(0),
 128   ftrigger(0),
 129   fPtPion(0),
 130   fInvMass(0),
 131   fRECOPtDStar(0),
 132   fDStar(0),
 133   fDiff(0),
 134   fDiffSideBand(0),
 135   fDStarMass(0),
 136   fPhi(0),
 137   fPhiBkg(0),
 138   fTrueDiff(0),
 139   fResZ(0),
 140   fResZBkg(0),
 141   fcharmpt(0),
 142   fdstarE(0),
 143   fEjet(0),
 144   fPhijet(0),
 145   fEtaJet(0),
 146   fdstarpt(0)
 147 {
 148   //
 149   // Constructor. Initialization of Inputs and Outputs
 150   //
 151   Info("AliAnalysisTaskSEDStarJets","Calling Constructor");
 152
 153   DefineOutput(1,TList::Class());
 154 }
 155
 156 //___________________________________________________________________________
 157 AliAnalysisTaskSEDStarJets& AliAnalysisTaskSEDStarJets::operator=(const AliAnalysisTaskSEDStarJets& c)
 158 {
 159   //
 160   // Assignment operator
 161   //
 162   if (this!=&c) {
 163     AliAnalysisTaskSE::operator=(c) ;
 164   }
 165
 166   return *this;
 167 }
 168
 169 //___________________________________________________________________________
 170 AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const AliAnalysisTaskSEDStarJets& c) :
 171   AliAnalysisTaskSE(c),
 172   fCountReco(c.fCountReco),
 173   fCountRecoAcc(c.fCountRecoAcc),
 174   fCountRecoITSClusters(c.fCountRecoITSClusters),
 175   fCountRecoPPR(c.fCountRecoPPR),
 176   fCountDStar(c.fCountDStar),
 177   fEvents(c.fEvents),
 178   fMinITSClusters(c.fMinITSClusters),
 179   fComputeD0(c.fComputeD0),
 180   fUseMCInfo(c.fUseMCInfo),
 181   ftopologicalCut(c.ftopologicalCut),
 182   fRequireNormalization(c.fRequireNormalization),
 183   fLorentzTrack1(c.fLorentzTrack1),
 184   fLorentzTrack2(c.fLorentzTrack2),
 185   fLorentzTrack3(c.fLorentzTrack3),
 186   fLorentzTrack4(c.fLorentzTrack4),
 187   fOutput(c.fOutput),
 188   fD0ptvsSoftPtSignal(c.fD0ptvsSoftPtSignal),
 189   fD0ptvsSoftPt(c.fD0ptvsSoftPt),
 190   ftrigger(c.ftrigger),
 191   fPtPion(c.fPtPion),
 192   fInvMass(c.fInvMass),
 193   fRECOPtDStar(c.fRECOPtDStar),
 194   fDStar(c.fDStar),
 195   fDiff(c.fDiff),
 196   fDiffSideBand(c.fDiffSideBand),
 197   fDStarMass(c.fDStarMass),
 198   fPhi(c.fPhi),
 199   fPhiBkg(c.fPhiBkg),
 200   fTrueDiff(c.fTrueDiff),
 201   fResZ(c.fResZ),
 202   fResZBkg(c.fResZBkg),
 203   fcharmpt(c.fcharmpt),
 204   fdstarE(c.fdstarE),
 205   fEjet(c.fEjet),
 206   fPhijet(c.fPhijet),
 207   fEtaJet(c.fEtaJet),
 208   fdstarpt(c.fdstarpt)
 209
 210 {
 211   //
 212   // Copy Constructor
 213   //
 214 }
 215
 216 //___________________________________________________________________________
 217 AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() {
 218   // destructor
 219
 220   Info("~AliAnalysisTaskSEDStarJets","Calling Destructor");
 221   if (fOutput) {
 222     delete fOutput;
 223     fOutput = 0;
 224   }
 225 }
 226
 227 //_________________________________________________
 228 void AliAnalysisTaskSEDStarJets::UserExec(Option_t *)
 229 {
 230   // user exec
 231   if (!fInputEvent) {
 232     Error("UserExec","NO EVENT FOUND!");
 233     return;
 234   }
 235
 236   // Load the event
 237   fEvents++;
 238   AliDebug(2,Form("Event %d",fEvents));
 239   if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
 240   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 241
 242   TClonesArray *arrayVerticesHF=0;
 243
 244   if(!aodEvent && AODEvent() && IsStandardAOD()) {
 245     // In case there is an AOD handler writing a standard AOD, use the AOD
 246     // event in memory rather than the input (ESD) event.
 247     aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 248     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 249     // have to taken from the AOD event hold by the AliAODExtension
 250     AliAODHandler* aodHandler = (AliAODHandler*)
 251       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 252     if(aodHandler->GetExtensions()) {
 253       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 254       AliAODEvent *aodFromExt = ext->GetAOD();
 255       arrayVerticesHF=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 256     }
 257   } else {
 258     arrayVerticesHF=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
 259   }
 260
 261   if (!arrayVerticesHF){
 262     AliInfo("Could not find array of HF vertices, skipping the event");
 263     return;
 264   }else AliDebug(2, Form("Found %d vertices",arrayVerticesHF->GetEntriesFast()));
 265
 266   // Simulate a jet triggered sample
 267   TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets();
 268   if(aodEvent->GetNJets()<=0) return;
 269   AliInfo("found a jet: processing D* in jet analysis");
 270
 271   // AOD primary vertex
 272   AliAODVertex *prVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 273   Double_t primaryPos[3];
 274   prVtx->GetXYZ(primaryPos);
 275   Bool_t vtxFlag = kTRUE;
 276   TString title=prVtx->GetTitle();
 277   if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
 278
 279   // counters for efficiencies
 280   Int_t icountReco = 0;
 281   Int_t icountRecoAcc = 0;
 282   Int_t icountRecoITSClusters = 0;
 283   Int_t icountRecoPPR = 0;
 284   Int_t fiDstar    = 0;
 285   Int_t fDStarD0   = 0;
 286
 287   // TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 288
 289   if(fUseMCInfo){
 290     TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 291     //loop on the MC event - some basic MC info on D*, D0 and soft pion
 292     if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 293
 294     for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 295       AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 296       if (!mcPart) {
 297         AliWarning("Particle not found in tree, skipping");
 298         continue;
 299       }
 300
 301       // charm pt
 302       if(TMath::Abs(mcPart->GetPdgCode())==4){
 303         fcharmpt->Fill(mcPart->Pt());
 304       }
 305
 306       // fill energy and pt for D* in acceptance with correct prongs
 307       Bool_t isOk = DstarInMC(mcPart,mcArray);
 308
 309       if (isOk){ //D*
 310         AliDebug(2, "Found a DStar in MC with correct prongs and in acceptance");
 311         fdstarE ->Fill(mcPart->E());
 312         fdstarpt->Fill(mcPart->Pt());
 313
 314         // check the MC-Acceptance level cuts
 315         // since standard CF functions are not applicable, using Kine Cuts on daughters
 316
 317         Int_t daughter0 = mcPart->GetDaughter(0);
 318         Int_t daughter1 = mcPart->GetDaughter(1);
 319
 320         AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
 321
 322         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
 323         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
 324
 325         Double_t eta0 = mcPartDaughter0->Eta();
 326         Double_t eta1 = mcPartDaughter1->Eta();
 327         Double_t y0   = mcPartDaughter0->Y();
 328         Double_t y1   = mcPartDaughter1->Y();
 329         Double_t pt0  = mcPartDaughter0->Pt();
 330         Double_t pt1  = mcPartDaughter1->Pt();
 331
 332         AliDebug(2, Form("D* Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
 333         AliDebug(2, Form("D* Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
 334
 335         Int_t daughD00 = 0;
 336         Int_t daughD01 = 0;
 337
 338         // D0 daughters - do not need to check D0-kpi, already done
 339
 340         if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){
 341           daughD00 = mcPartDaughter0->GetDaughter(0);
 342           daughD01 = mcPartDaughter0->GetDaughter(1);
 343         }else{
 344           daughD00 = mcPartDaughter1->GetDaughter(0);
 345           daughD01 = mcPartDaughter1->GetDaughter(1);
 346         }
 347
 348         AliAODMCParticle* mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00));
 349         AliAODMCParticle* mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01));
 350
 351         if (!mcD0PartDaughter0 || !mcD0PartDaughter1) {
 352           AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done...");
 353         }
 354
 355         // D0 daughters - needed for acceptance
 356         Double_t pD0pt0 =  mcD0PartDaughter0->Pt();
 357         Double_t pD0pt1 =  mcD0PartDaughter1->Pt();
 358         Double_t pD0eta0 = mcD0PartDaughter0->Eta();
 359         Double_t pD0eta1 = mcD0PartDaughter1->Eta();
 360
 361         // ACCEPTANCE REQUESTS ---------
 362
 363         // soft pion
 364         Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 > 0.05);
 365         // Do daughters
 366         Bool_t theD0daught0inAcceptance = (TMath::Abs(pD0eta0) <= 0.9 && pD0pt0 >= 0.1);
 367         Bool_t theD0daught1inAcceptance = (TMath::Abs(pD0eta1) <= 0.9 && pD0pt1 >= 0.1);
 368
 369         if (daught1inAcceptance && theD0daught0inAcceptance && theD0daught1inAcceptance) {
 370
 371           AliDebug(2, "Daughter particles in acceptance");
 372
 373           // check on the vertex
 374           if (vtxFlag){
 375             printf("Vertex cut passed 2\n");
 376             fDStarD0++;
 377             Bool_t refitFlag = kTRUE;
 378             for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
 379               AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
 380
 381               // refit only for D0 daughters
 382               if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) {
 383                 if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
 384                   refitFlag = kFALSE;
 385                 }
 386               }
 387             }
 388             if (refitFlag){
 389               printf("Refit cut passed\n");
 390             }
 391             else{
 392               AliDebug(3,"Refit cut not passed\n");
 393             }
 394           }
 395           else{
 396             AliDebug(3,"Vertex cut not passed\n");
 397           }
 398         }
 399         else{
 400           AliDebug(3,"Acceptance cut not passed\n");
 401         }
 402       }
 403     }
 404
 405     AliDebug(2, Form("Found %i MC particles that are D* in Kpipi and satisfy Acc cuts!!",fDStarD0));
 406     fCountDStar += fDStarD0;
 407
 408   } // End of MC
 409
 410   // Now perform the D* in jet reconstruction
 411
 412   // Normalization factor
 413   if(fRequireNormalization){
 414     ftrigger->Fill(1);
 415   }
 416
 417   Int_t nJets = 0; // for reco D0 countings
 418
 419   for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) { // main loop on jets
 420
 421     AliAODJet* jet = (AliAODJet*)arrayofJets->At(iJets);
 422
 423     //jets variables
 424     Double_t ejet   = jet->E();
 425     Double_t phiJet = jet->Phi();
 426     Double_t etaJet = jet->Eta();
 427
 428     // fill energy, eta and phi of the jet
 429     fEjet   ->Fill(ejet);
 430     fPhijet ->Fill(phiJet);
 431     fEtaJet ->Fill(etaJet);
 432
 433     nJets++;
 434
 435     // loop over the tracks to search for candidates soft pions
 436     for (Int_t i=0; i<aodEvent->GetNTracks(); i++){
 437
 438       AliAODTrack* aodTrack = aodEvent->GetTrack(i);
 439       Double_t vD0[4] = {0.,0.,0.,0.};
 440       Double_t vD0bar[4] ={0.,0.,0.,0.};
 441
 442       Int_t pCharge = aodTrack->Charge();
 443
 444       // few selections on soft pion
 445       Bool_t tPrimCand =  aodTrack->IsPrimaryCandidate(); // is it primary?
 446
 447       if(aodTrack->Pt()>= 5 || aodTrack->Pt()< 0.08 ) continue; //cut on soft pion pt VERY large
 448                                                                   //~ D*s of pt >80GeV with a soft pion of 5GeV!
 449       if(TMath::Abs(aodTrack->Eta())>0.9) continue;
 450
 451       // if D*+ analysis tha D0 and pi+  otherwise pi-
 452       if(fComputeD0 && pCharge!= 1 ) continue;
 453       if(!fComputeD0 && pCharge!= -1 ) continue;
 454
 455       if(tPrimCand && arrayVerticesHF->GetEntriesFast()>0){ // isPion and is Primary, no PID for now
 456
 457         // label to the candidate soft pion
 458         Int_t pLabel = -1;
 459         if(fUseMCInfo) pLabel = aodTrack->GetLabel();
 460
 461         // prepare the TLorentz vector for the pion
 462         Float_t pionMass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 463         fLorentzTrack3.SetPxPyPzE(aodTrack->Px(),aodTrack->Py(),aodTrack->Pz(),aodTrack->E(pionMass));
 464
 465         // search for the D0
 466         for (Int_t iVertex = 0; iVertex<arrayVerticesHF->GetEntriesFast(); iVertex++) {
 467
 468           Double_t invM      = 0;
 469           Double_t invMDStar = 0;
 470           Double_t dPhi      = 0;
 471
 472           AliAODRecoDecayHF2Prong* vtx = (AliAODRecoDecayHF2Prong*)arrayVerticesHF->At(iVertex);
 473
 474           Double_t pt = vtx->Pt();
 475
 476           // acceptance variables
 477           Bool_t acceptanceProng0 = (TMath::Abs(vtx->EtaProng(0))<= 0.9 && vtx->PtProng(0) >= 0.1);
 478           Bool_t acceptanceProng1 = (TMath::Abs(vtx->EtaProng(1))<= 0.9 && vtx->PtProng(1) >= 0.1);
 479           Bool_t acceptanceSoftPi = (TMath::Abs(aodTrack->Eta())<= 0.9 && aodTrack->Pt() >= 0.05);
 480
 481           Int_t pdgDgD0toKpi[2]={321,211};
 482
 483           Int_t mcLabel =-1;
 484
 485           Bool_t isDStar = kFALSE; // just to count
 486
 487           //switch of MC for real data
 488           if(fUseMCInfo){
 489             TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 490             mcLabel = vtx->MatchToMC(421, mcArray,2,pdgDgD0toKpi) ;   //MC D0
 491             // matching to MC D*
 492             if(mcLabel !=-1 && pLabel!=-1 && nJets ==1) { // count only once in case of multijets
 493
 494               // search for a D0 and a pi with mother and check it is a D*
 495               AliAODMCParticle* theMCpion = (AliAODMCParticle*)mcArray->At(pLabel);
 496               Int_t motherMCPion = theMCpion->GetMother();
 497               AliAODMCParticle* theMCD0 = (AliAODMCParticle*)mcArray->At(mcLabel);
 498               Int_t motherMCD0 = theMCD0->GetMother();
 499
 500               if(motherMCPion!=-1 && (motherMCD0 == motherMCPion)){
 501                 AliAODMCParticle* mcMother = (AliAODMCParticle*)mcArray->At(motherMCPion);
 502                 if(TMath::Abs(mcMother->GetPdgCode()) == 413 && TMath::Abs(theMCpion->GetPdgCode())==211){
 503                   isDStar = kTRUE;
 504                   fiDstar++;
 505                 }
 506               }
 507             }
 508           }
 509
 510
 511           if (acceptanceProng0 && acceptanceProng1 && acceptanceSoftPi) {
 512
 513             AliDebug(2,"D0 reco daughters in acceptance");
 514             if(isDStar && nJets==1) icountRecoAcc++;
 515
 516             // D0 daughters
 517             AliAODTrack *track0 = (AliAODTrack*)vtx->GetDaughter(0);
 518             AliAODTrack *track1 = (AliAODTrack*)vtx->GetDaughter(1);
 519
 520             // check for ITS refit (already required at the ESD filter level )
 521             Bool_t kRefitITS = kTRUE;
 522
 523             if((!(track0->GetStatus()&AliESDtrack::kITSrefit)|| (!(track1->GetStatus()&AliESDtrack::kITSrefit)))) {
 524               kRefitITS = kFALSE;
 525             }
 526
 527             Int_t ncls0=0,ncls1=0,ncls2=0;
 528
 529             for(Int_t l=0;l<6;l++) {
 530               if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
 531               if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
 532               if(TESTBIT(aodTrack->GetITSClusterMap(),l)) ncls2++;
 533             }
 534
 535             // clusters in ITS for D0 daugthers and soft pion
 536             if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>=4 && kRefitITS) {
 537
 538               if(isDStar && nJets==1) icountRecoITSClusters++;
 539
 540               // D0 cutting varibles
 541               Double_t cosThetaStar = 9999.;
 542               Double_t pTpi = 0.;
 543               Double_t pTK = 0.;
 544               Double_t d01 = 0;
 545               Double_t d00 = 0;
 546
 547               // D0, D0bar
 548               if (fComputeD0){
 549                 cosThetaStar = vtx->CosThetaStarD0();
 550                 pTpi = vtx->PtProng(0);
 551                 pTK =  vtx->PtProng(1);
 552                 d01 =   vtx->Getd0Prong(0)*1E4;
 553                 d00 =   vtx->Getd0Prong(1)*1E4;
 554               }else{
 555                 cosThetaStar = vtx->CosThetaStarD0bar();
 556                 pTpi = vtx->PtProng(1);
 557                 pTK =  vtx->PtProng(0);
 558                 d01 =   vtx->Getd0Prong(1)*1E4;
 559                 d00 =   vtx->Getd0Prong(0)*1E4;
 560               }
 561
 562               AliDebug(2,"D0 reco daughters in acceptance");
 563
 564               Double_t dca =   vtx->GetDCA()*1E4;
 565               Double_t d0d0 =  vtx->Prodd0d0()*1E8;
 566
 567               TVector3 mom(vtx->Px(),vtx->Py(),vtx->Pz());
 568               TVector3 flight((vtx->Xv())- primaryPos[0],(vtx->Yv())- primaryPos[1],(vtx->Zv())- primaryPos[2]);
 569               Double_t pta = mom.Angle(flight);
 570               Double_t cosPointingAngle = TMath::Cos(pta);
 571
 572               // Crstian Ivan D* cuts. Multidimensional optimization
 573               Double_t cuts[7] = {9999999., 1.1, 0., 9999999.,999999., 9999999., 0.};
 574
 575               if (pt <= 1){ // first bin not optimized
 576                 cuts[0] = 400;
 577                 cuts[1] = 0.8;
 578                 cuts[2] = 0.21;
 579                 cuts[3] = 500;
 580                 cuts[4] = 500;
 581                 cuts[5] = -20000;
 582                 cuts[6] = 0.6;
 583               }
 584               else if (pt > 1 && pt <= 2){
 585                 cuts[0] = 200;
 586                 cuts[1] = 0.7;
 587                 cuts[2] = 0.8;
 588                 cuts[3] = 210;
 589                 cuts[4] = 210;
 590                 cuts[5] = -20000;
 591                 cuts[6] = 0.9;
 592               }
 593               else if (pt > 2 && pt <= 3){
 594                 cuts[0] = 400;
 595                 cuts[1] = 0.8;
 596                 cuts[2] = 0.8;
 597                 cuts[3] = 420;
 598                 cuts[4] = 350;
 599                 cuts[5] = -8500;
 600                 cuts[6] = 0.9;
 601               }
 602               else if (pt > 3 && pt <= 5){
 603                 cuts[0] = 160;
 604                 cuts[1] = 1.0;
 605                 cuts[2] = 1.2;
 606                 cuts[3] = 420;
 607                 cuts[4] = 560;
 608                 cuts[5] = -8500;
 609                 cuts[6] = 0.9;
 610               }
 611               else if (pt > 5){
 612                 cuts[0] = 800;
 613                 cuts[1] = 1.0;
 614                 cuts[2] = 1.2;
 615                 cuts[3] = 700;
 616                 cuts[4] = 700;
 617                 cuts[5] = 10000;
 618                 cuts[6] = 0.9;
 619               }
 620               // apply D0 cuts
 621
 622               if (dca < cuts[0]
 623                   && TMath::Abs(cosThetaStar) < cuts[1]
 624                   && pTpi > cuts[2]
 625                   && pTK > cuts[2]
 626                   && TMath::Abs(d01) < cuts[3]
 627                   && TMath::Abs(d00) < cuts[4]
 628                   && d0d0 < cuts[5]
 629                   && cosPointingAngle > cuts[6]
 630                 ){
 631
 632                 if(fComputeD0){ // D0 from default
 633
 634                   if(vtx->ChargeProng(0)==-1){
 635                     fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,321) );
 636                     fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,211) );
 637                   }else{
 638                     fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,211) );
 639                     fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,321) );
 640                   }
 641
 642                   vD0[0] =  (fLorentzTrack1+fLorentzTrack2).Px();
 643                   vD0[1] =  (fLorentzTrack1+fLorentzTrack2).Py();
 644                   vD0[2] =  (fLorentzTrack1+fLorentzTrack2).Pz();
 645                   vD0[3] =  (fLorentzTrack1+fLorentzTrack2).E();
 646
 647                   fLorentzTrack4.SetPxPyPzE(vD0[0],vD0[1],vD0[2],vD0[3]);
 648
 649                 }else{ // D0bar analysis
 650
 651                   if(vtx->ChargeProng(0)==-1){
 652                     fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,211) );
 653                     fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,321) );
 654                   }else{
 655                     fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,321) );
 656                     fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,211) );
 657                   }
 658
 659                   vD0bar[0] = (fLorentzTrack1+fLorentzTrack2).Px();
 660                   vD0bar[1] = (fLorentzTrack1+fLorentzTrack2).Py();
 661                   vD0bar[2] = (fLorentzTrack1+fLorentzTrack2).Pz();
 662                   vD0bar[3] = (fLorentzTrack1+fLorentzTrack2).E();
 663
 664                   fLorentzTrack4.SetPxPyPzE(vD0bar[0],vD0bar[1],vD0bar[2],vD0bar[3]);
 665                 }
 666
 667                 // D0-D0bar related variables
 668
 669                 invM = GetInvariantMass(fLorentzTrack1,fLorentzTrack2);
 670                 if(nJets==1) fInvMass->Fill(invM);
 671
 672                 if(isDStar && nJets==1) icountRecoPPR++;
 673
 674                 //DStar invariant mass
 675                 invMDStar = GetInvariantMassDStar(fLorentzTrack3,fLorentzTrack4);
 676
 677                 //conversion from phi TLorentzVerctor to phi aliroot
 678                 Double_t kconvert = 0;
 679                 Double_t phiDStar = (fLorentzTrack3 + fLorentzTrack4).Phi();
 680                 kconvert = phiDStar;
 681                 if(phiDStar<0) kconvert = phiDStar + 2*(TMath::Pi());
 682                 phiDStar = kconvert;
 683
 684                 // dphi between jet and D*
 685                 dPhi = phiJet - phiDStar;
 686
 687                 //Just for plotting pourposal
 688                 if(dPhi<=-1.58) dPhi = TMath::Abs(dPhi);
 689                 if(dPhi>4.72){
 690                   dPhi = dPhi-2*(TMath::Pi());
 691                 }
 692
 693                 //Alternative cutting procedure
 694                 //the cut on cosThetaStar is to reduce near collinear
 695                 //background from jet fragmentation
 696
 697                 Bool_t tCut = EvaluateCutOnPiD0pt(vtx,aodTrack);
 698
 699                 if(ftopologicalCut && tCut) continue;
 700                 if(ftopologicalCut && TMath::Abs(cosThetaStar)>cuts[1]) continue;
 701
 702                 if(invM >= 1.829 && invM <= 1.901){ // ~4 sigma cut on D0 mass
 703
 704                   if(isDStar && nJets==1) {
 705                     fDStarMass->Fill(invMDStar);
 706                     fTrueDiff->Fill(invMDStar-invM);
 707                   }
 708
 709                   if(nJets==1) {  // avoid double counting
 710                     fDStar->Fill(invMDStar);
 711                     fDiff->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi)
 712                   }
 713
 714                   // now the dphi signal and the fragmentation function
 715                   if((invMDStar-invM)<=0.148 && (invMDStar-invM)>=0.142) {
 716
 717                     //fill candidates D* and soft pion reco pt
 718                     if(nJets==1) fPtPion->Fill(aodTrack->Pt());
 719                     if(nJets==1) fRECOPtDStar->Fill((fLorentzTrack3 + fLorentzTrack4).Pt());
 720                     fPhi ->Fill(dPhi);
 721
 722                     if(dPhi>=-0.5 && dPhi<=0.5){  // evaluate in the near side
 723                       Double_t dStarMom = (fLorentzTrack3 + fLorentzTrack4).P();
 724                       Double_t zFrag = (TMath::Cos(dPhi)*dStarMom)/ejet;
 725                       fResZ->Fill(TMath::Abs(zFrag));
 726                     }
 727                   }
 728                 }
 729
 730                 // evaluate side band background
 731                 SideBandBackground(invM, invMDStar, ejet, dPhi, nJets);
 732
 733                 invM      = 0;
 734                 invMDStar = 0;
 735
 736                } // end PPR cuts
 737             } // end ITS cluster
 738           } // end acceptance
 739         } // D0 cand
 740       } // softpion
 741     } // tracks
 742   } // jets
 743
 744   AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
 745
 746   fCountReco+= fiDstar;
 747   fCountRecoAcc+= icountRecoAcc;
 748   fCountRecoITSClusters+= icountRecoITSClusters;
 749   fCountRecoPPR+= icountRecoPPR;
 750
 751   PostData(1,fOutput);
 752 }
 753
 754 //________________________________________ terminate ___________________________
 755 void AliAnalysisTaskSEDStarJets::Terminate(Option_t*)
 756 {
 757   // The Terminate() function is the last function to be called during
 758   // a query. It always runs on the client, it can be used to present
 759   // the results graphically or save the results to file.
 760
 761   Info("Terminate","");
 762   AliAnalysisTaskSE::Terminate();
 763
 764   AliInfo(Form("Found %i of that MC D*->D0pi(D0->kpi) are in acceptance, in %d events",fCountDStar,fEvents));
 765
 766   AliInfo(Form("Found %i RECO particles after cuts that are DStar, in %d events",fCountReco,fEvents));
 767   AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
 768   AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents));
 769   AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and satisfy D* cuts, in %d events",fCountRecoPPR,fEvents));
 770
 771   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 772   if (!fOutput) {
 773     printf("ERROR: fOutput not available\n");
 774     return;
 775   }
 776
 777   fcharmpt      = dynamic_cast<TH1F*>(fOutput->FindObject("fcharmpt"));
 778   fdstarE       = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarE"));
 779   fdstarpt      = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarpt"));
 780   fDStarMass    = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
 781   fTrueDiff     = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
 782   fInvMass      = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
 783   fPtPion       = dynamic_cast<TH1F*>(fOutput->FindObject("fPtPion "));
 784   fDStar        = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar"));
 785   fDiff         = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff"));
 786   fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
 787   ftrigger      = dynamic_cast<TH1F*>(fOutput->FindObject("ftrigger"));
 788   fRECOPtDStar  = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtDStar"));
 789   fEjet         = dynamic_cast<TH1F*>(fOutput->FindObject("fEjet"));
 790   fPhijet       = dynamic_cast<TH1F*>(fOutput->FindObject("fPhijet"));
 791   fEtaJet       = dynamic_cast<TH1F*>(fOutput->FindObject("fEtaJet"));
 792   fPhi          = dynamic_cast<TH1F*>(fOutput->FindObject("fPhi"));
 793   fResZ         = dynamic_cast<TH1F*>(fOutput->FindObject("fResZ"));
 794   fResZBkg      = dynamic_cast<TH1F*>(fOutput->FindObject("fResZBkg"));
 795   fPhiBkg       = dynamic_cast<TH1F*>(fOutput->FindObject("fPhiBkg"));
 796   fD0ptvsSoftPtSignal = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPtSignal"));
 797   fD0ptvsSoftPt       = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPt"));
 798
 799 }
 800 //___________________________________________________________________________
 801
 802 void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() {
 803  // output
 804
 805   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 806
 807   //slot #1
 808   OpenFile(1);
 809   fOutput = new TList();
 810   fOutput->SetOwner();
 811   // define histograms
 812   DefineHistoFroAnalysis();
 813
 814   return;
 815 }
 816
 817 //_______________________________D0 invariant mass________________________________
 818
 819 Double_t AliAnalysisTaskSEDStarJets::GetInvariantMass(TLorentzVector LorentzTrack1, TLorentzVector LorentzTrack2){
 820
 821   return TMath::Abs((LorentzTrack1+LorentzTrack2).M());   // invariant mass of two tracks
 822 }
 823 //______________________________D* invariant mass_________________________________
 824
 825 Double_t AliAnalysisTaskSEDStarJets::GetInvariantMassDStar(TLorentzVector LorentzTrack4, TLorentzVector LorentzTrack3){
 826
 827   return TMath::Abs((LorentzTrack4+LorentzTrack3).M());   // invariant mass of two tracks
 828 }
 829 //_________________________________D* in MC _______________________________________
 830
 831 Bool_t  AliAnalysisTaskSEDStarJets::DstarInMC(AliAODMCParticle* const mcPart, TClonesArray* mcArray){
 832   // D* in MC
 833
 834   Bool_t isOk = kFALSE;
 835
 836   if(TMath::Abs(mcPart->GetPdgCode())!=413) return isOk;
 837
 838   // getting the daughters
 839   Int_t daughter0 = mcPart->GetDaughter(0);
 840   Int_t daughter1 = mcPart->GetDaughter(1);
 841
 842   AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
 843   if (daughter0 == 0 || daughter1 == 0) {
 844     AliDebug(2, "Error! the D* MC doesn't have correct daughters!!");
 845     return isOk;
 846   }
 847
 848   if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
 849     AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
 850     return isOk;
 851   }
 852
 853   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
 854   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
 855   if (!mcPartDaughter0 || !mcPartDaughter1) {
 856     AliWarning("At least one Daughter Particle not found in tree, skipping");
 857     return isOk;
 858   }
 859
 860   if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
 861         TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
 862       !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
 863         TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
 864     AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
 865     return isOk;
 866   }
 867
 868   Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
 869   Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
 870
 871   // getting vertex from daughters
 872   mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
 873   mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
 874
 875   // check if the secondary vertex is the same for both
 876   if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
 877     AliError("Daughters have different secondary vertex, skipping the track");
 878     return isOk;
 879   }
 880
 881   AliAODMCParticle* neutralDaugh = mcPartDaughter0;
 882
 883   if (!EvaluateIfD0toKpi(neutralDaugh,mcArray)) {
 884     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 885     return isOk;
 886   }
 887
 888   isOk = kTRUE;
 889
 890   return isOk;
 891
 892 }
 893
 894 Bool_t AliAnalysisTaskSEDStarJets::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray)const{
 895
 896   //
 897   // chack wether D0 is decaing into kpi
 898   //
 899
 900   Bool_t isHadronic = kFALSE;
 901
 902   Int_t daughterD00 = neutralDaugh->GetDaughter(0);
 903   Int_t daughterD01 = neutralDaugh->GetDaughter(1);
 904
 905   AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
 906   if (daughterD00 == 0 || daughterD01 == 0) {
 907     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 908     return isHadronic;
 909   }
 910
 911   if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
 912     AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
 913     return isHadronic;
 914   }
 915
 916   AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
 917   AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
 918   if (!mcPartDaughterD00 || !mcPartDaughterD01) {
 919     AliWarning("At least one Daughter Particle not found in tree, skipping");
 920     return isHadronic;
 921   }
 922
 923   if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
 924         TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
 925       !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
 926         TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
 927     AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
 928     return isHadronic;
 929   }
 930
 931   isHadronic = kTRUE;
 932
 933
 934   return isHadronic;
 935
 936 }
 937
 938
 939 //___________________________________ hiostograms _______________________________________
 940
 941 Bool_t  AliAnalysisTaskSEDStarJets::DefineHistoFroAnalysis(){
 942
 943   // Invariant mass related histograms
 944   fInvMass = new TH1F("invMass","Kpi invariant mass distribution",1500,.5,3.5);
 945   fInvMass->SetStats(kTRUE);
 946   fInvMass->GetXaxis()->SetTitle("GeV/c");
 947   fInvMass->GetYaxis()->SetTitle("Entries");
 948   fOutput->Add(fInvMass);
 949
 950   fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4);
 951   fDStar->SetStats(kTRUE);
 952   fDStar->GetXaxis()->SetTitle("GeV/c");
 953   fDStar->GetYaxis()->SetTitle("Entries");
 954   fOutput->Add(fDStar);
 955
 956   fDiff = new TH1F("Diff","M(kpipi)-M(kpi)",750,0.1,0.2);
 957   fDiff->SetStats(kTRUE);
 958   fDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV/c^2");
 959   fDiff->GetYaxis()->SetTitle("Entries");
 960   fOutput->Add(fDiff);
 961
 962   fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
 963   fDiffSideBand->SetStats(kTRUE);
 964   fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
 965   fDiffSideBand->GetYaxis()->SetTitle("Entries");
 966   fOutput->Add(fDiffSideBand);
 967
 968   fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
 969   fOutput->Add(fDStarMass);
 970
 971   fTrueDiff  = new TH1F("dstar","True Reco diff",750,0,0.2);
 972   fOutput->Add(fTrueDiff);
 973
 974   // trigger normalization
 975   ftrigger = new TH1F("Normalization","Normalization factor for correlations",1,0,10);
 976   ftrigger->SetStats(kTRUE);
 977   fOutput->Add(ftrigger);
 978
 979   //correlation fistograms
 980   fPhi = new TH1F("phi","Delta phi between Jet axis and DStar ",25,-1.57,4.72);
 981   fPhi->SetStats(kTRUE);
 982   fPhi->GetXaxis()->SetTitle("#Delta #phi (rad)");
 983   fPhi->GetYaxis()->SetTitle("Entries");
 984   fOutput->Add(fPhi);
 985
 986   fPhiBkg = new TH1F("phiBkg","Delta phi between Jet axis and DStar background ",25,-1.57,4.72);
 987   fPhiBkg->SetStats(kTRUE);
 988   fPhiBkg->GetXaxis()->SetTitle("#Delta #phi (rad)");
 989   fPhiBkg->GetYaxis()->SetTitle("Entries");
 990   fOutput->Add(fPhiBkg);
 991
 992   fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",600,0,15);
 993   fRECOPtDStar->SetStats(kTRUE);
 994   fRECOPtDStar->SetLineColor(2);
 995   fRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
 996   fRECOPtDStar->GetYaxis()->SetTitle("Entries");
 997   fOutput->Add(fRECOPtDStar);
 998
 999   fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,5);
1000   fPtPion->SetStats(kTRUE);
1001   fPtPion->GetXaxis()->SetTitle("GeV/c");
1002   fPtPion->GetYaxis()->SetTitle("Entries");
1003   fOutput->Add(fPtPion);
1004
1005   fcharmpt = new TH1F("charmpt","MC Charm pt distribution",    10000,0,5000);
1006   fdstarE  = new TH1F("dstare", "MC D* energy distribution",   10000,0,5000);
1007   fdstarpt = new TH1F("dstarpt","MC D* pt distribution",       10000,0,5000);
1008   fOutput->Add(fcharmpt);
1009   fOutput->Add(fdstarE);
1010   fOutput->Add(fdstarpt);
1011
1012   // jet related fistograms
1013   fEjet      = new TH1F("ejet",  "UA1 algorithm jet energy distribution",1000,0,500);
1014   fPhijet    = new TH1F("Phijet","UA1 algorithm jet #phi distribution",  200,-7,7);
1015   fEtaJet    = new TH1F("Etajet","UA1 algorithm jet #eta distribution",  200,-7,7);
1016   fOutput->Add(fEjet);
1017   fOutput->Add(fPhijet);
1018   fOutput->Add(fEtaJet);
1019
1020   fResZ      = new TH1F("FragFunc","Fragmentation function ",50,0,1);
1021   fResZBkg   = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1);
1022   fOutput->Add(fResZ);
1023   fOutput->Add(fResZBkg);
1024
1025   fD0ptvsSoftPt       = new TH2F("D0piRec","Candidates (background + signal)",100,0,3,100,0,15);
1026   fD0ptvsSoftPtSignal = new TH2F("D0PiSignal","Signal",100,0,3,100,0,15);
1027   fOutput->Add(fD0ptvsSoftPt);
1028   fOutput->Add(fD0ptvsSoftPtSignal);
1029
1030   return kTRUE;
1031
1032 }
1033
1034 //______________________________Phase space cut alternative to PPR _______________________________________
1035
1036 Bool_t AliAnalysisTaskSEDStarJets::EvaluateCutOnPiD0pt(AliAODRecoDecayHF2Prong* const vtx,  AliAODTrack* const aodTrack) {
1037
1038   // The soft pion pt and DO pt are strongly correlated. It can be shown that ~ 95% of the signal in constrained
1039   // into a narrow band defined by  10 < D0pt/softPt < 18. This cut can be used with a relaxed CosThetaStar cut
1040   // to reconstruct the D*.
1041
1042   Double_t softPt   = 0;
1043   Double_t d0ptReco = 0;
1044
1045   softPt = aodTrack->Pt();
1046   d0ptReco = vtx->Pt();
1047   fD0ptvsSoftPt->Fill(softPt,d0ptReco);
1048
1049   if(softPt>0){
1050     Double_t ratio = d0ptReco/softPt;
1051     if( ratio <=10. || ratio>=18. ) return kFALSE;
1052   }
1053
1054   fD0ptvsSoftPtSignal->Fill(softPt,d0ptReco);
1055   return kTRUE;
1056 }
1057
1058 //______________________________ side band background for D*___________________________________
1059
1060 void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t ejet, Double_t dPhi, Int_t nJets){
1061
1062   //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1063   // (expected detector resolution) on the left and right frm the D0 mass. Each band
1064   //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
1065
1066   if((invM>=1.763 && invM<=1.811) || (invM>=1.919 && invM<=1.963)){
1067
1068     if(nJets==1)fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background
1069
1070     if ((invMDStar-invM)<=0.148 && (invMDStar-invM)>=0.142) {
1071       fPhiBkg->Fill(dPhi);
1072
1073       if(dPhi>=-0.5 && dPhi<=0.5){  // evaluate in the near side
1074
1075         Double_t dStarMomBkg = (fLorentzTrack3 + fLorentzTrack4).P();
1076         Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/ejet;
1077         fResZBkg->Fill(TMath::Abs(zFragBkg));
1078
1079       }
1080     }
1081   }
1082 }