Updates from Chiara for flavor tagging on reco and MC truth level (Chiara)
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 Nov 2013 11:58:52 +0000 (11:58 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 Nov 2013 11:58:52 +0000 (11:58 +0000)
PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx
PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.h
PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.cxx
PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.h
PWGJE/FlavourJetTasks/macros/AddTaskFlavourJetCorrelations.C
PWGJE/FlavourJetTasks/macros/AddTaskSEDmesonsFilterCJ.C
PWGJE/FlavourJetTasks/macros/AddTasksFlavourJet.C
PWGJE/FlavourJetTasks/macros/AnalysisTrainCorrJetsLocal.C

index 788a121..906c920 100644 (file)
-/**************************************************************************
- * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-//
-//
-//  Analysis Taks for reconstructed particle correlation 
-//  (first implementation done for D mesons) with jets 
-//  (use the so called Emcal framework)
-//
-//-----------------------------------------------------------------------
-// Authors:
-// C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
-// A. Grelli (Utrecht University) a.grelli@uu.nl
-// X. Zhang (LBNL)  XMZhang@lbl.gov
-//-----------------------------------------------------------------------
-
-#include <TDatabasePDG.h>
-#include <TParticle.h>
-#include <TVector3.h>
-#include "TROOT.h"
-#include <TH3F.h>
-
-#include "AliAnalysisTaskFlavourJetCorrelations.h"
-#include "AliAODMCHeader.h"
-#include "AliAODHandler.h"
-#include "AliAnalysisManager.h"
-#include "AliLog.h"
-#include "AliEmcalJet.h"
-#include "AliAODRecoDecay.h"
-#include "AliAODRecoCascadeHF.h"
-#include "AliAODRecoDecayHF2Prong.h"
-#include "AliESDtrack.h"
-#include "AliAODMCParticle.h"
-#include "AliPicoTrack.h"
-#include "AliRDHFCutsD0toKpi.h"
-#include "AliRDHFCutsDStartoKpipi.h"
-#include "AliJetContainer.h"
-
-ClassImp(AliAnalysisTaskFlavourJetCorrelations)
-
-//__________________________________________________________________________
-AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
-  AliAnalysisTaskEmcalJet("",kFALSE),
-  fUseMCInfo(kTRUE), 
-  fCandidateType(),
-  fPDGmother(),
-  fNProngs(),
-  fPDGdaughters(),
-  fBranchName(),
-  fOutput(0),
-  fCuts(0),
-  fMinMass(),
-  fMaxMass(),  
-  fJetArrName(0),
-  fCandArrName(0),
-  fLeadingJetOnly(kFALSE),
-  fJetRadius(0)
-{
-  //
-  // Default ctor
-  //
-}
-//___________________________________________________________________________
-AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype, TString jetArrName) :
-  AliAnalysisTaskEmcalJet(name,kFALSE),
-  fUseMCInfo(kTRUE),
-  fCandidateType(),
-  fPDGmother(),
-  fNProngs(),
-  fPDGdaughters(),
-  fBranchName(),
-  fOutput(0),
-  fCuts(0),
-  fMinMass(),
-  fMaxMass(),  
-  fJetArrName(0),
-  fCandArrName(0),
-  fLeadingJetOnly(kFALSE),
-  fJetRadius(0)
-{
-  //
-  // Constructor. Initialization of Inputs and Outputs
-  //
-  Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
-  fCuts=cuts;
-  fCandidateType=candtype;
-  const Int_t nptbins=fCuts->GetNPtBins();
-  Float_t defaultSigmaD013[13]={0.012, 0.012, 0.012, 0.015, 0.015,0.018,0.018,0.020,0.020,0.030,0.030,0.037,0.040};
-
-  switch(fCandidateType){
-  case 0:
-    fPDGmother=421;
-    fNProngs=2;
-    fPDGdaughters[0]=211;//pi 
-    fPDGdaughters[1]=321;//K
-    fPDGdaughters[2]=0; //empty
-    fPDGdaughters[3]=0; //empty
-    fBranchName="D0toKpi";
-    fCandArrName="D0";
-    break;
-  case 1: 
-    fPDGmother=413;
-    fNProngs=3;
-    fPDGdaughters[1]=211;//pi soft
-    fPDGdaughters[0]=421;//D0
-    fPDGdaughters[2]=211;//pi fromD0
-    fPDGdaughters[3]=321; // K from D0
-    fBranchName="Dstar";
-    fCandArrName="Dstar";
-    //fSigmaD0=new Float_t[nptbins];
-    if(nptbins<=13){
-      for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
-    } else {
-      AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
-    }
-    break;
-  default:
-    printf("%d not accepted!!\n",fCandidateType);
-    break;
-  }
-
-  if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
-  if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
-  fJetArrName=jetArrName;
-  Printf("Jet read is %s",fJetArrName.Data());
-  
-  DefineOutput(1,TList::Class()); // histos
-  DefineOutput(2,AliRDHFCuts::Class()); // my cuts
-}
-//___________________________________________________________________________
-AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
-  //
-  // destructor
-  //
-
-  Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");  
-    delete fOutput;
-    delete fCuts;
-    
-}
-
-//___________________________________________________________
-void AliAnalysisTaskFlavourJetCorrelations::Init(){
-  //
-  // Initialization
-  //
-  if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
-  switch(fCandidateType){
-  case 0:
-    {
-      AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
-      copyfCuts->SetName("AnalysisCutsDzero");
-      // Post the data
-      PostData(2,copyfCuts);
-    }
-    break;
-  case 1:
-    {
-      AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
-      copyfCuts->SetName("AnalysisCutsDStar");
-      // Post the cuts
-      PostData(2,copyfCuts);
-    }
-    break;
-  default:
-    return;
-  }
-
-  return;
-}
-
-//___________________________________________________________________________
-
-void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() { 
- // output 
-  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
-  fOutput = new TList();
-  fOutput->SetOwner();
-  fOutput->SetName("fOutputDjet");
-  // define histograms
-  DefineHistoForAnalysis();
-  PostData(1,fOutput);
-
-  return;
-}
-
-//_________________________________________________
-void AliAnalysisTaskFlavourJetCorrelations::UserExec(Option_t *)
-{
-  // user exec
-
-  // Load the event
-  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
-  TClonesArray *arrayDStartoD0pi=0;
-  TClonesArray *trackArr = 0;
-  //TClonesArray *clusArr = 0;  
-  TClonesArray *jetArr = 0;
-  TClonesArray *candidatesArr = 0;
-//TClonesArray *isselArr = 0;
-
-  if (!aodEvent && AODEvent() && IsStandardAOD()) {
-
-    // In case there is an AOD handler writing a standard AOD, use the AOD 
-    // event in memory rather than the input (ESD) event.    
-    aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
-
-    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
-    // have to taken from the AOD event hold by the AliAODExtension
-    AliAODHandler* aodHandler = (AliAODHandler*) 
-      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
-    if(aodHandler->GetExtensions()) {
-      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
-      AliAODEvent *aodFromExt = ext->GetAOD();
-      arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
-    }
-  } else {
-    arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
-  }
-
-  if (!arrayDStartoD0pi) {
-    AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
-//  return;
-  } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   
-
-  TClonesArray* mcArray = 0x0;
-  if (fUseMCInfo) {
-    mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
-    if (!mcArray) {
-      printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
-      return;
-    }
-  }
-
-  //retrieve jets
-  trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
-  //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
-  jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
-  fJetRadius=GetJetContainer()->GetJetRadius();
-
-  if(!trackArr){
-    AliInfo(Form("Could not find the track array\n"));
-    return;
-  }
-
-  if(!jetArr){
-    Printf("JET array not found");
-    return;
-  }
-
-  //retrieve reconstructed particles selected
-  /*
-  TString listname="AliAnalysisTaskSEDmesonsForJetCorrelations";
-  TList *l=dynamic_cast<TList*>(InputEvent()->FindListObject(listname));
-  TClonesArray *cla=dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(listname));
-  // if(l){
-  //   l->ls();
-  // } else{
-  //   Printf("list not found!!!!!!!!!!!");
-  //   return;
-  // } 
-  if(!cla){
-    Printf("cla not found!!!!!!!!!!!");
-    return;
-  } else {
-    cla->ls();
-  }
-  */
-
-  TString arrname="fCandidateArray";
-  candidatesArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s%s",arrname.Data(),fCandArrName.Data())));
-  if (!candidatesArr) {
-    Printf("%s%s array not found",arrname.Data(),fCandArrName.Data());
-    InputEvent()->GetList()->ls();
-    return;
-  }
-
-  /*
-  arrname="fIsSelectedArray";
-  isselArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s%s",arrname.Data(),fCandArrName.Data())));
-  if(!isselArr){
-    Printf("%s%s array not found",arrname.Data(),fCandArrName.Data());
-    InputEvent()->ls();
-    return;
-  }
-  */
-
-  //Histograms
-  TH1I* hstat=(TH1I*)fOutput->FindObject("hstat");
-  TH1F* hPtJetTrks=(TH1F*)fOutput->FindObject("hPtJetTrks");
-  TH1F* hPhiJetTrks=(TH1F*)fOutput->FindObject("hPhiJetTrks");
-  TH1F* hEtaJetTrks=(TH1F*)fOutput->FindObject("hEtaJetTrks");
-  TH1F* hEjetTrks=(TH1F*)fOutput->FindObject("hEjetTrks");
-  TH1F* hPtJet=(TH1F*)fOutput->FindObject("hPtJet");
-  TH1F* hPhiJet=(TH1F*)fOutput->FindObject("hPhiJet");
-  TH1F* hEtaJet=(TH1F*)fOutput->FindObject("hEtaJet");
-  TH1F* hEjet=(TH1F*)fOutput->FindObject("hEjet");
-  TH1F* hNtrArr=(TH1F*)fOutput->FindObject("hNtrArr");
-  TH1F* hNJetPerEv=(TH1F*)fOutput->FindObject("hNJetPerEv");
-  TH1F* hdeltaRJetTracks=((TH1F*)fOutput->FindObject("hdeltaRJetTracks"));
-
-  hstat->Fill(0);
-
-  // fix for temporary bug in ESDfilter 
-  // the AODs with null vertex pointer didn't pass the PhysSel
-  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
-
-  //Event selection
-  Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
-  TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
-  if(!iseventselected) return;
-  
-  hstat->Fill(1);
-  
-  //trigger on jets
-  Int_t njets=jetArr->GetEntriesFast();
-  if(njets == 0) return;
-
-  const Int_t nD=arrayDStartoD0pi->GetEntriesFast();
-  hstat->Fill(2,nD);
-  
-  // counters for efficiencies
-  Int_t icountReco = 0;
-  
-  //D* and D0 prongs needed to MatchToMC method
-  // Int_t pdgDgDStartoD0pi[2]={421,211};
-  // Int_t pdgDgD0toKpi[2]={321,211};
-  //D0 from D0 bar
-  // we start with jets
-  Double_t ejet   = 0;
-  Double_t phiJet = 0;
-  Double_t etaJet = 0;
-  Double_t ptjet = 0;
-  Double_t leadingJet =0;
-
-  Int_t ntrarr=trackArr->GetEntriesFast();
-  hNtrArr->Fill(ntrarr);
-  
-  for(Int_t i=0;i<ntrarr;i++){
-    AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));
-    //reusing histograms
-    hPtJetTrks->Fill(jtrack->Pt());
-    hPhiJetTrks->Fill(jtrack->Phi());
-    hEtaJetTrks->Fill(jtrack->Eta());
-    hEjetTrks->Fill(jtrack->E());
-  }
-
-  
-  //option to use only the leading jet
-  if(fLeadingJetOnly){
-    for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {    
-      AliEmcalJet* jetL = (AliEmcalJet*)jetArr->At(iJetsL);
-      ptjet   = jetL->Pt();
-      if(ptjet>leadingJet ) leadingJet = ptjet;
-
-    }
-  }
-
-  Int_t cntjet=0;
-  //loop over jets and consider only the leading jet in the event
-  for (Int_t iJets = 0; iJets<njets; iJets++) {    
-    AliEmcalJet* jet = (AliEmcalJet*)jetArr->At(iJets);
-    if(!AcceptJet(jet)) {
-       hstat->Fill(5);
-       continue;
-    }
-    std::vector<int> DmesonInJetLabels(10);
-    //Int_t iD=0;
-    //jets variables
-    ejet   = jet->E();
-    phiJet = jet->Phi();
-    etaJet = jet->Eta();
-    ptjet = jet->Pt();
-    
-    // choose the leading jet
-    if(fLeadingJetOnly && (ejet<leadingJet)) continue;
-    //used for normalization
-    hstat->Fill(3);
-    cntjet++;
-    // fill energy, eta and phi of the jet
-    hEjet   ->Fill(ejet);
-    hPhiJet ->Fill(phiJet);
-    hEtaJet ->Fill(etaJet);
-    hPtJet  ->Fill(ptjet);
-    
-    //loop on jet particles
-    Int_t ntrjet=  jet->GetNumberOfTracks();    
-    for(Int_t itrk=0;itrk<ntrjet;itrk++){
-      
-      AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);     
-      hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
-      
-      // check MC in for the traks within the jet, look at D mesons
-      // within the jet area
-      //if(fUseMCInfo) FillMCDJetInfo(jetTrk,jet,mcArray,ptjet);
-      
-    }//end loop on jet tracks
-    
-    //retrieve charm candidates selected
-    Int_t candidates = candidatesArr->GetEntriesFast();
-    for(Int_t ic = 0; ic < candidates; ic++) {
-      // D* candidates
-      AliAODRecoDecayHF* charm = 0x0;
-      charm=(AliAODRecoDecayHF*)candidatesArr->At(ic);
-      
-      FillHistogramsRecoJetCorr(charm, jet);      
-
-    }
-  } // end of jet loop 
-
-  hNJetPerEv->Fill(cntjet);
-
-  AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
-  PostData(1,fOutput);
-
-}
-
-//________________________________________ terminate ___________________________
-void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
-{    
-  // The Terminate() function is the last function to be called during
-  // a query. It always runs on the client, it can be used to present
-  // the results graphically or save the results to file.
-  
-  Info("Terminate"," terminate");
-  AliAnalysisTaskSE::Terminate();
-
-  fOutput = dynamic_cast<TList*> (GetOutputData(1));
-  if (!fOutput) {     
-    printf("ERROR: fOutput not available\n");
-    return;
-  }
-}
-
-//_________________________________________________________________
-void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
-  Float_t mass=0;
-  Int_t abspdg=TMath::Abs(pdg);
-  
-  mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
-  // compute the Delta mass for the D*
-   if(fCandidateType==kDstartoKpipi){
-     Float_t mass1=0;
-     mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
-     mass = mass-mass1;
-  }
-   std::cout<<"mass ---------------"<<std::endl;
-  fMinMass = mass-range;
-  fMaxMass = mass+range;
-  
-  AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
-  if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
-}
-//_________________________________________________________________
-void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
-  if(uplimit>lowlimit)
-    {
-      fMinMass = lowlimit;
-      fMaxMass = uplimit;
-    }
-  else{
-    printf("Error! Lower limit larger than upper limit!\n");
-    fMinMass = uplimit - uplimit*0.2;
-    fMaxMass = uplimit;
-  }
-}
-
-//__________________________________________________________________
-Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
-  if(nptbins>30) {
-    AliInfo("Maximum number of bins allowed is 30!");
-    return kFALSE;
-  }
-  if(!width) return kFALSE;
-  for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
-  return kTRUE;
-}
-
-//__________________________________________________________________
-Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
-  if(!part || !jet){
-    printf("Particle or jet do not exist!\n");
-    return -99;
-  }
-  Double_t p[3],pj[3];
-  Bool_t okpp=part->PxPyPz(p);
-  Bool_t okpj=jet->PxPyPz(pj);
-  if(!okpp || !okpj){
-    printf("Problems getting momenta\n");
-    return -999;
-  }
-  Double_t pjet=jet->P();
-  Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);
-  return z;
-}
-//___________________________________ histograms _______________________________________
-
-Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
-  // Statistics 
-  TH1I* hstat=new TH1I("hstat","Statistics",6,-0.5,5.5);
-  hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
-  hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
-  hstat->GetXaxis()->SetBinLabel(3,"N cand sel cuts");
-  hstat->GetXaxis()->SetBinLabel(4,"N jets");
-  hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
-  hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
-  // if(fUseMCInfo) {
-  //   hstat->GetXaxis()->SetBinLabel(7,"N D");
-  //   hstat->GetXaxis()->SetBinLabel(8,"N D in jet");
-
-  // }
-
-  hstat->SetNdivisions(1);
-  fOutput->Add(hstat);
-
-  const Int_t nbinsmass=200;
-
-  
-  if(fCandidateType==kDstartoKpipi) 
-    {
-      // TH2F *hDiff = new TH2F("hDiff","M(kpipi)-M(kpi)",500,fMinMass,fMaxMass,100, 0.,30.);
-      // hDiff->SetStats(kTRUE);
-      // hDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV");
-      // hDiff->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-      // fOutput->Add(hDiff);
-  
-      TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,100, 0.,30.);
-      hDiffSideBand->SetStats(kTRUE);
-      hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
-      hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-      fOutput->Add(hDiffSideBand); 
-      //correlation histograms
-      // fPhi = new TH1F("phi","Delta phi between Jet axis and DStar ",25,-1.57,4.72);
-      // fPhi->SetStats(kTRUE);
-      // fPhi->GetXaxis()->SetTitle("#Delta #phi (rad)");
-      // fPhi->GetYaxis()->SetTitle("Entries");
-      // fOutput->Add(fPhi);
-
-      // fDphiD0Dstar = new TH1F("phiD0Dstar","Delta phi between D0 and DStar ",1000,-6.5,6.5);
-      // fOutput->Add(fDphiD0Dstar);
-
-      // fPhiBkg = new TH1F("phiBkg","Delta phi between Jet axis and DStar background ",25,-1.57,4.72);
-      // fPhiBkg->SetStats(kTRUE);
-      // fPhiBkg->GetXaxis()->SetTitle("#Delta #phi (rad)");
-      // fPhiBkg->GetYaxis()->SetTitle("Entries");
-      // fOutput->Add(fPhiBkg);
-
-      // TH1F* hRECOPtDStar = new TH1F("hRECODStar","RECO DStar pt distribution",30,0,30);
-      // hRECOPtDStar->SetStats(kTRUE);
-      // hRECOPtDStar->SetLineColor(2);
-      // hRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
-      // hRECOPtDStar->GetYaxis()->SetTitle("Entries");
-      // fOutput->Add(hRECOPtDStar);
-  
-      // TH1F* hRECOPtBkg = new TH1F("hRECOptBkg","RECO pt distribution side bands",30,0,30);
-      // hRECOPtBkg->SetStats(kTRUE);
-      // hRECOPtBkg->SetLineColor(2);
-      // hRECOPtBkg->GetXaxis()->SetTitle("GeV/c");
-      // hRECOPtBkg->GetYaxis()->SetTitle("Entries");
-      // fOutput->Add(hRECOPtBkg);
-
-      TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
-      hPtPion->SetStats(kTRUE);
-      hPtPion->GetXaxis()->SetTitle("GeV/c");
-      hPtPion->GetYaxis()->SetTitle("Entries");
-      fOutput->Add(hPtPion);
-
-    }
-
-  // jet related fistograms
-  
-  TH1F* hEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
-  TH1F* hPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  200,0,6.30);
-  TH1F* hEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  100,-1.5,1.5);
-  TH1F* hPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",500,0,200);
-  
-  TH1F* hEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
-  TH1F* hPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  200,0,6.30);
-  TH1F* hEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta",  100,-1.5,1.5);
-  TH1F* hPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",500,0,200);
-
-  TH2F* hPtJetWithD=new TH2F("hPtJetWithD","D-Jet Pt distribution; p_{T} (GeV/c);inv mass (GeV)",500,0,200,nbinsmass,fMinMass,fMaxMass);
-  TH1F* hPtJetWithDsb=new TH1F("hPtJetWithDsb","D(background)-Jet Pt distribution; p_{T} (GeV/c)",500,0,200);
-  TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
-
-  TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
-  TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
-
-  fOutput->Add(hEjetTrks);
-  fOutput->Add(hPhiJetTrks);
-  fOutput->Add(hEtaJetTrks);
-  fOutput->Add(hPtJetTrks);
-  fOutput->Add(hEjet);
-  fOutput->Add(hPhiJet);
-  fOutput->Add(hEtaJet);
-  fOutput->Add(hPtJet);
-  fOutput->Add(hPtJetWithD);
-  fOutput->Add(hPtJetWithDsb);
-  fOutput->Add(hdeltaRJetTracks);
-  fOutput->Add(hNtrArr);
-  fOutput->Add(hNJetPerEv);
-
-  /*
-  if(fUseMCInfo){
-    fhMomjetpartPdg=new TH1F("fhMomjetpartPdg","Jet particles' mothers distribution;PDG code;Counts;",1100,-550,550);
-    fOutput->Add(fhMomjetpartPdg);
-    fptDinjetallvsptjMC=new TH2F("fptDinjetallvsptjMC","p_{t} distribution of D within a jet (all DeltaR) vs p_{t}^{jet};p_{t}^{D};p_{t}^{jet}",100, 0.,30.,500,0.,200.);
-    fOutput->Add(fptDinjetallvsptjMC);
-    fptJwithDMC=new TH1F("fptJwithDMC","p_{t}^{jet} with a D meson (all DeltaR);p_{t}^{jet};Counts",500,0.,200.);
-    fOutput->Add(fptJwithDMC);
-
-    fptDinjet04vsptjMC=new TH2F("fptDinjet04vsptjMC","p_{t} distribution of D within a jet (DeltaR 0.4) vs p_{t}^{jet};p_{t}^{D};p_{t}^{jet}",100, 0.,30.,500,0.,200.);
-    fOutput->Add(fptDinjet04vsptjMC);
-    TH1F* hdeltaRDMC=new TH1F("hdeltaRDMC","#Delta R for MC tagged D mesons; #Delta R_{D}^{MC}",200, 0.,10.);
-    fOutput->Add(hdeltaRDMC);
-  }
-  */
-  TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
-  fOutput->Add(hDeltaRD);
-
-  TH3F* hdeltaPhiDja=new TH3F("hdeltaPhiDja", "#Delta#phi D-jet (jet p_{T} > threshold)",nbinsmass,fMinMass,fMaxMass,100, 0.,30.,50,(-1)*TMath::Pi()/2.,3./2.*TMath::Pi());
-  hdeltaPhiDja->GetXaxis()->SetTitle("mass (GeV)");
-  hdeltaPhiDja->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-  hdeltaPhiDja->GetYaxis()->SetTitle("#Delta#phi (rad)");
-  // TH3F* hdeltaPhiDjl=new TH3F("hdeltaPhiDjl", Form("#Delta#phi D-jet (jet p_{T} < %.0f (GeV/c))",fJetPtThr),500,fMinMass,fMaxMass,100, 0.,30.,50,(-1)*TMath::Pi()/2.,3./2.*TMath::Pi());
-  // hdeltaPhiDjl->GetXaxis()->SetTitle("mass (GeV)");
-  // hdeltaPhiDjl->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-  // hdeltaPhiDjl->GetYaxis()->SetTitle("#Delta#phi (rad)");
-  // TH3F* hdeltaPhiDjh=new TH3F("hdeltaPhiDjh", Form("#Delta#phi D-jet (jet p_{T} > %.0f (GeV/c))",fJetPtThr),500,fMinMass,fMaxMass,100, 0.,30.,50,(-1)*TMath::Pi()/2.,3./2.*TMath::Pi());
-  // hdeltaPhiDjh->GetXaxis()->SetTitle("mass (GeV/c)");
-  // hdeltaPhiDjh->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-  // hdeltaPhiDjh->GetYaxis()->SetTitle("#Delta#phi (rad)");
-  fOutput->Add(hdeltaPhiDja);
-  // fOutput->Add(hdeltaPhiDjl);
-  // fOutput->Add(hdeltaPhiDjh);
-
-  //background (side bands for the Dstar and like sign for D0)
-
-  TH3F* hdeltaPhiDjaB=new TH3F("hdeltaPhiDjaB", "#Delta#phi D-jet (all jet p_{T})",nbinsmass,fMinMass,fMaxMass,100, 0.,30.,50,(-1)*TMath::Pi()/2.,3./2.*TMath::Pi());
-  hdeltaPhiDjaB->GetXaxis()->SetTitle("mass (GeV)");
-  hdeltaPhiDjaB->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-  hdeltaPhiDjaB->GetYaxis()->SetTitle("#Delta#phi (rad)");
-  // TH3F* hdeltaPhiDjlB=new TH3F("hdeltaPhiDjlB", Form("#Delta#phi D-jet (jet p_{T} < %.0f (GeV/c))",fJetPtThr),1500,fMinMass,fMaxMass,100, 0.,30.,50,(-1)*TMath::Pi()/2.,3./2.*TMath::Pi());
-  // hdeltaPhiDjlB->GetXaxis()->SetTitle("mass (GeV)");
-  // hdeltaPhiDjlB->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-  // hdeltaPhiDjlB->GetYaxis()->SetTitle("#Delta#phi (rad)");
-  // TH3F* hdeltaPhiDjhB=new TH3F("hdeltaPhiDjhB", Form("#Delta#phi D-jet (jet p_{T} > %.0f (GeV/c))",fJetPtThr),1500,fMinMass,fMaxMass,100, 0.,30.,50,(-1)*TMath::Pi()/2.,3./2.*TMath::Pi());
-  // hdeltaPhiDjhB->GetXaxis()->SetTitle("mass (GeV)");
-  // hdeltaPhiDjhB->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-  // hdeltaPhiDjhB->GetYaxis()->SetTitle("#Delta#phi (rad)");
-  fOutput->Add(hdeltaPhiDjaB);
-  // fOutput->Add(hdeltaPhiDjlB);
-  // fOutput->Add(hdeltaPhiDjhB);
-
-  TH2F* hInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,100,0.,50.);
-  hInvMassptD->SetStats(kTRUE);
-  hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
-  hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-
-  fOutput->Add(hInvMassptD);
-  // fMasspjDeltaR=new TH3F("fMasspjDeltaR","Mass vs p^{jet} vs #Delta R;Mass (Gev/c);p^{jet}(Gev);#Delta R",1500,fMinMass,fMaxMass,100,0.,50.,100,0.,1.);
-  // fMasspjDeltaR->SetStats(kFALSE);
-  // fOutput->Add(fMasspjDeltaR);
-
-  TH3F* hzptD=new TH3F("hzptD",Form("Fragmentation function (DeltaR < %f)",fJetRadius),100,0.,1.2,nbinsmass,fMinMass,fMaxMass,100,0,50);
-  hzptD->SetStats(kTRUE);
-  hzptD->GetXaxis()->SetTitle("z=p_{D} #cdot p_{j}/|p_{j}|^{2}");
-  hzptD->GetYaxis()->SetTitle("mass (GeV)");
-  hzptD->GetZaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-  fOutput->Add(hzptD);
-
-  TH3F* hzptDB=new TH3F("hzptDB",Form("Fragmentation function (DeltaR < %f) - Side Bands",fJetRadius),100,0.,1.2,nbinsmass,fMinMass,fMaxMass,100,0.,50.);
-  hzptDB->SetStats(kTRUE);
-  hzptDB->GetXaxis()->SetTitle("z=p_{D} #cdot p_{j}/|p_{j}|^{2}");
-  hzptDB->GetYaxis()->SetTitle("mass (GeV)");
-  hzptDB->GetZaxis()->SetTitle("p_{t}^{D} (GeV/c)");
-  fOutput->Add(hzptDB);
-  //TH1F* hResZ      = new TH1F("hResZ","Fragmentation function ",50,0,1);
-  //  TH1F* hResZBkg   = new TH1F("hResZBkg","Fragmentation function background",50,0,1);  
-
-  //fOutput->Add(hResZ);
-  //fOutput->Add(hResZBkg);
-  PostData(1, fOutput);
-  
-  return kTRUE; 
-}
-
-void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliAODRecoDecayHF* candidate, AliEmcalJet *jet){
-
-  Double_t ptD=candidate->Pt();
-  Double_t ptjet=jet->Pt();
-  Double_t deltaR=DeltaR(candidate,jet);
-  Double_t deltaphi = jet->Phi()-candidate->Phi();
-  if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());
-  if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());
-  Double_t z=Z(candidate,jet);
-
-  TH1F* hDeltaRD=(TH1F*)fOutput->FindObject("hDeltaRD");
-  hDeltaRD->Fill(deltaR); 
-  TH1I* hstat=(TH1I*)fOutput->FindObject("hstat");
-  hstat->Fill(4);
-
-  if(fCandidateType==kD0toKpi) {
-
-    FillHistogramsD0JetCorr(candidate,deltaphi,z,ptD,ptjet,deltaR, AODEvent());
-
-  }
-
-  if(fCandidateType==kDstartoKpipi) {
-    AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
-    FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);
-
-  }
-
-}
-
-void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){
-  Double_t masses[2]={0.,0.};
-  Int_t pdgdaughtersD0[2]={211,321};//pi,K 
-  Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
-
-  masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
-  masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
-  
-  TH2F* hPtJetWithD=(TH2F*)fOutput->FindObject("hPtJetWithD");
-
-  Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
-  if(isselected==1 || isselected==3) {
-     
-     if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0]);
-
-     FillMassHistograms(masses[0],dPhi,z, ptD, deltaR);
-  }
-  if(isselected>=2) {
-     if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1]);
-
-     FillMassHistograms(masses[1],dPhi,z, ptD, deltaR);
-  }
-
-}
-
-void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){
-  AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
-  Double_t deltamass= dstar->DeltaInvMass(); 
-  Double_t massD0= dstar->InvMassD0();
-
-  TH1F* hPtPion=(TH1F*)fOutput->FindObject("hPtPion");
-  hPtPion->Fill(softpi->Pt());
-  
-  TH2F* hPtJetWithD=(TH2F*)fOutput->FindObject("hPtJetWithD");
-  if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass);
-  
-  FillMassHistograms(deltamass,dPhi,z, ptD, deltaR);
-  // evaluate side band background
-  TH2F* hDiffSideBand=(TH2F*)fOutput->FindObject("hDiffSideBand");
-
-  TH3F* hdeltaPhiDjaB=(TH3F*)fOutput->FindObject("hdeltaPhiDjaB");
-
-  TH3F* hzptDB=(TH3F*)fOutput->FindObject("hzptDB");
-
-  TH1F* hPtJetWithDsb=(TH1F*)fOutput->FindObject("hPtJetWithDsb");
-
-  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
-
-  Int_t bin = fCuts->PtBin(ptD);
-  Float_t fourSigmal= mPDGD0-3.5*fSigmaD0[bin] , sixSigmal= mPDGD0-5.*fSigmaD0[bin];
-  Float_t fourSigmar= mPDGD0+3.5*fSigmaD0[bin] , sixSigmar= mPDGD0+5.*fSigmaD0[bin];
-
-  if((massD0>sixSigmal && massD0<=fourSigmal) || (massD0>fourSigmar && massD0<=sixSigmar)){  
-    hDiffSideBand->Fill(deltamass,ptD); // M(Kpipi)-M(Kpi) side band background    
-    hdeltaPhiDjaB->Fill(deltamass,ptD,dPhi);
-
-    if(deltaR<fJetRadius){  // evaluate in the near side       
-      hzptDB->Fill(z,deltamass,ptD);
-      hPtJetWithDsb->Fill(ptj);
-    }
-    
-  }  //  SideBandBackground(dstar, dPhi, z, ptD, deltaR);
-}
-
-void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t dphi, Double_t z,Double_t ptD, Double_t deltaR){
-  TH3F* hdeltaPhiDja=((TH3F*)fOutput->FindObject("hdeltaPhiDja"));
-  hdeltaPhiDja->Fill(mass,ptD,dphi);
-
-  if(deltaR<fJetRadius) {
-    TH3F* hzptD=(TH3F*)fOutput->FindObject("hzptD");
-    hzptD->Fill(z,mass,ptD);
-
-    TH2F* hInvMassptD=(TH2F*)fOutput->FindObject("hInvMassptD");
-    hInvMassptD->Fill(mass,ptD);
-  }
-}
-//______________________________ side band background for D*___________________________________
-
-// void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t deltaR){
-
-//   //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas 
-//   // (expected detector resolution) on the left and right frm the D0 mass. Each band
-//   //  has a width of ~5 sigmas. Two band needed  for opening angle considerations   
-//   TH2F* hDiffSideBand=(TH2F*)fOutput->FindObject("hDiffSideBand");
-
-//   TH3F* hdeltaPhiDjaB=(TH3F*)fOutput->FindObject("hdeltaPhiDjaB");
-
-//   TH3F* hzptDB=(TH3F*)fOutput->FindObject("hzptDB");
-
-//   Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
-
-//   Int_t bin = fCuts->PtBin(candDstar->Pt());
-//   Float_t fourSigmal= mPDGD0-3.5*fSigmaD0[bin] , sixSigmal= mPDGD0-5.*fSigmaD0[bin];
-//   Float_t fourSigmar= mPDGD0+3.5*fSigmaD0[bin] , sixSigmar= mPDGD0+5.*fSigmaD0[bin];
-
-//   Double_t invM=candDstar->InvMassD0(),  deltaM=candDstar->DeltaInvMass(); //invMDstar=candDstar->InvMassDstarKpipi(),
-//   Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
-//   Double_t ptD=candDstar->Pt();
-//   //Double_t ptjet=jet->Pt();
-//   Double_t dPhi=jet->Phi()-candDstar->Phi();
-//   Double_t deltaR=DeltaR(candDstar,jet);
-//   if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
-//   if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
-
-//   if((invM>sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)){  
-//     hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background    
-//     hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
-
-//     if(deltaR<fJetRadius){  // evaluate in the near side    
-//       hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);    
-//     }
-
-//   }
-// }
-
-Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
-  //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
-
-  if(!p1 || !p2) return -1;
-  Double_t phi1=p1->Phi(),eta1=p1->Eta();
-  Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
-
-  Double_t dPhi=phi1-phi2;
-  if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
-  if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
-
-  Double_t dEta=eta1-eta2;
-  Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
-  return deltaR;
-}
-
-/*
-//_____________________________________________________________________
-
-Bool_t AliAnalysisTaskFlavourJetCorrelations::IsD(Int_t pdg) const{
-  Int_t abspdg=TMath::Abs(pdg);
-  if(abspdg>400 && abspdg<500) return kTRUE;
-  else return kFALSE;
-}
-
-//_____________________________________________________________________
-
-Bool_t AliAnalysisTaskFlavourJetCorrelations::IsD(Int_t pdg,Int_t abspdgD) const{
-  Int_t abspdg=TMath::Abs(pdg);
-  if(abspdg==abspdgD) return kTRUE;
-  else return kFALSE;
-}
-
-//_____________________________________________________________________
-
-Bool_t AliAnalysisTaskFlavourJetCorrelations::PartFromC(AliMCParticle* mother) const{
-  Int_t pdgmoth=mother->PdgCode();
-  if(TMath::Abs(pdgmoth)==4) return kTRUE;
-  else return kFALSE;
-}
-
-Int_t AliAnalysisTaskFlavourJetCorrelations::GetFirstMother(Int_t labpart, TClonesArray *mcArray) const{
-  AliAODMCParticle* partMC=(AliAODMCParticle*)mcArray->UncheckedAt(labpart);
-  if (!partMC) return -2;
-  Int_t labmom=labpart;
-  Printf("Starting from %d",labmom);
-  while(labmom>-1){
-    labmom=labpart;
-    partMC=(AliAODMCParticle*)mcArray->UncheckedAt(labmom);
-    if (!partMC) return -2;
-    labmom= partMC->GetMother();
-    Printf("Lab mom %d",labmom);
-  }
-  Printf("Return labpart %d", labpart);
-  return labpart;
-}
-
-
-// -------------------------------------- check the PDG -----------------------------------------
-
-Int_t AliAnalysisTaskFlavourJetCorrelations::FindPDGInFamily(Int_t labpart,Int_t pdgcode, TClonesArray *mcArray) const{
-
-  //return the label of the particle which is a "pdgcode" in the family
-  Printf("FindPDGInFamily label %d, pdg %d, mcarray %p",labpart,pdgcode,mcArray);
-  AliAODMCParticle* partMC=(AliAODMCParticle*)mcArray->UncheckedAt(labpart);
-  if (!partMC) return -2;
-  Int_t labmom=labpart;
-  Printf("Starting from %d",labmom);
-  while(labmom>-1){
-
-    partMC=(AliAODMCParticle*)mcArray->UncheckedAt(labmom);
-    if (!partMC) return -2;
-    Int_t pdgmom=partMC->GetPdgCode();
-    if(pdgmom==pdgcode) return labmom;
-    labmom= partMC->GetMother();
-    Printf("Lab mom %d",labmom);
-  }
-
-  return -1;
-
-}
-
-// ------------------------- check on MC the distance between D meson and jet ----------------------
-
-Bool_t AliAnalysisTaskFlavourJetCorrelations::FillMCDJetInfo(AliPicoTrack *jetTrk, AliEmcalJet* jet, TClonesArray *mcArray, Double_t ptjet){
-  
-  Bool_t foundD = kFALSE;
-  vector<int> DmesonInJetLabels(10);
-  
-  Int_t jtlabel=jetTrk->GetLabel();
-  if(jtlabel<=0) return foundD;
-  AliAODMCParticle* jetMCpart=(AliAODMCParticle*)mcArray->UncheckedAt(jtlabel);
-  if(!jetMCpart) return foundD;
-  printf("AliMCParticle %d, %p\n",1,jetMCpart);
-  
-  Int_t labDmeson=FindPDGInFamily(jtlabel,fPDGmother, mcArray);
-  if(labDmeson>0){
-    AliAODMCParticle *partDmeson=(AliAODMCParticle*)mcArray->UncheckedAt(labDmeson);
-    fhMomjetpartPdg->Fill(partDmeson->GetPdgCode());
-    
-    //tmp
-    Int_t momjetpartlabel=labDmeson;
-    
-    Int_t iD=5;
-    Bool_t exists=kFALSE;
-    for(Int_t k=0;k<iD;k++){
-      if(momjetpartlabel==DmesonInJetLabels[k]) {//mother already found
-       exists=kTRUE;
-       break;
-      }
-    }
-    if(!exists) {
-      DmesonInJetLabels[iD]=momjetpartlabel;
-      AliDebug(2,Form("D meson number %d found: label %d\n",iD,DmesonInJetLabels[iD]));
-      hstat->Fill(6);
-      iD++;
-      foundD=kTRUE;
-      
-    }
-  }
-  
-  if(fUseMCInfo && foundD) {
-    fptJwithDMC->Fill(ptjet); //filled only once per jet, not per each D meson
-    Int_t iD=5;
-    // loop over the D within the jet 
-    for(Int_t kD=0;kD<iD;kD++){
-      
-      AliAODMCParticle* momjetpart=(AliAODMCParticle*)mcArray->At(DmesonInJetLabels[kD]);
-      Double_t ptD=momjetpart->Pt(),etaD=momjetpart->Eta(), phiD=momjetpart->Phi();
-      fptDinjetallvsptjMC->Fill(ptD,ptjet);
-      
-      Double_t deltaRD=DeltaR(jet,momjetpart);
-      
-      ((TH1F*)fOutput->FindObject("hdeltaRDMC"))->Fill(deltaRD); //Delta R of D mesons (MC)
-      
-      Double_t z=Z(momjetpart,jet);
-      
-      if(deltaRD<fJetRadius) {
-       hstat->Fill(7);
-       //comment if you prefer to ask for DeltaR of the daughters < fJetRadius and uncomment below
-       fptDinjet04vsptjMC->Fill(ptD,ptjet);
-       fzptDptj->Fill(z,ptD,ptjet);
-      }
-      
-      
-    }//end loop on MC D
-    
-    return foundD;
-    
-  }
-}
-*/  
+/**************************************************************************\r
+* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+*                                                                        *\r
+* Author: The ALICE Off-line Project.                                    *\r
+* Contributors are mentioned in the code where appropriate.              *\r
+*                                                                        *\r
+* Permission to use, copy, modify and distribute this software and its   *\r
+* documentation strictly for non-commercial purposes is hereby granted   *\r
+* without fee, provided that the above copyright notice appears in all   *\r
+* copies and that both the copyright notice and this permission notice   *\r
+* appear in the supporting documentation. The authors make no claims     *\r
+* about the suitability of this software for any purpose. It is          *\r
+* provided "as is" without express or implied warranty.                  *\r
+**************************************************************************/\r
+//\r
+//\r
+//  Analysis Taks for reconstructed particle correlation \r
+//  (first implementation done for D mesons) with jets \r
+//  (use the so called Emcal framework)\r
+//\r
+//-----------------------------------------------------------------------\r
+// Authors:\r
+// C. Bianchin (Utrecht University) chiara.bianchin@cern.ch\r
+// A. Grelli (Utrecht University) a.grelli@uu.nl\r
+// X. Zhang (LBNL)  XMZhang@lbl.gov\r
+//-----------------------------------------------------------------------\r
+\r
+#include <TDatabasePDG.h>\r
+#include <TParticle.h>\r
+#include "TROOT.h"\r
+#include <TH3F.h>\r
+\r
+#include "AliAnalysisTaskFlavourJetCorrelations.h"\r
+#include "AliAODMCHeader.h"\r
+#include "AliAODHandler.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliLog.h"\r
+#include "AliEmcalJet.h"\r
+#include "AliJetContainer.h"\r
+#include "AliAODRecoDecay.h"\r
+#include "AliAODRecoCascadeHF.h"\r
+#include "AliAODRecoDecayHF2Prong.h"\r
+#include "AliESDtrack.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliPicoTrack.h"\r
+#include "AliRDHFCutsD0toKpi.h"\r
+#include "AliRDHFCutsDStartoKpipi.h"\r
+\r
+ClassImp(AliAnalysisTaskFlavourJetCorrelations)\r
+\r
+\r
+//_______________________________________________________________________________\r
+\r
+AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :\r
+AliAnalysisTaskEmcalJet("",kFALSE),\r
+fUseMCInfo(kTRUE), \r
+fUseReco(kTRUE),\r
+fCandidateType(),\r
+fPDGmother(),\r
+fNProngs(),\r
+fPDGdaughters(),\r
+fBranchName(),\r
+fmyOutput(0),\r
+fCuts(0),\r
+fMinMass(),\r
+fMaxMass(),  \r
+fJetArrName(0),\r
+fCandArrName(0),\r
+fLeadingJetOnly(kFALSE),\r
+fJetRadius(0)\r
+{\r
+   //\r
+   // Default ctor\r
+   //\r
+   //SetMakeGeneralHistograms(kTRUE);\r
+   \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :\r
+AliAnalysisTaskEmcalJet(name,kFALSE),\r
+fUseMCInfo(kTRUE),\r
+fUseReco(kTRUE),  \r
+fCandidateType(),\r
+fPDGmother(),\r
+fNProngs(),\r
+fPDGdaughters(),\r
+fBranchName(),\r
+fmyOutput(0),\r
+fCuts(0),\r
+fMinMass(),\r
+fMaxMass(),  \r
+fJetArrName(0),\r
+fCandArrName(0),\r
+fLeadingJetOnly(kFALSE),\r
+fJetRadius(0)\r
+{\r
+   //\r
+   // Constructor. Initialization of Inputs and Outputs\r
+   //\r
+   \r
+   Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");\r
+   fCuts=cuts;\r
+   fCandidateType=candtype;\r
+   const Int_t nptbins=fCuts->GetNPtBins();\r
+   Float_t defaultSigmaD013[13]={0.012, 0.012, 0.012, 0.015, 0.015,0.018,0.018,0.020,0.020,0.030,0.030,0.037,0.040};\r
+   \r
+   switch(fCandidateType){\r
+   case 0:\r
+      fPDGmother=421;\r
+      fNProngs=2;\r
+      fPDGdaughters[0]=211;//pi \r
+      fPDGdaughters[1]=321;//K\r
+      fPDGdaughters[2]=0; //empty\r
+      fPDGdaughters[3]=0; //empty\r
+      fBranchName="D0toKpi";\r
+      fCandArrName="D0";\r
+      break;\r
+   case 1: \r
+      fPDGmother=413;\r
+      fNProngs=3;\r
+      fPDGdaughters[1]=211;//pi soft\r
+      fPDGdaughters[0]=421;//D0\r
+      fPDGdaughters[2]=211;//pi fromD0\r
+      fPDGdaughters[3]=321; // K from D0\r
+      fBranchName="Dstar";\r
+      fCandArrName="Dstar";\r
+      \r
+      if(nptbins<=13){\r
+        for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];\r
+      } else {\r
+        AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));\r
+      }\r
+      break;\r
+   default:\r
+      printf("%d not accepted!!\n",fCandidateType);\r
+      break;\r
+   }\r
+   \r
+   if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);\r
+   if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);\r
+   \r
+   DefineOutput(1,TList::Class()); // histos\r
+   DefineOutput(2,AliRDHFCuts::Class()); // my cuts\r
+   \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {\r
+   //\r
+   // destructor\r
+   //\r
+   \r
+   Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");  \r
+   \r
+   delete fmyOutput;\r
+   delete fCuts;\r
+   \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::Init(){\r
+   //\r
+   // Initialization\r
+   //\r
+   \r
+   if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");\r
+   switch(fCandidateType){\r
+   case 0:\r
+      {\r
+        AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));\r
+        copyfCuts->SetName("AnalysisCutsDzero");\r
+        // Post the data\r
+        PostData(2,copyfCuts);\r
+      }\r
+      break;\r
+   case 1:\r
+      {\r
+        AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));\r
+        copyfCuts->SetName("AnalysisCutsDStar");\r
+        // Post the cuts\r
+        PostData(2,copyfCuts);\r
+      }\r
+      break;\r
+   default:\r
+      return;\r
+   }\r
+   \r
+   return;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() { \r
+   // output \r
+   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());\r
+   fmyOutput = new TList();\r
+   fmyOutput->SetOwner();\r
+   fmyOutput->SetName("pippo");\r
+   // define histograms\r
+   DefineHistoForAnalysis();\r
+   PostData(1,fmyOutput);\r
+   \r
+   return;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::UserExec(Option_t *)\r
+{\r
+   // user exec\r
+   if (!fInitialized){\r
+      AliAnalysisTaskEmcalJet::ExecOnce();\r
+   }\r
+   \r
+   // Load the event\r
+   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);\r
+   \r
+   TClonesArray *arrayDStartoD0pi=0;\r
+   TClonesArray *trackArr = 0;\r
+   TClonesArray *candidatesArr = 0;\r
+   TClonesArray *sbcandArr = 0;\r
+   \r
+   if (!aodEvent && AODEvent() && IsStandardAOD()) {\r
+      \r
+      // In case there is an AOD handler writing a standard AOD, use the AOD \r
+      // event in memory rather than the input (ESD) event.    \r
+      aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());\r
+      \r
+      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r
+      // have to taken from the AOD event hold by the AliAODExtension\r
+      AliAODHandler* aodHandler = (AliAODHandler*) \r
+      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r
+      if(aodHandler->GetExtensions()) {\r
+        AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r
+        AliAODEvent *aodFromExt = ext->GetAOD();\r
+        arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());\r
+      }\r
+   } else {\r
+      arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());\r
+   }\r
+   \r
+   if (!arrayDStartoD0pi) {\r
+      AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));\r
+      //  return;\r
+   } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   \r
+   \r
+   TClonesArray* mcArray = 0x0;\r
+   if (fUseMCInfo) {\r
+      mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));\r
+      if (!mcArray) {\r
+        printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");\r
+        return;\r
+      }\r
+   }\r
+   \r
+   //retrieve jets\r
+   trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));\r
+   //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));\r
+   //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));\r
+   fJetRadius=GetJetContainer(0)->GetJetRadius();\r
+   \r
+   if(!trackArr){\r
+      AliInfo(Form("Could not find the track array\n"));\r
+      return;\r
+   }\r
+   \r
+   \r
+   TString arrname="fCandidateArray";\r
+   TString candarrname=Form("%s%s%s",arrname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");\r
+   candidatesArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(candarrname));\r
+   if (!candidatesArr) {\r
+      Printf("%s array not found",candarrname.Data());\r
+      InputEvent()->GetList()->ls();\r
+      return;\r
+   }\r
+   //Printf("ncandidates found %d",candidatesArr->GetEntriesFast());\r
+   \r
+   TString arrSBname="fSideBandArray";\r
+   TString sbarrname=Form("%s%s%s",arrSBname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");\r
+   sbcandArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sbarrname));\r
+   if (fCandidateType==1 && !sbcandArr) {\r
+      Printf("%s array not found",sbarrname.Data());\r
+      InputEvent()->GetList()->ls();\r
+      return;\r
+   }\r
+   //Printf("nSBCand or Bkg found %d",sbcandArr->GetEntriesFast());\r
+   \r
+   \r
+   //Histograms\r
+   TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat");\r
+   TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks");\r
+   TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks");\r
+   TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks");\r
+   TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks");\r
+   TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet");\r
+   TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet");\r
+   TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet");\r
+   TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet");\r
+   TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr");\r
+   TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv");\r
+   TH1F* hdeltaRJetTracks=((TH1F*)fmyOutput->FindObject("hdeltaRJetTracks"));\r
+   \r
+   hstat->Fill(0);\r
+   \r
+   // fix for temporary bug in ESDfilter \r
+   // the AODs with null vertex pointer didn't pass the PhysSel\r
+   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;\r
+   \r
+   //Event selection\r
+   Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);\r
+   TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();\r
+   if(!iseventselected) return;\r
+   \r
+   hstat->Fill(1);\r
+   \r
+   //trigger on jets\r
+   \r
+   Int_t njets=GetJetContainer()->GetNJets();\r
+   if(njets == 0) return;\r
+\r
+   //retrieve charm candidates selected\r
+   Int_t candidates = candidatesArr->GetEntriesFast();\r
+   \r
+   // we start with jets\r
+   Double_t ejet   = 0;\r
+   Double_t phiJet = 0;\r
+   Double_t etaJet = 0;\r
+   Double_t ptjet = 0;\r
+   Double_t leadingJet =0;\r
+   \r
+   Int_t ntrarr=trackArr->GetEntriesFast();\r
+   hNtrArr->Fill(ntrarr);\r
+   \r
+   for(Int_t i=0;i<ntrarr;i++){\r
+      AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));\r
+      //reusing histograms\r
+      hPtJetTrks->Fill(jtrack->Pt());\r
+      hPhiJetTrks->Fill(jtrack->Phi());\r
+      hEtaJetTrks->Fill(jtrack->Eta());\r
+      hEjetTrks->Fill(jtrack->E());\r
+   }\r
+   \r
+   \r
+   //option to use only the leading jet\r
+   if(fLeadingJetOnly){\r
+      for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {    \r
+        AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);\r
+        ptjet   = jetL->Pt();\r
+        if(ptjet>leadingJet ) leadingJet = ptjet;\r
+        \r
+      }\r
+   }\r
+   \r
+   Int_t cntjet=0;\r
+   //loop over jets and consider only the leading jet in the event\r
+   for (Int_t iJets = 0; iJets<njets; iJets++) {    \r
+      //Printf("Jet N %d",iJets);\r
+      AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);\r
+      //Printf("Corr task Accept Jet");\r
+      if(!AcceptJet(jet)) {\r
+        hstat->Fill(5);\r
+        continue;\r
+      }\r
+      \r
+      //jets variables\r
+      ejet   = jet->E();\r
+      phiJet = jet->Phi();\r
+      etaJet = jet->Eta();\r
+      ptjet = jet->Pt();\r
+      \r
+      // choose the leading jet\r
+      if(fLeadingJetOnly && (ejet<leadingJet)) continue;\r
+      //used for normalization\r
+      hstat->Fill(3);\r
+      cntjet++;\r
+      // fill energy, eta and phi of the jet\r
+      hEjet   ->Fill(ejet);\r
+      hPhiJet ->Fill(phiJet);\r
+      hEtaJet ->Fill(etaJet);\r
+      hPtJet  ->Fill(ptjet);\r
+      \r
+      //loop on jet particles\r
+      Int_t ntrjet=  jet->GetNumberOfTracks();    \r
+      for(Int_t itrk=0;itrk<ntrjet;itrk++){\r
+        \r
+        AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);     \r
+        hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));\r
+        \r
+      }//end loop on jet tracks\r
+      \r
+      //Printf("N candidates %d ", candidates);\r
+      for(Int_t ic = 0; ic < candidates; ic++) {\r
+        \r
+        // D* candidates\r
+        AliVParticle* charm=0x0;\r
+        charm=(AliVParticle*)candidatesArr->At(ic);\r
+        if(!charm) continue;\r
+        hstat->Fill(2);\r
+        \r
+        FlagFlavour(charm, jet);\r
+        if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);\r
+        \r
+        FillHistogramsRecoJetCorr(charm, jet);\r
+        \r
+      }\r
+      //retrieve side band background candidates for Dstar\r
+      Int_t nsbcand = 0; if(sbcandArr) nsbcand=sbcandArr->GetEntriesFast();\r
+      \r
+      for(Int_t ib=0;ib<nsbcand;ib++){\r
+        if(fCandidateType==kDstartoKpipi && !fUseMCInfo){\r
+           AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)sbcandArr->At(ib);\r
+           if(!sbcand) continue;\r
+           SideBandBackground(sbcand,jet);\r
+        }\r
+        if(fUseMCInfo){\r
+           AliAODRecoDecayHF* charmbg = 0x0;\r
+           charmbg=(AliAODRecoDecayHF*)candidatesArr->At(ib);\r
+           if(!charmbg) continue;\r
+           MCBackground(charmbg,jet);\r
+        }\r
+      }\r
+   } // end of jet loop\r
+   \r
+   hNJetPerEv->Fill(cntjet);\r
+   \r
+   PostData(1,fmyOutput);\r
+   \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)\r
+{    \r
+   // The Terminate() function is the last function to be called during\r
+   // a query. It always runs on the client, it can be used to present\r
+   // the results graphically or save the results to file.\r
+   \r
+   Info("Terminate"," terminate");\r
+   AliAnalysisTaskSE::Terminate();\r
+   \r
+   fmyOutput = dynamic_cast<TList*> (GetOutputData(1));\r
+   if (!fmyOutput) {     \r
+      printf("ERROR: fmyOutput not available\n");\r
+      return;\r
+   }\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){\r
+   Float_t mass=0;\r
+   Int_t abspdg=TMath::Abs(pdg);\r
+   \r
+   mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();\r
+   // compute the Delta mass for the D*\r
+   if(fCandidateType==kDstartoKpipi){\r
+      Float_t mass1=0;\r
+      mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+      mass = mass-mass1;\r
+   }\r
+   \r
+   fMinMass = mass-range;\r
+   fMaxMass = mass+range;\r
+   \r
+   AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));\r
+   if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){\r
+   if(uplimit>lowlimit)\r
+   {\r
+      fMinMass = lowlimit;\r
+      fMaxMass = uplimit;\r
+   }\r
+   else{\r
+      printf("Error! Lower limit larger than upper limit!\n");\r
+      fMinMass = uplimit - uplimit*0.2;\r
+      fMaxMass = uplimit;\r
+   }\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){\r
+   if(nptbins>30) {\r
+      AliInfo("Maximum number of bins allowed is 30!");\r
+      return kFALSE;\r
+   }\r
+   if(!width) return kFALSE;\r
+   for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];\r
+   return kTRUE;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{\r
+   if(!part || !jet){\r
+      printf("Particle or jet do not exist!\n");\r
+      return -99;\r
+   }\r
+   Double_t p[3],pj[3];\r
+   Bool_t okpp=part->PxPyPz(p);\r
+   Bool_t okpj=jet->PxPyPz(pj);\r
+   if(!okpp || !okpj){\r
+      printf("Problems getting momenta\n");\r
+      return -999;\r
+   }\r
+   Double_t pjet=jet->P();\r
+   Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);\r
+   return z;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){\r
+   \r
+   // Statistics \r
+   TH1I* hstat=new TH1I("hstat","Statistics",6,-0.5,5.5);\r
+   hstat->GetXaxis()->SetBinLabel(1,"N ev anal");\r
+   hstat->GetXaxis()->SetBinLabel(2,"N ev sel");\r
+   hstat->GetXaxis()->SetBinLabel(3,"N cand sel cuts");\r
+   hstat->GetXaxis()->SetBinLabel(4,"N jets");\r
+   hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");\r
+   hstat->GetXaxis()->SetBinLabel(6,"N jet rej");\r
+   hstat->SetNdivisions(1);\r
+   fmyOutput->Add(hstat);\r
+   \r
+   const Int_t nbinsmass=200;\r
+   \r
+   \r
+   if(fCandidateType==kDstartoKpipi) \r
+   {\r
+      \r
+      TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,100, 0.,30.);\r
+      hDiffSideBand->SetStats(kTRUE);\r
+      hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");\r
+      hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");\r
+      fmyOutput->Add(hDiffSideBand); \r
+      \r
+      \r
+      TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);\r
+      hPtPion->SetStats(kTRUE);\r
+      hPtPion->GetXaxis()->SetTitle("GeV/c");\r
+      hPtPion->GetYaxis()->SetTitle("Entries");\r
+      fmyOutput->Add(hPtPion);\r
+      \r
+   }\r
+   \r
+   // jet related fistograms\r
+   \r
+   TH1F* hEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);\r
+   TH1F* hPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  200,0,6.30);\r
+   TH1F* hEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  100,-1.5,1.5);\r
+   TH1F* hPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",500,0,200);\r
+   \r
+   TH1F* hEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);\r
+   TH1F* hPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  200,0,6.30);\r
+   TH1F* hEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta",  100,-1.5,1.5);\r
+   TH1F* hPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",500,0,200);\r
+   \r
+   TH3F* hPtJetWithD=new TH3F("hPtJetWithD","D-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2})",500,0,200,nbinsmass,fMinMass,fMaxMass,100, 0.,30.);\r
+   //for the MC this histogram is filled with the real background\r
+   TH3F* hPtJetWithDsb=new TH3F("hPtJetWithDsb","D(background)-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2});p_{T}^{D} (GeV/c)",500,0,200,nbinsmass,fMinMass,fMaxMass,100, 0.,30.);\r
+   \r
+   TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);\r
+   \r
+   TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);\r
+   TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);\r
+   \r
+   fmyOutput->Add(hEjetTrks);\r
+   fmyOutput->Add(hPhiJetTrks);\r
+   fmyOutput->Add(hEtaJetTrks);\r
+   fmyOutput->Add(hPtJetTrks);\r
+   fmyOutput->Add(hEjet);\r
+   fmyOutput->Add(hPhiJet);\r
+   fmyOutput->Add(hEtaJet);\r
+   fmyOutput->Add(hPtJet);\r
+   fmyOutput->Add(hPtJetWithD);\r
+   fmyOutput->Add(hPtJetWithDsb);\r
+   fmyOutput->Add(hdeltaRJetTracks);\r
+   fmyOutput->Add(hNtrArr);\r
+   fmyOutput->Add(hNJetPerEv);\r
+   \r
+   TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);\r
+   fmyOutput->Add(hDeltaRD);\r
+   \r
+   //background (side bands for the Dstar and like sign for D0)\r
+   TH2F* hInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,100,0.,50.);\r
+   hInvMassptD->SetStats(kTRUE);\r
+   hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");\r
+   hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");\r
+   \r
+   fmyOutput->Add(hInvMassptD);\r
+   \r
+   if(fUseMCInfo){\r
+      TH2F* hInvMassptDbg = new TH2F("hInvMassptDbg",Form("Background D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,100,0.,50.);\r
+      hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));\r
+      hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");\r
+      fmyOutput->Add(hInvMassptDbg);\r
+      \r
+   }\r
+   PostData(1, fmyOutput);\r
+   \r
+   return kTRUE; \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){\r
+   \r
+   Double_t ptD=candidate->Pt();\r
+   Double_t ptjet=jet->Pt();\r
+   Double_t deltaR=DeltaR(candidate,jet);\r
+   Double_t deltaphi = jet->Phi()-candidate->Phi();\r
+   if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());\r
+   if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());\r
+   Double_t z=Z(candidate,jet);\r
+   \r
+   TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD");\r
+   hDeltaRD->Fill(deltaR); \r
+   if(fUseReco){\r
+      if(fCandidateType==kD0toKpi) {\r
+        AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;\r
+        FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent());\r
+        \r
+      }\r
+      \r
+      if(fCandidateType==kDstartoKpipi) {\r
+        AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;\r
+        FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);\r
+        \r
+      }\r
+   } else{ //generated level\r
+      //AliInfo("Non reco");\r
+      FillHistogramsMCGenDJetCorr(ptD,ptjet,deltaR);\r
+   }\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){\r
+\r
+  //dPhi and z not used at the moment,but will be (re)added\r
+\r
+   Double_t masses[2]={0.,0.};\r
+   Int_t pdgdaughtersD0[2]={211,321};//pi,K \r
+   Int_t pdgdaughtersD0bar[2]={321,211};//K,pi \r
+   \r
+   masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0\r
+   masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar\r
+   \r
+   TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");\r
+   \r
+   Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);\r
+   if(isselected==1 || isselected==3) {\r
+      \r
+      if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);\r
+      \r
+      FillMassHistograms(masses[0], ptD, deltaR);\r
+   }\r
+   if(isselected>=2) {\r
+      if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);\r
+      \r
+      FillMassHistograms(masses[1], ptD, deltaR);\r
+   }\r
+   \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){\r
+  //dPhi and z not used at the moment,but will be (re)added\r
+\r
+   AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();\r
+   Double_t deltamass= dstar->DeltaInvMass(); \r
+   //Double_t massD0= dstar->InvMassD0();\r
+   \r
+   \r
+   TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion");\r
+   hPtPion->Fill(softpi->Pt());\r
+   \r
+   TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");\r
+   if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);\r
+   \r
+   FillMassHistograms(deltamass, ptD, deltaR);\r
+   \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t ptD,Double_t ptjet,Double_t deltaR){\r
+   \r
+   Double_t pdgmass=0;\r
+   TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");\r
+   if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+   if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+   \r
+   if(deltaR<fJetRadius) hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet\r
+   \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){\r
+   \r
+   if(deltaR<fJetRadius) {\r
+      TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD");\r
+      hInvMassptD->Fill(mass,ptD);\r
+   }\r
+}\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){\r
+   Double_t deltaR=DeltaR(charm, jet);\r
+   AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;\r
+   if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;\r
+   if (deltaR<fJetRadius) jet->AddFlavourTag(tag);\r
+   \r
+   return;\r
+   \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){\r
+   \r
+   //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas \r
+   // (expected detector resolution) on the left and right frm the D0 mass. Each band\r
+   //  has a width of ~5 sigmas. Two band needed  for opening angle considerations   \r
+   TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand");\r
+   TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");\r
+   \r
+   Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+   \r
+   Int_t bin = fCuts->PtBin(candDstar->Pt());\r
+   Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];\r
+   Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];\r
+   \r
+   Double_t invM=candDstar->InvMassD0(),  deltaM=candDstar->DeltaInvMass(); \r
+   //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);\r
+   Double_t ptD=candDstar->Pt();\r
+   Double_t ptjet=jet->Pt();\r
+   Double_t dPhi=jet->Phi()-candDstar->Phi();\r
+   Double_t deltaR=DeltaR(candDstar,jet);\r
+   if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();\r
+   if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();\r
+   \r
+   //fill the background histogram with the side bands or when looking at MC with the real background\r
+   if((invM>sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)){  \r
+      hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background    \r
+      //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);\r
+      \r
+      if(deltaR<fJetRadius){  // evaluate in the near side     \r
+        //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);\r
+        hPtJetWithDsb->Fill(ptjet,deltaM,ptD);\r
+      }\r
+   }\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){\r
+   \r
+   //need mass, deltaR, pt jet, ptD\r
+   //two cases: D0 and Dstar\r
+   \r
+   Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());\r
+   if(!isselected) return;\r
+   \r
+   Double_t ptD=candbg->Pt();\r
+   Double_t ptjet=jet->Pt();\r
+   Double_t deltaR=DeltaR(candbg,jet);\r
+   \r
+   TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg");\r
+   TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");\r
+   \r
+   \r
+   if(fCandidateType==kDstartoKpipi){\r
+      AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;\r
+      Double_t deltaM=dstarbg->DeltaInvMass();\r
+      hInvMassptDbg->Fill(deltaM,ptD);\r
+      if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);\r
+   }\r
+   \r
+   if(fCandidateType==kD0toKpi){\r
+      Double_t masses[2]={0.,0.};\r
+      Int_t pdgdaughtersD0[2]={211,321};//pi,K \r
+      Int_t pdgdaughtersD0bar[2]={321,211};//K,pi \r
+      \r
+      masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0\r
+      masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar\r
+      \r
+      \r
+      if(isselected==1 || isselected==3) {\r
+        if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);\r
+        hInvMassptDbg->Fill(masses[0],ptD);\r
+      }\r
+      if(isselected>=2) {\r
+        if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);\r
+        hInvMassptDbg->Fill(masses[1],ptD);\r
+      }\r
+      \r
+      \r
+   }\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {\r
+   //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)\r
+   \r
+   if(!p1 || !p2) return -1;\r
+   Double_t phi1=p1->Phi(),eta1=p1->Eta();\r
+   Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;\r
+   \r
+   Double_t dPhi=phi1-phi2;\r
+   if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());\r
+   if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());\r
+   \r
+   Double_t dEta=eta1-eta2;\r
+   Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );\r
+   return deltaR;\r
+   \r
+}\r
+\r
+\r
index ee90a2a..00cb429 100644 (file)
-#ifndef ALIANALYSISTASKFLAVOURJETCORRELATIONS_H
-#define ALIANALYSISTASKFLAVOURJETCORRELATIONS_H
-/**************************************************************************
- * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-//-----------------------------------------------------------------------
-// Author : A. Grelli,  Utrecht University
-//          C. Bianchin, Utrecht University
-//          X. Zhang, LBNL
-//-----------------------------------------------------------------------
-
-
-#include <TH2F.h>
-#include "AliAODEvent.h"
-#include "AliPicoTrack.h"
-#include "AliAnalysisTaskEmcalJet.h"
-
-class TH3F;
-class TParticle ;
-class TClonesArray ;
-class AliMCParticle;
-class AliAODMCParticle;
-class AliRDHFCuts;
-class AliEmcalJet;
-class AliAODRecoDecayHF;
-class AliAODRecoCascadeHF;
-class AliAODEvent;
-
-class AliAnalysisTaskFlavourJetCorrelations : public AliAnalysisTaskEmcalJet 
-{
-
- public:
-
-  enum ECandidateType{ kD0toKpi, kDstartoKpipi };
-
-  AliAnalysisTaskFlavourJetCorrelations();
-  AliAnalysisTaskFlavourJetCorrelations(const Char_t* name,AliRDHFCuts* cuts, ECandidateType candtype, TString jetArrName);
-  virtual ~AliAnalysisTaskFlavourJetCorrelations();
-
-  virtual void     UserCreateOutputObjects();
-  virtual void     UserExec(Option_t *option);
-  virtual void     Terminate(Option_t *);
-  virtual void     Init();
-  virtual void     LocalInit() {Init();}
-
-  // inizializations
-  Bool_t DefineHistoForAnalysis();  
-
-  // set MC usage
-  void   SetMC(Bool_t theMCon) {fUseMCInfo = theMCon;}
-  Bool_t GetMC() const {return fUseMCInfo;}
-  
-  void SetMassLimits(Double_t range, Int_t pdg);
-  void SetMassLimits(Double_t lowlimit, Double_t uplimit);
-
-  //jet reconstruction algorithm
-  void SetJetArrayName(TString jetArrName) {fJetArrName=jetArrName;};
-  TString GetJetArrayName() const {return fJetArrName;};
-
-  // trigger on jet events
-  void SetTriggerOnLeadingJet(Bool_t triggerOnLeadingJet) {fLeadingJetOnly=triggerOnLeadingJet;};
-  Bool_t GetTriggerOnLeadingJet() const {return fLeadingJetOnly;}
-
-
-  // Array of D0 width for the Dstar
-  Bool_t SetD0WidthForDStar(Int_t nptbins,Float_t* width);
-
-  //Bool_t   FillMCDJetInfo(AliPicoTrack *jetTrk,AliEmcalJet* jet, TClonesArray *mcArray,Double_t ptjet);
-  void FillHistogramsRecoJetCorr(AliAODRecoDecayHF* candidate, AliEmcalJet *jet);
-  void FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR, AliAODEvent* aodEvent);
-
-  void FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR);
-
-  void FillMassHistograms(Double_t mass,Double_t dphi, Double_t z,Double_t ptD, Double_t deltaR);
- private:
-  
-  AliAnalysisTaskFlavourJetCorrelations(const AliAnalysisTaskFlavourJetCorrelations &source);
-  AliAnalysisTaskFlavourJetCorrelations& operator=(const AliAnalysisTaskFlavourJetCorrelations& source); 
-
-  Double_t Z(AliVParticle* part,AliEmcalJet* jet) const;
-  Float_t DeltaR(AliVParticle *p1, AliVParticle *p2) const;
-
-  /*
-  Bool_t IsD(Int_t pdg) const;
-  Bool_t IsD(Int_t pdg,Int_t abspdgD) const;
-  Bool_t PartFromC(AliMCParticle* mother) const;
-  Int_t GetFirstMother(Int_t lab,TClonesArray* mcarr) const; 
-  Int_t FindPDGInFamily(Int_t labpart,Int_t pdgcode, TClonesArray *mcArray) const;
-  */
-
-  Bool_t fUseMCInfo;             //  Use MC info
-  Int_t  fCandidateType;         // Dstar or D0
-  Int_t  fPDGmother;             // PDG code of D meson
-  Int_t  fNProngs;               // number of prong of the decay channel  
-  Int_t  fPDGdaughters[4];       // PDG codes of daughters
-  Float_t fSigmaD0[30];          //
-  TString fBranchName;           // AOD branch name
-  TList *fOutput;                //! user output
-  AliRDHFCuts *fCuts;            // Cuts 
-
-  Double_t fMinMass;             // mass lower limit histogram
-  Double_t fMaxMass;             // mass upper limit histogram
-
-  TString  fJetArrName;          // name of the jet array, taken from the task running the jet finder
-  TString fCandArrName;          // string which correspond to the candidate type
-  Bool_t fLeadingJetOnly;        // use only the leading jet in the event to make the correlations
-  Double_t fJetRadius;           // jet radius (filled from the JetContainer)
-
-  ClassDef(AliAnalysisTaskFlavourJetCorrelations,2); // class for charm-jet correlations
-};
-
-#endif
+#ifndef ALIANALYSISTASKFLAVOURJETCORRELATIONS_H\r
+#define ALIANALYSISTASKFLAVOURJETCORRELATIONS_H\r
+/**************************************************************************\r
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+//-----------------------------------------------------------------------\r
+// Author : A. Grelli,  Utrecht University\r
+//          C. Bianchin, Utrecht University\r
+//          X. Zhang, LBNL\r
+//-----------------------------------------------------------------------\r
+\r
+\r
+#include <TH2F.h>\r
+#include "AliAODEvent.h"\r
+#include "AliPicoTrack.h"\r
+#include "AliAnalysisTaskEmcalJet.h"\r
+\r
+class TH3F;\r
+class TParticle ;\r
+class TClonesArray ;\r
+class AliMCParticle;\r
+class AliAODMCParticle;\r
+class AliRDHFCuts;\r
+class AliEmcalJet;\r
+class AliAODRecoDecayHF;\r
+class AliAODRecoCascadeHF;\r
+class AliAODEvent;\r
+\r
+class AliAnalysisTaskFlavourJetCorrelations : public AliAnalysisTaskEmcalJet \r
+{\r
+\r
+ public:\r
+\r
+  enum ECandidateType{ kD0toKpi, kDstartoKpipi };\r
+\r
+  AliAnalysisTaskFlavourJetCorrelations();\r
+  AliAnalysisTaskFlavourJetCorrelations(const Char_t* name,AliRDHFCuts* cuts, ECandidateType candtype);\r
+  virtual ~AliAnalysisTaskFlavourJetCorrelations();\r
+\r
+  virtual void     UserCreateOutputObjects();\r
+  virtual void     UserExec(Option_t *option);\r
+  virtual void     Terminate(Option_t *);\r
+  virtual void     Init();\r
+  virtual void     LocalInit() {Init();}\r
+\r
+  // inizializations\r
+  Bool_t DefineHistoForAnalysis();  \r
+\r
+  // set MC usage\r
+  void   SetMC(Bool_t theMCon) {fUseMCInfo = theMCon;}\r
+  Bool_t GetMC() const {return fUseMCInfo;}\r
+  // set usage of reconstructed tracks\r
+  void   SetUseReco(Bool_t reco) {fUseReco=reco;}\r
+  Bool_t GetUseReco() {return fUseReco;}\r
+  \r
+  \r
+  void SetMassLimits(Double_t range, Int_t pdg);\r
+  void SetMassLimits(Double_t lowlimit, Double_t uplimit);\r
+\r
+  //jet reconstruction algorithm\r
+  void SetJetArrayName(TString jetArrName) {fJetArrName=jetArrName;};\r
+  TString GetJetArrayName() const {return fJetArrName;};\r
+\r
+  // trigger on jet events\r
+  void SetTriggerOnLeadingJet(Bool_t triggerOnLeadingJet) {fLeadingJetOnly=triggerOnLeadingJet;};\r
+  Bool_t GetTriggerOnLeadingJet() const {return fLeadingJetOnly;}\r
+\r
+\r
+  // Array of D0 width for the Dstar\r
+  Bool_t SetD0WidthForDStar(Int_t nptbins,Float_t* width);\r
+\r
+  //Bool_t   FillMCDJetInfo(AliPicoTrack *jetTrk,AliEmcalJet* jet, TClonesArray *mcArray,Double_t ptjet);\r
+  void FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet);\r
+  void FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR, AliAODEvent* aodEvent);\r
+\r
+  void FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR);\r
+  void FillHistogramsMCGenDJetCorr(Double_t ptD,Double_t ptjet,Double_t deltaR);\r
+  void SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet);\r
+  void MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet);\r
+  void FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR);\r
+  void FlagFlavour(AliVParticle* charm, AliEmcalJet* jet);\r
+  \r
+ private:\r
+  \r
+  AliAnalysisTaskFlavourJetCorrelations(const AliAnalysisTaskFlavourJetCorrelations &source);\r
+  AliAnalysisTaskFlavourJetCorrelations& operator=(const AliAnalysisTaskFlavourJetCorrelations& source); \r
+\r
+  Double_t Z(AliVParticle* part,AliEmcalJet* jet) const;\r
+  Float_t DeltaR(AliVParticle *p1, AliVParticle *p2) const;\r
+\r
+\r
+  Bool_t fUseMCInfo;             //  Use MC info\r
+  Bool_t fUseReco;               // use reconstructed tracks when running on MC\r
+  Int_t  fCandidateType;         // Dstar or D0\r
+  Int_t  fPDGmother;             // PDG code of D meson\r
+  Int_t  fNProngs;               // number of prong of the decay channel  \r
+  Int_t  fPDGdaughters[4];       // PDG codes of daughters\r
+  Float_t fSigmaD0[30];          //\r
+  TString fBranchName;           // AOD branch name\r
+  TList *fmyOutput;                //! user output\r
+  AliRDHFCuts *fCuts;            // Cuts \r
+\r
+  Double_t fMinMass;             // mass lower limit histogram\r
+  Double_t fMaxMass;             // mass upper limit histogram\r
+\r
+  TString  fJetArrName;          // name of the jet array, taken from the task running the jet finder\r
+  TString fCandArrName;          // string which correspond to the candidate type\r
+  Bool_t fLeadingJetOnly;        // use only the leading jet in the event to make the correlations\r
+  Double_t fJetRadius;           // jet radius (filled from the JetContainer)\r
+\r
+  ClassDef(AliAnalysisTaskFlavourJetCorrelations,2); // class for charm-jet correlations\r
+};\r
+\r
+#endif\r
index 69474e1..19f66d9 100644 (file)
-/**************************************************************************
- * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-//
-//
-// Task for selecting D mesons to be used as an input for D-Jet correlations
-//
-//-----------------------------------------------------------------------
-// Authors:
-// C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
-// A.Grelli (Utrecht University) a.grelli@uu.nl
-// Xiaoming Zhang (LBNL)  XMZhang@lbl.gov
-//-----------------------------------------------------------------------
-
-#include <TDatabasePDG.h>
-#include <TParticle.h>
-#include <TVector3.h>
-#include "TROOT.h"
-#include <TH3F.h>
-
-#include "AliRDHFCutsDStartoKpipi.h"
-#include "AliRDHFCutsD0toKpi.h"
-#include "AliAODMCHeader.h"
-#include "AliAODHandler.h"
-#include "AliAnalysisManager.h"
-#include "AliLog.h"
-#include "AliAODVertex.h"
-#include "AliAODRecoDecay.h"
-#include "AliAODRecoCascadeHF.h"
-#include "AliAODRecoDecayHF2Prong.h"
-#include "AliESDtrack.h"
-#include "AliAODMCParticle.h"
-#include "AliAnalysisTaskSEDmesonsFilterCJ.h"
-
-ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)
-
-//_____________________________________________________________________________
-AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
-  AliAnalysisTaskSE(),
-  fUseMCInfo(kFALSE),
-  fCandidateType(0),
-  fCandidateName(""),
-  fPDGmother(0),
-  fNProngs(0),
-  fBranchName(""),
-  fOutput(0),
-//fOutputCandidates(0),
-  fCuts(0),
-  fMinMass(0.),
-  fMaxMass(0.),
-  fCandidateArray(0)
-//fIsSelectedArray(0)
-{
-//
-// Default constructor
-//
-
-  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
-  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
-}
-
-//_____________________________________________________________________________
-AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
-  AliAnalysisTaskSE(name),
-  fUseMCInfo(kFALSE),
-  fCandidateType(candtype),
-  fCandidateName(""),
-  fPDGmother(0),
-  fNProngs(0),
-//fPDGdaughters(),
-  fBranchName(""),
-  fOutput(0),
-//fOutputCandidates(0),
-  fCuts(cuts),
-  fMinMass(0.),
-  fMaxMass(0.),
-  fCandidateArray(0)
-//fIsSelectedArray(0)
-{
-//
-// Constructor. Initialization of Inputs and Outputs
-//
-  Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
-
-  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
-  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
-
-  const Int_t nptbins = fCuts->GetNPtBins();
-  Float_t defaultSigmaD013[13] = { 0.012, 0.012, 0.012, 0.015, 0.015, 0.018, 0.018, 0.020, 0.020, 0.030, 0.030, 0.037, 0.040 };
-
-  switch (fCandidateType) {
-  case 0 :
-    fCandidateName = "D0";
-    fPDGmother = 421;
-    fNProngs = 2;
-    fPDGdaughters[0] = 211;  // pi 
-    fPDGdaughters[1] = 321;  // K
-    fPDGdaughters[2] = 0;    // empty
-    fPDGdaughters[3] = 0;    // empty
-    fBranchName = "D0toKpi";
-    break;
-  case 1 :
-    fCandidateName = "Dstar";
-    fPDGmother = 413;
-    fNProngs = 3;
-    fPDGdaughters[1] = 211; // pi soft
-    fPDGdaughters[0] = 421; // D0
-    fPDGdaughters[2] = 211; // pi fromD0
-    fPDGdaughters[3] = 321; // K from D0
-    fBranchName = "Dstar";
-
-    if (nptbins<=13) {
-      for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
-    } else {
-      AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
-    }
-    break;
-  default :
-    printf("%d not accepted!!\n",fCandidateType);
-    break;
-  }
-
-  if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
-  if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
-
-  DefineOutput(1, TList::Class());       // histos
-  DefineOutput(2, AliRDHFCuts::Class()); // my cuts
-}
-
-//_____________________________________________________________________________
-AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
-{
-//
-// destructor
-//
-
-  Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");  
-  if (fOutput) { delete fOutput; fOutput = 0; }
-  if (fCuts)   { delete fCuts;   fCuts   = 0; }
-  if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }
-//if (fIsSelectedArray) { delete fIsSelectedArray; fIsSelectedArray = 0; }
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEDmesonsFilterCJ::Init()
-{
-//
-// Initialization
-//
-
-  if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
-
-  switch (fCandidateType) {
-    case 0: {
-              AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
-              copyfCutsDzero->SetName("AnalysisCutsDzero");
-              PostData(2, copyfCutsDzero);  // Post the data
-            } break;
-    case 1: {
-              AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
-              copyfCutsDstar->SetName("AnalysisCutsDStar");
-              PostData(2, copyfCutsDstar); // Post the cuts
-            } break;
-    default: return;
-  }
-
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
-{ 
-//
-// output 
-//
-
-  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
-  
-//OpenFile(1);
-  fOutput = new TList(); fOutput->SetOwner();
-  DefineHistoForAnalysis(); // define histograms
-//fOutputCandidates= new TList();
-//fOutputCandidates->SetOwner();
-
-  fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
-  fCandidateArray->SetName(Form("fCandidateArray%s",fCandidateName.Data()));
-
-//fOutputCandidates->Add(fCandidateArray);
-//fIsSelectedArray = new TClonesArray("Int_t",0);
-//fIsSelectedArray->SetName(Form("fIsSelectedArray%s",suffix.Data()));
-//fOutputCandidates->Add(fIsSelectedArray);
-
-  PostData(1, fOutput);
-//PostData(3, fOutputCandidates);
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *)
-{
-//
-// user exec
-//
-
-  // add cadidate branch
-  fCandidateArray->Delete();
-  if (!(InputEvent()->FindListObject(Form("fCandidateArray%s",fCandidateName.Data())))) InputEvent()->AddObject(fCandidateArray);
-
-  // Load the event
-  AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
-  TClonesArray *arrayDStartoD0pi = 0;
-  if (!aodEvent && AODEvent() && IsStandardAOD()) {
-
-    // In case there is an AOD handler writing a standard AOD, use the AOD 
-    // event in memory rather than the input (ESD) event.    
-    aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
-
-    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
-    // have to taken from the AOD event hold by the AliAODExtension
-    AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
-    if(aodHandler->GetExtensions()) {
-      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
-      AliAODEvent *aodFromExt = ext->GetAOD();
-      arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
-    }
-  } else {
-    arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
-  }
-
-  if (!arrayDStartoD0pi) {
-    AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
-    return;
-  } else {
-    AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   
-  }
-  
-  TClonesArray* mcArray = 0x0;
-  if (fUseMCInfo) {
-    mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
-    if (!mcArray) {
-      printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
-      return;
-    }
-  }
-
-  //Histograms
-  TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
-  TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
-
-  TH1F* hPtPion=0x0;
-  if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
-  hstat->Fill(0);
-
-  // fix for temporary bug in ESDfilter 
-  // the AODs with null vertex pointer didn't pass the PhysSel
-  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
-
-  //Event selection
-  Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
-//TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
-  if (!iseventselected) return;
-  hstat->Fill(1);
-
-  const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
-  hstat->Fill(2, nD);
-
-//Int_t icountReco = 0;  // counters for efficiencies
-
-  //D* and D0 prongs needed to MatchToMC method
-  Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
-  Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi
-
-  //D0 from D0 bar
-  Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K 
-  Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi 
-
-  Int_t iCand =0;
-  Int_t isSelected = 0;
-  AliAODRecoDecayHF *charmCand = 0;
-
-  Int_t mcLabel = -9999;
-  Int_t pdgMeson = 413;
-  if (fCandidateType==kD0toKpi) pdgMeson = 421;
-
-  for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
-    charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
-    if (!charmCand) return;
-
-
-    if (fUseMCInfo) { // Look in MC, try to simulate the z
-      if (fCandidateType==kDstartoKpipi) {
-       AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
-       mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
-      }
-
-      if (fCandidateType==kD0toKpi) mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
-
-      if (mcLabel<=0) continue;               
-    }
-
-    // region of interest + cuts
-    if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(),charmCand->Y(pdgMeson))) continue;    
-
-    isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
-    if (!isSelected) continue;
-
-    new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
-//  new ((*fIsSelectedArray)[iCand]) Int_t(isSelected);
-    iCand++;
-
-    Double_t masses[2];
-    Double_t ptD = charmCand->Pt();
-    if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
-
-      //softpion from D* decay
-      AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
-      AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor();  
-      hstat->Fill(3);
-
-      // select D* in the D0 window.
-      // In the cut object window is loose to allow side bands
-      Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
-
-      // retrieve the corresponding pt bin for the candidate
-      // and set the expected D0 width (x3)
-//    static const Int_t n = fCuts->GetNPtBins();
-      Int_t bin = fCuts->PtBin(ptD);
-      if(bin<0 || bin>=fCuts->GetNPtBins()) {
-       AliError(Form("Pt %.3f out of bounds",ptD));
-       continue;
-      }
-
-      AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
-      if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {  
-       masses[0] = temp->DeltaInvMass(); //D*
-       masses[1] = 0.; //dummy for D*
-
-       //D*  delta mass
-       hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
-
-       // D* pt and soft pion pt for good candidates           
-       Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
-       Double_t invmassDelta = temp->DeltaInvMass();
-       if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
-      }
-    } //Dstar specific
-
-    if (fCandidateType==kD0toKpi) { //D0->Kpi
-
-      //needed quantities
-      masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
-      masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
-      hstat->Fill(3);
-
-      // mass vs pt
-      if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
-      if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
-    } //D0 specific
-
-    charmCand = 0;  
-  } // end of D cand loop
-
-//AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
-
-//PostData(1,fOutput);
-//PostData(3,fOutputCandidates);
-
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
-{
-// The Terminate() function is the last function to be called during
-// a query. It always runs on the client, it can be used to present
-// the results graphically or save the results to file.
-  
-  Info("Terminate"," terminate");
-  AliAnalysisTaskSE::Terminate();
-
-  fOutput = dynamic_cast<TList*>(GetOutputData(1));
-  if (!fOutput) {
-    printf("ERROR: fOutput not available\n");
-    return;
-  }
-
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
-{
-//
-// AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
-//
-
-  Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
-
-  // compute the Delta mass for the D*
-  if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
-//cout << "mass ---------------" << endl;
-
-  fMinMass = mass - range;
-  fMaxMass = mass + range;
-
-  AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
-  if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
-
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
-{
-//
-// AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
-//
-
-  if (uplimit>lowlimit) {
-    fMinMass = lowlimit;
-    fMaxMass = uplimit;
-  } else {
-    printf("Error! Lower limit larger than upper limit!\n");
-    fMinMass = uplimit - uplimit*0.2;
-    fMaxMass = uplimit;
-  }
-
-  return;
-}
-
-//_____________________________________________________________________________
-Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
-{
-//
-// AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
-//
-
-  if (nptbins>30) {
-    AliInfo("Maximum number of bins allowed is 30!");
-    return kFALSE;
-  }
-
-  if (!width) return kFALSE;
-  for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
-
-  return kTRUE;
-}
-
-
-//_____________________________________________________________________________
-Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
-{
-//
-// AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
-//
-
-  // Statistics 
-  TH1I* hstat = new TH1I("hstat","Statistics",4,-0.5,3.5);
-  hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
-  hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
-  hstat->GetXaxis()->SetBinLabel(3, "N cand");
-  hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
-/*if(fUseMCInfo) {
-    hstat->GetXaxis()->SetBinLabel(7,"N D");
-    hstat->GetXaxis()->SetBinLabel(8,"N D in jet");
-  }*/
-  hstat->SetNdivisions(1);
-  fOutput->Add(hstat);
-
-  // Invariant mass related histograms
-  const Int_t nbinsmass = 200;
-  TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
-  hInvMass->SetStats(kTRUE);
-  hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
-  hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
-  fOutput->Add(hInvMass);
-  
-  if (fCandidateType==kDstartoKpipi) {
-    TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
-    hPtPion->SetStats(kTRUE);
-    hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    hPtPion->GetYaxis()->SetTitle("entries");
-    fOutput->Add(hPtPion);
-  }
-
-  return kTRUE; 
-}
+/**************************************************************************\r
+* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+*                                                                        *\r
+* Author: The ALICE Off-line Project.                                    *\r
+* Contributors are mentioned in the code where appropriate.              *\r
+*                                                                        *\r
+* Permission to use, copy, modify and distribute this software and its   *\r
+* documentation strictly for non-commercial purposes is hereby granted   *\r
+* without fee, provided that the above copyright notice appears in all   *\r
+* copies and that both the copyright notice and this permission notice   *\r
+* appear in the supporting documentation. The authors make no claims     *\r
+* about the suitability of this software for any purpose. It is          *\r
+* provided "as is" without express or implied warranty.                  *\r
+**************************************************************************/\r
+//\r
+//\r
+// Task for selecting D mesons to be used as an input for D-Jet correlations\r
+//\r
+//-----------------------------------------------------------------------\r
+// Authors:\r
+// C. Bianchin (Utrecht University) chiara.bianchin@cern.ch\r
+// A.Grelli (Utrecht University) a.grelli@uu.nl\r
+// Xiaoming Zhang (LBNL)  XMZhang@lbl.gov\r
+//-----------------------------------------------------------------------\r
+\r
+#include <TDatabasePDG.h>\r
+#include <TParticle.h>\r
+#include <TVector3.h>\r
+#include "TROOT.h"\r
+#include <TH3F.h>\r
+\r
+#include "AliRDHFCutsDStartoKpipi.h"\r
+#include "AliRDHFCutsD0toKpi.h"\r
+#include "AliAODMCHeader.h"\r
+#include "AliAODHandler.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliLog.h"\r
+#include "AliAODVertex.h"\r
+#include "AliAODRecoDecay.h"\r
+#include "AliAODRecoCascadeHF.h"\r
+#include "AliAODRecoDecayHF2Prong.h"\r
+#include "AliESDtrack.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliAnalysisTaskSEDmesonsFilterCJ.h"\r
+\r
+ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)\r
+\r
+//_______________________________________________________________________________\r
+\r
+AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :\r
+AliAnalysisTaskSE(),\r
+fUseMCInfo(kFALSE),\r
+fUseReco(kTRUE),\r
+fCandidateType(0),\r
+fCandidateName(""),\r
+fPDGmother(0),\r
+fNProngs(0),\r
+fBranchName(""),\r
+fOutput(0),\r
+fCuts(0),\r
+fMinMass(0.),\r
+fMaxMass(0.),\r
+fCandidateArray(0)\r
+\r
+{\r
+   //\r
+   // Default constructor\r
+   //\r
+   \r
+   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;\r
+   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :\r
+AliAnalysisTaskSE(name),\r
+fUseMCInfo(kFALSE),\r
+fUseReco(kTRUE),\r
+fCandidateType(candtype),\r
+fCandidateName(""),\r
+fPDGmother(0),\r
+fNProngs(0),\r
+fBranchName(""),\r
+fOutput(0),\r
+fCuts(cuts),\r
+fMinMass(0.),\r
+fMaxMass(0.),\r
+fCandidateArray(0)\r
+{\r
+   //\r
+   // Constructor. Initialization of Inputs and Outputs\r
+   //\r
+   \r
+   Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");\r
+   \r
+   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;\r
+   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;\r
+   \r
+   const Int_t nptbins = fCuts->GetNPtBins();\r
+   Float_t defaultSigmaD013[13] = { 0.012, 0.012, 0.012, 0.015, 0.015, 0.018, 0.018, 0.020, 0.020, 0.030, 0.030, 0.037, 0.040 };\r
+   \r
+   switch (fCandidateType) {\r
+   case 0 :\r
+      fCandidateName = "D0";\r
+      fPDGmother = 421;\r
+      fNProngs = 2;\r
+      fPDGdaughters[0] = 211;  // pi \r
+      fPDGdaughters[1] = 321;  // K\r
+      fPDGdaughters[2] = 0;    // empty\r
+      fPDGdaughters[3] = 0;    // empty\r
+      fBranchName = "D0toKpi";\r
+      break;\r
+   case 1 :\r
+      fCandidateName = "Dstar";\r
+      fPDGmother = 413;\r
+      fNProngs = 3;\r
+      fPDGdaughters[1] = 211; // pi soft\r
+      fPDGdaughters[0] = 421; // D0\r
+      fPDGdaughters[2] = 211; // pi fromD0\r
+      fPDGdaughters[3] = 321; // K from D0\r
+      fBranchName = "Dstar";\r
+      \r
+      if (nptbins<=13) {\r
+        for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];\r
+      } else {\r
+        AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));\r
+      }\r
+      break;\r
+   default :\r
+      printf("%d not accepted!!\n",fCandidateType);\r
+      break;\r
+   }\r
+   \r
+   if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);\r
+   if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);\r
+   \r
+   DefineOutput(1, TList::Class());       // histos\r
+   DefineOutput(2, AliRDHFCuts::Class()); // my cuts\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()\r
+{\r
+   //\r
+   // destructor\r
+   //\r
+   \r
+   Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");  \r
+   \r
+   if (fOutput) { delete fOutput; fOutput = 0; }\r
+   if (fCuts)   { delete fCuts;   fCuts   = 0; }\r
+   if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }\r
+   \r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskSEDmesonsFilterCJ::Init()\r
+{\r
+   //\r
+   // Initialization\r
+   //\r
+   \r
+   if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");\r
+   \r
+   switch (fCandidateType) {\r
+   case 0: {\r
+        AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));\r
+        copyfCutsDzero->SetName("AnalysisCutsDzero");\r
+        PostData(2, copyfCutsDzero);  // Post the data\r
+   } break;\r
+case 1: {\r
+      AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));\r
+      copyfCutsDstar->SetName("AnalysisCutsDStar");\r
+      PostData(2, copyfCutsDstar); // Post the cuts\r
+} break;\r
+default: return;\r
+   }\r
+   \r
+   return;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()\r
+{ \r
+   //\r
+   // output \r
+   //\r
+   \r
+   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());\r
+   \r
+   fOutput = new TList(); fOutput->SetOwner();\r
+   DefineHistoForAnalysis(); // define histograms\r
+   \r
+   fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);\r
+   fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));\r
+   \r
+   if (fCandidateType==kDstartoKpipi){\r
+      fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0); //this is for the DStar only!\r
+      fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));\r
+   }\r
+   \r
+   PostData(1, fOutput);\r
+   return;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){\r
+   //\r
+   // user exec\r
+   //\r
+   \r
+   // add cadidate branch\r
+   fCandidateArray->Delete();\r
+   if (!(InputEvent()->FindListObject(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen")))) InputEvent()->AddObject(fCandidateArray);\r
+   if (fCandidateType==kDstartoKpipi){\r
+      fSideBandArray->Delete();\r
+      if (!(InputEvent()->FindListObject(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen")))) InputEvent()->AddObject(fSideBandArray);\r
+   }\r
+   //Printf("Arr names %s, %s",fCandidateArray->GetName(),fSideBandArray->GetName());\r
+   // Load the event\r
+   AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);\r
+   \r
+   TClonesArray *arrayDStartoD0pi = 0;\r
+   if (!aodEvent && AODEvent() && IsStandardAOD()) {\r
+      \r
+      // In case there is an AOD handler writing a standard AOD, use the AOD \r
+      // event in memory rather than the input (ESD) event.    \r
+      aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());\r
+      \r
+      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r
+      // have to taken from the AOD event hold by the AliAODExtension\r
+      AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r
+      if(aodHandler->GetExtensions()) {\r
+        AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r
+        AliAODEvent *aodFromExt = ext->GetAOD();\r
+        arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());\r
+      }\r
+   } else {\r
+      arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());\r
+   }\r
+   \r
+   if (!arrayDStartoD0pi) {\r
+      AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));\r
+      return;\r
+   } else {\r
+      AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   \r
+   }\r
+   \r
+   TClonesArray* mcArray = 0x0;\r
+   if (fUseMCInfo) {\r
+      mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));\r
+      if (!mcArray) {\r
+        printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");\r
+        return;\r
+      }\r
+   }\r
+   \r
+   //Histograms\r
+   TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");\r
+   TH1F* hnSBCandEv=(TH1F*)fOutput->FindObject("hnSBCandEv");\r
+   TH1F* hnCandEv=(TH1F*)fOutput->FindObject("hnCandEv");\r
+   TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");\r
+   \r
+   TH1F* hPtPion=0x0;\r
+   if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");\r
+   hstat->Fill(0);\r
+   \r
+   // fix for temporary bug in ESDfilter \r
+   // the AODs with null vertex pointer didn't pass the PhysSel\r
+   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;\r
+   \r
+   //Event selection\r
+   Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);\r
+   //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();\r
+   if (!iseventselected) return;\r
+   hstat->Fill(1);\r
+   \r
+   const Int_t nD = arrayDStartoD0pi->GetEntriesFast();\r
+   if(fUseReco) hstat->Fill(2, nD);\r
+   \r
+   //D* and D0 prongs needed to MatchToMC method\r
+   Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi\r
+   Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi\r
+   \r
+   //D0 from D0 bar\r
+   Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K \r
+   Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi \r
+   \r
+   Int_t iCand =0;\r
+   Int_t iSBCand=0;\r
+   Int_t isSelected = 0;\r
+   AliAODRecoDecayHF *charmCand = 0;\r
+   AliAODMCParticle *charmPart = 0;\r
+   Bool_t isMCBkg=kFALSE;\r
+   \r
+   Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+   \r
+   Int_t mcLabel = -9999;\r
+   Int_t pdgMeson = 413;\r
+   if (fCandidateType==kD0toKpi) pdgMeson = 421;\r
+   \r
+   for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates\r
+      charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates\r
+      if (!charmCand) continue;\r
+      \r
+      \r
+      if (fUseMCInfo) { // Look in MC, try to simulate the z\r
+        if (fCandidateType==kDstartoKpipi) {\r
+           AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;\r
+           mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);\r
+        }\r
+        \r
+        if (fCandidateType==kD0toKpi) \r
+           mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);\r
+        \r
+        if (mcLabel<=0) isMCBkg=kTRUE;\r
+        else hstat->Fill(2);\r
+        if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);\r
+      }\r
+      \r
+      Double_t ptD = charmCand->Pt();\r
+      \r
+      // region of interest + cuts\r
+      if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;    \r
+      \r
+      if(!fUseMCInfo && fCandidateType==kDstartoKpipi){\r
+        //select by track cuts the side band candidates (don't want mass cut)\r
+        isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent); \r
+        if (!isSelected) continue;\r
+        //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])\r
+        Int_t bin = fCuts->PtBin(ptD);\r
+        if(bin<0 || bin>=fCuts->GetNPtBins()) {\r
+           AliError(Form("Pt %.3f out of bounds",ptD));\r
+           continue;\r
+        }\r
+        AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;\r
+        //if data and Dstar from D0 side band\r
+        if (((temp->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/||  ((temp->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (temp->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){   \r
+           \r
+           new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*temp);\r
+           iSBCand++;\r
+        }\r
+      }\r
+      //candidate selected by cuts and PID\r
+      isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected\r
+      if (!isSelected) continue;\r
+      \r
+      //for data and MC signal fill fCandidateArray\r
+      if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){\r
+        // for data or MC with the requirement fUseReco fill with candidates\r
+        if(fUseReco) {\r
+           new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);\r
+           //Printf("Filling reco");\r
+           hstat->Fill(3);\r
+        }\r
+        // for MC with requirement particle level fill with AliAODMCParticle\r
+        else if (fUseMCInfo) {\r
+           new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);\r
+           //Printf("Filling gen");\r
+           hstat->Fill(3);\r
+        }\r
+        \r
+        iCand++;       \r
+      }\r
+      //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)\r
+      else if(fUseReco){\r
+        new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);\r
+        iSBCand++;\r
+      }\r
+      \r
+      \r
+      Double_t masses[2];\r
+      if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi\r
+        \r
+        //softpion from D* decay\r
+        AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;\r
+        AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor();  \r
+        \r
+        // select D* in the D0 window.\r
+        // In the cut object window is loose to allow for side bands\r
+        \r
+        \r
+        // retrieve the corresponding pt bin for the candidate\r
+        // and set the expected D0 width (x3)\r
+        //    static const Int_t n = fCuts->GetNPtBins();\r
+        Int_t bin = fCuts->PtBin(ptD);\r
+        if(bin<0 || bin>=fCuts->GetNPtBins()) {\r
+           AliError(Form("Pt %.3f out of bounds",ptD));\r
+           continue;\r
+        }\r
+        \r
+        AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));\r
+        //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value\r
+        if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {        \r
+           masses[0] = temp->DeltaInvMass(); //D*\r
+           masses[1] = 0.; //dummy for D*\r
+           \r
+           //D*  delta mass\r
+           hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins\r
+           \r
+           // D* pt and soft pion pt for good candidates               \r
+           Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();\r
+           Double_t invmassDelta = temp->DeltaInvMass();\r
+           if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());\r
+        }\r
+      } //Dstar specific\r
+      \r
+      if (fCandidateType==kD0toKpi) { //D0->Kpi\r
+        \r
+        //needed quantities\r
+        masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0\r
+        masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar\r
+        hstat->Fill(3);\r
+        \r
+        // mass vs pt\r
+        if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);\r
+        if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);\r
+      } //D0 specific\r
+      \r
+      charmCand = 0;\r
+      isMCBkg=kFALSE;\r
+   } // end of D cand loop\r
+   \r
+   hnCandEv->Fill(fCandidateArray->GetEntriesFast());\r
+   if (fCandidateType==kDstartoKpipi || fUseMCInfo) {\r
+      Int_t nsbcand=fSideBandArray->GetEntriesFast();\r
+      hstat->Fill(4,nsbcand);\r
+      hnSBCandEv->Fill(nsbcand);\r
+   }\r
+   \r
+   return;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)\r
+{\r
+   // The Terminate() function is the last function to be called during\r
+   // a query. It always runs on the client, it can be used to present\r
+   // the results graphically or save the results to file.\r
+   \r
+   Info("Terminate"," terminate");\r
+   AliAnalysisTaskSE::Terminate();\r
+   \r
+   fOutput = dynamic_cast<TList*>(GetOutputData(1));\r
+   if (!fOutput) {\r
+      printf("ERROR: fOutput not available\n");\r
+      return;\r
+   }\r
+   \r
+   return;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)\r
+{\r
+   //\r
+   // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits\r
+   //\r
+   \r
+   Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();\r
+   \r
+   // compute the Delta mass for the D*\r
+   if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+   \r
+   \r
+   fMinMass = mass - range;\r
+   fMaxMass = mass + range;\r
+   \r
+   AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));\r
+   if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");\r
+   \r
+   return;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)\r
+{\r
+   //\r
+   // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits\r
+   //\r
+   \r
+   if (uplimit>lowlimit) {\r
+      fMinMass = lowlimit;\r
+      fMaxMass = uplimit;\r
+   } else {\r
+      printf("Error! Lower limit larger than upper limit!\n");\r
+      fMinMass = uplimit - uplimit*0.2;\r
+      fMaxMass = uplimit;\r
+   }\r
+   \r
+   return;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)\r
+{\r
+   //\r
+   // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar\r
+   //\r
+   \r
+   if (nptbins>30) {\r
+      AliInfo("Maximum number of bins allowed is 30!");\r
+      return kFALSE;\r
+   }\r
+   \r
+   if (!width) return kFALSE;\r
+   for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];\r
+   \r
+   return kTRUE;\r
+}\r
+\r
+//_______________________________________________________________________________\r
+\r
+Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()\r
+{\r
+   //\r
+   // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis\r
+   //\r
+   \r
+   // Statistics \r
+   TH1I* hstat = new TH1I("hstat","Statistics",5,-0.5,4.5);\r
+   hstat->GetXaxis()->SetBinLabel(1, "N ev anal");\r
+   hstat->GetXaxis()->SetBinLabel(2, "N ev sel");\r
+   if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");\r
+   else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");\r
+   hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");\r
+   if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");  \r
+   hstat->SetNdivisions(1);\r
+   fOutput->Add(hstat);\r
+   \r
+   TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);\r
+   fOutput->Add(hnCandEv);\r
+   \r
+   // Invariant mass related histograms\r
+   const Int_t nbinsmass = 200;\r
+   TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);\r
+   hInvMass->SetStats(kTRUE);\r
+   hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");\r
+   hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");\r
+   fOutput->Add(hInvMass);\r
+   \r
+   if (fCandidateType==kDstartoKpipi) {\r
+      TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);\r
+      fOutput->Add(hnSBCandEv);\r
+      \r
+      TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);\r
+      hPtPion->SetStats(kTRUE);\r
+      hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");\r
+      hPtPion->GetYaxis()->SetTitle("entries");\r
+      fOutput->Add(hPtPion);\r
+   }\r
+   \r
+   return kTRUE; \r
+}\r
index eaae3cc..3370f78 100644 (file)
-#ifndef ALIANALYSISTASKSEDMESONSFILTERCJ_H
-#define ALIANALYSISTASKSEDMESONSFILTERCJ_H
-/**************************************************************************
- * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-//-----------------------------------------------------------------------
-// Author : A. Grelli,  Utrecht University
-//          C. Bianchin, Utrecht University
-//          X. Zhang, LBNL
-//-----------------------------------------------------------------------
-
-
-#include <TH2F.h>
-#include "AliAODEvent.h"
-#include "AliPicoTrack.h"
-#include "AliAnalysisTaskSE.h"
-
-class TH3F;
-class TString;
-class TParticle ;
-class TClonesArray ;
-class AliMCParticle;
-class AliAODMCParticle;
-class AliRDHFCuts;
-class AliAODRecoCascadeHF;
-
-class AliAnalysisTaskSEDmesonsFilterCJ : public AliAnalysisTaskSE 
-{
-
- public :
-
-  enum ECandidateType{ kD0toKpi, kDstartoKpipi };
-  
-  AliAnalysisTaskSEDmesonsFilterCJ();
-  AliAnalysisTaskSEDmesonsFilterCJ(const Char_t* name,AliRDHFCuts* cuts,ECandidateType candtype);
-  virtual ~AliAnalysisTaskSEDmesonsFilterCJ();
-
-  virtual void     UserCreateOutputObjects();
-  virtual void     UserExec(Option_t *option);
-  virtual void     Terminate(Option_t *);
-  virtual void     Init();
-  virtual void     LocalInit() { Init(); }
-
-  // inizializations
-  Bool_t DefineHistoForAnalysis();
-
-  // set MC usage
-  void   SetMC(Bool_t theMCon) { fUseMCInfo = theMCon; }
-  Bool_t GetMC() const { return fUseMCInfo; }
-  
-  void SetMassLimits(Double_t range, Int_t pdg);
-  void SetMassLimits(Double_t lowlimit, Double_t uplimit);
-
-  // Array of D0 width for the Dstar
-  Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width);
-
- private :
-  
-  AliAnalysisTaskSEDmesonsFilterCJ(const AliAnalysisTaskSEDmesonsFilterCJ &source);
-  AliAnalysisTaskSEDmesonsFilterCJ& operator=(const AliAnalysisTaskSEDmesonsFilterCJ& source); 
-
-  Bool_t fUseMCInfo;               //  Use MC info
-
-  UInt_t  fCandidateType;          //  Dstar or D0
-  TString fCandidateName;          //  Dstar or D0
-
-  Int_t fPDGmother;                //  PDG code of D meson
-  Int_t fNProngs;                  //  number of prong of the decay channel  
-  Int_t fPDGdaughters[4];          //  PDG codes of daughters
-  Float_t fSigmaD0[30];            //  D0 sigma for Dstar
-
-  TString fBranchName;             //  AOD branch name
-  TList  *fOutput;                 //! user output
-//TList *fOutputCandidates;        //! output of array of candidates (kExchange)
-
-  AliRDHFCuts *fCuts;              // Cuts 
-  Double_t fMinMass;               //  mass lower limit histogram
-  Double_t fMaxMass;               //  mass upper limit histogram
-
-  TClonesArray *fCandidateArray;   //! contains candidates selected by AliRDHFCuts
-//TClonesArray *fIsSelectedArray;  //! contains result of IsSelected for candidates which pass the cuts (needed for D0)
-
-  ClassDef(AliAnalysisTaskSEDmesonsFilterCJ,1); // class for charm-jet correlations
-};
-
-#endif
+#ifndef ALIANALYSISTASKSEDMESONSFILTERCJ_H\r
+#define ALIANALYSISTASKSEDMESONSFILTERCJ_H\r
+/**************************************************************************\r
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+//-----------------------------------------------------------------------\r
+// Author : A. Grelli,  Utrecht University\r
+//          C. Bianchin, Utrecht University\r
+//          X. Zhang, LBNL\r
+//-----------------------------------------------------------------------\r
+\r
+\r
+#include <TH2F.h>\r
+#include "AliAODEvent.h"\r
+#include "AliPicoTrack.h"\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class TH3F;\r
+class TString;\r
+class TParticle ;\r
+class TClonesArray ;\r
+class AliMCParticle;\r
+class AliAODMCParticle;\r
+class AliRDHFCuts;\r
+class AliAODRecoCascadeHF;\r
+\r
+class AliAnalysisTaskSEDmesonsFilterCJ : public AliAnalysisTaskSE \r
+{\r
+\r
+ public :\r
+\r
+  enum ECandidateType{ kD0toKpi, kDstartoKpipi };\r
+  \r
+  AliAnalysisTaskSEDmesonsFilterCJ();\r
+  AliAnalysisTaskSEDmesonsFilterCJ(const Char_t* name,AliRDHFCuts* cuts,ECandidateType candtype);\r
+  virtual ~AliAnalysisTaskSEDmesonsFilterCJ();\r
+\r
+  virtual void     UserCreateOutputObjects();\r
+  virtual void     UserExec(Option_t *option);\r
+  virtual void     Terminate(Option_t *);\r
+  virtual void     Init();\r
+  virtual void     LocalInit() { Init(); }\r
+\r
+  // inizializations\r
+  Bool_t DefineHistoForAnalysis();\r
+\r
+  // set MC usage\r
+  void   SetMC(Bool_t theMCon) { fUseMCInfo = theMCon; }\r
+  Bool_t GetMC() const { return fUseMCInfo; }\r
+  \r
+  // set usage of generated or reconstucted quantities (relevant for MC)\r
+  void SetUseReco(Bool_t useReco=kTRUE) { fUseReco= useReco;}\r
+  Bool_t GetUseReco() const {return fUseReco;}\r
\r
+  void SetMassLimits(Double_t range, Int_t pdg);\r
+  void SetMassLimits(Double_t lowlimit, Double_t uplimit);\r
+\r
+  // Array of D0 width for the Dstar\r
+  Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width);\r
+\r
+ private :\r
+  \r
+  AliAnalysisTaskSEDmesonsFilterCJ(const AliAnalysisTaskSEDmesonsFilterCJ &source);\r
+  AliAnalysisTaskSEDmesonsFilterCJ& operator=(const AliAnalysisTaskSEDmesonsFilterCJ& source); \r
+\r
+  Bool_t fUseMCInfo;               //  Use MC info\r
+  Bool_t fUseReco;                 // use reconstructed tracks when running on MC\r
+\r
+  UInt_t  fCandidateType;          //  Dstar or D0\r
+  TString fCandidateName;          //  Dstar or D0\r
+\r
+  Int_t fPDGmother;                //  PDG code of D meson\r
+  Int_t fNProngs;                  //  number of prong of the decay channel  \r
+  Int_t fPDGdaughters[4];          //  PDG codes of daughters\r
+  Float_t fSigmaD0[30];            //  D0 sigma for Dstar\r
+\r
+  TString fBranchName;             //  AOD branch name\r
+  TList  *fOutput;                 //! user output\r
+\r
+  AliRDHFCuts *fCuts;              // Cuts \r
+  Double_t fMinMass;               //  mass lower limit histogram\r
+  Double_t fMaxMass;               //  mass upper limit histogram\r
+\r
+  TClonesArray *fCandidateArray;   //! contains candidates selected by AliRDHFCuts\r
+  TClonesArray *fSideBandArray;    //! contains candidates selected by AliRDHFCuts::IsSelected(kTracks), to be used for side bands (DStar case only!!)\r
+\r
+  ClassDef(AliAnalysisTaskSEDmesonsFilterCJ,2); // class for charm-jet correlations\r
+};\r
+\r
+#endif\r
index cba20c0..0958e23 100644 (file)
-AliAnalysisTaskFlavourJetCorrelations *AddTaskFlavourJetCorrelations(
-  AliAnalysisTaskFlavourJetCorrelations::ECandidateType cand = AliAnalysisTaskFlavourJetCorrelations::kDstartoKpipi,
-  TString filename = "DStartoKpipiCuts.root",
-  Bool_t theMCon = kFALSE,
-  TString jetArrname = "",
-  TString suffix = "",
-  Bool_t triggerOnLeadingJet = kFALSE,
-  Float_t R = 0.4,
-  Float_t jptcut = 10.,
-  Int_t acceptance = 1 /*1= 0-2pi, 2=emcal cut*/,
-  Double_t areaCut = 0.)
-{
-  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
-  if (!mgr) {
-    ::Error("AddTaskFlavourJetCorrelations::AddTaskFlavourJetCorrelations", "No analysis manager to connect to.");
-    return NULL;
-  } 
-
-  Bool_t useStdC = kFALSE;
-  TFile* filecuts=TFile::Open(filename);
-  if (!filecuts || (filecuts && !filecuts->IsOpen())) {
-    cout<<"Input file not found: use std cuts"<<endl;
-    useStdC = kTRUE;
-  }
-
-  AliRDHFCuts *analysiscuts = 0x0;
-  switch (cand) {
-  case 0 :
-    if (useStdC) {
-      analysiscuts = new AliRDHFCutsD0toKpi();
-      analysiscuts->SetStandardCutsPP2010();
-    } else
-      analysiscuts = (AliRDHFCutsD0toKpi*)filecuts->Get("D0toKpiCuts");
-    break;
-  case 1 :
-    if(useStdC) {
-      analysiscuts = new AliRDHFCutsDStartoKpipi();
-      analysiscuts->SetStandardCutsPP2010();
-    } else
-      analysiscuts = (AliRDHFCutsDStartoKpipi*)filecuts->Get("DStartoKpipiCuts");
-    analysiscuts->SetName("DStartoKpipiCuts");
-    break;
-  }
-  
-  if (!analysiscuts) { // mm let's see if everything is ok
-    AliFatal("Specific AliRDHFCuts not found");
-    return;
-  }
-
-  printf("CREATE TASK\n"); //CREATE THE TASK
-
-  // create the task
-  AliAnalysisTaskFlavourJetCorrelations *task = new AliAnalysisTaskFlavourJetCorrelations("AnaTaskFlavourJetCorrelations", analysiscuts, cand, jetArrname);
-
-  //D meson settings
-  task->SetMC(theMCon);
-  task->SetTriggerOnLeadingJet(triggerOnLeadingJet);
-
-  //jet settings
-  task->SetJetRadius(R);
-  task->SetJetPtCut(jptcut);
-
-  Float_t etaCov=0.9;
-  if (acceptance==2) etaCov=0.7; //EMCal
-
-  Float_t minEta = -etaCov+R;
-  Float_t maxEta =  etaCov-R;
-  if (acceptance) task->SetJetEtaLimits(minEta, maxEta);
-
-  Float_t minPhi = 0.;
-  Float_t maxPhi = 2.*TMath::Pi();
-  if (acceptance==2) { /*80-180 degree*/ }
-  if (acceptance) task->SetJetPhiLimits(minPhi, maxPhi);
-
-  //Float_t area=0.6*TMath::Pi()*R*R;
-  task->SetJetAreaCut(areaCut);
-  
-  mgr->AddTask(task);
-
-  // Create and connect containers for input/output
-  TString outputfile = AliAnalysisManager::GetCommonFileName();
-  outputfile += ":PWG3_D2H_DEmcalJet";
-  outputfile += suffix;
-
-  TString candname="DStar"; 
-  if(cand==0)  candname="D0";
-
-  TString nameContainer1="hCor";
-  TString nameContainer2="cutsJ";
-
-  nameContainer1 += candname;
-  nameContainer2 += candname;
-  
-  nameContainer1 += suffix;
-  nameContainer2 += suffix;
-  // ------ input data ------
-  AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();
-
-  // ----- output data -----
-  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(nameContainer1, TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
-  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(nameContainer2, AliRDHFCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());
-
-  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
-  mgr->ConnectOutput(task,1,coutput1);
-  mgr->ConnectOutput(task,2,coutput2);
-
-
-  Printf("Input and Output connected to the manager");
-  return task ;
-}
+AliAnalysisTaskFlavourJetCorrelations *AddTaskFlavourJetCorrelations(\r
+  AliAnalysisTaskFlavourJetCorrelations::ECandidateType cand = AliAnalysisTaskFlavourJetCorrelations::kDstartoKpipi,\r
+  TString filename = "DStartoKpipiCuts.root",\r
+  Bool_t theMCon = kFALSE,\r
+  Bool_t reco = kTRUE /*must be true if theMCon is false*/,\r
+  TString jetArrname = "",\r
+  TString suffix = "",\r
+  Bool_t triggerOnLeadingJet = kFALSE,\r
+  Int_t leadingHadType = 0 /*0 = charged, 1 = neutral, 2 = both*/,\r
+  Float_t R = 0.4,\r
+  Float_t jptcut = 10.,\r
+  const char *cutType = "TPC",\r
+  Double_t percjetareacut = 1.)\r
+{\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error("AddTaskFlavourJetCorrelations::AddTaskFlavourJetCorrelations", "No analysis manager to connect to.");\r
+    return NULL;\r
+  } \r
+\r
+  Bool_t useStdC = kFALSE;\r
+  TFile* filecuts=TFile::Open(filename);\r
+  if (!filecuts || (filecuts && !filecuts->IsOpen())) {\r
+    cout<<"Input file not found: use std cuts"<<endl;\r
+    useStdC = kTRUE;\r
+  }\r
+\r
+  AliRDHFCuts *analysiscuts = 0x0;\r
+  switch (cand) {\r
+  case 0 :\r
+    if (useStdC) {\r
+      analysiscuts = new AliRDHFCutsD0toKpi();\r
+      analysiscuts->SetStandardCutsPP2010();\r
+    } else\r
+      analysiscuts = (AliRDHFCutsD0toKpi*)filecuts->Get("D0toKpiCuts");\r
+    break;\r
+  case 1 :\r
+    if(useStdC) {\r
+      analysiscuts = new AliRDHFCutsDStartoKpipi();\r
+      analysiscuts->SetStandardCutsPP2010();\r
+    } else\r
+      analysiscuts = (AliRDHFCutsDStartoKpipi*)filecuts->Get("DStartoKpipiCuts");\r
+    analysiscuts->SetName("DStartoKpipiCuts");\r
+    break;\r
+  }\r
+  \r
+  if (!analysiscuts) { // mm let's see if everything is ok\r
+    AliFatal("Specific AliRDHFCuts not found");\r
+    return;\r
+  }\r
+\r
+  printf("CREATE TASK\n"); //CREATE THE TASK\r
+\r
+  // create the task\r
+  AliAnalysisTaskFlavourJetCorrelations *task = new AliAnalysisTaskFlavourJetCorrelations("AnaTaskFlavourJetCorrelations", \r
+     analysiscuts, cand);\r
+  task->SetJetsName(jetArrname);\r
+  task->SetMC(theMCon);\r
+  task->SetUseReco(reco);\r
+  task->SetTriggerOnLeadingJet(triggerOnLeadingJet);\r
+\r
+  mgr->AddTask(task);\r
+\r
+  if(theMCon) {\r
+     suffix+="MC";\r
+     if(reco) suffix+="rec";  \r
+  }\r
+\r
+  // Create and connect containers for input/output\r
+  TString outputfile = AliAnalysisManager::GetCommonFileName();\r
+  outputfile += ":PWG3_D2H_DEmcalJet";\r
+  outputfile += suffix;\r
+\r
+  TString candname="DStar"; \r
+  if(cand==0)  candname="D0";\r
+\r
+  TString nameContainer1="hCor";\r
+  TString nameContainer2="cutsJ";\r
+\r
+  nameContainer1 += candname;\r
+  nameContainer2 += candname;\r
+  \r
+  nameContainer1 += suffix;\r
+  nameContainer2 += suffix;\r
+  // ------ input data ------\r
+  AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();\r
+\r
+  // ----- output data -----\r
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(nameContainer1, TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());\r
+  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(nameContainer2, AliRDHFCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());\r
+\r
+  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());\r
+  mgr->ConnectOutput(task,1,coutput1);\r
+  mgr->ConnectOutput(task,2,coutput2);\r
+\r
+\r
+  Printf("Input and Output connected to the manager");\r
+  return task ;\r
+}\r
index e1efec5..43e63b9 100644 (file)
@@ -1,87 +1,94 @@
-AliAnalysisTaskSEDmesonsFilterCJ *AddTaskSEDmesonsFilterCJ(
-  AliAnalysisTaskSEDmesonsFilterCJ::ECandidateType cand = AliAnalysisTaskSEDmesonsFilterCJ::kDstartoKpipi,
-  TString filename = "DStartoKpipiCuts.root",
-  Bool_t theMCon = kFALSE,
-  TString suffix = "")
-{
-  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
-  if (!mgr) {
-    ::Error("AddTaskSEDmesonsFilterCJ", "No analysis manager to connect to.");
-    return NULL;
-  } 
-
-  Bool_t useStdC = kFALSE;
-  TFile* filecuts=TFile::Open(filename);
-  if(!filecuts || (filecuts && !filecuts->IsOpen())) {
-    cout<<"Input file not found: use std cuts"<<endl;
-    useStdC = kTRUE;
-  }
-
-  AliRDHFCuts *analysiscuts=0x0;
-  switch (cand) {
-  case 0 :
-    if(useStdC) {
-      analysiscuts = new AliRDHFCutsD0toKpi();
-      analysiscuts->SetStandardCutsPP2010();
-    } else
-      analysiscuts = (AliRDHFCutsD0toKpi*)filecuts->Get("D0toKpiCuts");
-    break;
-  case 1 :
-    if(useStdC) {
-      analysiscuts = new AliRDHFCutsDStartoKpipi();
-      analysiscuts->SetStandardCutsPP2010();
-    } else
-      analysiscuts = (AliRDHFCutsDStartoKpipi*)filecuts->Get("DStartoKpipiCuts");
-    analysiscuts->SetName("DStartoKpipiCuts");
-    break;
-  }
-  
-  if (!analysiscuts) { // mm let's see if everything is ok
-    AliFatal("Specific AliRDHFCuts not found");
-    return;
-  } 
-
-  printf("CREATE TASK\n"); //CREATE THE TASK
-
-  // create the task
-  AliAnalysisTaskSEDmesonsFilterCJ *task = new AliAnalysisTaskSEDmesonsFilterCJ("AnaTaskSEDmesonsFilterCJ",analysiscuts,cand);
-  task->SetMC(theMCon); //D meson settings
-  mgr->AddTask(task);
-
-  TString candname="DStar"; 
-  if(cand==0)  candname="D0";
-  
-  // Create and connect containers for input/output
-  TString outputfile = AliAnalysisManager::GetCommonFileName();
-  outputfile += ":PWG3_D2H_DmesonsForJetCorrelations";
-  outputfile += suffix;
-
-  TString nameContainer0="histograms";
-  TString nameContainer1="cuts";
-  //TString nameContainer2="DcandidatesSel";
-
-  nameContainer0 += candname;
-  nameContainer1 += candname;
-  
-  nameContainer0 += suffix;
-  nameContainer1 += suffix;
-//nameContainer2 += suffix;
-
-  // ------ input data ------
-  AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();
-  
-  // ----- output data -----
-  
-  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(nameContainer0, TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
-  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(nameContainer1, AliRDHFCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());
-//AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(nameContainer2, TList::Class(),AliAnalysisManager::kExchangeContainer, outputfile.Data());
-  
-  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
-  mgr->ConnectOutput(task,1,coutput1);
-  mgr->ConnectOutput(task,2,coutput2);
-//mgr->ConnectOutput(task,3,coutput3);
-
-  Printf("Input and Output connected to the manager");
-  return task ;
-}
-
+AliAnalysisTaskSEDmesonsFilterCJ *AddTaskSEDmesonsFilterCJ(\r
+  AliAnalysisTaskSEDmesonsFilterCJ::ECandidateType cand = AliAnalysisTaskSEDmesonsFilterCJ::kDstartoKpipi,\r
+  TString filename = "DStartoKpipiCuts.root",\r
+  Bool_t theMCon = kFALSE,\r
+  Bool_t reco = kTRUE /*must be true if theMCon is false*/,\r
+  TString suffix = "")\r
+{\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error("AddTaskSEDmesonsFilterCJ", "No analysis manager to connect to.");\r
+    return NULL;\r
+  } \r
+\r
+  Bool_t useStdC = kFALSE;\r
+  TFile* filecuts=TFile::Open(filename);\r
+  if(!filecuts || (filecuts && !filecuts->IsOpen())) {\r
+    cout<<"Input file not found: use std cuts"<<endl;\r
+    useStdC = kTRUE;\r
+  }\r
+\r
+  AliRDHFCuts *analysiscuts=0x0;\r
+  switch (cand) {\r
+  case 0 :\r
+    if(useStdC) {\r
+      analysiscuts = new AliRDHFCutsD0toKpi();\r
+      analysiscuts->SetStandardCutsPP2010();\r
+    } else\r
+      analysiscuts = (AliRDHFCutsD0toKpi*)filecuts->Get("D0toKpiCuts");\r
+    break;\r
+  case 1 :\r
+    if(useStdC) {\r
+      analysiscuts = new AliRDHFCutsDStartoKpipi();\r
+      analysiscuts->SetStandardCutsPP2010();\r
+    } else\r
+      analysiscuts = (AliRDHFCutsDStartoKpipi*)filecuts->Get("DStartoKpipiCuts");\r
+    analysiscuts->SetName("DStartoKpipiCuts");\r
+    break;\r
+  }\r
+  \r
+  if (!analysiscuts) { // mm let's see if everything is ok\r
+    AliFatal("Specific AliRDHFCuts not found");\r
+    return;\r
+  } \r
+\r
+  printf("CREATE TASK\n"); //CREATE THE TASK\r
+\r
+  // create the task\r
+  AliAnalysisTaskSEDmesonsFilterCJ *task = new AliAnalysisTaskSEDmesonsFilterCJ("AnaTaskSEDmesonsFilterCJ",analysiscuts,cand);\r
+  if(!theMCon) reco=kTRUE;\r
+  task->SetMC(theMCon); //D meson settings\r
+  task->SetUseReco(reco);\r
+  mgr->AddTask(task);\r
+  if(theMCon) {\r
+     suffix+="MC";\r
+     if(reco) suffix+="rec";  \r
+  }\r
+  \r
+  TString candname="DStar"; \r
+  if(cand==0)  candname="D0";\r
+  \r
+  // Create and connect containers for input/output\r
+  TString outputfile = AliAnalysisManager::GetCommonFileName();\r
+  outputfile += ":PWG3_D2H_DmesonsForJetCorrelations";\r
+  outputfile += suffix;\r
+\r
+  TString nameContainer0="histograms";\r
+  TString nameContainer1="cuts";\r
+  //TString nameContainer2="DcandidatesSel";\r
+\r
+  nameContainer0 += candname;\r
+  nameContainer1 += candname;\r
+  \r
+  nameContainer0 += suffix;\r
+  nameContainer1 += suffix;\r
+//nameContainer2 += suffix;\r
+\r
+  // ------ input data ------\r
+  AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();\r
+  \r
+  // ----- output data -----\r
+  \r
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(nameContainer0, TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());\r
+  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(nameContainer1, AliRDHFCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());\r
+//AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(nameContainer2, TList::Class(),AliAnalysisManager::kExchangeContainer, outputfile.Data());\r
+  \r
+  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());\r
+  mgr->ConnectOutput(task,1,coutput1);\r
+  mgr->ConnectOutput(task,2,coutput2);\r
+//mgr->ConnectOutput(task,3,coutput3);\r
+\r
+  Printf("Input and Output connected to the manager");\r
+  return task ;\r
+}\r
+\r
index 5f84f74..c38cab7 100644 (file)
 void AddTasksFlavourJet(const Int_t iCandType = 1 /*0 = D0, 1=Dstar...*/,
-                       const TString sCutFile = "cutsHF/D0toKpiCutsppRecVtxNoPileupRejNoEMCAL.root",
-                        const Double_t dJetPtCut   = 1.,
-                        const Double_t dJetAreaCut = 0.,
-                        const Int_t iAccCut = 1,
-                        const TString sRunPeriod = "LHC10b",
-                        const Int_t    uBeamType = 0,
-                        const UInt_t uTriggerMask = AliVEvent::kMB, /*for jets; the D mesons trigger is defined in the cut object*/
-                        const Bool_t bIsMC = kFALSE,
-                        TString sText=""/*completes the name of the candidate task lists*/
-)
+   const TString sCutFile = "cutsHF/D0toKpiCutsppRecVtxNoPileupRejNoEMCAL.root",
+   const Double_t dJetPtCut   = 1.,
+   const Double_t dJetAreaCut = 0.,
+   const char *acctype = "TPC",
+   const TString sRunPeriod = "LHC10b",
+   const Int_t    uBeamType = 0,
+   const UInt_t uTriggerMask = AliVEvent::kMB, /*for jets; the D mesons trigger is defined in the cut object*/
+   const Bool_t bIsMC = kFALSE,
+   const Bool_t bIsReco = kFALSE,
+   const Bool_t bIsMap = kFALSE,
+   TString sText=""/*completes the name of the candidate task lists*/
+   )
 {
-  const TString sInputTrk  = "tracks";
-  const TString sUsedTrks  = "PicoTracks";
-  const TString sUsedClus  = "";
-
-  const Int_t iJetAlgo = 1;
-  const Int_t iJetType = 1;
-  /*
-  const Int_t    nRadius = 3;
-  const Double_t aRadius[] = {  0.2,   0.4,   0.6  };
-  const TString  sRadius[] = { "R02", "R04", "R06" };
-  */
-  const Int_t    nRadius = 1;
-  const Double_t aRadius[] = {  0.4  };
-  const TString  sRadius[] = { "R04" };
-
-//=============================================================================
-
-  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
-
-  if (!mgr) {
-    ::Error("AddTasksFlavourJet.C::AddTasksFlavourJet", "No analysis manager to connect to.");
-    return;
-  }
-
-  TString type = mgr->GetInputEventHandler()->GetDataType();
-  if (!type.Contains("ESD") && !type.Contains("AOD")) {
-    ::Error("AddTasksFlavourJet.C::AddTasksFlavourJet", "Task manager to have an ESD or AOD input handler.");
-    return;
-  }
-
-  if (!mgr->GetInputEventHandler()) {
-    ::Error("AddTasksFlavourJet.C::AddTasksFlavourJet", "This task requires an input event handler");
-    return;
-  }
-//=============================================================================
-
-  UInt_t uAnaType = (((iJetType==0) ||     (iJetType==2)) ? 1 : 0);
-  Int_t  iLeading =  ((iJetType==0) ? 3 : ((iJetType==1)  ? 0 : 1));
-
-  //D mesons -- PID
-  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
-  AliAnalysisTaskSE *taskRespPID = AddTaskPIDResponse(bIsMC);
-
-  // -- D meson selection
-  gROOT->LoadMacro("AliAnalysisTaskSEDmesonsFilterCJ.cxx++g");
-  gROOT->LoadMacro("AddTaskSEDmesonsFilterCJ.C");
-  AliAnalysisTaskSEDmesonsFilterCJ *taskDmesonsFilter = AddTaskSEDmesonsFilterCJ(iCandType,sCutFile,bIsMC,sText);
-
-  // EMCal framework
-  // -- Physics selection task
-  gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalPhysicsSelection.C");
-  AliPhysicsSelectionTask *physSelTask = AddTaskEmcalPhysicsSelection(kTRUE, kTRUE, uTriggerMask, 5, 5, 10, kTRUE, -1, -1, -1, -1);
-
-  if (!physSelTask) {
-    cout << "no physSelTask"; 
-    return; 
-  }
-
-  // -- 
-  gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalSetup.C");
-  AliEmcalSetupTask *taskSetupEMCal = AddTaskEmcalSetup();
-  taskSetupEMCal->SetGeoPath("$ALICE_ROOT/OADB/EMCAL");
-  taskSetupEMCal->SelectCollisionCandidates(uTriggerMask);
-
-  // Jet preparation
-  gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskJetPreparation.C");
-  AddTaskJetPreparation(type,bIsMC,sRunPeriod);
-
-  gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalPicoTrackMaker.C");
-  AliEmcalPicoTrackMaker *taskPicoTrack = AddTaskEmcalPicoTrackMaker(sUsedTrks.Data(),sInputTrk.Data(),sRunPeriod.Data());
-  taskPicoTrack->SelectCollisionCandidates(uTriggerMask);
-
-  gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJet.C");
-  gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJetSample.C");
-  gROOT->LoadMacro("AliAnalysisTaskFlavourJetCorrelations.cxx++g");
-  gROOT->LoadMacro("AddTaskFlavourJetCorrelations.C");
-
-  for (Int_t i=0; i<nRadius; i++) {
-    AliEmcalJetTask *taskFJ = AddTaskEmcalJet(sUsedTrks.Data(),sUsedClus.Data(),iJetAlgo,aRadius[i],iJetType);
-    taskFJ->SelectCollisionCandidates(uTriggerMask);
-
-    AliAnalysisTaskFlavourJetCorrelations *taskDmesonCJ = AddTaskFlavourJetCorrelations(iCandType,
-                                                                                        sCutFile,
-                                                                                        bIsMC,
-                                                                                        taskFJ->GetName(),
-                                                                                        Form("JetR%s",sRadius[i].Data()),
-                                                                                        iLeading,
-                                                                                        aRadius[i],
-                                                                                        dJetPtCut,
-                                                                                        iAccCut,
-                                                                                        dJetAreaCut);
-
-    taskDmesonCJ->SetName(Form("AliAnalysisTaskSEEmcalJetDmesonsCJ_%s",sRadius[i].Data()));
-   taskDmesonCJ->SetForceBeamType(uBeamType);
-    taskDmesonCJ->SetAnaType(uAnaType);
-    taskDmesonCJ->SetLeadingHadronType(iLeading);
-
-//  taskDmesonCJ->SelectCollisionCandidates(uTriggerMask);
-  }
-
-  return;
+   const TString sInputTrkMC  = "MCParticlesSelected";
+   const TString sInputTrkRec  = "tracks";
+   const TString sUsedTrks  = "PicoTracks";
+   const TString sUsedClus  = "";
+   TString sInputTrk = bIsReco ? sInputTrkRec : sInputTrkMC;
+   const Int_t iJetAlgo = 1;
+   const Int_t iJetType = 1;
+   /*
+   const Int_t    nRadius = 3;
+   const Double_t aRadius[] = {  0.2,   0.4,   0.6  };
+   const TString  sRadius[] = { "R02", "R04", "R06" };
+   */
+   const Int_t    nRadius = 1;
+   const Double_t aRadius[] = {  0.4  };
+   const TString  sRadius[] = { "R04" };
+   
+   //=============================================================================
+   
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   
+   if (!mgr) {
+      ::Error("AddTasksFlavourJet.C::AddTasksFlavourJet", "No analysis manager to connect to.");
+      return;
+   }
+   
+   TString type = mgr->GetInputEventHandler()->GetDataType();
+   if (!type.Contains("ESD") && !type.Contains("AOD")) {
+      ::Error("AddTasksFlavourJet.C::AddTasksFlavourJet", "Task manager to have an ESD or AOD input handler.");
+      return;
+   }
+   
+   if (!mgr->GetInputEventHandler()) {
+      ::Error("AddTasksFlavourJet.C::AddTasksFlavourJet", "This task requires an input event handler");
+      return;
+   }
+   //=============================================================================
+   
+   UInt_t uAnaType = (((iJetType==0) ||     (iJetType==2)) ? 1 : 0);
+   Int_t  iLeading =  ((iJetType==0) ? 3 : ((iJetType==1)  ? 0 : 1));
+   Int_t leadHadType=0; /* 0=charged, 1=neutral, 2=both*/
+   //D mesons -- PID
+   gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
+   AliAnalysisTaskSE *taskRespPID = AddTaskPIDResponse(bIsMC);
+   
+   // -- D meson selection
+  
+   gROOT->LoadMacro("$ALICE_ROOT/PWGJE/FlavourJetTasks/macros/AddTaskSEDmesonsFilterCJ.C");
+   AliAnalysisTaskSEDmesonsFilterCJ *taskDmesonsFilter = AddTaskSEDmesonsFilterCJ(iCandType,sCutFile,bIsMC,bIsReco,sText);
+   if(bIsMap) {
+      AliAnalysisTaskSEDmesonsFilterCJ *taskMCDmesonsFilter = AddTaskSEDmesonsFilterCJ(iCandType,sCutFile,bIsMC,kFALSE,sText);
+   }
+   // EMCal framework
+   // -- Physics selection task
+   gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalPhysicsSelection.C");
+   AliPhysicsSelectionTask *physSelTask = AddTaskEmcalPhysicsSelection(kTRUE, kTRUE, uTriggerMask, 5, 5, 10, kTRUE, -1, -1, -1, -1);
+   
+   if (!physSelTask) {
+      cout << "no physSelTask"; 
+      return; 
+   }
+   
+   // -- 
+   gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalSetup.C");
+   AliEmcalSetupTask *taskSetupEMCal = AddTaskEmcalSetup();
+   taskSetupEMCal->SetGeoPath("$ALICE_ROOT/OADB/EMCAL");
+   taskSetupEMCal->SelectCollisionCandidates(uTriggerMask);
+   
+   // Jet preparation
+   //gROOT->LoadMacro("/data/Work/jets/testEMCalJetFramework/code/v4/AddTaskJetPreparation.C");
+   gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskJetPreparation.C");
+   AddTaskJetPreparation(type,sRunPeriod,"PicoTracks",bIsMC ? "MCParticlesSelected" : "",/*next 7 emcal default settings*/"","",2.,0.,0.03,0.015,0.15,uTriggerMask, kFALSE /*track cluster*/,kFALSE /*do histos*/,kTRUE /*make pico tracks*/,kFALSE /*make emcal trigger*/,kFALSE /*is emcal train*/);
+   
+   gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalPicoTrackMaker.C");
+   AliEmcalPicoTrackMaker *taskPicoTrack = AddTaskEmcalPicoTrackMaker(sUsedTrks.Data(),sInputTrk.Data(),sRunPeriod.Data());
+   taskPicoTrack->SelectCollisionCandidates(uTriggerMask);
+   
+   gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJet.C");
+   //gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJetSample.C");
+   
+   gROOT->LoadMacro("$ALICE_ROOT/PWGJE/FlavourJetTasks/macros/AddTaskFlavourJetCorrelations.C");
+   gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskJetResponseMaker.C");
+   
+   for (Int_t i=0; i<nRadius; i++) {
+      //jet reconstruction for correlation
+      //AliEmcalJetTask *taskFJ = AddTaskEmcalJet(bIsMC&&!bIsReco ? sInputTrkMC.Data() : sUsedTrks.Data(),sUsedClus.Data(),iJetAlgo,aRadius[i],iJetType);
+      //I think like follows it's fine: sUsedTrks contains the picotracks of MC or reco tracks according to IsReco
+      AliEmcalJetTask *taskFJ = AddTaskEmcalJet(sUsedTrks.Data(),sUsedClus.Data(),iJetAlgo,aRadius[i],iJetType);
+      
+      taskFJ->SelectCollisionCandidates(uTriggerMask);
+      
+            
+     // AliAnalysisTaskEmcalJetSample *taskjetsample=AddTaskEmcalJetSample(sUsedTrks.Data(),"",taskFJ->GetName(),
+       // "Rho",aRadius[i],dJetPtCut,dJetAreaCut,AliAnalysisTaskEmcalJet::kTPC,iLeading,
+       // "AliAnalysisTaskEmcalJetSample");
+      //taskjetsample->SelectCollisionCandidates(uTriggerMask);
+      
+      //correlation with D meson
+      
+      AliAnalysisTaskFlavourJetCorrelations *taskDmesonCJ = AddTaskFlavourJetCorrelations(
+        iCandType,
+        sCutFile,
+        bIsMC,
+        bIsReco,
+        taskFJ->GetName(),
+        Form("JetR%s",sRadius[i].Data()),
+        iLeading,
+        leadHadType,
+        aRadius[i],
+        dJetPtCut,
+        acctype
+       /*percjetareacut=1.*/);
+      
+      taskDmesonCJ->SetName(Form("AliAnalysisTaskSEEmcalJetDmesonsCJ_%s",sRadius[i].Data()));
+      taskDmesonCJ->SetForceBeamType(uBeamType);
+      taskDmesonCJ->SetAnaType(uAnaType);
+      taskDmesonCJ->SetLeadingHadronType(iLeading);
+      
+      AliEmcalJetTask *taskMCJ;
+      //jet reconstruction for correction map
+      if(bIsMap){ 
+        taskMCJ = AddTaskEmcalJet(sInputTrkMC.Data(),sUsedClus.Data(),iJetAlgo,aRadius[i],
+           iJetType);
+        
+        AliAnalysisTaskFlavourJetCorrelations *taskMCDmesonCJ = AddTaskFlavourJetCorrelations(
+           iCandType,
+           sCutFile,
+           kTRUE,
+           kFALSE,
+           taskMCJ->GetName(),
+           Form("JetR%s",sRadius[i].Data()),
+           iLeading,
+           leadHadType,
+           aRadius[i],
+           dJetPtCut,
+           acctype
+           /*percjetareacut=1.*/);
+        
+        taskMCDmesonCJ->SetName(Form("AliAnalysisTaskSEEmcalJetMCDmesonsCJ_%s",sRadius[i].Data()));
+        taskMCDmesonCJ->SetForceBeamType(uBeamType);
+        taskMCDmesonCJ->SetAnaType(uAnaType);
+        taskMCDmesonCJ->SetLeadingHadronType(iLeading);
+      }
+      //  taskDmesonCJ->SelectCollisionCandidates(uTriggerMask);
+      
+      // definition of correction map
+      Int_t tag=0;
+      if(iCandType == 0) tag=AliEmcalJet::kD0;
+      if(iCandType == 1) tag=AliEmcalJet::kDStar;
+      Printf("************** tag = %d", tag&0x1);
+      AliJetResponseMaker* taskResp=AddTaskJetResponseMaker(
+        sUsedTrks.Data(),sUsedClus.Data(),taskFJ->GetName(),"",aRadius[i],
+        sInputTrkMC.Data(),"",taskMCJ->GetName(),"",aRadius[i],dJetPtCut,dJetAreaCut,5,0,AliJetResponseMaker::kGeometrical, 0.25,0.25,"TPC",-999,-999,-999,"AliJetResponseMaker", kFALSE, 0, -10,10, tag );
+      taskResp->SetMinJetMCPt(0); //added to bypass a return not needed (feature of PrepareJetTask)
+      //taskResp->SetHistoType(1);
+   }
+   
+   return;
 }
index 944d34f..2090575 100644 (file)
@@ -30,27 +30,28 @@ AliAnalysisGrid* CreateAlienHandler(const char* uniqueName, const char* gridDir,
 void AnalysisTrainCorrJetsLocal (
          const char*    dataType            = "AOD",                       // set the analysis type, AOD, ESD or sESD
          Bool_t         useGrid             = kTRUE,                      // local or grid
-        TString localfilename = "file_aodlhc10d.txt",
+        TString localfilename = "/data/Work/jets/testEMCalJetFramework/ptJdistrAug2nd/setMCtasks/map/inheritDev/files_LHC10f7a.txt",
          const char*    gridMode            = "test",                      // set the grid run mode (can be "full", "test", "offline", "submit" or "terminate")
-         const char*    pattern             = "*ESDs/pass2/AOD135/*AliAOD.root", //"*/*/AliAOD.root" "*ESDs/pass1/AOD106/*AliAOD.root",    // file pattern (here one can specify subdirs like passX etc.) (used on grid)
-         const char*    gridDir             = "/alice/data/2010/LHC10d",   // /alice/data/2011/LHC11d /alice/sim/2012/LHC12f2b   dir on alien, where the files live (used on grid)
-         const char*    runNumbers          = /*126437*/" 126432 126425 126424 126422 126409 126408 126407 126406 126405 126404 126403 126359 126352 126351 126350 126285 126284 126283 126168 126167 126160 126158 126097 126090 126088 126082 126081 126078 126073 126008 126007 126004 125855 125851 125850 125849 125848 125847 125844 125843 125842 125633 125632 125630 125296 125134 125101 125100 125097 125085 125023 124751 122375 122374",             // considered run numbers (used on grid) /*LHC12g 188359 188362, LHC11a 146860 146859*/ /*LHC12f2b 158285 159582 */
+         const char*    pattern             = "*AOD136/*AliAOD.root", //"*/*/AliAOD.root" "*ESDs/pass1/AOD106/*AliAOD.root",    // file pattern (here one can specify subdirs like passX etc.) (used on grid)
+         const char*    gridDir             = "/alice/sim/LHC10f7a",   // /alice/data/2011/LHC11d /alice/sim/2012/LHC12f2b   dir on alien, where the files live (used on grid)
+         const char*    runNumbers          = /*114931 */" 115186 115193 115393 115401 115414 116102 116288 116402 116403 116562 116571 116574 116643 116645 117048 117050 117052 117053 117054 117059 117060 117063 117065 117077 117086 117092 117099 117109 117112 117116 117220 117222 119159 119161 119163 119841 119842 119844 119845 119846 119849 119853 119856 119859 119862 120067 120069 120072 120073 120076 120079 120244 120503 120504 120505 120616 120617 120671 120741 120750 120758 120820 120821 120822 120823 120824 120825 120829 122374 122375 124187 124191 124355 124358 124362 124367 124378 124380 124381 124604 124605 124606 124607 124608 124702 124746 124750 124751 125023 125085 125097 125100 125101 125133 125134 125139 125140 125156 125186 125296 125628 125630 125632 125633 125842 125843 125844 125847 125848 125849 125850 125851 125855 126004 126007 126008 126073 126078 126081 126082 126088 126090 126097 126158 126160 126167 126168 126283 126284 126285 126350 126351 126352 126359 126403 126404 126405 126406 126407 126408 126409 126422 126424 126425 126432 126437 127719 127724 127729 127730 127814 127819 127930 127940 128263 128778 128913 129536 129599 129639 129641 129654 129659 129666 129667 129723 129725 129726 129729 129735 129736 129738 129742 129744 129959 129960 129961 129962 129966 129983 130149 130151 130157 130158 130172 130178 130179 130342 130343 130354 130356 130358 130360 130375 130479 130480 130481 130517 130519 130520 130524 130526 130601 130608 130620 130621 130623 130628 130696 130704 130793 130795 130798 130799 130834 130840 130842 130844 130847 130848",             // considered run numbers (used on grid) /*LHC12g 188359 188362, LHC11a 146860 146859*/ /*LHC12f2b 158285 159582 */
         const Int_t nrunspermaster= 100,
-         UInt_t         numLocalFiles       = 3,                          // number of files analyzed locally  
-         const char*    runPeriod           = "LHC10d",                    // set the run period (used on grid)
-         const char*    uniqueName          = "DJetNewCode",     // sets base string for the name of the task on the grid
+         UInt_t         numLocalFiles       = 10,                          // number of files analyzed locally  
+         const char*    runPeriod           = "LHC10f7a",                    // set the run period (used on grid)
+         const char*    uniqueName          = "DJetNewCodeMCSandBchJ",     // sets base string for the name of the task on the grid
          UInt_t         pSel                = AliVEvent::kAny,             // used event selection for every task except for the analysis tasks
          Bool_t         useTender           = kFALSE,                      // trigger, if tender task should be used
-         Bool_t         isMC                = kFALSE,                      // trigger, if MC handler should be used
-
+         Bool_t         isMC                = kTRUE,                      // trigger, if MC handler should be used
+        Bool_t         isReco                = kTRUE,
+        Bool_t         isMap                = kTRUE,
          // Here you have to specify additional code files you want to use but that are not in aliroot
-         const char*    addCXXs             = "AliAnalysisTaskSEDmesonsFilterCJ.cxx AliAnalysisTaskFlavourJetCorrelations.cxx",
-         const char*    addHs               = "AliAnalysisTaskSEDmesonsFilterCJ.h AliAnalysisTaskFlavourJetCorrelations.h",
+         const char*    addCXXs             = "",
+         const char*    addHs               = "",
 
          // These two settings depend on the dataset and your quotas on the AliEN services
          Int_t          maxFilesPerWorker   = 20,
          Int_t          workerTTL           = 86000,
-        Int_t          nfiletestmode       = 1
+        Int_t          nfiletestmode       = 50
          )
 {
 
@@ -111,21 +112,24 @@ void AnalysisTrainCorrJetsLocal (
   if(!useGrid)
     cout << "Using " << localFiles.Data() << " as input file list.\n";
 
-  gROOT->LoadMacro("AddTasksFlavourJet.C");
+  gROOT->LoadMacro("$ALICE_ROOT/PWGJE/FlavourJetTasks/macros/AddTasksFlavourJet.C");
   //List of arguments:
 
-// const Int_t iCandtype=1 /*0 = D0, 1=Dstar...*/,
-// const TString sCutFile = "cutsHF/D0toKpiCutsppRecVtxNoPileupRejNoEMCAL.root",
-// const Double_t dJetPtCut   = 1.,
-// const Double_t dJetAreaCut = 0.,
-// const Int_t iAccCut = 1,
-// const TString sRunPeriod = "LHC10b",
-// const Int_t    uBeamType   = 0,
-// const UInt_t uTriggerMask = AliVEvent::kMB, /*for jets and phys sel; the D mesons trigger is defined in the cut object*/
-// const Bool_t bIsMC=kFALSE,
-// TString sText=""/*completes the name of the candidate task lists*/
-
-  AddTasksFlavourJet(1,"/data/Work/jets/testEMCalJetFramework/CutFilesMB/DStartoKpipiCuts_new.root",10.,0.,1,runPeriod,0,pSel,isMC,"");
+  // const Int_t iCandType = 1 /*0 = D0, 1=Dstar...*/,
+  //  const TString sCutFile = "cutsHF/D0toKpiCutsppRecVtxNoPileupRejNoEMCAL.root",
+  //  const Double_t dJetPtCut   = 1.,
+  //  const Double_t dJetAreaCut = 0.,
+  //  const char *acctype = "TPC",
+  //  const TString sRunPeriod = "LHC10b",
+  //  const Int_t    uBeamType = 0,
+  //  const UInt_t uTriggerMask = AliVEvent::kMB, /*for jets; the D mesons trigger is defined in the cut object*/
+  //  const Bool_t bIsMC = kFALSE,
+  //  const Bool_t bIsReco = kFALSE,
+  //  const Bool_t bIsMap = kFALSE,
+  //  TString sText=""/*completes the name of the candidate task lists*/
+
+  //AddTasksFlavourJet(1,"/data/Work/jets/testEMCalJetFramework/CutFilesMB/DStartoKpipiCuts_new.root",10.,0.,"TPC",runPeriod,0,pSel,isMC,isReco,isMap,"");
+  AddTasksFlavourJet(1,"/data/Work/jets/testEMCalJetFramework/CutFilesMB/DStartoKpipiCuts.root",5.,0.,"TPC",runPeriod,0,pSel,isMC,isReco,isMap,"");
 
   // Set the physics selection for all given tasks
   TObjArray *toptasks = mgr->GetTasks();
@@ -221,6 +225,7 @@ void LoadLibs()
   gSystem->Load("libSTAT.so");
   gSystem->Load("libTRDrec.so");
   gSystem->Load("libHMPIDbase.so");
+  gSystem->Load("libPWGTools.so");
   gSystem->Load("libPWGPP.so");
   gSystem->Load("libPWGHFbase");
   gSystem->Load("libPWGDQdielectron");
@@ -253,11 +258,13 @@ void LoadLibs()
   gSystem->Load("libJETAN");
   gSystem->Load("libFASTJETAN");
   gSystem->Load("libPWGJEEMCALJetTasks");
+  gSystem->Load("libPWGJEFlavourJetTasks");
 
 
   // include paths
   gSystem->AddIncludePath("-Wno-deprecated");
   gSystem->AddIncludePath("-I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/EMCAL -I$ALICE_ROOT/PWG/EMCAL -I$ALICE_ROOT/PWGJE/EMCALJetTasks");
+  gSystem->AddIncludePath(" -I$ALICE_ROOT/PWGHF/ -I$ALICE_ROOT/PWGHF/base -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWGJE/FlavourJetTasks");
   gSystem->AddIncludePath("-I$ALICE_ROOT/PWGDQ/dielectron -I$ALICE_ROOT/PWGHF/hfe");
   gSystem->AddIncludePath("-I$ALICE_ROOT/JETAN -I$ALICE_ROOT/JETAN/fastjet");
 }
@@ -279,7 +286,7 @@ AliAnalysisGrid* CreateAlienHandler(const char* uniqueName, const char* gridDir,
   }else tmpName +="_20130412_122423";
   */
   TString tmpAdditionalLibs("");
-  tmpAdditionalLibs = Form("libTree.so libVMC.so libGeom.so libGui.so libXMLParser.so libMinuit.so libMinuit2.so libProof.so libPhysics.so libSTEERBase.so libESD.so libAOD.so libOADB.so libANALYSIS.so libCDB.so libRAWDatabase.so libSTEER.so libANALYSISalice.so libCORRFW.so libTOFbase.so libRAWDatabase.so libRAWDatarec.so libTPCbase.so libTPCrec.so libITSbase.so libITSrec.so libTRDbase.so libTENDER.so libSTAT.so libTRDrec.so libHMPIDbase.so libPWGPP.so libPWGHFbase.so libPWGDQdielectron.so libPWGHFhfe.so libPWGflowBase.so libPWGflowTasks.so libPWGHFvertexingHF.so libEMCALUtils.so libPHOSUtils.so libPWGCaloTrackCorrBase.so libEMCALraw.so libEMCALbase.so libEMCALrec.so libTRDbase.so libVZERObase.so libVZEROrec.so libTENDER.so libTENDERSupplies.so libPWGEMCAL.so libPWGGAEMCALTasks.so libPWGTools.so libPWGCFCorrelationsBase.so libPWGCFCorrelationsDPhi.so  libCGAL.so libJETAN.so libfastjet.so libSISConePlugin.so libFASTJETAN.so libPWGJE.so libPWGmuon.so libPWGJEEMCALJetTasks.so %s %s",additionalCode.Data(),additionalHeaders.Data());
+  tmpAdditionalLibs = Form("libTree.so libVMC.so libGeom.so libGui.so libXMLParser.so libMinuit.so libMinuit2.so libProof.so libPhysics.so libSTEERBase.so libESD.so libAOD.so libOADB.so libANALYSIS.so libCDB.so libRAWDatabase.so libSTEER.so libANALYSISalice.so libCORRFW.so libTOFbase.so libRAWDatabase.so libRAWDatarec.so libTPCbase.so libTPCrec.so libITSbase.so libITSrec.so libTRDbase.so libTENDER.so libSTAT.so libTRDrec.so libHMPIDbase.so libPWGTools.so libPWGPP.so libPWGHFbase.so libPWGDQdielectron.so libPWGHFhfe.so libPWGflowBase.so libPWGflowTasks.so libPWGHFvertexingHF.so libEMCALUtils.so libPHOSUtils.so libPWGCaloTrackCorrBase.so libEMCALraw.so libEMCALbase.so libEMCALrec.so libTRDbase.so libVZERObase.so libVZEROrec.so libTENDER.so libTENDERSupplies.so libPWGEMCAL.so libPWGGAEMCALTasks.so libPWGTools.so libPWGCFCorrelationsBase.so libPWGCFCorrelationsDPhi.so  libCGAL.so libJETAN.so libfastjet.so libSISConePlugin.so libFASTJETAN.so libPWGJE.so libPWGmuon.so libPWGJEEMCALJetTasks.so libPWGJEFlavourJetTasks.so %s %s",additionalCode.Data(),additionalHeaders.Data());
 
 
   TString macroName("");
@@ -295,8 +302,8 @@ AliAnalysisGrid* CreateAlienHandler(const char* uniqueName, const char* gridDir,
 
   // Here you can set the (Ali)ROOT version you want to use
   plugin->SetAPIVersion("V1.1x");
-  plugin->SetROOTVersion("v5-34-05");
-  plugin->SetAliROOTVersion("v5-04-38-AN");
+  plugin->SetROOTVersion("v5-34-08");
+  plugin->SetAliROOTVersion("v5-05-11-AN");
 
   plugin->SetGridDataDir(gridDir); // e.g. "/alice/sim/LHC10a6"
   plugin->SetDataPattern(pattern); //dir structure in run directory
@@ -312,7 +319,7 @@ AliAnalysisGrid* CreateAlienHandler(const char* uniqueName, const char* gridDir,
 
   plugin->SetAnalysisSource(additionalCode.Data());
   plugin->SetAdditionalLibs(tmpAdditionalLibs.Data());
-  plugin->AddIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/AOD -I$ALICE_ROOT/ANALYSIS  -I$ALICE_ROOT/OADB -I$ALICE_ROOT/PWGHF -I$ALICE_ROOT/PWGHF/base -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWG/FLOW/Base -I$ALICE_ROOT/PWG/FLOW/Tasks  -I$ALICE_ROOT/PWGJE  -I$ALICE_ROOT/JETAN -I$ALICE_ROOT/PWGJE/EMCALJetTasks -g");
+  plugin->AddIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/AOD -I$ALICE_ROOT/ANALYSIS  -I$ALICE_ROOT/OADB -I$ALICE_ROOT/PWGHF -I$ALICE_ROOT/PWGHF/base -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWG/FLOW/Base -I$ALICE_ROOT/PWG/FLOW/Tasks  -I$ALICE_ROOT/PWGJE  -I$ALICE_ROOT/JETAN -I$ALICE_ROOT/PWGJE/EMCALJetTasks -I$ALICE_ROOT/PWGJE/FlavourJetTasks -g");
   plugin->AddExternalPackage("boost::v1_43_0");
   plugin->AddExternalPackage("cgal::v3.6");
   plugin->AddExternalPackage("fastjet::v2.4.2");