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