X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSGammaJet.cxx;h=5be2356ef2b3ad10f68a29bad94e870fc7b163c4;hb=0f20e40cb80a7bae46238ecb48ce772aee16518a;hp=c724b7c2074e310c4d446c6e906f2d13a6ac1176;hpb=1305bc01348e596eb849da88ed7cd5bbad2a3cd7;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSGammaJet.cxx b/PHOS/AliPHOSGammaJet.cxx index c724b7c2074..5be2356ef2b 100644 --- a/PHOS/AliPHOSGammaJet.cxx +++ b/PHOS/AliPHOSGammaJet.cxx @@ -14,9 +14,34 @@ **************************************************************************/ /* $Id$ */ +/* History of cvs commits: + * + * $Log$ + * Revision 1.11 2006/01/31 20:30:52 hristov + * Including TFile.h + * + * Revision 1.10 2006/01/23 18:04:08 hristov + * Removing meaningless const + * + * Revision 1.9 2006/01/12 16:23:26 schutz + * ESD is properly read with methods of macros/AliReadESD.C copied in it + * + * Revision 1.8 2005/12/20 07:08:32 schutz + * corrected error in call AliReadESD + * + * Revision 1.6 2005/05/28 14:19:04 schutz + * Compilation warnings fixed by T.P. + * + */ + //_________________________________________________________________________ // Class for the analysis of gamma-jet correlations -// +// Basically it seaches for a prompt photon in the PHOS acceptance, +// if so we construct a jet around the highest pt particle in the opposite +// side in azimuth, inside the TPC and EMCAL acceptances. First the leading +// particle and then the jet have to fullfill several conditions +// (energy, direction ..) to be accepted. Then the fragmentation function +// of this jet is constructed // //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN) ////////////////////////////////////////////////////////////////////////////// @@ -24,24 +49,22 @@ // --- ROOT system --- -#include "TString.h" -#include "TFile.h" -#include "TLorentzVector.h" -#include "TList.h" -#include "TTree.h" -#include "TNtuple.h" -#include "TParticle.h" -#include "TCanvas.h" -#include "TPaveLabel.h" -#include "TPad.h" -#include "TArrayC.h" -#include "AliRun.h" +#include +#include +#include +#include +#include +#include + #include "AliPHOSGammaJet.h" #include "AliPHOSGetter.h" #include "AliPHOSGeometry.h" -#include "TF1.h" +#include "AliPHOSFastGlobalReconstruction.h" +#include "AliESD.h" +#include "AliESDtrack.h" +#include "AliESDCaloCluster.h" #include "Riostream.h" -#include "../PYTHIA6/AliGenPythia.h" + ClassImp(AliPHOSGammaJet) @@ -51,9 +74,10 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask() // ctor fAngleMaxParam.Set(4) ; fAngleMaxParam.Reset(0.); - fAnyConeOrPt = kFALSE ; + fAnyConeOrPt = 0 ; fEtaCut = 0. ; - fFastRec = 0 ; + fESDdata = 0 ; + fFastRec = 0 ; fInvMassMaxCut = 0. ; fInvMassMinCut = 0. ; fJetRatioMaxCut = 0. ; @@ -61,11 +85,11 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask() fJetTPCRatioMaxCut = 0. ; fJetTPCRatioMinCut = 0. ; fMinDistance = 0. ; - fNEvent = 0 ; - fNCone = 0 ; - fNPt = 0 ; - fOnlyCharged = kFALSE ; - fOptFast = kFALSE ; + fNEvent = 0 ; + fNCone = 0 ; + fNPt = 0 ; + fOnlyCharged = 0 ; + fOptFast = 0 ; fPhiMaxCut = 0. ; fPhiMinCut = 0. ; fPtCut = 0. ; @@ -75,8 +99,13 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask() fRatioMinCut = 0. ; fCone = 0 ; fPtThreshold = 0 ; - fTPCCutsLikeEMCAL = kFALSE ; - fSelect = kFALSE ; + fTPCCutsLikeEMCAL = 0 ; + fSelect = 0 ; + + fDirName = "" ; + fESDTree = "" ; + fPattern = "" ; + for(Int_t i = 0; i<10; i++){ fCones[i] = 0.0 ; fNameCones[i] = "" ; @@ -100,7 +129,7 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask() } } - fOption = "" ; + fOptionGJ = "" ; fOutputFile = new TFile(gDirectory->GetName()) ; fInputFileName = gDirectory->GetName() ; fOutputFileName = gDirectory->GetName() ; @@ -130,7 +159,7 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask() AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) : TTask("GammaJet","Analysis of gamma-jet correlations") { - // ctor + // ctor fInputFileName = inputfilename; fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName); AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ; @@ -145,11 +174,12 @@ AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) // cpy ctor fAngleMaxParam = gj.fAngleMaxParam; fAnyConeOrPt = gj.fAnyConeOrPt; + fESDdata = gj.fESDdata; fEtaCut = gj.fEtaCut ; fInvMassMaxCut = gj.fInvMassMaxCut ; fInvMassMinCut = gj.fInvMassMinCut ; fFastRec = gj.fFastRec ; - fOption = gj.fOption ; + fOptionGJ = gj.fOptionGJ ; fMinDistance = gj.fMinDistance ; fOptFast = gj.fOptFast ; fOnlyCharged = gj.fOnlyCharged ; @@ -183,6 +213,10 @@ AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) fCone = gj.fCone ; fPtThreshold = gj.fPtThreshold ; + fDirName = gj.fDirName ; + fESDTree = gj.fESDTree ; + fPattern = gj.fPattern ; + SetName (gj.GetName()) ; SetTitle(gj.GetTitle()) ; @@ -220,8 +254,7 @@ AliPHOSGammaJet::~AliPHOSGammaJet() void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, TClonesArray * plCh, TClonesArray * plNe, - TClonesArray * plNePHOS, - const AliPHOSGeometry * geom ) + TClonesArray * plNePHOS) { // List of particles copied to a root file. @@ -242,11 +275,11 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, Char_t fi[100]; sprintf(fi,"bgrd/%d/list.root",iEvent); TFile * f = TFile::Open(fi) ; - cout<GetName()<GetName()<Get("T"); - TBranch *branch = T->GetBranch("primaries"); + TTree *t = (TTree*) f->Get("T"); + TBranch *branch = t->GetBranch("primaries"); branch->SetAddress(&particle); Int_t index = particleList->GetEntries() ; @@ -257,10 +290,14 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, Int_t iParticle = 0 ; Int_t m = 0; Double_t x = 0., z = 0.; - // cout<<"bkg entries "<GetEntries()<GetEntries()<PHOSGeometry() ; + if(!fOptFast){ - for (iParticle=0 ; iParticle < T->GetEntries() ; iParticle++) { - T->GetEvent(iParticle) ; + for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) { + t->GetEvent(iParticle) ; m = 0 ; x = 0. ; @@ -280,7 +317,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, (dynamic_cast(plCh->At(indexCh)))->SetStatusCode(0) ; indexCh++ ; - if(strstr(fOption,"deb all")||strstr(fOption,"deb")){ + if(strstr(fOptionGJ,"deb all")||strstr(fOptionGJ,"deb")){ dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); } @@ -290,7 +327,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, if(m != 0) {//Is in PHOS - if(strstr(fOption,"deb all")|| strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); @@ -302,7 +339,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, if((particle->Phi()>fPhiEMCALCut[0]) && (particle->Phi()(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); @@ -328,8 +365,8 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, // fFastRec = new AliPHOSFastGlobalReconstruction(fHIJINGFileName); // fFastRec->FastReconstruction(iiEvent); - for (iParticle=0 ; iParticle < T->GetEntries() ; iParticle++) { - T->GetEvent(iParticle) ; + for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) { + t->GetEvent(iParticle) ; m = 0 ; x = 0. ; z = 0. ; @@ -397,7 +434,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, //Check if decay photons are too close for PHOS cellDistance = angle*460; //cm if (cellDistance < fMinDistance) { - if(strstr(fOption,"deb all")|| strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); @@ -423,7 +460,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, indexNe++ ; } if(m != 0){ - if(strstr(fOption,"deb all")|| strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ; @@ -452,7 +489,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x); if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()(fListHistos->FindObject("PtSpectra"))-> Fill(photon1->Pt()); new((*particleList)[index]) TParticle(*photon1) ; @@ -465,7 +502,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, p1 = kTRUE; } if(m != 0){ - if(strstr(fOption,"deb all") || strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(photon1->Pt()); new((*plNePHOS)[indexNePHOS]) TParticle(*photon1) ; @@ -486,7 +523,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x); if(photon2->Phi()>fPhiEMCALCut[0] && photon2->Phi()(fListHistos->FindObject("PtSpectra"))-> Fill(photon2->Pt()); new((*particleList)[index]) TParticle(*photon2) ; @@ -498,7 +535,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, indexNe++ ; } if(m != 0){ - if(strstr(fOption,"deb all") || strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(photon2->Pt()); new((*plNePHOS)[indexNePHOS]) TParticle(*photon2) ; @@ -529,14 +566,14 @@ Double_t AliPHOSGammaJet::CalculateJetRatioLimit(const Double_t ptg, const Double_t *x) { //Info("CalculateLimit","x1 %f, x2%f",x[0],x[1]); - Double_t Epp = par[0] + par[1] * ptg ; - Double_t Spp = par[2] + par[3] * ptg ; + Double_t epp = par[0] + par[1] * ptg ; + Double_t spp = par[2] + par[3] * ptg ; Double_t f = x[0] + x[1] * ptg ; - Double_t Epb = Epp + par[4] ; - Double_t Spb = TMath::Sqrt(Spp*Spp+ par[5]*par[5]) ; - Double_t rat = (Epb - Spb * f) / ptg ; - //Info("CalculateLimit","Epp %f, Spp %f, f %f", Epp, Spp, f); - //Info("CalculateLimit","Epb %f, Spb %f, rat %f", Epb, Spb, rat); + Double_t epb = epp + par[4] ; + Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ; + Double_t rat = (epb - spb * f) / ptg ; + //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f); + //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat); return rat ; } @@ -545,13 +582,14 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, TClonesArray * particleList, TClonesArray * plCh, TClonesArray * plNe, - TClonesArray * plNePHOS, - const AliPHOSGeometry * geom ) + TClonesArray * plNePHOS) { //Info("CreateParticleList","Inside"); - AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ; + AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ; + const AliPHOSGeometry * geom = gime->PHOSGeometry() ; gime->Event(iEvent, "X") ; + Int_t index = particleList->GetEntries() ; Int_t indexCh = plCh->GetEntries() ; Int_t indexNe = plNe->GetEntries() ; @@ -579,7 +617,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, ->GetParticle(particle->GetPdgCode())->Charge(); if((charge != 0) && (particle->Pt() > fChargedPtCut)){ - if(strstr(fOption,"deb all")|| strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); new((*plCh)[indexCh++]) TParticle(*particle) ; @@ -589,7 +627,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x); if(m != 0) {//Is in PHOS - if(strstr(fOption,"deb all")|| strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); @@ -598,7 +636,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, if((particle->Phi()>fPhiEMCALCut[0]) && (particle->Phi()(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); new((*plNe)[indexNe++]) TParticle(*particle) ; @@ -693,7 +731,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, if(part.Pt() > fNeutralPtCut){ if(particle->Phi()>fPhiEMCALCut[0] && particle->Phi()(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); @@ -707,7 +745,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, indexNe++; }//InEMCAL if(m != 0){ - if(strstr(fOption,"deb all")|| strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(particle->Pt()); new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ; @@ -734,7 +772,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x); if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()(fListHistos->FindObject("PtSpectra"))-> Fill(photon1->Pt()); new((*plNe)[indexNe++]) TParticle(*photon1) ; @@ -743,7 +781,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, p1 = kTRUE; } if(m != 0){ - if(strstr(fOption,"deb all") || strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(photon1->Pt()); new((*plNePHOS)[indexNePHOS++]) TParticle(*photon1) ; @@ -761,14 +799,14 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x); if(photon2->Phi()>fPhiEMCALCut[0] && photon2->Phi()(fListHistos->FindObject("PtSpectra"))-> Fill(photon2->Pt()); new((*plNe)[indexNe++]) TParticle(*photon2) ; new((*particleList)[index++]) TParticle(*photon2) ; } if(m != 0){ - if(strstr(fOption,"deb all") || strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb")) dynamic_cast(fListHistos->FindObject("PtSpectra"))-> Fill(photon2->Pt()); new((*plNePHOS)[indexNePHOS++]) TParticle(*photon2) ; @@ -791,6 +829,134 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, }//OptFast //gime->Delete() ; } +//____________________________________________________________________________ +void AliPHOSGammaJet::CreateParticleListFromESD(TClonesArray * pl, + TClonesArray * plCh, + TClonesArray * plNe, + TClonesArray * plNePHOS, + const AliESD * esd){ + + //Create a list of particles from the ESD. These particles have been measured + //by the Central Tracking system (TPC), PHOS and EMCAL + //(EMCAL not available for the moment). + //Info("CreateParticleListFromESD","Inside"); + + Int_t index = pl->GetEntries() ; + Int_t npar = 0 ; + Float_t *pid = new Float_t[AliPID::kSPECIESN]; + + //########### PHOS ############## + //Info("CreateParticleListFromESD","Fill ESD PHOS list"); + Int_t begphos = esd->GetFirstPHOSCluster(); + Int_t endphos = esd->GetFirstPHOSCluster() + + esd->GetNumberOfPHOSClusters() ; + Int_t indexNePHOS = plNePHOS->GetEntries() ; + if(strstr(fOptionGJ,"deb all")) + Info("CreateParticleListFromESD","PHOS: first particle %d, last particle %d", + begphos,endphos); + + for (npar = begphos; npar < endphos; npar++) {//////////////PHOS track loop + AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd + + //Create a TParticle to fill the particle list + + Float_t en = clus->GetClusterEnergy() ; + Float_t *p = new Float_t(); + clus->GetGlobalPosition(p) ; + TVector3 pos(p[0],p[1],p[2]) ; + Double_t phi = pos.Phi(); + Double_t theta= pos.Theta(); + Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);; + Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta); + Double_t pz = en*TMath::Cos(theta); + + TParticle * particle = new TParticle() ; + particle->SetMomentum(px,py,pz,en) ; + + //Select only photons + + pid=clus->GetPid(); + //cout<<"pid "< 0.75) + new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ; + } + + //########### TPC ##################### + //Info("CreateParticleListFromESD","Fill ESD TPC list"); + Int_t begtpc = 0 ; + Int_t endtpc = esd->GetNumberOfTracks() ; + Int_t indexCh = plCh->GetEntries() ; + if(strstr(fOptionGJ,"deb all")) + Info("CreateParticleListFromESD","TPC: first particle %d, last particle %d", + begtpc,endtpc); + + for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop + AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd + + Double_t en = track ->GetTPCsignal() ; + TVector3 mom = track->P3() ; + Double_t px = mom.Px(); + Double_t py = mom.Py(); + Double_t pz = mom.Pz(); //Check with TPC people if this is correct. + + //cout<<"TPC signal "<SetMomentum(px,py,pz,en) ; + + new((*plCh)[indexCh++]) TParticle(*particle) ; + new((*pl)[index++]) TParticle(*particle) ; + + } + + //################ EMCAL ############## + Double_t v[3] ; //vertex ; + esd->GetVertex()->GetXYZ(v) ; + //##########Uncomment when ESD for EMCAL works ########## + //Info("CreateParticleListFromESD","Fill ESD EMCAL list"); + + Int_t begem = esd->GetFirstEMCALCluster(); + Int_t endem = esd->GetFirstEMCALCluster() + + esd->GetNumberOfEMCALClusters() ; + Int_t indexNe = plNe->GetEntries() ; + if(strstr(fOptionGJ,"deb all")) + Info("CreateParticleListFromESD","EMCAL: first particle %d, last particle %d", + begem,endem); + + for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop + AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd + + Float_t en = clus->GetClusterEnergy() ; + Float_t *p = new Float_t(); + clus->GetGlobalPosition(p) ; + TVector3 pos(p[0],p[1],p[2]) ; + Double_t phi = pos.Phi(); + Double_t theta= pos.Theta(); + Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);; + Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta); + Double_t pz = en*TMath::Cos(theta); + //cout<<"EMCAL signal "<SetMomentum(px,py,pz,en) ; + + Int_t pdg = 0; + // //Uncomment if PID IS WORKING, photon and pi0 idenitification. + // //if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen. + // //pdg = 22; + // //else if( pid[AliPID::kPi0] > 0.75) + // //pdg = 111; + pdg = 22; //No PID, assume all photons + TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, + px, py, pz, en, v[0], v[1], v[2], 0); + + new((*plNe)[indexNe++]) TParticle(*particle) ; + new((*pl)[index++]) TParticle(*particle) ; + } + + // Info("CreateParticleListFromESD","End Inside"); + +} @@ -798,12 +964,26 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent, void AliPHOSGammaJet::Exec(Option_t *option) { // does the job - fOption = option; + fOptionGJ = option; MakeHistos() ; + + AliESD * esd = 0; + TChain * t = 0 ; + + if(fESDdata){ + // Create chain of esd trees + const UInt_t kNevent = static_cast(GetNEvent()) ; + t = ReadESD(kNevent, fDirName, fESDTree, fPattern) ; + if(!t) { + AliError("Could not create the TChain") ; + //break ; + } + + // ESD + t->SetBranchAddress("ESD",&esd); // point to the container esd where to put the event from the esdTree - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; - const AliPHOSGeometry * geom = gime->PHOSGeometry() ; - + } + // AliGenPythia* pyth = (AliGenPythia*) gAlice->Generator(); // pyth->Init(); @@ -814,7 +994,7 @@ void AliPHOSGammaJet::Exec(Option_t *option) TClonesArray * plNePHOS = new TClonesArray("TParticle",1000); for (Int_t iEvent = 0 ; iEvent < fNEvent ; iEvent++) { - if(strstr(fOption,"deb")||strstr(fOption,"deb all")) + if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all")) Info("Exec", "Event %d", iEvent) ; fRan.SetSeed(0); @@ -825,45 +1005,56 @@ void AliPHOSGammaJet::Exec(Option_t *option) TLorentzVector jet (0,0,0,0); TLorentzVector jettpc(0,0,0,0); - - CreateParticleList(iEvent, particleList, plCh,plNe,plNePHOS, geom ); -// TLorentzVector pyjet(0,0,0,0); + if(fESDdata){ -// Int_t nJ, nJT; -// Float_t jets[4][10]; -// pyth->SetJetReconstructionMode(1); -// pyth->LoadEvent(); -// pyth->GetJets(nJ, nJT, jets); - -// Float_t pxJ = jets[0][0]; -// Float_t pyJ = jets[1][0]; -// Float_t pzJ = jets[2][0]; -// Float_t eJ = jets[3][0]; -// pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ; - -// if(nJT > 1){ -// //Info("Exec",">>>>>>>>>>Number of jets !!!! %d",nJT); -// for (Int_t iJ = 1; iJ < nJT; iJ++) { -// Float_t pxJ = jets[0][iJ]; -// Float_t pyJ = jets[1][iJ]; -// Float_t pzJ = jets[2][iJ]; -// Float_t eJ = jets[3][iJ]; -// pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ; -// //Info("Exec",">>>>>Pythia Jet: %d, Phi %f, Eta %f, Pt %f", -// // iJ,pyjet.Phi(),pyjet.Eta(),pyjet.Pt()); -// } + Int_t iNbytes = t->GetEntry(iEvent); // store event in esd + //cout<<"nbytes "<SetJetReconstructionMode(1); + // pyth->LoadEvent(); + // pyth->GetJets(nJ, nJT, jets); + + // Float_t pxJ = jets[0][0]; + // Float_t pyJ = jets[1][0]; + // Float_t pzJ = jets[2][0]; + // Float_t eJ = jets[3][0]; + // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ; -// } + // if(nJT > 1){ + // //Info("Exec",">>>>>>>>>>Number of jets !!!! %d",nJT); + // for (Int_t iJ = 1; iJ < nJT; iJ++) { + // Float_t pxJ = jets[0][iJ]; + // Float_t pyJ = jets[1][iJ]; + // Float_t pzJ = jets[2][iJ]; + // Float_t eJ = jets[3][iJ]; + // pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ; + // //Info("Exec",">>>>>Pythia Jet: %d, Phi %f, Eta %f, Pt %f", + // // iJ,pyjet.Phi(),pyjet.Eta(),pyjet.Pt()); + // } + + // } + + if(fHIJING) + AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS); + } - if(fHIJING) - AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS, geom); - - - Bool_t IsInPHOS = kFALSE ; - GetGammaJet(plNePHOS, ptg, phig, etag, IsInPHOS) ; + Bool_t iIsInPHOS = kFALSE ; + GetGammaJet(plNePHOS, ptg, phig, etag, iIsInPHOS) ; - if(IsInPHOS){ + if(iIsInPHOS){ //Info("Exec"," In PHOS") ; dynamic_cast(fListHistos->FindObject("NGamma"))->Fill(ptg); @@ -871,7 +1062,7 @@ void AliPHOSGammaJet::Exec(Option_t *option) ->Fill(ptg,phig); dynamic_cast(fListHistos->FindObject("EtaGamma")) ->Fill(ptg,etag); - if(strstr(fOption,"deb")||strstr(fOption,"deb all")) + if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all")) Info("Exec", "Gamma: pt %f, phi %f, eta %f", ptg, phig,etag) ; @@ -912,7 +1103,7 @@ void AliPHOSGammaJet::Exec(Option_t *option) ->Fill(ptg,phig-phich); dynamic_cast(fListHistos->FindObject("DeltaEtaCharge")) ->Fill(ptg,etag-etach); - if(strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb")) Info("Exec"," Charged Leading") ; } if((ptpi > ptch) && insidepi){ @@ -927,18 +1118,18 @@ void AliPHOSGammaJet::Exec(Option_t *option) dynamic_cast(fListHistos->FindObject("DeltaEtaPi0")) ->Fill(ptg,etag-etapi); - if(ptpi > 0. && strstr(fOption,"deb")) + if(ptpi > 0. && strstr(fOptionGJ,"deb")) Info("Exec"," Pi0 Leading") ; } - if(strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb")) Info("Exec","Leading pt %f, phi %f",ptl,phil); if(insidech || insidepi){ if(!fAnyConeOrPt){ MakeJet(particleList, ptg, phig, ptl, phil, etal, "", jet); - if(strstr(fOption,"deb")){ + if(strstr(fOptionGJ,"deb")){ // Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f", // pyjet.Phi(),pyjet.Eta(),pyjet.Pt()); Info("Exec","TPC+EMCAL Jet: Phi %f, Eta %f, Pt %f", @@ -958,7 +1149,7 @@ void AliPHOSGammaJet::Exec(Option_t *option) //TPC if(fOnlyCharged && ptch > 0.) { - if(strstr(fOption,"deb")) + if(strstr(fOptionGJ,"deb")) Info("Exec","Leading TPC pt %f, phi %f",ptch,phich); dynamic_cast(fListHistos->FindObject("TPCRatio")) @@ -972,7 +1163,7 @@ void AliPHOSGammaJet::Exec(Option_t *option) MakeJet(plCh, ptg, phig, ptch, phich, etach, "TPC",jettpc); - if(strstr(fOption,"deb")){ + if(strstr(fOptionGJ,"deb")){ // Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f", // pyjet.Phi(),pyjet.Eta(),pyjet.Pt()); Info("Exec","TPC Jet: Phi %f, Eta %f, Pt %f", @@ -1011,7 +1202,8 @@ void AliPHOSGammaJet::Exec(Option_t *option) void AliPHOSGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg, TString conf, TString type) { - + //Fill jet fragmentation histograms if !fAnyCone, + //only for fCone and fPtThres TParticle * particle = 0 ; Int_t ipr = -1 ; Float_t charge = 0; @@ -1042,7 +1234,8 @@ void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg, TString conf, TString type, TString cone, TString ptcut) { - + //Fill jet fragmentation histograms if fAnyCone, + //for several cones and pt thresholds TParticle *particle = 0; Int_t ipr=-1; Float_t charge = 0; @@ -1072,9 +1265,9 @@ void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg, //____________________________________________________________________________ void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt, - Double_t &phi, Double_t &eta, Bool_t &Is) - + Double_t &phi, Double_t &eta, Bool_t &Is) const { + //Search for the prompt photon in PHOS with pt > fPtCut pt = -10.; eta = -10.; phi = -10.; @@ -1082,15 +1275,12 @@ void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt, for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){ TParticle * particle = dynamic_cast(pl->At(ipr)) ; - if( (particle->Pt() > fPtCut ) && ( particle->GetStatusCode() == 1) - && ( particle->GetPdgCode() == 22 || particle->GetPdgCode() == 111) ){ - - if (particle->Pt() > pt) { - pt = particle->Pt(); - phi = particle->Phi() ; - eta = particle->Eta() ; - Is = kTRUE; - } + if((particle->Pt() > fPtCut) && (particle->Pt() > pt)){ + + pt = particle->Pt(); + phi = particle->Phi() ; + eta = particle->Eta() ; + Is = kTRUE; } } } @@ -1098,9 +1288,10 @@ void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt, //____________________________________________________________________________ void AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl, Double_t ptg, Double_t phig, - Double_t &pt, Double_t &eta, Double_t &phi) + Double_t &pt, Double_t &eta, Double_t &phi) const { - + //Search for the charged particle with highest with + //Phi=Phi_gamma-Pi and pT=0.1E_gamma pt = -100.; eta = -100; phi = -100; @@ -1129,8 +1320,11 @@ void AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl, //____________________________________________________________________________ void AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl, Double_t ptg, Double_t phig, - Double_t &pt, Double_t &eta, Double_t &phi) + Double_t &pt, Double_t &eta, Double_t &phi) { + + //Search for the neutral pion with highest with + //Phi=Phi_gamma-Pi and pT=0.1E_gamma pt = -100.; eta = -100.; phi = -100.; @@ -1209,8 +1403,7 @@ void AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl, while ( (particle = (TParticle*)next()) ) { iPrimary++; - // if(particle->GetStatusCode() == 1){ - + ksPdg = particle->GetPdgCode(); ptl = particle->Pt(); if(ksPdg == 111){ //2 gamma @@ -1314,6 +1507,9 @@ void AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl, //____________________________________________________________________________ void AliPHOSGammaJet::InitParameters() { + + //Initialize the parameters of the analysis. + fAngleMaxParam.Set(4) ; fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4}; fAngleMaxParam.AddAt(-0.25,1) ; @@ -1321,7 +1517,7 @@ void AliPHOSGammaJet::InitParameters() fAngleMaxParam.AddAt(-2e-4,3) ; fAnyConeOrPt = kFALSE ; fOutputFileName = "GammaJet.root" ; - fOption = ""; + fOptionGJ = ""; fHIJINGFileName = "galice.root" ; fHIJING = kFALSE ; fMinDistance = 3.6 ; @@ -1330,6 +1526,7 @@ void AliPHOSGammaJet::InitParameters() fInvMassMinCut = 0.12 ; fOnlyCharged = kFALSE ; fOptFast = kFALSE ; + fESDdata = kTRUE ; fPhiEMCALCut[0] = 60. *TMath::Pi()/180.; fPhiEMCALCut[1] = 180.*TMath::Pi()/180.; fPhiMaxCut = 3.4 ; @@ -1348,6 +1545,10 @@ void AliPHOSGammaJet::InitParameters() fJetTPCRatioMinCut = 0.3 ; fSelect = kFALSE ; + fDirName = "./" ; + fESDTree = "esdTree" ; + fPattern = "." ; + //Cut depending on gamma energy fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applyed @@ -1408,7 +1609,9 @@ void AliPHOSGammaJet::InitParameters() } //__________________________________________________________________________- -Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e){ +Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) { + //Check if the opening angle of the candidate pairs is inside + //our selection windowd Bool_t result = kFALSE; Double_t mpi0 = 0.1349766; Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e) @@ -1427,6 +1630,7 @@ Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e){ //__________________________________________________________________________- Bool_t AliPHOSGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj, const TString type ){ + //Check if the energy of the reconstructed jet is within an energy window Double_t par[6]; Double_t xmax[2]; @@ -2184,7 +2388,9 @@ void AliPHOSGammaJet::MakeJet(TClonesArray * pl, Double_t ptl, Double_t phil, Double_t etal, TString conf, TLorentzVector & jet) { - + //Fill the jet with the particles around the leading particle with + //R=fCone and pt_th = fPtThres. Calculate the energy of the jet and + //check if we select it. Fill jet histograms Float_t ptcut = 0. ; if(fHIJING){ if(ptg > fPtJetSelectionCut) ptcut = 2. ; @@ -2239,7 +2445,7 @@ void AliPHOSGammaJet::MakeJet(TClonesArray * pl, ptjet = jet.Pt(); ptbkg = bkg.Pt(); - if(strstr(fOption,"deb") || strstr(fOption,"deb all")) + if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all")) Info("MakeJet","Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg); @@ -2259,7 +2465,7 @@ void AliPHOSGammaJet::MakeJet(TClonesArray * pl, if(IsJetSelected(ptg,ptjet,conf) || fSelect){ - if(strstr(fOption,"deb") || strstr(fOption,"deb all")) + if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all")) Info("MakeJet","JetSelected"); dynamic_cast(fListHistos->FindObject("N"+conf+"Jet"))-> Fill(ptg); @@ -2278,6 +2484,10 @@ void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg, Double_t phig, Double_t ptl, Double_t phil, Double_t etal, TString conf){ + + //Fill the jet with the particles around the leading particle with + //R=fCone(i) and pt_th = fPtThres(i). Calculate the energy of the jet and + //check if we select it. Fill jet i histograms TClonesArray * jetList = new TClonesArray("TParticle",1000); TClonesArray * bkgList = new TClonesArray("TParticle",1000); @@ -2337,7 +2547,7 @@ void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg, //Fill histograms if(ptjet > 0.) { - if(strstr(fOption,"deb")){ + if(strstr(fOptionGJ,"deb")){ Info("MakeJetAnyPt","cone %f, ptcut %f",fCones[icone],fPtThres[ipt]); Info("MakeJetAnyPt","pT: Gamma %f, Jet %f, Bkg %f",ptg,ptjet,ptbkg); } @@ -2387,7 +2597,7 @@ void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg, //____________________________________________________________________________ void AliPHOSGammaJet::MakePhoton(TLorentzVector & particle) { - + //Fast reconstruction for photons Double_t energy = particle.E() ; Double_t modenergy = MakeEnergy(energy) ; //Info("MakePhoton","Energy %f, Modif %f",energy,modenergy); @@ -2475,6 +2685,7 @@ void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0, //____________________________________________________________________________ void AliPHOSGammaJet::Plot(TString what, Option_t * option) const { + //Plot some relevant histograms of the analysis TH2F * h = dynamic_cast(fOutputFile->Get(what)); if(h){ h->Draw(); @@ -2520,12 +2731,13 @@ void AliPHOSGammaJet::Plot(TString what, Option_t * option) const } //____________________________________________________________________________ -void AliPHOSGammaJet::Print(char * opt) +void AliPHOSGammaJet::Print(const Option_t * opt) const { + + //Print some relevant parameters set for the analysis if(! opt) return; - // Print Info("Print", "%s %s", GetName(), GetTitle() ) ; printf("Eta cut : %f\n", fEtaCut) ; @@ -2543,11 +2755,73 @@ void AliPHOSGammaJet::Print(char * opt) } +//__________________________________________________________________________ +TChain * AliPHOSGammaJet::ReadESDfromdisk(const UInt_t eventsToRead, + const TString dirName, + const TString esdTreeName, + const char * pattern) +{ + // Reads ESDs from Disk + TChain * rv = 0; + + AliInfo( Form("\nReading files in %s \nESD tree name is %s \nReading %d events", + dirName.Data(), esdTreeName.Data(), eventsToRead) ) ; + + // create a TChain of all the files + TChain * cESDTree = new TChain(esdTreeName) ; + + // read from the directory file until the require number of events are collected + void * from = gSystem->OpenDirectory(dirName) ; + if (!from) { + AliError( Form("Directory %s does not exist") ) ; + rv = 0 ; + } + else{ // reading file names from directory + const char * subdir ; + // search all subdirectories witch matching pattern + while( (subdir = gSystem->GetDirEntry(from)) && + (cESDTree->GetEntries() < eventsToRead)) { + if ( strstr(subdir, pattern) != 0 ) { + char file[200] ; + sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir); + AliInfo( Form("Adding %s\n", file) ); + cESDTree->Add(file) ; + } + } // while file + + AliInfo( Form(" %d events read", cESDTree->GetEntriesFast()) ) ; + rv = cESDTree ; + + } // reading file names from directory + return rv ; +} + +//__________________________________________________________________________ +TChain * AliPHOSGammaJet::ReadESD(const UInt_t eventsToRead, + const TString dirName, + const TString esdTreeName, + const char * pattern) +{ + // Read AliESDs files and return a Chain of events + + if ( dirName == "" ) { + AliError("Give the name of the DIR where to find files") ; + return 0 ; + } + if ( esdTreeName == "" ) + return ReadESDfromdisk(eventsToRead, dirName) ; + else if ( strcmp(pattern, "") == 0 ) + return ReadESDfromdisk(eventsToRead, dirName, esdTreeName) ; + else + return ReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ; +} + //___________________________________________________________________ void AliPHOSGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone, Double_t eta, Double_t phi) { - + + //Check if the particle is inside the cone defined by the leading particle b = kFALSE; if(phi > TMath::TwoPi())