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