+++ /dev/null
-/**************************************************************************
- * 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. *
- **************************************************************************/
-//
-//
-// Base class for DStar in Jets Analysis
-//
-//-----------------------------------------------------------------------
-// Author A.Grelli
-// ERC-QGP Utrecht University - a.grelli@uu.nl
-//-----------------------------------------------------------------------
-
-#include <TDatabasePDG.h>
-#include <TParticle.h>
-#include <TVector3.h>
-#include "TROOT.h"
-
-#include "AliAnalysisTaskSEDStarJets.h"
-#include "AliRDHFCutsDStartoKpipi.h"
-#include "AliStack.h"
-#include "AliMCEvent.h"
-#include "AliAODMCHeader.h"
-#include "AliAODHandler.h"
-#include "AliAnalysisManager.h"
-#include "AliLog.h"
-#include "AliAODVertex.h"
-#include "AliAODJet.h"
-#include "AliAODRecoDecay.h"
-#include "AliAODRecoCascadeHF.h"
-#include "AliAODRecoDecayHF2Prong.h"
-#include "AliESDtrack.h"
-#include "AliAODMCParticle.h"
-
-ClassImp(AliAnalysisTaskSEDStarJets)
-
-//__________________________________________________________________________
-AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets() :
- AliAnalysisTaskSE(),
- fEvents(0),
- fchargeFrCorr(0),
- fUseMCInfo(kTRUE),
- fRequireNormalization(kTRUE),
- fOutput(0),
- fCuts(0),
- ftrigger(0),
- fPtPion(0),
- fInvMass(0),
- fRECOPtDStar(0),
- fRECOPtBkg(0),
- fDStar(0),
- fDiff(0),
- fDiffSideBand(0),
- fDStarMass(0),
- fPhi(0),
- fPhiBkg(0),
- fTrueDiff(0),
- fResZ(0),
- fResZBkg(0),
- fEjet(0),
- fPhijet(0),
- fEtaJet(0),
- theMCFF(0),
- fDphiD0Dstar(0),
- fPtJet(0)
-{
- //
- // Default ctor
- //
-}
-//___________________________________________________________________________
-AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
- AliAnalysisTaskSE(name),
- fEvents(0),
- fchargeFrCorr(0),
- fUseMCInfo(kTRUE),
- fRequireNormalization(kTRUE),
- fOutput(0),
- fCuts(0),
- ftrigger(0),
- fPtPion(0),
- fInvMass(0),
- fRECOPtDStar(0),
- fRECOPtBkg(0),
- fDStar(0),
- fDiff(0),
- fDiffSideBand(0),
- fDStarMass(0),
- fPhi(0),
- fPhiBkg(0),
- fTrueDiff(0),
- fResZ(0),
- fResZBkg(0),
- fEjet(0),
- fPhijet(0),
- fEtaJet(0),
- theMCFF(0),
- fDphiD0Dstar(0),
- fPtJet(0)
-{
- //
- // Constructor. Initialization of Inputs and Outputs
- //
- fCuts=cuts;
- Info("AliAnalysisTaskSEDStarJets","Calling Constructor");
-
- DefineOutput(1,TList::Class()); // histos
- DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
-}
-//___________________________________________________________________________
-AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() {
- //
- // destructor
- //
-
- Info("~AliAnalysisTaskSEDStarJets","Calling Destructor");
- if (fOutput) {
- delete fOutput;
- fOutput = 0;
- }
-
- if (fCuts) {
- delete fCuts;
- fCuts = 0;
- }
-
- if (ftrigger) { delete ftrigger; ftrigger = 0;}
- if (fPtPion) { delete fPtPion; fPtPion = 0;}
- if (fInvMass) { delete fInvMass; fInvMass = 0;}
- if (fRECOPtDStar) { delete fRECOPtDStar; fRECOPtDStar = 0;}
- if (fRECOPtBkg) { delete fRECOPtBkg; fRECOPtBkg = 0;}
- if (fDStar) { delete fDStar; fDStar = 0;}
- if (fDiff) { delete fDiff; fDiff = 0;}
- if (fDiffSideBand) { delete fDiffSideBand; fDiffSideBand = 0;}
- if (fDStarMass) { delete fDStarMass; fDStarMass = 0;}
- if (fPhi) { delete fPhi; fPhi = 0;}
- if (fPhiBkg) { delete fPhiBkg; fPhiBkg = 0;}
- if (fTrueDiff){ delete fTrueDiff; fTrueDiff = 0;}
- if (fResZ) { delete fResZ; fResZ = 0;}
- if (fResZBkg) { delete fResZBkg; fResZBkg = 0;}
- if (fEjet) { delete fEjet; fEjet = 0;}
- if (fPhijet) { delete fPhijet; fPhijet = 0;}
- if (fEtaJet) { delete fEtaJet; fEtaJet = 0;}
- if (theMCFF) { delete theMCFF; theMCFF = 0;}
- if (fDphiD0Dstar) { delete fDphiD0Dstar; fDphiD0Dstar = 0;}
- if (fPtJet) { delete fPtJet; fPtJet = 0;}
-}
-
-//___________________________________________________________
-void AliAnalysisTaskSEDStarJets::Init(){
- //
- // Initialization
- //
- if(fDebug > 1) printf("AnalysisTaskSEDStarJets::Init() \n");
- AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
- // Post the cuts
- PostData(2,copyfCuts);
-
- return;
-}
-
-//_________________________________________________
-void AliAnalysisTaskSEDStarJets::UserExec(Option_t *)
-{
- // user exec
- if (!fInputEvent) {
- Error("UserExec","NO EVENT FOUND!");
- return;
- }
-
- fEvents++;
- AliInfo(Form("Event %d",fEvents));
- if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
-
- // 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("Dstar");
- }
- } else {
- arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
- }
-
- if (!arrayDStartoD0pi){
- AliInfo("Could not find array of HF vertices, skipping the event");
- return;
- }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
-
- // 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;
-
- // Simulate a jet triggered sample
- TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets();
- if(aodEvent->GetNJets()<=0) return;
-
- // counters for efficiencies
- Int_t icountReco = 0;
-
- // Normalization factor
- if(fRequireNormalization){
- ftrigger->Fill(1);
- }
-
- //D* and D0 prongs needed to MatchToMC method
- Int_t pdgDgDStartoD0pi[2]={421,211};
- Int_t pdgDgD0toKpi[2]={321,211};
-
- Double_t max =0;
- Double_t ejet = 0;
- Double_t phiJet = 0;
- Double_t etaJet = 0;
- Double_t ptjet = 0;
-
- //loop over jets and consider only the leading jey in the event
- for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) {
- AliAODJet* jet = (AliAODJet*)arrayofJets->At(iJets);
-
- //jets variables
- ejet = jet->E();
-
- if(ejet>max){
- max = jet->E();
- phiJet = jet->Phi();
- etaJet = jet->Eta();
- ptjet = jet->Pt();
- }
-
- // fill energy, eta and phi of the jet
- fEjet ->Fill(ejet);
- fPhijet ->Fill(phiJet);
- fEtaJet ->Fill(etaJet);
- fPtJet->Fill(ptjet);
- }
-
- //loop over D* candidates
- for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
-
- // D* candidates
- AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
- AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
-
- Double_t finvM =-999;
- Double_t finvMDStar =-999;
- Double_t dPhi =-999;
- Bool_t isDStar =0;
- Int_t mcLabel = -9999;
-
- // find associated MC particle for D* ->D0toKpi
- if(fUseMCInfo){
- TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
- if(!mcArray) {
- printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
- return;
- }
- mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
-
- if(mcLabel>=0) isDStar = 1;
- if(mcLabel>0){
- Double_t zMC =-999;
- AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
- //fragmentation function in mc
- zMC= FillMCFF(mcPart,mcArray,mcLabel);
- if(zMC>0) theMCFF->Fill(zMC);
- }
- }
-
- // soft pion
- AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
- Double_t pt = dstarD0pi->Pt();
-
- // track quality cuts
- Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
- if(!isTkSelected) continue;
-
- // region of interest + topological cuts + PID
- if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
- Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
- if (!isSelected) continue;
-
- // fill histos
- finvM = dstarD0pi->InvMassD0();
- fInvMass->Fill(finvM);
-
- //DStar invariant mass
- finvMDStar = dstarD0pi->InvMassDstarKpipi();
-
- Double_t EGjet = 0;
- Double_t dStarMom = dstarD0pi->P();
- Double_t phiDStar = dstarD0pi->Phi();
- Double_t phiD0 = theD0particle->Phi();
- //check suggested by Federico
- Double_t dPhiD0Dstar = phiD0 - phiDStar;
-
- dPhi = phiJet - phiDStar;
-
- //plot right dphi
- if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
- if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
-
- Double_t corrFactorCharge = (ejet/100)*20;
- EGjet = ejet + corrFactorCharge;
-
- // fill D* candidates
- Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
- if(finvM >= (mPDGD0-0.05) && finvM <=(mPDGD0+0.05)){ // ~3 sigma (sigma=17MeV, conservative)
-
- if(isDStar == 1) {
- fDphiD0Dstar->Fill(dPhiD0Dstar);
- fDStarMass->Fill(finvMDStar);
- fTrueDiff->Fill(finvMDStar-finvM);
- }
- if(isDStar == 0) fDphiD0Dstar->Fill(dPhiD0Dstar); // angle between D0 and D*
-
- fDStar->Fill(finvMDStar);
- fDiff ->Fill(finvMDStar-finvM);
-
- Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
- Double_t invmassDelta = dstarD0pi->DeltaInvMass();
-
- // now the dphi signal and the fragmentation function
- if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0019){
- //fill candidates D* and soft pion reco pt
-
- fRECOPtDStar->Fill(pt);
- fPtPion->Fill(track2->Pt());
-
- fPhi ->Fill(dPhi);
-
- Double_t jetCone = 0.4;
- if(dPhi>=-jetCone && dPhi<=jetCone){ // evaluate in the near side inside UA1 radius
- Double_t zFrag = (TMath::Cos(dPhi)*dStarMom)/EGjet;
- fResZ->Fill(TMath::Abs(zFrag));
- }
- }
- }
- // evaluate side band background
- SideBandBackground(finvM, finvMDStar, dStarMom, EGjet, dPhi);
-
- } // D* cand
-
-AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
-
-PostData(1,fOutput);
-}
-
-//________________________________________ terminate ___________________________
-void AliAnalysisTaskSEDStarJets::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;
- }
-
- fDStarMass = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
- fTrueDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
- fInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
- fPtPion = dynamic_cast<TH1F*>(fOutput->FindObject("fPtPion "));
- fDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar"));
- fDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff"));
- fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
- ftrigger = dynamic_cast<TH1F*>(fOutput->FindObject("ftrigger"));
- fRECOPtDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtDStar"));
- fRECOPtBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtBkg"));
- fEjet = dynamic_cast<TH1F*>(fOutput->FindObject("fEjet"));
- fPhijet = dynamic_cast<TH1F*>(fOutput->FindObject("fPhijet"));
- fEtaJet = dynamic_cast<TH1F*>(fOutput->FindObject("fEtaJet"));
- fPhi = dynamic_cast<TH1F*>(fOutput->FindObject("fPhi"));
- fResZ = dynamic_cast<TH1F*>(fOutput->FindObject("fResZ"));
- fResZBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fResZBkg"));
- fPhiBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fPhiBkg"));
- theMCFF = dynamic_cast<TH1F*>(fOutput->FindObject("theMCFF"));
- fDphiD0Dstar = dynamic_cast<TH1F*>(fOutput->FindObject("fDphiD0Dstar"));
- fPtJet = dynamic_cast<TH1F*>(fOutput->FindObject("fPtJet"));
-
-}
-//___________________________________________________________________________
-
-void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() {
- // output
- Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
-
- //slot #1
- OpenFile(1);
- fOutput = new TList();
- fOutput->SetOwner();
- // define histograms
- DefineHistoFroAnalysis();
-
- return;
-}
-//___________________________________ hiostograms _______________________________________
-
-Bool_t AliAnalysisTaskSEDStarJets::DefineHistoFroAnalysis(){
-
- // Invariant mass related histograms
- fInvMass = new TH1F("invMass","Kpi invariant mass distribution",1500,.5,3.5);
- fInvMass->SetStats(kTRUE);
- fInvMass->GetXaxis()->SetTitle("GeV/c");
- fInvMass->GetYaxis()->SetTitle("Entries");
- fOutput->Add(fInvMass);
-
- fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4);
- fDStar->SetStats(kTRUE);
- fDStar->GetXaxis()->SetTitle("GeV/c");
- fDStar->GetYaxis()->SetTitle("Entries");
- fOutput->Add(fDStar);
-
- fDiff = new TH1F("Diff","M(kpipi)-M(kpi)",750,0.1,0.2);
- fDiff->SetStats(kTRUE);
- fDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV/c^2");
- fDiff->GetYaxis()->SetTitle("Entries");
- fOutput->Add(fDiff);
-
- fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
- fDiffSideBand->SetStats(kTRUE);
- fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
- fDiffSideBand->GetYaxis()->SetTitle("Entries");
- fOutput->Add(fDiffSideBand);
-
- fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
- fOutput->Add(fDStarMass);
-
- fTrueDiff = new TH1F("dstar","True Reco diff",750,0,0.2);
- fOutput->Add(fTrueDiff);
-
- // trigger normalization
- ftrigger = new TH1F("Normalization","Normalization factor for correlations",1,0,10);
- ftrigger->SetStats(kTRUE);
- fOutput->Add(ftrigger);
-
- //correlation fistograms
- 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);
-
- fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",30,0,30);
- fRECOPtDStar->SetStats(kTRUE);
- fRECOPtDStar->SetLineColor(2);
- fRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
- fRECOPtDStar->GetYaxis()->SetTitle("Entries");
- fOutput->Add(fRECOPtDStar);
-
- fRECOPtBkg = new TH1F("RECOptBkg","RECO pt distribution side bands",30,0,30);
- fRECOPtBkg->SetStats(kTRUE);
- fRECOPtBkg->SetLineColor(2);
- fRECOPtBkg->GetXaxis()->SetTitle("GeV/c");
- fRECOPtBkg->GetYaxis()->SetTitle("Entries");
- fOutput->Add(fRECOPtBkg);
-
- fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,10);
- fPtPion->SetStats(kTRUE);
- fPtPion->GetXaxis()->SetTitle("GeV/c");
- fPtPion->GetYaxis()->SetTitle("Entries");
- fOutput->Add(fPtPion);
-
- // jet related fistograms
- fEjet = new TH1F("ejet", "UA1 algorithm jet energy distribution",1000,0,500);
- fPhijet = new TH1F("Phijet","UA1 algorithm jet #phi distribution", 200,-7,7);
- fEtaJet = new TH1F("Etajet","UA1 algorithm jet #eta distribution", 200,-7,7);
- fPtJet = new TH1F("PtJet", "UA1 algorithm jet Pt distribution",1000,0,500);
- fOutput->Add(fEjet);
- fOutput->Add(fPhijet);
- fOutput->Add(fEtaJet);
- fOutput->Add(fPtJet);
-
- theMCFF = new TH1F("FragFuncMC","Fragmentation function in MC for FC ",100,0,10);
- fResZ = new TH1F("FragFunc","Fragmentation function ",50,0,1);
- fResZBkg = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1);
- fOutput->Add(theMCFF);
- fOutput->Add(fResZ);
- fOutput->Add(fResZBkg);
-
- return kTRUE;
-}
-
-//______________________________ side band background for D*___________________________________
-
-void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t dStarMomBkg, Double_t EGjet, Double_t dPhi){
-
- // 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
-
- if((invM>=1.7 && invM<=1.8) || (invM>=1.92 && invM<=2.19)){
- fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background
- if ((invMDStar-invM)<=0.14732 && (invMDStar-invM)>=0.14352) {
- fPhiBkg->Fill(dPhi);
- fRECOPtBkg->Fill(dStarMomBkg);
- if(dPhi>=-0.4 && dPhi<=0.4){ // evaluate in the near side
- Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/EGjet;
- fResZBkg->Fill(TMath::Abs(zFragBkg));
- }
- }
- }
-}
-
-//_____________________________________________________________________________________________-
-double AliAnalysisTaskSEDStarJets::FillMCFF(AliAODMCParticle* mcPart, TClonesArray* mcArray, Int_t mcLabel){
- //
- // GS from MC
- // UA1 jet algorithm reproduced in MC
- //
- Double_t zMC2 =-999;
-
- Double_t leading =0;
- Double_t PartE = 0;
- Double_t PhiLeading = -999;
- Double_t EtaLeading = -999;
- Double_t PtLeading = -999;
- Int_t counter =-999;
-
- if (!mcPart) return zMC2;
- if (!mcArray) return zMC2;
-
- //find leading particle
- for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
- AliAODMCParticle* Part = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
- if (!Part) {
- AliWarning("MC Particle not found in tree, skipping");
- continue;
- }
-
- // remove quarks and the leading particle (it will be counted later)
- if(iPart == mcLabel) continue;
- if(iPart <= 8) continue;
-
- //remove resonances not directly detected in detector
- Int_t PDGCode = Part->GetPdgCode();
-
- // be sure the particle reach the detector
- Double_t x = Part->Xv();
- Double_t y = Part->Yv();
- Double_t z = Part->Zv();
-
- if(TMath::Abs(PDGCode)== 2212 && x<3 && y<3) continue;
- if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue;
- if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 && TMath::Abs(PDGCode)!=11 && TMath::Abs(PDGCode)!=13 && TMath::Abs(PDGCode)!=2212) continue;
-
- Int_t daug0 = -999;
- Double_t xd =-999;
- Double_t yd =-999;
- Double_t zd =-999;
-
- daug0 = Part->GetDaughter(0);
-
- if(daug0>=0){
- AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
- if(tdaug){
- xd = tdaug->Xv();
- yd = tdaug->Yv();
- zd = tdaug->Zv();
- }
- }
- if(TMath::Abs(xd)<3 || TMath::Abs(yd)<3) continue;
-
- Bool_t AliceAcc = (TMath::Abs(Part->Eta()) <= 0.9);
- if(!AliceAcc) continue;
-
- PartE = Part->E();
-
- if(PartE>leading){
- leading = Part->E();
- PhiLeading = Part->Phi();
- EtaLeading = Part->Eta();
- PtLeading = Part->Pt();
- counter = iPart;
- }
-
- }
-
- Double_t jetEnergy = 0;
-
- //reconstruct the jet
- for (Int_t iiPart=0; iiPart<mcArray->GetEntriesFast(); iiPart++) {
- AliAODMCParticle* tPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iiPart));
- if (!tPart) {
- AliWarning("MC Particle not found in tree, skipping");
- continue;
- }
- // remove quarks and the leading particle (it will be counted later)
- if(iiPart == counter) continue; // do not count again the leading particle
- if(iiPart == mcLabel) continue;
- if(iiPart <= 8) continue;
-
- //remove resonances not directly detected in detector
- Int_t PDGCode = tPart->GetPdgCode();
-
- // be sure the particle reach the detector
- Double_t x = tPart->Xv();
- Double_t y = tPart->Yv();
- Double_t z = tPart->Zv();
-
- if(TMath::Abs(PDGCode)== 2212 && (x<3 && y<3)) continue;
- if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue; // has to be generated at least in the silicon tracker or beam pipe
- if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 && TMath::Abs(PDGCode)!=11 && TMath::Abs(PDGCode)!=13 && TMath::Abs(PDGCode)!=2212) continue;
-
-
- Int_t daug0 = -999;
- Double_t xd =-999;
- Double_t yd =-999;
- Double_t zd =-999;
-
- daug0 = tPart->GetDaughter(0);
-
- if(daug0>=0){
- AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
- if(tdaug){
- xd = tdaug->Xv();
- yd = tdaug->Yv();
- zd = tdaug->Zv();
- }
- }
- if(TMath::Abs(xd)<3 && TMath::Abs(yd)<3) continue;
- //remove particles not in ALICE acceptance
- if(tPart->Pt()<0.07) continue;
- Bool_t AliceAcc = (TMath::Abs(tPart->Eta()) <= 0.9);
- if(!AliceAcc) continue;
-
- Double_t EtaMCp = tPart->Eta();
- Double_t PhiMCp = tPart->Phi();
-
- Double_t DphiMClead = PhiLeading-PhiMCp;
-
- if(DphiMClead<=-(TMath::Pi())/2) DphiMClead = DphiMClead+2*(TMath::Pi());
- if(DphiMClead>(3*(TMath::Pi()))/2) DphiMClead = DphiMClead-2*(TMath::Pi());
-
- Double_t deta = (EtaLeading-EtaMCp);
- //cone radius
- Double_t radius = TMath::Sqrt((DphiMClead*DphiMClead)+(deta*deta));
-
- if(radius>0.4) continue; // in the jet cone
- if(tPart->Charge()==0) continue; // only charged fraction
-
- jetEnergy = jetEnergy+(tPart->E());
- }
-
- jetEnergy = jetEnergy + leading;
-
- // delta phi D*, jet axis
- Double_t dPhi = PhiLeading - (mcPart->Phi());
-
- //plot right dphi
- if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
- if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
-
- zMC2 = (TMath::Cos(dPhi)*(mcPart->P()))/jetEnergy;
-
- return zMC2;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff