PWGHF/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 //-----------------------------------------------------------------------
  20 //                         Author A.Grelli
  21 //              ERC-QGP Utrecht University - a.grelli@uu.nl
  22 //-----------------------------------------------------------------------
  23
  24 #include <TDatabasePDG.h>
  25 #include <TParticle.h>
  26 #include <TVector3.h>
  27 #include "TROOT.h"
  28
  29 #include "AliAnalysisTaskSEDStarJets.h"
  30 #include "AliRDHFCutsDStartoKpipi.h"
  31 #include "AliStack.h"
  32 #include "AliMCEvent.h"
  33 #include "AliAODMCHeader.h"
  34 #include "AliAODHandler.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliLog.h"
  37 #include "AliAODVertex.h"
  38 #include "AliAODJet.h"
  39 #include "AliAODRecoDecay.h"
  40 #include "AliAODRecoCascadeHF.h"
  41 #include "AliAODRecoDecayHF2Prong.h"
  42 #include "AliESDtrack.h"
  43 #include "AliAODMCParticle.h"
  44
  45 ClassImp(AliAnalysisTaskSEDStarJets)
  46
  47 //__________________________________________________________________________
  48 AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets() :
  49   AliAnalysisTaskSE(),
  50   fEvents(0),
  51   fchargeFrCorr(0),
  52   fUseMCInfo(kTRUE),
  53   fRequireNormalization(kTRUE),
  54   fOutput(0),
  55   fCuts(0),
  56   ftrigger(0),
  57   fPtPion(0),
  58   fInvMass(0),
  59   fRECOPtDStar(0),
  60   fRECOPtBkg(0),
  61   fDStar(0),
  62   fDiff(0),
  63   fDiffSideBand(0),
  64   fDStarMass(0),
  65   fPhi(0),
  66   fPhiBkg(0),
  67   fTrueDiff(0),
  68   fResZ(0),
  69   fResZBkg(0),
  70   fEjet(0),
  71   fPhijet(0),
  72   fEtaJet(0),
  73   theMCFF(0),
  74   fDphiD0Dstar(0),
  75   fPtJet(0)
  76 {
  77   //
  78   // Default ctor
  79   //
  80 }
  81 //___________________________________________________________________________
  82 AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
  83   AliAnalysisTaskSE(name),
  84   fEvents(0),
  85   fchargeFrCorr(0),
  86   fUseMCInfo(kTRUE),
  87   fRequireNormalization(kTRUE),
  88   fOutput(0),
  89   fCuts(0),
  90   ftrigger(0),
  91   fPtPion(0),
  92   fInvMass(0),
  93   fRECOPtDStar(0),
  94   fRECOPtBkg(0),
  95   fDStar(0),
  96   fDiff(0),
  97   fDiffSideBand(0),
  98   fDStarMass(0),
  99   fPhi(0),
 100   fPhiBkg(0),
 101   fTrueDiff(0),
 102   fResZ(0),
 103   fResZBkg(0),
 104   fEjet(0),
 105   fPhijet(0),
 106   fEtaJet(0),
 107   theMCFF(0),
 108   fDphiD0Dstar(0),
 109   fPtJet(0)
 110 {
 111   //
 112   // Constructor. Initialization of Inputs and Outputs
 113   //
 114   fCuts=cuts;
 115   Info("AliAnalysisTaskSEDStarJets","Calling Constructor");
 116
 117   DefineOutput(1,TList::Class()); // histos
 118   DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
 119 }
 120 //___________________________________________________________________________
 121 AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() {
 122   //
 123   // destructor
 124   //
 125
 126   Info("~AliAnalysisTaskSEDStarJets","Calling Destructor");
 127   if (fOutput) {
 128     delete fOutput;
 129     fOutput = 0;
 130   }
 131
 132   if (fCuts) {
 133     delete fCuts;
 134     fCuts = 0;
 135   }
 136
 137   if (ftrigger) { delete ftrigger; ftrigger = 0;}
 138   if (fPtPion)  { delete fPtPion;  fPtPion = 0;}
 139   if (fInvMass) { delete fInvMass; fInvMass = 0;}
 140   if (fRECOPtDStar) { delete fRECOPtDStar; fRECOPtDStar = 0;}
 141   if (fRECOPtBkg)   { delete fRECOPtBkg; fRECOPtBkg = 0;}
 142   if (fDStar) { delete fDStar; fDStar = 0;}
 143   if (fDiff)  { delete fDiff; fDiff = 0;}
 144   if (fDiffSideBand) { delete fDiffSideBand; fDiffSideBand = 0;}
 145   if (fDStarMass)    { delete fDStarMass; fDStarMass = 0;}
 146   if (fPhi)     { delete fPhi; fPhi = 0;}
 147   if (fPhiBkg)  { delete fPhiBkg; fPhiBkg = 0;}
 148   if (fTrueDiff){ delete fTrueDiff; fTrueDiff = 0;}
 149   if (fResZ)    { delete fResZ;  fResZ = 0;}
 150   if (fResZBkg) { delete fResZBkg; fResZBkg = 0;}
 151   if (fEjet)    { delete fEjet; fEjet = 0;}
 152   if (fPhijet)  { delete fPhijet; fPhijet = 0;}
 153   if (fEtaJet)  { delete fEtaJet; fEtaJet = 0;}
 154   if (theMCFF)  { delete theMCFF; theMCFF = 0;}
 155   if (fDphiD0Dstar) { delete fDphiD0Dstar; fDphiD0Dstar = 0;}
 156   if (fPtJet) { delete fPtJet; fPtJet = 0;}
 157 }
 158
 159 //___________________________________________________________
 160 void AliAnalysisTaskSEDStarJets::Init(){
 161   //
 162   // Initialization
 163   //
 164   if(fDebug > 1) printf("AnalysisTaskSEDStarJets::Init() \n");
 165   AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
 166   // Post the cuts
 167   PostData(2,copyfCuts);
 168
 169   return;
 170 }
 171
 172 //_________________________________________________
 173 void AliAnalysisTaskSEDStarJets::UserExec(Option_t *)
 174 {
 175   // user exec
 176   if (!fInputEvent) {
 177     Error("UserExec","NO EVENT FOUND!");
 178     return;
 179   }
 180
 181   fEvents++;
 182   AliInfo(Form("Event %d",fEvents));
 183   if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
 184
 185   // Load the event
 186   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 187
 188   TClonesArray *arrayDStartoD0pi=0;
 189
 190   if(!aodEvent && AODEvent() && IsStandardAOD()) {
 191     // In case there is an AOD handler writing a standard AOD, use the AOD
 192     // event in memory rather than the input (ESD) event.
 193     aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 194     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 195     // have to taken from the AOD event hold by the AliAODExtension
 196     AliAODHandler* aodHandler = (AliAODHandler*)
 197       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 198     if(aodHandler->GetExtensions()) {
 199       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 200       AliAODEvent *aodFromExt = ext->GetAOD();
 201       arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 202     }
 203   } else {
 204     arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 205   }
 206
 207   if (!arrayDStartoD0pi){
 208     AliInfo("Could not find array of HF vertices, skipping the event");
 209     return;
 210   }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 211
 212   // fix for temporary bug in ESDfilter
 213   // the AODs with null vertex pointer didn't pass the PhysSel
 214   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 215
 216   // Simulate a jet triggered sample
 217   TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets();
 218   if(aodEvent->GetNJets()<=0) return;
 219
 220   // counters for efficiencies
 221   Int_t icountReco = 0;
 222
 223   // Normalization factor
 224   if(fRequireNormalization){
 225     ftrigger->Fill(1);
 226   }
 227
 228   //D* and D0 prongs needed to MatchToMC method
 229   Int_t pdgDgDStartoD0pi[2]={421,211};
 230   Int_t pdgDgD0toKpi[2]={321,211};
 231
 232   Double_t max =0;
 233   Double_t ejet   = 0;
 234   Double_t phiJet = 0;
 235   Double_t etaJet = 0;
 236   Double_t ptjet = 0;
 237
 238   //loop over jets and consider only the leading jey in the event
 239   for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) {
 240     AliAODJet* jet = (AliAODJet*)arrayofJets->At(iJets);
 241
 242     //jets variables
 243     ejet   = jet->E();
 244
 245     if(ejet>max){
 246       max = jet->E();
 247       phiJet = jet->Phi();
 248       etaJet = jet->Eta();
 249       ptjet = jet->Pt();
 250     }
 251
 252     // fill energy, eta and phi of the jet
 253     fEjet   ->Fill(ejet);
 254     fPhijet ->Fill(phiJet);
 255     fEtaJet ->Fill(etaJet);
 256     fPtJet->Fill(ptjet);
 257   }
 258
 259   //loop over D* candidates
 260   for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
 261
 262     // D* candidates
 263     AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
 264     AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 265
 266     Double_t finvM =-999;
 267     Double_t finvMDStar =-999;
 268     Double_t dPhi =-999;
 269     Bool_t isDStar =0;
 270     Int_t mcLabel = -9999;
 271
 272     // find associated MC particle for D* ->D0toKpi
 273     if(fUseMCInfo){
 274       TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 275       if(!mcArray) {
 276         printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
 277         return;
 278       }
 279       mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
 280
 281       if(mcLabel>=0) isDStar = 1;
 282       if(mcLabel>0){
 283         Double_t zMC =-999;
 284         AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
 285       //fragmentation function in mc
 286         zMC= FillMCFF(mcPart,mcArray,mcLabel);
 287         if(zMC>0) theMCFF->Fill(zMC);
 288       }
 289     }
 290
 291     // soft pion
 292     AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
 293     Double_t pt = dstarD0pi->Pt();
 294
 295     // track quality cuts
 296     Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
 297     if(!isTkSelected) continue;
 298
 299     // region of interest + topological cuts + PID
 300     if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
 301     Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
 302     if (!isSelected) continue;
 303
 304     // fill histos
 305     finvM = dstarD0pi->InvMassD0();
 306     fInvMass->Fill(finvM);
 307
 308     //DStar invariant mass
 309     finvMDStar = dstarD0pi->InvMassDstarKpipi();
 310
 311     Double_t EGjet = 0;
 312     Double_t dStarMom = dstarD0pi->P();
 313     Double_t phiDStar = dstarD0pi->Phi();
 314     Double_t phiD0    = theD0particle->Phi();
 315     //check suggested by Federico
 316     Double_t dPhiD0Dstar = phiD0 - phiDStar;
 317
 318     dPhi = phiJet - phiDStar;
 319
 320     //plot right dphi
 321     if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
 322     if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
 323
 324     Double_t corrFactorCharge = (ejet/100)*20;
 325     EGjet =  ejet + corrFactorCharge;
 326
 327     // fill D* candidates
 328     Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 329     if(finvM >= (mPDGD0-0.05) && finvM <=(mPDGD0+0.05)){ // ~3 sigma (sigma=17MeV, conservative)
 330
 331       if(isDStar == 1) {
 332         fDphiD0Dstar->Fill(dPhiD0Dstar);
 333         fDStarMass->Fill(finvMDStar);
 334         fTrueDiff->Fill(finvMDStar-finvM);
 335       }
 336       if(isDStar == 0) fDphiD0Dstar->Fill(dPhiD0Dstar); // angle between D0 and D*
 337
 338       fDStar->Fill(finvMDStar);
 339       fDiff ->Fill(finvMDStar-finvM);
 340
 341       Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
 342       Double_t invmassDelta = dstarD0pi->DeltaInvMass();
 343
 344       // now the dphi signal and the fragmentation function
 345       if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0019){
 346         //fill candidates D* and soft pion reco pt
 347
 348         fRECOPtDStar->Fill(pt);
 349         fPtPion->Fill(track2->Pt());
 350
 351         fPhi ->Fill(dPhi);
 352
 353         Double_t jetCone = 0.4;
 354         if(dPhi>=-jetCone && dPhi<=jetCone){  // evaluate in the near side inside UA1 radius
 355           Double_t zFrag = (TMath::Cos(dPhi)*dStarMom)/EGjet;
 356           fResZ->Fill(TMath::Abs(zFrag));
 357         }
 358       }
 359     }
 360     // evaluate side band background
 361     SideBandBackground(finvM, finvMDStar, dStarMom, EGjet, dPhi);
 362
 363   } // D* cand
 364
 365 AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
 366
 367 PostData(1,fOutput);
 368 }
 369
 370 //________________________________________ terminate ___________________________
 371 void AliAnalysisTaskSEDStarJets::Terminate(Option_t*)
 372 {
 373   // The Terminate() function is the last function to be called during
 374   // a query. It always runs on the client, it can be used to present
 375   // the results graphically or save the results to file.
 376
 377   Info("Terminate"," terminate");
 378   AliAnalysisTaskSE::Terminate();
 379
 380   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 381   if (!fOutput) {
 382     printf("ERROR: fOutput not available\n");
 383     return;
 384   }
 385
 386   fDStarMass    = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
 387   fTrueDiff     = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
 388   fInvMass      = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
 389   fPtPion       = dynamic_cast<TH1F*>(fOutput->FindObject("fPtPion "));
 390   fDStar        = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar"));
 391   fDiff         = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff"));
 392   fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
 393   ftrigger      = dynamic_cast<TH1F*>(fOutput->FindObject("ftrigger"));
 394   fRECOPtDStar  = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtDStar"));
 395   fRECOPtBkg    = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtBkg"));
 396   fEjet         = dynamic_cast<TH1F*>(fOutput->FindObject("fEjet"));
 397   fPhijet       = dynamic_cast<TH1F*>(fOutput->FindObject("fPhijet"));
 398   fEtaJet       = dynamic_cast<TH1F*>(fOutput->FindObject("fEtaJet"));
 399   fPhi          = dynamic_cast<TH1F*>(fOutput->FindObject("fPhi"));
 400   fResZ         = dynamic_cast<TH1F*>(fOutput->FindObject("fResZ"));
 401   fResZBkg      = dynamic_cast<TH1F*>(fOutput->FindObject("fResZBkg"));
 402   fPhiBkg       = dynamic_cast<TH1F*>(fOutput->FindObject("fPhiBkg"));
 403   theMCFF       = dynamic_cast<TH1F*>(fOutput->FindObject("theMCFF"));
 404   fDphiD0Dstar  = dynamic_cast<TH1F*>(fOutput->FindObject("fDphiD0Dstar"));
 405   fPtJet  = dynamic_cast<TH1F*>(fOutput->FindObject("fPtJet"));
 406
 407 }
 408 //___________________________________________________________________________
 409
 410 void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() {
 411  // output
 412   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 413
 414   //slot #1
 415   OpenFile(1);
 416   fOutput = new TList();
 417   fOutput->SetOwner();
 418   // define histograms
 419   DefineHistoFroAnalysis();
 420
 421   return;
 422 }
 423 //___________________________________ hiostograms _______________________________________
 424
 425 Bool_t  AliAnalysisTaskSEDStarJets::DefineHistoFroAnalysis(){
 426
 427   // Invariant mass related histograms
 428   fInvMass = new TH1F("invMass","Kpi invariant mass distribution",1500,.5,3.5);
 429   fInvMass->SetStats(kTRUE);
 430   fInvMass->GetXaxis()->SetTitle("GeV/c");
 431   fInvMass->GetYaxis()->SetTitle("Entries");
 432   fOutput->Add(fInvMass);
 433
 434   fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4);
 435   fDStar->SetStats(kTRUE);
 436   fDStar->GetXaxis()->SetTitle("GeV/c");
 437   fDStar->GetYaxis()->SetTitle("Entries");
 438   fOutput->Add(fDStar);
 439
 440   fDiff = new TH1F("Diff","M(kpipi)-M(kpi)",750,0.1,0.2);
 441   fDiff->SetStats(kTRUE);
 442   fDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV/c^2");
 443   fDiff->GetYaxis()->SetTitle("Entries");
 444   fOutput->Add(fDiff);
 445
 446   fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
 447   fDiffSideBand->SetStats(kTRUE);
 448   fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
 449   fDiffSideBand->GetYaxis()->SetTitle("Entries");
 450   fOutput->Add(fDiffSideBand);
 451
 452   fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
 453   fOutput->Add(fDStarMass);
 454
 455   fTrueDiff  = new TH1F("dstar","True Reco diff",750,0,0.2);
 456   fOutput->Add(fTrueDiff);
 457
 458   // trigger normalization
 459   ftrigger = new TH1F("Normalization","Normalization factor for correlations",1,0,10);
 460   ftrigger->SetStats(kTRUE);
 461   fOutput->Add(ftrigger);
 462
 463   //correlation fistograms
 464   fPhi = new TH1F("phi","Delta phi between Jet axis and DStar ",25,-1.57,4.72);
 465   fPhi->SetStats(kTRUE);
 466   fPhi->GetXaxis()->SetTitle("#Delta #phi (rad)");
 467   fPhi->GetYaxis()->SetTitle("Entries");
 468   fOutput->Add(fPhi);
 469
 470   fDphiD0Dstar = new TH1F("phiD0Dstar","Delta phi between D0 and DStar ",1000,-6.5,6.5);
 471   fOutput->Add(fDphiD0Dstar);
 472
 473   fPhiBkg = new TH1F("phiBkg","Delta phi between Jet axis and DStar background ",25,-1.57,4.72);
 474   fPhiBkg->SetStats(kTRUE);
 475   fPhiBkg->GetXaxis()->SetTitle("#Delta #phi (rad)");
 476   fPhiBkg->GetYaxis()->SetTitle("Entries");
 477   fOutput->Add(fPhiBkg);
 478
 479   fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",30,0,30);
 480   fRECOPtDStar->SetStats(kTRUE);
 481   fRECOPtDStar->SetLineColor(2);
 482   fRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
 483   fRECOPtDStar->GetYaxis()->SetTitle("Entries");
 484   fOutput->Add(fRECOPtDStar);
 485
 486   fRECOPtBkg = new TH1F("RECOptBkg","RECO pt distribution side bands",30,0,30);
 487   fRECOPtBkg->SetStats(kTRUE);
 488   fRECOPtBkg->SetLineColor(2);
 489   fRECOPtBkg->GetXaxis()->SetTitle("GeV/c");
 490   fRECOPtBkg->GetYaxis()->SetTitle("Entries");
 491   fOutput->Add(fRECOPtBkg);
 492
 493   fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,10);
 494   fPtPion->SetStats(kTRUE);
 495   fPtPion->GetXaxis()->SetTitle("GeV/c");
 496   fPtPion->GetYaxis()->SetTitle("Entries");
 497   fOutput->Add(fPtPion);
 498
 499   // jet related fistograms
 500   fEjet      = new TH1F("ejet",  "UA1 algorithm jet energy distribution",1000,0,500);
 501   fPhijet    = new TH1F("Phijet","UA1 algorithm jet #phi distribution",  200,-7,7);
 502   fEtaJet    = new TH1F("Etajet","UA1 algorithm jet #eta distribution",  200,-7,7);
 503   fPtJet      = new TH1F("PtJet",  "UA1 algorithm jet Pt distribution",1000,0,500);
 504   fOutput->Add(fEjet);
 505   fOutput->Add(fPhijet);
 506   fOutput->Add(fEtaJet);
 507   fOutput->Add(fPtJet);
 508
 509   theMCFF    = new TH1F("FragFuncMC","Fragmentation function in MC for FC ",100,0,10);
 510   fResZ      = new TH1F("FragFunc","Fragmentation function ",50,0,1);
 511   fResZBkg   = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1);
 512   fOutput->Add(theMCFF);
 513   fOutput->Add(fResZ);
 514   fOutput->Add(fResZBkg);
 515
 516   return kTRUE;
 517 }
 518
 519 //______________________________ side band background for D*___________________________________
 520
 521 void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t dStarMomBkg, Double_t EGjet, Double_t dPhi){
 522
 523   //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
 524   // (expected detector resolution) on the left and right frm the D0 mass. Each band
 525   //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
 526
 527   if((invM>=1.7 && invM<=1.8) || (invM>=1.92 && invM<=2.19)){
 528     fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background
 529     if ((invMDStar-invM)<=0.14732 && (invMDStar-invM)>=0.14352) {
 530       fPhiBkg->Fill(dPhi);
 531       fRECOPtBkg->Fill(dStarMomBkg);
 532       if(dPhi>=-0.4 && dPhi<=0.4){  // evaluate in the near side
 533         Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/EGjet;
 534         fResZBkg->Fill(TMath::Abs(zFragBkg));
 535       }
 536     }
 537   }
 538 }
 539
 540 //_____________________________________________________________________________________________-
 541 double AliAnalysisTaskSEDStarJets::FillMCFF(AliAODMCParticle* mcPart, TClonesArray* mcArray, Int_t mcLabel){
 542   //
 543   // GS from MC
 544   // UA1 jet algorithm reproduced in MC
 545   //
 546   Double_t zMC2 =-999;
 547
 548   Double_t leading =0;
 549   Double_t PartE = 0;
 550   Double_t PhiLeading = -999;
 551   Double_t EtaLeading = -999;
 552   Double_t PtLeading = -999;
 553   Int_t counter =-999;
 554
 555   if (!mcPart) return zMC2;
 556   if (!mcArray) return zMC2;
 557
 558   //find leading particle
 559   for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 560     AliAODMCParticle* Part = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 561     if (!Part) {
 562       AliWarning("MC Particle not found in tree, skipping");
 563       continue;
 564     }
 565
 566     // remove quarks and the leading particle (it will be counted later)
 567     if(iPart == mcLabel) continue;
 568     if(iPart <= 8) continue;
 569
 570     //remove resonances not directly detected in detector
 571     Int_t PDGCode = Part->GetPdgCode();
 572
 573     // be sure the particle reach the detector
 574     Double_t x = Part->Xv();
 575     Double_t y = Part->Yv();
 576     Double_t z = Part->Zv();
 577
 578     if(TMath::Abs(PDGCode)== 2212 && x<3 && y<3) continue;
 579     if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue;
 580     if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 &&  TMath::Abs(PDGCode)!=11 &&  TMath::Abs(PDGCode)!=13 &&  TMath::Abs(PDGCode)!=2212) continue;
 581
 582     Int_t daug0 = -999;
 583     Double_t xd =-999;
 584     Double_t yd =-999;
 585     Double_t zd =-999;
 586
 587     daug0 = Part->GetDaughter(0);
 588
 589     if(daug0>=0){
 590       AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
 591       if(tdaug){
 592       xd = tdaug->Xv();
 593       yd = tdaug->Yv();
 594       zd = tdaug->Zv();
 595       }
 596     }
 597     if(TMath::Abs(xd)<3 || TMath::Abs(yd)<3) continue;
 598
 599     Bool_t AliceAcc = (TMath::Abs(Part->Eta()) <= 0.9);
 600     if(!AliceAcc) continue;
 601
 602     PartE  = Part->E();
 603
 604     if(PartE>leading){
 605       leading = Part->E();
 606       PhiLeading = Part->Phi();
 607       EtaLeading = Part->Eta();
 608       PtLeading = Part->Pt();
 609       counter = iPart;
 610     }
 611
 612   }
 613
 614   Double_t jetEnergy = 0;
 615
 616   //reconstruct the jet
 617   for (Int_t iiPart=0; iiPart<mcArray->GetEntriesFast(); iiPart++) {
 618     AliAODMCParticle* tPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iiPart));
 619     if (!tPart) {
 620       AliWarning("MC Particle not found in tree, skipping");
 621       continue;
 622     }
 623     // remove quarks and the leading particle (it will be counted later)
 624     if(iiPart == counter) continue; // do not count again the leading particle
 625     if(iiPart == mcLabel) continue;
 626     if(iiPart <= 8) continue;
 627
 628     //remove resonances not directly detected in detector
 629     Int_t PDGCode = tPart->GetPdgCode();
 630
 631     // be sure the particle reach the detector
 632     Double_t x = tPart->Xv();
 633     Double_t y = tPart->Yv();
 634     Double_t z = tPart->Zv();
 635
 636     if(TMath::Abs(PDGCode)== 2212 && (x<3 && y<3)) continue;
 637     if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue; // has to be generated at least in the silicon tracker or beam pipe
 638     if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 &&  TMath::Abs(PDGCode)!=11 &&  TMath::Abs(PDGCode)!=13 &&  TMath::Abs(PDGCode)!=2212) continue;
 639
 640
 641     Int_t daug0 = -999;
 642     Double_t xd =-999;
 643     Double_t yd =-999;
 644     Double_t zd =-999;
 645
 646     daug0 = tPart->GetDaughter(0);
 647
 648     if(daug0>=0){
 649       AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
 650       if(tdaug){
 651       xd = tdaug->Xv();
 652       yd = tdaug->Yv();
 653       zd = tdaug->Zv();
 654       }
 655     }
 656     if(TMath::Abs(xd)<3 && TMath::Abs(yd)<3) continue;
 657     //remove particles not in ALICE acceptance
 658     if(tPart->Pt()<0.07) continue;
 659     Bool_t AliceAcc = (TMath::Abs(tPart->Eta()) <= 0.9);
 660     if(!AliceAcc) continue;
 661
 662     Double_t EtaMCp = tPart->Eta();
 663     Double_t PhiMCp = tPart->Phi();
 664
 665     Double_t DphiMClead = PhiLeading-PhiMCp;
 666
 667     if(DphiMClead<=-(TMath::Pi())/2) DphiMClead = DphiMClead+2*(TMath::Pi());
 668     if(DphiMClead>(3*(TMath::Pi()))/2) DphiMClead = DphiMClead-2*(TMath::Pi());
 669
 670     Double_t deta = (EtaLeading-EtaMCp);
 671     //cone radius
 672     Double_t radius = TMath::Sqrt((DphiMClead*DphiMClead)+(deta*deta));
 673
 674     if(radius>0.4) continue; // in the jet cone
 675     if(tPart->Charge()==0) continue; // only charged fraction
 676
 677     jetEnergy = jetEnergy+(tPart->E());
 678   }
 679
 680   jetEnergy = jetEnergy + leading;
 681
 682   // delta phi D*, jet axis
 683   Double_t dPhi = PhiLeading - (mcPart->Phi());
 684
 685   //plot right dphi
 686   if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
 687   if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
 688
 689   zMC2 = (TMath::Cos(dPhi)*(mcPart->P()))/jetEnergy;
 690
 691   return zMC2;
 692 }