X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PHOS%2FAliPHOSGammaJet.cxx;h=b04234927419d0e88325bcd9be7fd072c5b48602;hp=04f26a5d7c54449e8ea1d25495e10816214703f7;hb=bfbd5665ee83ec5e12fe28ed27bff3182595f115;hpb=e957fea89c462908a01dcfa4f611e99a98080071 diff --git a/PHOS/AliPHOSGammaJet.cxx b/PHOS/AliPHOSGammaJet.cxx index 04f26a5d7c5..b0423492741 100644 --- a/PHOS/AliPHOSGammaJet.cxx +++ b/PHOS/AliPHOSGammaJet.cxx @@ -14,124 +14,2635 @@ **************************************************************************/ /* $Id$ */ +/* History of cvs commits: + * + * $Log$ + * 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) ////////////////////////////////////////////////////////////////////////////// // --- ROOT system --- -#include "TRandom.h" -#include "TLorentzVector.h" -#include "TList.h" -#include "TParticle.h" +#include +#include +#include +#include +#include +#include + #include "AliPHOSGammaJet.h" #include "AliPHOSGetter.h" +#include "AliPHOSGeometry.h" +#include "AliPHOSFastGlobalReconstruction.h" +#include "AliESD.h" +#include "AliESDtrack.h" +#include "Riostream.h" + ClassImp(AliPHOSGammaJet) //____________________________________________________________________________ -AliPHOSGammaJet::AliPHOSGammaJet() { +AliPHOSGammaJet::AliPHOSGammaJet() : TTask() +{ // ctor + fAngleMaxParam.Set(4) ; + fAngleMaxParam.Reset(0.); + fAnyConeOrPt = 0 ; + fEtaCut = 0. ; + fESDdata = 0 ; + fFastRec = 0 ; + fInvMassMaxCut = 0. ; + fInvMassMinCut = 0. ; + fJetRatioMaxCut = 0. ; + fJetRatioMinCut = 0. ; + fJetTPCRatioMaxCut = 0. ; + fJetTPCRatioMinCut = 0. ; + fMinDistance = 0. ; + fNEvent = 0 ; + fNCone = 0 ; + fNPt = 0 ; + fOnlyCharged = 0 ; + fOptFast = 0 ; + fPhiMaxCut = 0. ; + fPhiMinCut = 0. ; + fPtCut = 0. ; + fNeutralPtCut = 0. ; + fChargedPtCut = 0. ; + fRatioMaxCut = 0. ; + fRatioMinCut = 0. ; + fCone = 0 ; + fPtThreshold = 0 ; + fTPCCutsLikeEMCAL = 0 ; + fSelect = 0 ; + + fDirName = "" ; + fESDTree = "" ; + fPattern = "" ; + + for(Int_t i = 0; i<10; i++){ + fCones[i] = 0.0 ; + fNameCones[i] = "" ; + fPtThres[i] = 0.0 ; + fPtJetSelectionCut = 0.0 ; + fNamePtThres[i] = "" ; + if( i < 6 ){ + fJetXMin1[i] = 0.0 ; + fJetXMin2[i] = 0.0 ; + fJetXMax1[i] = 0.0 ; + fJetXMax2[i] = 0.0 ; + fBkgMean[i] = 0.0 ; + fBkgRMS[i] = 0.0 ; + if( i < 2 ){ + fJetE1[i] = 0.0 ; + fJetE2[i] = 0.0 ; + fJetSigma1[i] = 0.0 ; + fJetSigma2[i] = 0.0 ; + fPhiEMCALCut[i] = 0.0 ; + } + } + } + + fOptionGJ = "" ; + fOutputFile = new TFile(gDirectory->GetName()) ; + fInputFileName = gDirectory->GetName() ; + fOutputFileName = gDirectory->GetName() ; + fHIJINGFileName = gDirectory->GetName() ; + fHIJING = 0 ; + fPosParaA = 0. ; + fPosParaB = 0. ; + fRan = 0 ; + fResPara1 = 0. ; + fResPara2 = 0. ; + fResPara3 = 0. ; + + fListHistos = new TObjArray(100) ; + TList * list = gDirectory->GetListOfKeys() ; + TIter next(list) ; + TH2F * h = 0 ; + Int_t index ; + for (index = 0 ; index < list->GetSize()-1 ; index++) { + //-1 to avoid GammaJet Task + h = dynamic_cast(gDirectory->Get(list->At(index)->GetName())) ; + fListHistos->Add(h) ; + } + List() ; } //____________________________________________________________________________ -AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) { - // ctor - AliPHOSGetter::Instance(inputfilename) ; +AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) : + TTask("GammaJet","Analysis of gamma-jet correlations") +{ + // ctor + fInputFileName = inputfilename; + fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName); + AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ; + fNEvent = gime->MaxEvent(); + InitParameters(); + fListHistos = 0 ; } //____________________________________________________________________________ -AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) { +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 ; + fOptionGJ = gj.fOptionGJ ; + fMinDistance = gj.fMinDistance ; + fOptFast = gj.fOptFast ; + fOnlyCharged = gj.fOnlyCharged ; + fOutputFile = gj.fOutputFile ; + fInputFileName = gj.fInputFileName ; + fOutputFileName = gj.fOutputFileName ; + fHIJINGFileName = gj.fHIJINGFileName ; + fHIJING = gj.fHIJING ; + fRatioMaxCut = gj.fRatioMaxCut ; + fRatioMinCut = gj.fRatioMinCut ; + fJetRatioMaxCut = gj.fJetRatioMaxCut ; + fJetRatioMinCut = gj.fJetRatioMinCut ; + fJetTPCRatioMaxCut = gj.fJetRatioMaxCut ; + fJetTPCRatioMinCut = gj.fJetRatioMinCut ; + fNEvent = gj.fNEvent ; + fNCone = gj.fNCone ; + fNPt = gj.fNPt ; + fResPara1 = gj.fResPara1 ; + fResPara2 = gj.fResPara2 ; + fResPara3 = gj.fResPara3 ; + fPtCut = gj.fPtCut ; + fNeutralPtCut = gj.fNeutralPtCut ; + fChargedPtCut = gj.fChargedPtCut ; + fPtJetSelectionCut = gj.fPtJetSelectionCut ; + fPhiMaxCut = gj.fPhiMaxCut ; + fPhiMinCut = gj.fPhiMinCut ; + fPosParaA = gj.fPosParaA ; + fPosParaB = gj.fPosParaB ; + fSelect = gj.fSelect ; + fTPCCutsLikeEMCAL = gj.fTPCCutsLikeEMCAL ; + fCone = gj.fCone ; + fPtThreshold = gj.fPtThreshold ; + + fDirName = gj.fDirName ; + fESDTree = gj.fESDTree ; + fPattern = gj.fPattern ; + + SetName (gj.GetName()) ; + SetTitle(gj.GetTitle()) ; + + for(Int_t i = 0; i<10; i++){ + fCones[i] = gj.fCones[i] ; + fNameCones[i] = gj.fNameCones[i] ; + fPtThres[i] = gj.fPtThres[i] ; + fNamePtThres[i] = gj.fNamePtThres[i] ; + if( i < 6 ){ + fJetXMin1[i] = gj.fJetXMin1[i] ; + fJetXMin2[i] = gj.fJetXMin2[i] ; + fJetXMax1[i] = gj.fJetXMax1[i] ; + fJetXMax2[i] = gj.fJetXMax2[i] ; + fBkgMean[i] = gj.fBkgMean[i] ; + fBkgRMS[i] = gj.fBkgRMS[i] ; + if( i < 2 ){ + fJetE1[i] = gj.fJetE1[i] ; + fJetE2[i] = gj.fJetE2[i] ; + fJetSigma1[i] = gj.fJetSigma1[i] ; + fJetSigma2[i] = gj.fJetSigma2[i] ; + fPhiEMCALCut[i] = gj.fPhiEMCALCut[i] ; + } + } + } } //____________________________________________________________________________ -AliPHOSGammaJet::~AliPHOSGammaJet() { +AliPHOSGammaJet::~AliPHOSGammaJet() +{ + fOutputFile->Close() ; + } //____________________________________________________________________________ -void AliPHOSGammaJet::Exec(Option_t * opt) +void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, + TClonesArray * plCh, + TClonesArray * plNe, + TClonesArray * plNePHOS) { - // does the job - - if (! opt) - return ; + + // List of particles copied to a root file. + +// Char_t sb[] = "bgrd/"; +// // cout<GetName()<GetName()<Get("T"); + TBranch *branch = t->GetBranch("primaries"); + branch->SetAddress(&particle); + + Int_t index = particleList->GetEntries() ; + Int_t indexNe = plNe->GetEntries() ; + Int_t indexCh = plCh->GetEntries() ; + Int_t indexNePHOS = plNePHOS->GetEntries() ; + Double_t charge = 0.; + Int_t iParticle = 0 ; + Int_t m = 0; + Double_t x = 0., z = 0.; + // cout<<"bkg entries "<GetEntries()<MaxEvent() ; - Int_t iEvent = 0 ; - for ( iEvent = 0 ; iEvent < nEvent1 ; iEvent++) { - if (iEvent <= 100 || iEvent%10 == 0) - Info("Exec", "Event %d", iEvent) ; - TParticle * particle = 0 ; - //-----------------Fill list with particles-------------------- + const AliPHOSGeometry * geom = gime->PHOSGeometry() ; + + if(!fOptFast){ + for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) { + t->GetEvent(iParticle) ; + + m = 0 ; + x = 0. ; + z = 0. ; + + charge = TDatabasePDG::Instance() + ->GetParticle(particle->GetPdgCode())->Charge(); + + if((charge != 0) && (particle->Pt() > fChargedPtCut)){ + if(TMath::Abs(particle->Eta())(particleList->At(index))) + ->SetStatusCode(0) ; + index++ ; + + new((*plCh)[indexCh]) TParticle(*particle) ; + (dynamic_cast(plCh->At(indexCh)))->SetStatusCode(0) ; + indexCh++ ; + + if(strstr(fOptionGJ,"deb all")||strstr(fOptionGJ,"deb")){ + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + } + } + else if((charge == 0) && (particle->Pt() > fNeutralPtCut) ){ + geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x); + + if(m != 0) + {//Is in PHOS + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + + new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ; + (dynamic_cast(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ; + indexNePHOS++ ; + } + + if((particle->Phi()>fPhiEMCALCut[0]) && + (particle->Phi()(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + + new((*particleList)[index]) TParticle(*particle) ; + (dynamic_cast(particleList->At(index))) + ->SetStatusCode(0) ; + index++ ; + new((*plNe)[indexNe]) TParticle(*particle) ; + (dynamic_cast(plNe->At(indexNe)))->SetStatusCode(0) ; + indexNe++ ; + } + } + } + } + } //No OptFast + else{ + + Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass(); TLorentzVector pPi0, pGamma1, pGamma2 ; - Double_t angle = 0 ; - TList particleList ; - Int_t n = -1; - gime->Event(iEvent, "X") ; - Int_t nparticles = gime->NPrimaries() ; - Int_t iParticle=0 ; - for (iParticle=0 ; iParticle < nparticles ; iParticle++) { - //Keep original partons - particle = gime->Primary(iParticle) ; - if((particle->GetStatusCode()== 21)){ - particleList.Add(particle); - n++; + Double_t angle = 0, cellDistance = 0.; + Bool_t p1 = kFALSE; + + // fFastRec = new AliPHOSFastGlobalReconstruction(fHIJINGFileName); + // fFastRec->FastReconstruction(iiEvent); + + for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) { + t->GetEvent(iParticle) ; + m = 0 ; + x = 0. ; + z = 0. ; + charge = TDatabasePDG::Instance() + ->GetParticle(particle->GetPdgCode())->Charge(); + + if((charge != 0) && (particle->Pt() > fChargedPtCut) + && (TMath::Abs(particle->Eta())(particleList->At(index))) + ->SetStatusCode(0) ; + index++ ; + + new((*plCh)[indexCh]) TParticle(*particle) ; + (dynamic_cast(plCh->At(indexCh)))->SetStatusCode(0) ; + indexCh++ ; } - //Keep Stable particles within eta range - Float_t etacut = 0. ; - if((particle->GetStatusCode() == 1)&& - (TMath::Abs(particle->Eta())GetPdgCode() != 111){ - particleList.Add(particle); - n++; - } - //Decay Pi0 and keep it with different status name - //Keep decay photons - if(particle->GetPdgCode() == 111) { - //cout<<"Pi0 "<Momentum(pPi0); - Pi0Decay(particle->GetMass(),pPi0,pGamma1,pGamma2,angle); - TParticle * photon1 = new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(), - pGamma1.Pz(),pGamma1.E(),0,0,0,0); - TParticle * photon2 = new TParticle(22,1,0,0,0,0,pGamma2.Px(),pGamma2.Py(), - pGamma2.Pz(),pGamma2.E(),0,0,0,0); - photon1->SetWeight(1); - photon2->SetWeight(2); - particle->SetStatusCode(2); - particleList.Add(particle); - particleList.Add(photon1); - particleList.Add(photon2); - //photon1->Print(); - //photon2->Print(); - }//if pi0 - }//final particle etacut + else if(charge == 0){ + geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x); + if((particle->GetPdgCode() != 111) && (particle->Pt() > fNeutralPtCut) && + (TMath::Abs(particle->Eta())Px(),particle->Py(), + particle->Pz(),particle->Energy()); + MakePhoton(part); + if(part.Pt() > fNeutralPtCut){ + + if(particle->Phi()>fPhiEMCALCut[0] && + particle->Phi()(particleList->At(index))) + ->SetStatusCode(0) ; + (dynamic_cast(particleList->At(index))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + index++ ; + + + new((*plNe)[indexNe]) TParticle(*particle) ; + (dynamic_cast(plNe->At(indexNe))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + (dynamic_cast(plNe->At(indexNe)))->SetStatusCode(0) ; + indexNe++ ; + } + if(m != 0) + { + new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ; + (dynamic_cast(plNePHOS->At(indexNePHOS))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + (dynamic_cast(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ; + indexNePHOS++ ; + } + } + } + if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) && + (TMath::Abs(particle->Eta())Px(),particle->Py(),particle->Pz(), + particle->Energy()); + + //Decay + + Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle); + //Check if decay photons are too close for PHOS + cellDistance = angle*460; //cm + if (cellDistance < fMinDistance) { + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + + //Pi0 inside phi EMCAL acceptance + + TLorentzVector part(particle->Px(),particle->Py(), + particle->Pz(),particle->Energy()); + MakePhoton(part); + if(part.Pt() > fNeutralPtCut){ + if(particle->Phi()>fPhiEMCALCut[0] && + particle->Phi()(particleList->At(index)))->SetStatusCode(0) ; + (dynamic_cast(particleList->At(index))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + index++ ; + + new((*plNe)[indexNe]) TParticle(*particle) ; + (dynamic_cast(plNe->At(indexNe))) ->SetStatusCode(0) ; + (dynamic_cast(plNe->At(indexNe))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + indexNe++ ; + } + if(m != 0){ + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ; + (dynamic_cast(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ; + (dynamic_cast(plNePHOS->At(indexNePHOS))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + indexNePHOS++; + }//In PHOS + } + }// if cell(fListHistos->FindObject("AnglePair")) + ->Fill(pPi0.E(),angle); + + p1 = kFALSE; + if(pGamma1.Pt() > 0. && TMath::Abs(pGamma1.Eta()) fNeutralPtCut ){ + + TParticle * photon1 = + new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(), + pGamma1.Pz(),pGamma1.E(),0,0,0,0); + 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) ; + (dynamic_cast(particleList->At(index)))->SetStatusCode(0) ; + index++ ; + + new((*plNe)[indexNe]) TParticle(*photon1) ; + (dynamic_cast(plNe->At(indexNe))) ->SetStatusCode(0) ; + indexNe++ ; + p1 = kTRUE; + } + if(m != 0){ + if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(photon1->Pt()); + new((*plNePHOS)[indexNePHOS]) TParticle(*photon1) ; + (dynamic_cast(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ; + indexNePHOS++; + p1 = kTRUE; + }//Is in PHOS + } + } + if(pGamma2.Pt() > 0. && TMath::Abs(pGamma2.Eta()) fNeutralPtCut){ + + TParticle * photon2 = + new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(), + pGamma2.Pz(),pGamma2.E(),0,0,0,0); + 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) ; + (dynamic_cast(particleList->At(index)))->SetStatusCode(0) ; + index++ ; + + new((*plNe)[indexNe]) TParticle(*photon2) ; + (dynamic_cast(plNe->At(indexNe))) ->SetStatusCode(0) ; + indexNe++ ; + } + if(m != 0){ + if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(photon2->Pt()); + new((*plNePHOS)[indexNePHOS]) TParticle(*photon2) ; + (dynamic_cast(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ; + indexNePHOS++; + } + if(p1){ + // e = (pGamma1+pGamma2).E(); + // if(IsAngleInWindow(angle,e)) + dynamic_cast + (fListHistos->FindObject("AnglePairAccepted"))-> + Fill(pPi0.E(),angle); + } + } + }//photon2 in acceptance + }//if angle > mindist + }//if pi0 + } + }//for (iParticlePHOSGeometry() ; + gime->Event(iEvent, "X") ; + + + Int_t index = particleList->GetEntries() ; + Int_t indexCh = plCh->GetEntries() ; + Int_t indexNe = plNe->GetEntries() ; + Int_t indexNePHOS = plNePHOS->GetEntries() ; + Int_t iParticle = 0 ; + Double_t charge = 0.; + Int_t m = 0; + Double_t x = 0., z = 0.; + if(!fOptFast){ + + for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) { + const TParticle * particle = gime->Primary(iParticle) ; + + m = 0 ; + x = 0. ; + z = 0. ; + + //Keep Stable particles within eta range + if((particle->GetStatusCode() == 1) && + (particle->Pt() > 0)){ + if(TMath::Abs(particle->Eta())GetParticle(particle->GetPdgCode())->Charge(); + if((charge != 0) && (particle->Pt() > fChargedPtCut)){ + + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + new((*plCh)[indexCh++]) TParticle(*particle) ; + new((*particleList)[index++]) TParticle(*particle) ; + } + else if((charge == 0) && (particle->Pt() > fNeutralPtCut)){ + geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x); + if(m != 0) + {//Is in PHOS + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + + new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ; + } + if((particle->Phi()>fPhiEMCALCut[0]) && + (particle->Phi()(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + new((*plNe)[indexNe++]) TParticle(*particle) ; + new((*particleList)[index++]) TParticle(*particle) ; + } + } + } + }//final particle, etacut }//for (iParticleGetParticle(111)->Mass(); + TLorentzVector pPi0, pGamma1, pGamma2 ; + Double_t angle = 0, cellDistance = 0.; + + fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName); + fFastRec->FastReconstruction(iEvent); + + for (iParticle=0 ; iParticle < gime->NPrimaries() ; iParticle++) { + const TParticle * particle = gime->Primary(iParticle) ; + m = 0 ; + x = 0. ; + z = 0. ; + //Keep Stable particles within eta range + if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){ + + //Fill lists + + charge = TDatabasePDG::Instance() + ->GetParticle(particle->GetPdgCode())->Charge(); + if((charge != 0) && (particle->Pt() > fChargedPtCut) && (TMath::Abs(particle->Eta())ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x); + if((particle->GetPdgCode() != 111) && particle->Pt() > 0 && + (TMath::Abs(particle->Eta())Px(),particle->Py(), + particle->Pz(),particle->Energy()); + + MakePhoton(part); + + if(part.Pt() > fNeutralPtCut){ + if(particle->Phi()>fPhiEMCALCut[0] && + particle->Phi()(particleList->At(index))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + index++ ; + + new((*plNe)[indexNe]) TParticle(*particle) ; + (dynamic_cast(plNe->At(indexNe))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + indexNe++ ; + } + if(m != 0) + { + new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ; + (dynamic_cast(plNePHOS->At(indexNePHOS))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + indexNePHOS++ ; + } + }// Small pt + } //No Pi0 + if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) && + (TMath::Abs(particle->Eta())Px(),particle->Py(),particle->Pz(), + particle->Energy()); + + //Decay + + Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle); + //Check if decay photons are too close for PHOS + cellDistance = angle*460; //cm + + if (cellDistance < fMinDistance) { + + //Pi0 inside phi EMCAL acceptance + + + TLorentzVector part(particle->Px(),particle->Py(), + particle->Pz(),particle->Energy()); + MakePhoton(part); + + if(part.Pt() > fNeutralPtCut){ + if(particle->Phi()>fPhiEMCALCut[0] && + particle->Phi()(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + + new((*plNe)[indexNe]) TParticle(*particle) ; + (dynamic_cast(plNe->At(indexNe))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + new((*particleList)[index]) TParticle(*particle) ; + (dynamic_cast(particleList->At(index))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + index++; + indexNe++; + }//InEMCAL + if(m != 0){ + if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(particle->Pt()); + new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ; + (dynamic_cast(plNePHOS->At(indexNePHOS))) + ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E()); + indexNePHOS++; + }//In PHOS + }//Small Pt + }// if cell(fListHistos->FindObject("AnglePair")) + ->Fill(pPi0.E(),angle); + + Bool_t p1 = kFALSE; + + if(pGamma1.Pt() > 0 && TMath::Abs(pGamma1.Eta()) fNeutralPtCut){ + TParticle * photon1 = + new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(), + pGamma1.Pz(),pGamma1.E(),0,0,0,0); + 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) ; + new((*particleList)[index++]) TParticle(*photon1) ; + //photon1->Print(); + p1 = kTRUE; + } + if(m != 0){ + if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(photon1->Pt()); + new((*plNePHOS)[indexNePHOS++]) TParticle(*photon1) ; + p1 = kTRUE; + } + } + } + if(pGamma2.Pt() > 0 && TMath::Abs(pGamma2.Eta()) fNeutralPtCut ){ + TParticle * photon2 = + new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(), + pGamma2.Pz(),pGamma2.E(),0,0,0,0); + 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(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb")) + dynamic_cast(fListHistos->FindObject("PtSpectra"))-> + Fill(photon2->Pt()); + new((*plNePHOS)[indexNePHOS++]) TParticle(*photon2) ; + } + + if(p1){ + // Float_t e = (pGamma1+pGamma2).E(); + // if(IsAngleInWindow(angle,e)) + dynamic_cast + (fListHistos->FindObject("AnglePairAccepted"))-> + Fill(pPi0.E(),angle); + } + } + }//photon2 in acceptance + }//if angle > mindist + }//if pi0 + }//If neutral + }//final particle, etacut + }//for (iParticleDelete() ; +} +//____________________________________________________________________________ +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 ; + Double_t pid[AliPID::kSPECIESN]; + + //########### PHOS ############## + //Info("CreateParticleListFromESD","Fill ESD PHOS list"); + Int_t begphos = esd->GetFirstPHOSParticle(); + Int_t endphos = esd->GetFirstPHOSParticle() + + esd->GetNumberOfPHOSParticles() ; + 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 + AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd + + //Create a TParticle to fill the particle list + + Double_t en = track->GetPHOSsignal() ; + Double_t * p = new Double_t(); + track->GetPHOSposition(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 + + track->GetPHOSpid(pid); + //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->GetFirstEMCALParticle(); + Int_t endem = esd->GetFirstEMCALParticle() + + esd->GetNumberOfEMCALParticles() ; + 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 + AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd + + Double_t en = track->GetEMCALsignal() ; + Double_t *p = new Double_t(); + track->GetEMCALposition(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"); + +} + + + +//____________________________________________________________________________ +void AliPHOSGammaJet::Exec(Option_t *option) +{ + // does the job + 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 + + } + + +// AliGenPythia* pyth = (AliGenPythia*) gAlice->Generator(); +// pyth->Init(); + + TClonesArray * particleList = new TClonesArray("TParticle",1000); + TClonesArray * plCh = new TClonesArray("TParticle",1000); + TClonesArray * plNe = new TClonesArray("TParticle",1000); + TClonesArray * plNePHOS = new TClonesArray("TParticle",1000); + + for (Int_t iEvent = 0 ; iEvent < fNEvent ; iEvent++) { + if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all")) + Info("Exec", "Event %d", iEvent) ; + + fRan.SetSeed(0); + + Double_t phig = 0., phil = 0., phich = 0 , phipi = 0; + Double_t etag = 0., etal = 0., etach = 0., etapi = 0. ; + Double_t ptg = 0., ptl = 0., ptch = 0., ptpi = 0.; + + TLorentzVector jet (0,0,0,0); + TLorentzVector jettpc(0,0,0,0); + + if(fESDdata){ + + 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); + } + + Bool_t iIsInPHOS = kFALSE ; + GetGammaJet(plNePHOS, ptg, phig, etag, iIsInPHOS) ; + + if(iIsInPHOS){ + + //Info("Exec"," In PHOS") ; + dynamic_cast(fListHistos->FindObject("NGamma"))->Fill(ptg); + dynamic_cast(fListHistos->FindObject("PhiGamma")) + ->Fill(ptg,phig); + dynamic_cast(fListHistos->FindObject("EtaGamma")) + ->Fill(ptg,etag); + if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all")) + Info("Exec", "Gamma: pt %f, phi %f, eta %f", ptg, + phig,etag) ; + +// cout<<"n charged "<GetEntries()<GetEntries()<GetEntries()<GetEntries()<GetEntries()<GetEntries()< fPhiEMCALCut[0] && + (phich + fCone) < fPhiEMCALCut[1]){ + insidech = kTRUE ; + } + Bool_t insidepi = kFALSE ; + if((phipi - fCone) > fPhiEMCALCut[0] && + (phipi + fCone) < fPhiEMCALCut[1]){ + insidepi = kTRUE ; + } + + if ((ptch > 0 || ptpi > 0)){ + if((ptch > ptpi) && insidech){ + phil = phich ; + etal = etach ; + ptl = ptch ; + dynamic_cast(fListHistos->FindObject("ChargeRatio")) + ->Fill(ptg,ptch/ptg); + dynamic_cast(fListHistos->FindObject("DeltaPhiCharge")) + ->Fill(ptg,phig-phich); + dynamic_cast(fListHistos->FindObject("DeltaEtaCharge")) + ->Fill(ptg,etag-etach); + if(strstr(fOptionGJ,"deb")) + Info("Exec"," Charged Leading") ; + } + if((ptpi > ptch) && insidepi){ + phil = phipi ; + etal = etapi ; + ptl = ptpi ; + + dynamic_cast(fListHistos->FindObject("Pi0Ratio")) + ->Fill(ptg,ptpi/ptg); + dynamic_cast(fListHistos->FindObject("DeltaPhiPi0")) + ->Fill(ptg,phig-phipi); + dynamic_cast(fListHistos->FindObject("DeltaEtaPi0")) + ->Fill(ptg,etag-etapi); + + if(ptpi > 0. && strstr(fOptionGJ,"deb")) + Info("Exec"," Pi0 Leading") ; + } + + 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(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", + jet.Phi(),jet.Eta(),jet.Pt()); + } +// dynamic_cast(fListHistos->FindObject("DeltaPhiJet")) +// ->Fill(ptg,pyjet.Phi()-jet.Phi()); +// dynamic_cast(fListHistos->FindObject("DeltaEtaJet")) +// ->Fill(ptg,pyjet.Eta()-jet.Eta()); +// dynamic_cast(fListHistos->FindObject("DeltaPtJet")) +// ->Fill(ptg,pyjet.Pt()-jet.Pt()); + } + else + MakeJetAnyConeOrPt(particleList, ptg, phig, ptl, phil, etal, ""); + } + + //TPC + if(fOnlyCharged && ptch > 0.) + { + if(strstr(fOptionGJ,"deb")) + Info("Exec","Leading TPC pt %f, phi %f",ptch,phich); + + dynamic_cast(fListHistos->FindObject("TPCRatio")) + ->Fill(ptg,ptch/ptg); + dynamic_cast(fListHistos->FindObject("DeltaPhiTPC")) + ->Fill(ptg,phig-phich); + dynamic_cast(fListHistos->FindObject("DeltaEtaTPC")) + ->Fill(ptg,etag-etach); + + if(!fAnyConeOrPt){ + + MakeJet(plCh, ptg, phig, ptch, phich, etach, "TPC",jettpc); + + 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", + jettpc.Phi(),jettpc.Eta(),jettpc.Pt()); + } +// dynamic_cast(fListHistos->FindObject("DeltaPhiTPCJet")) +// ->Fill(ptg,pyjet.Phi()-jettpc.Phi()); +// dynamic_cast(fListHistos->FindObject("DeltaEtaTPCJet")) +// ->Fill(ptg,pyjet.Eta()-jettpc.Eta()); +// dynamic_cast(fListHistos->FindObject("DeltaPtTPCJet")) +// ->Fill(ptg,pyjet.Pt()-jettpc.Pt()); + } + else + MakeJetAnyConeOrPt(plCh, ptg, phig, ptch, phich, etach, "TPC"); + + } + } + } + + particleList->Delete() ; + plCh->Delete() ; + plNe->Delete() ; + plNePHOS->Delete() ; }//loop: events + + delete plNe ; + delete plCh ; + delete particleList ; + + fOutputFile->Write() ; + fOutputFile->cd(); + this->Write(); } +//____________________________________________________________________________ +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; + + TIter next(pl) ; + while ( (particle = dynamic_cast(next())) ) { + ipr++ ; + Double_t pt = particle->Pt(); + + charge = TDatabasePDG::Instance() + ->GetParticle(particle->GetPdgCode())->Charge(); + if(charge != 0){//Only jet Charged particles + dynamic_cast + (fListHistos->FindObject(type+conf+"Fragment")) + ->Fill(ptg,pt/ptg); + dynamic_cast + (fListHistos->FindObject(type+conf+"PtDist")) + ->Fill(ptg,pt); + } + } + if(type == "Bkg"){ + dynamic_cast + (fListHistos->FindObject("NBkg"+conf))->Fill(ipr); + } +} +//____________________________________________________________________________ +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; + + TIter next(pl) ; + while ( (particle = dynamic_cast(next())) ) { + ipr++; + Double_t pt = particle->Pt(); + charge = TDatabasePDG::Instance() + ->GetParticle(particle->GetPdgCode())->Charge(); + if(charge != 0){//Only jet Charged particles + dynamic_cast + (fListHistos->FindObject(type+conf+"FragmentCone"+cone+"Pt"+ptcut)) + ->Fill(ptg,pt/ptg); + dynamic_cast + (fListHistos->FindObject(type+conf+"PtDistCone"+cone+"Pt"+ptcut)) + ->Fill(ptg,pt); + } + }//while + + if(type == "Bkg"){ + dynamic_cast + (fListHistos->FindObject("NBkg"+conf+"Cone"+cone+"Pt"+ptcut)) + ->Fill(ipr); + } +} + +//____________________________________________________________________________ +void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt, + 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.; + + for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){ + TParticle * particle = dynamic_cast(pl->At(ipr)) ; + + if((particle->Pt() > fPtCut) && (particle->Pt() > pt)){ + + pt = particle->Pt(); + phi = particle->Phi() ; + eta = particle->Eta() ; + Is = kTRUE; + } + } +} + +//____________________________________________________________________________ +void AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl, + Double_t ptg, Double_t phig, + 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; + + for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){ + + TParticle * particle = dynamic_cast(pl->At(ipr)) ; + + Double_t ptl = particle->Pt(); + Double_t rat = ptl/ptg ; + Double_t phil = particle->Phi() ; + + if(((phig-phil)> fPhiMinCut) && ((phig-phil) fRatioMinCut) && (rat < fRatioMaxCut) && (ptl > pt)) { + eta = particle->Eta() ; + phi = phil ; + pt = ptl ; + //printf("GetLeadingCharge: %f %f %f %f \n", pt, eta, phi,rat) ; + } + } + //printf("GetLeadingCharge: %f %f %f \n", pt, eta, phi) ; + +} + + +//____________________________________________________________________________ +void AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl, + Double_t ptg, Double_t phig, + 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.; + Double_t ptl = -100.; + Double_t rat = -100.; + Double_t phil = -100. ; + + TIter next(pl); + TParticle * particle = 0; + Float_t ef = 0; + if(!fOptFast){ + Float_t e = 0; + while ( (particle = (TParticle*)next()) ) { + if( particle->GetPdgCode() == 111){ + ptl = particle->Pt(); + rat = ptl/ptg ; + phil = particle->Phi() ; + e = particle->Energy(); + dynamic_cast + (fListHistos->FindObject("AnglePairNoCut"))-> + Fill(e,0.1); + dynamic_cast + (fListHistos->FindObject("InvMassPairNoCut"))-> + Fill(ptg,1.35); + + if(((phig-phil)> fPhiMinCut) && ((phig-phil) fRatioMinCut) && (rat < fRatioMaxCut)) { + + dynamic_cast + (fListHistos->FindObject("AnglePairLeadingCut"))-> + Fill(e,0.1); + dynamic_cast + (fListHistos->FindObject("InvMassPairLeadingCut"))-> + Fill(ptg,1.35); + + dynamic_cast + (fListHistos->FindObject("AnglePairAngleCut"))-> + Fill(e,0.15); + dynamic_cast + (fListHistos->FindObject("InvMassPairAngleCut"))-> + Fill(ptg,1.36); + + dynamic_cast + (fListHistos->FindObject("InvMassPairAllCut"))-> + Fill(ptg,0.27); + dynamic_cast + (fListHistos->FindObject("AnglePairAllCut"))-> + Fill(e,1.34); + + + if(ptl > pt){ + eta = particle->Eta() ; + phi = phil ; + pt = ptl ; + ef = e; + //printf("GetLeadingPi0: %f %f %f %f %f \n", pt, eta, phi, rat, ptg) ; + } + } + } + + dynamic_cast + (fListHistos->FindObject("InvMassPairLeading"))-> + Fill(ptg,1.35); + dynamic_cast + (fListHistos->FindObject("AnglePairLeading"))-> + Fill(ef,0.1); + } + }//No fOptfast + else{ + Int_t iPrimary = -1; + TLorentzVector gammai,gammaj; + Double_t angle = 0., e = 0., invmass = 0.; + Double_t anglef = 0., ef = 0., invmassf = 0.; + Int_t ksPdg = 0; + Int_t jPrimary=-1; + + while ( (particle = (TParticle*)next()) ) { + iPrimary++; + + ksPdg = particle->GetPdgCode(); + ptl = particle->Pt(); + if(ksPdg == 111){ //2 gamma + rat = ptl/ptg ; + phil = particle->Phi() ; + if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) && + ((phig-phil)>fPhiMinCut)&&((phig-phil)Eta() ; + phi = phil ; + pt = ptl ; + }// cuts + }// pdg = 111 + if(ksPdg == 22){//1 gamma + particle->Momentum(gammai); + jPrimary=-1; + TIter next2(pl); + while ( (particle = (TParticle*)next2()) ) { + jPrimary++; + if(jPrimary>iPrimary){ + ksPdg = particle->GetPdgCode(); + if(ksPdg == 22){ + particle->Momentum(gammaj); + //Info("GetLeadingPi0","Egammai %f, Egammaj %f", + //gammai.Pt(),gammaj.Pt()); + + ptl = (gammai+gammaj).Pt(); + phil = (gammai+gammaj).Phi(); + if(phil < 0) + phil+=TMath::TwoPi(); + rat = ptl/ptg ; + invmass = (gammai+gammaj).M(); + angle = gammaj.Angle(gammai.Vect()); + //Info("GetLeadingPi0","Angle %f", angle); + e = (gammai+gammaj).E(); + + dynamic_cast + (fListHistos->FindObject("AnglePairNoCut"))-> + Fill(e,angle); + dynamic_cast + (fListHistos->FindObject("InvMassPairNoCut"))-> + Fill(ptg,invmass); + + if((rat > fRatioMinCut) && (rat < fRatioMaxCut) && + ((phig-phil)>fPhiMinCut)&&((phig-phil) + (fListHistos->FindObject("AnglePairLeadingCut"))-> + Fill(e,angle); + dynamic_cast + (fListHistos->FindObject("InvMassPairLeadingCut"))-> + Fill(ptg,invmass); + + if(IsAngleInWindow(angle,e)){ + dynamic_cast + (fListHistos->FindObject("AnglePairAngleCut"))-> + Fill(e,angle); + dynamic_cast + (fListHistos->FindObject("InvMassPairAngleCut"))-> + Fill(ptg,invmass); + + //Info("GetLeadingPi0","InvMass %f", invmass); + if((invmass>fInvMassMinCut) && (invmass + (fListHistos->FindObject("InvMassPairAllCut"))-> + Fill(ptg,invmass); + dynamic_cast + (fListHistos->FindObject("AnglePairAllCut"))-> + Fill(e,angle); + if(ptl > pt ){ + pt = ptl; + eta = particle->Eta() ; + phi = phil ; + ef = e ; + anglef = angle ; + invmassf = invmass ; + + } + }//cuts + }//(invmass>0.125) && (invmass<0.145) + }//gammaj.Angle(gammai.Vect())<0.04 + }//if pdg = 22 + }//iprimary 0.){ + dynamic_cast + (fListHistos->FindObject("InvMassPairLeading"))-> + Fill(ptg,invmassf); + dynamic_cast + (fListHistos->FindObject("AnglePairLeading"))-> + Fill(ef,anglef); + } + }//fOptFast + // printf("GetLeadingPi0: %f %f %f \n", pt, eta, phi) ; +} + + +//____________________________________________________________________________ +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) ; + fAngleMaxParam.AddAt(0.025,2) ; + fAngleMaxParam.AddAt(-2e-4,3) ; + fAnyConeOrPt = kFALSE ; + fOutputFileName = "GammaJet.root" ; + fOptionGJ = ""; + fHIJINGFileName = "galice.root" ; + fHIJING = kFALSE ; + fMinDistance = 3.6 ; + fEtaCut = 0.7 ; + fInvMassMaxCut = 0.15 ; + fInvMassMinCut = 0.12 ; + fOnlyCharged = kFALSE ; + fOptFast = kFALSE ; + fESDdata = kTRUE ; + fPhiEMCALCut[0] = 60. *TMath::Pi()/180.; + fPhiEMCALCut[1] = 180.*TMath::Pi()/180.; + fPhiMaxCut = 3.4 ; + fPhiMinCut = 2.9 ; + fPtCut = 10. ; + fNeutralPtCut = 0.5 ; + fChargedPtCut = 0.5 ; + fTPCCutsLikeEMCAL = kFALSE ; + //Jet selection parameters + //Fixed cut (old) + fRatioMaxCut = 1.0 ; + fRatioMinCut = 0.1 ; + fJetRatioMaxCut = 1.2 ; + fJetRatioMinCut = 0.8 ; + fJetTPCRatioMaxCut = 1.2 ; + fJetTPCRatioMinCut = 0.3 ; + fSelect = kFALSE ; + + fDirName = "./" ; + fESDTree = "esdTree" ; + fPattern = "." ; + + //Cut depending on gamma energy + + fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applyed + //Reconstructed jet energy dependence parameters + //e_jet = a1+e_gamma b2. + //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3 + fJetE1[0] = -5.75; fJetE1[1] = -4.1; + fJetE2[0] = 1.005; fJetE2[1] = 1.05; + + //Reconstructed sigma of jet energy dependence parameters + //s_jet = a1+e_gamma b2. + //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3 + fJetSigma1[0] = 2.65; fJetSigma1[1] = 2.75; + fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033; + + //Background mean energy and RMS + //Index 0-> No BKG; Index 1-> BKG > 2 GeV; + //Index 2-> (low pt jets)BKG > 0.5 GeV; + //Index > 2, same for TPC conf + fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5; + fBkgMean[3] = 0.; fBkgMean[4] = 6.4; fBkgMean[5] = 48.6; + fBkgRMS[0] = 0.; fBkgRMS[1] = 7.5; fBkgRMS[2] = 22.0; + fBkgRMS[3] = 0.; fBkgRMS[4] = 5.4; fBkgRMS[5] = 13.2; + + //Factor x of min/max = E -+ x * sigma. Obtained after selecting the + //limits for monoenergetic jets. + //Index 0-> No BKG; Index 1-> BKG > 2 GeV; + //Index 2-> (low pt jets) BKG > 0.5 GeV; + //Index > 2, same for TPC conf + + fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ; + fJetXMin1[3] =-2.0 ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1 ; + fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034; + fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ; + fJetXMax1[0] =-3.8 ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6 ; + fJetXMax1[3] =-2.7 ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7 ; + fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016; + fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027; + + + //Photon fast reconstruction + fResPara1 = 0.0255 ; // GeV + fResPara2 = 0.0272 ; + fResPara3 = 0.0129 ; + + fPosParaA = 0.096 ; // cm + fPosParaB = 0.229 ; + + //Different cones and pt thresholds to construct the jet + + fCone = 0.3 ; + fPtThreshold = 0. ; + fNCone = 1 ; + fNPt = 1 ; + fCones[0] = 0.3 ; fNameCones[0] = "03" ; + fPtThres[0] = 0.5 ; fNamePtThres[0] = "05" ; + +} + +//__________________________________________________________________________- +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) + +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e; + Double_t arg = (e*e-2*mpi0*mpi0)/(e*e); + Double_t min = 100. ; + if(arg>0.) + min = TMath::ACos(arg); + + if((angle=min)) + result = kTRUE; + + return result; +} + +//__________________________________________________________________________- +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]; + Double_t xmin[2]; + + Int_t iTPC = 0; + + if(type == "TPC" && !fTPCCutsLikeEMCAL){ + iTPC = 3 ;//If(fTPCCutsLikeEMCAL) take jet energy cuts like EMCAL + } + + + if(!fHIJING){ + //Phythia alone, jets with pt_th > 0.2, r = 0.3 + par[0] = fJetE1[0]; par[1] = fJetE2[0]; + //Energy of the jet peak + //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit + par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0]; + //Sigma of the jet peak + //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit + par[4] = fBkgMean[0 + iTPC]; par[5] = fBkgRMS[0 + iTPC]; + //Parameters reserved for HIJING bkg. + xmax[0] = fJetXMax1[0 + iTPC]; xmax[1] = fJetXMax2[0 + iTPC]; + xmin[0] = fJetXMin1[0 + iTPC]; xmin[1] = fJetXMin2[0 + iTPC]; + //Factor that multiplies sigma to obtain the best limits, + //by observation, of mono jet ratios (ptjet/ptg) + //X_jet = fJetX1[0]+fJetX2[0]*e_gamma + + } + else{ + if(ptg > fPtJetSelectionCut){ + //Phythia +HIJING with pt_th > 2 GeV/c, r = 0.3 + par[0] = fJetE1[0]; par[1] = fJetE2[0]; + //Energy of the jet peak, same as in pp + //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit + par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0]; + //Sigma of the jet peak, same as in pp + //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit + par[4] = fBkgMean[1 + iTPC]; par[5] = fBkgRMS[1 + iTPC]; + //Mean value and RMS of HIJING Bkg + xmax[0] = fJetXMax1[1 + iTPC]; xmax[1] = fJetXMax2[1 + iTPC]; + xmin[0] = fJetXMin1[1 + iTPC]; xmin[1] = fJetXMin2[1 + iTPC]; + //Factor that multiplies sigma to obtain the best limits, + //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg, + //pt_th > 2 GeV, r = 0.3 + //X_jet = fJetX1[0]+fJetX2[0]*e_gamma + + } + else{ + //Phythia + HIJING with pt_th > 0.5 GeV/c, r = 0.3 + par[0] = fJetE1[1]; par[1] = fJetE2[1]; + //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3 + //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit + par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1]; + //Sigma of the jet peak, pt_th > 2 GeV/c, r = 0.3 + //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit + par[4] = fBkgMean[2 + iTPC]; par[5] = fBkgRMS[2 + iTPC]; + //Mean value and RMS of HIJING Bkg in a 0.3 cone, pt > 2 GeV. + xmax[0] = fJetXMax1[2 + iTPC]; xmax[1] = fJetXMax2[2 + iTPC]; + xmin[0] = fJetXMin1[2 + iTPC]; xmin[1] = fJetXMin2[2 + iTPC]; + //Factor that multiplies sigma to obtain the best limits, + //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg, + //pt_th > 2 GeV, r = 0.3 + //X_jet = fJetX1[0]+fJetX2[0]*e_gamma + + }//If low pt jet in bkg + }//if Bkg + + //Calculate minimum and maximum limits of the jet ratio. + Double_t min = CalculateJetRatioLimit(ptg, par, xmin); + Double_t max = CalculateJetRatioLimit(ptg, par, xmax); + + //Info("IsJetSeleted","%s : Limits min %f, max %f, ptg / ptj %f", + // type.Data(),min,max,ptj/ptg); + if(( min < ptj/ptg ) && ( max > ptj/ptg)) + return kTRUE; + else + return kFALSE; + +} + +//____________________________________________________________________________ +void AliPHOSGammaJet::List() const +{ + // List the histos + + Info("List", "%d histograms found", fListHistos->GetEntries() ) ; + TIter next(fListHistos) ; + TH2F * h = 0 ; + while ( (h = dynamic_cast(next())) ) + Info("List", "%s", h->GetName()) ; +} + +//____________________________________________________________________________ +Double_t AliPHOSGammaJet::MakeEnergy(const Double_t energy) +{ + // Smears the energy according to the energy dependent energy resolution. + // A gaussian distribution is assumed + + Double_t sigma = SigmaE(energy) ; + return fRan.Gaus(energy, sigma) ; + + +} +//____________________________________________________________________________ +void AliPHOSGammaJet::MakeHistos() +{ + // Create histograms to be saved in output file and + // stores them in a TObjectArray + + fOutputFile = new TFile(fOutputFileName, "recreate") ; + + fListHistos = new TObjArray(10000) ; + + // Histos gamma pt vs leading pt + + TH1F * hPtSpectra = new TH1F + ("PtSpectra","p_{T i} vs p_{T #gamma}",200,0,200); + hPtSpectra->SetXTitle("p_{T} (GeV/c)"); + fListHistos->Add(hPtSpectra) ; + + //Histos ratio charged leading pt / gamma pt vs pt + + TH2F * hChargeRatio = new TH2F + ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}", + 120,0,120,120,0,1); + hChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}"); + hChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hChargeRatio) ; + + TH2F * hTPCRatio = new TH2F + ("TPCRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}", + 120,0,120,120,0,1); + hTPCRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}"); + hTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hTPCRatio) ; + + + TH2F * hPi0Ratio = new TH2F + ("Pi0Ratio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}", + 120,0,120,120,0,1); + hPi0Ratio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}"); + hPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hPi0Ratio) ; + + TH2F * hPhiGamma = new TH2F + ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7); + hPhiGamma->SetYTitle("#phi"); + hPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hPhiGamma) ; + + TH2F * hEtaGamma = new TH2F + ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8); + hEtaGamma->SetYTitle("#eta"); + hEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hEtaGamma) ; + +// //Jet reconstruction check +// TH2F * hDeltaPhiJet = new TH2F +// ("DeltaPhiJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}", +// 200,0,120,200,-1.5,1.5); +// hDeltaPhiJet->SetYTitle("#Delta #phi"); +// hDeltaPhiJet->SetXTitle("p_{T #gamma} (GeV/c)"); +// fListHistos->Add(hDeltaPhiJet) ; + +// TH2F * hDeltaPhiTPCJet = new TH2F +// ("DeltaPhiTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}", +// 200,0,120,200,-1.5,1.5); +// hDeltaPhiTPCJet->SetYTitle("#Delta #phi"); +// hDeltaPhiTPCJet->SetXTitle("p_{T #gamma} (GeV/c)"); +// fListHistos->Add(hDeltaPhiTPCJet) ; + +// TH2F * hDeltaEtaJet = new TH2F +// ("DeltaEtaJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}", +// 200,0,120,200,-1.5,1.5); +// hDeltaEtaJet->SetYTitle("#Delta #phi"); +// hDeltaEtaJet->SetXTitle("p_{T #gamma} (GeV/c)"); +// fListHistos->Add(hDeltaEtaJet) ; + +// TH2F * hDeltaEtaTPCJet = new TH2F +// ("DeltaEtaTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}", +// 200,0,120,200,-1.5,1.5); +// hDeltaEtaTPCJet->SetYTitle("#Delta #phi"); +// hDeltaEtaTPCJet->SetXTitle("p_{T #gamma} (GeV/c)"); +// fListHistos->Add(hDeltaEtaTPCJet) ; + +// TH2F * hDeltaPtJet = new TH2F +// ("DeltaPtJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}", +// 200,0,120,200,0.,100.); +// hDeltaPtJet->SetYTitle("#Delta #phi"); +// hDeltaPtJet->SetXTitle("p_{T #gamma} (GeV/c)"); +// fListHistos->Add(hDeltaPtJet) ; + +// TH2F * hDeltaPtTPCJet = new TH2F +// ("DeltaPtTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}", +// 200,0,120,200,0.,100.); +// hDeltaPtTPCJet->SetYTitle("#Delta #phi"); +// hDeltaPtTPCJet->SetXTitle("p_{T #gamma} (GeV/c)"); +// fListHistos->Add(hDeltaPtTPCJet) ; + + // + TH2F * hDeltaPhiCharge = new TH2F + ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}", + 200,0,120,200,0,6.4); + hDeltaPhiCharge->SetYTitle("#Delta #phi"); + hDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hDeltaPhiCharge) ; + + TH2F * hDeltaPhiTPC = new TH2F + ("DeltaPhiTPC","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}", + 200,0,120,200,0,6.4); + hDeltaPhiTPC->SetYTitle("#Delta #phi"); + hDeltaPhiTPC->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hDeltaPhiTPC) ; + TH2F * hDeltaPhiPi0 = new TH2F + ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}", + 200,0,120,200,0,6.4); + hDeltaPhiPi0->SetYTitle("#Delta #phi"); + hDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hDeltaPhiPi0) ; + + TH2F * hDeltaEtaCharge = new TH2F + ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}", + 200,0,120,200,-2,2); + hDeltaEtaCharge->SetYTitle("#Delta #eta"); + hDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hDeltaEtaCharge) ; + + TH2F * hDeltaEtaTPC = new TH2F + ("DeltaEtaTPC","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}", + 200,0,120,200,-2,2); + hDeltaEtaTPC->SetYTitle("#Delta #eta"); + hDeltaEtaTPC->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hDeltaEtaTPC) ; + + TH2F * hDeltaEtaPi0 = new TH2F + ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}", + 200,0,120,200,-2,2); + hDeltaEtaPi0->SetYTitle("#Delta #eta"); + hDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hDeltaEtaPi0) ; + + if(fOptFast){ + + TH2F * hAnglePair = new TH2F + ("AnglePair", + "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}", + 200,0,50,200,0,0.2); + hAnglePair->SetYTitle("Angle (rad)"); + hAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)"); + fListHistos->Add(hAnglePair) ; + + TH2F * hAnglePairAccepted = new TH2F + ("AnglePairAccepted", + "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window", + 200,0,50,200,0,0.2); + hAnglePairAccepted->SetYTitle("Angle (rad)"); + hAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)"); + fListHistos->Add(hAnglePairAccepted) ; + + TH2F * hAnglePairNoCut = new TH2F + ("AnglePairNoCut", + "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2); + hAnglePairNoCut->SetYTitle("Angle (rad)"); + hAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)"); + fListHistos->Add(hAnglePairNoCut) ; + + TH2F * hAnglePairLeadingCut = new TH2F + ("AnglePairLeadingCut", + "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}", + 200,0,50,200,0,0.2); + hAnglePairLeadingCut->SetYTitle("Angle (rad)"); + hAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)"); + fListHistos->Add(hAnglePairLeadingCut) ; + + TH2F * hAnglePairAngleCut = new TH2F + ("AnglePairAngleCut", + "Angle between all #gamma pair (angle + leading cut) vs p_{T #pi^{0}}" + ,200,0,50,200,0,0.2); + hAnglePairAngleCut->SetYTitle("Angle (rad)"); + hAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)"); + fListHistos->Add(hAnglePairAngleCut) ; + + TH2F * hAnglePairAllCut = new TH2F + ("AnglePairAllCut", + "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T #pi^{0}}" + ,200,0,50,200,0,0.2); + hAnglePairAllCut->SetYTitle("Angle (rad)"); + hAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)"); + fListHistos->Add(hAnglePairAllCut) ; + + TH2F * hAnglePairLeading = new TH2F + ("AnglePairLeading", + "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}", + 200,0,50,200,0,0.2); + hAnglePairLeading->SetYTitle("Angle (rad)"); + hAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)"); + fListHistos->Add(hAnglePairLeading) ; + + + TH2F * hInvMassPairNoCut = new TH2F + ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}", + 120,0,120,360,0,0.5); + hInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})"); + hInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hInvMassPairNoCut) ; + + TH2F * hInvMassPairLeadingCut = new TH2F + ("InvMassPairLeadingCut", + "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}", + 120,0,120,360,0,0.5); + hInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})"); + hInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hInvMassPairLeadingCut) ; + + TH2F * hInvMassPairAngleCut = new TH2F + ("InvMassPairAngleCut", + "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}", + 120,0,120,360,0,0.5); + hInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})"); + hInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hInvMassPairAngleCut) ; + + + TH2F * hInvMassPairAllCut = new TH2F + ("InvMassPairAllCut", + "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}", + 120,0,120,360,0,0.5); + hInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})"); + hInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hInvMassPairAllCut) ; + + TH2F * hInvMassPairLeading = new TH2F + ("InvMassPairLeading", + "Invariant Mass of #gamma pair selected vs p_{T #gamma}", + 120,0,120,360,0,0.5); + hInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})"); + hInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hInvMassPairLeading) ; + } + + //Count + + TH1F * hNGamma = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120); + hNGamma->SetYTitle("N"); + hNGamma->SetXTitle("p_{T #gamma}(GeV/c)"); + fListHistos->Add(hNGamma) ; + + TH1F * hNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000); + hNBkg->SetYTitle("counts"); + hNBkg->SetXTitle("N"); + fListHistos->Add(hNBkg) ; + + TH2F * hNLeading = new TH2F + ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120); + hNLeading->SetYTitle("p_{T charge} (GeV/c)"); + hNLeading->SetXTitle("p_{T #gamma}(GeV/c)"); + fListHistos->Add(hNLeading) ; + + + TH1F * hN = new TH1F("NJet","Accepted jets",240,0,120); + hN->SetYTitle("N"); + hN->SetXTitle("p_{T #gamma}(GeV/c)"); + fListHistos->Add(hN) ; + + + //Ratios and Pt dist of reconstructed (not selected) jets + //Jet + TH2F * hJetRatio = new TH2F + ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}", + 240,0,120,200,0,10); + hJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); + hJetRatio->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetRatio) ; + + + TH2F * hJetPt = new TH2F + ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200); + hJetPt->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); + hJetPt->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetPt) ; + + + //Bkg + + TH2F * hBkgRatio = new TH2F + ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}", + 240,0,120,200,0,10); + hBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}"); + hBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgRatio) ; + + + TH2F * hBkgPt = new TH2F + ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200); + hBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}"); + hBkgPt->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgPt) ; + + + //Jet Distributions + + TH2F * hJetFragment = + new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}", + 240,0.,120.,1000,0.,1.2); + hJetFragment->SetYTitle("x_{T}"); + hJetFragment->SetXTitle("p_{T #gamma}"); + fListHistos->Add(hJetFragment) ; + + TH2F * hBkgFragment = new TH2F + ("BkgFragment","x = p_{T i charged}/p_{T #gamma}", + 240,0.,120.,1000,0.,1.2); + hBkgFragment->SetYTitle("x_{T}"); + hBkgFragment->SetXTitle("p_{T #gamma}"); + fListHistos->Add(hBkgFragment) ; + + TH2F * hJetPtDist = + new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); + hJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetPtDist) ; + + TH2F * hBkgPtDist = new TH2F + ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); + hBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgPtDist) ; + + + if(fOnlyCharged){ + //Counts + TH1F * hNBkgTPC = new TH1F + ("NBkgTPC","TPC bkg multiplicity ",9000,0,9000); + hNBkgTPC->SetYTitle("counts"); + hNBkgTPC->SetXTitle("N"); + fListHistos->Add(hNBkgTPC) ; + + TH2F * hNTPCLeading = new TH2F + ("NTPCLeading","Accepted TPC jet leading",240,0,120,240,0,120); + hNTPCLeading->SetYTitle("p_{T charge} (GeV/c)"); + hNTPCLeading->SetXTitle("p_{T #gamma}(GeV/c)"); + fListHistos->Add(hNTPCLeading) ; + + TH1F * hNTPC = new TH1F("NTPCJet","Number of TPC jets",240,0,120); + hNTPC->SetYTitle("N"); + hNTPC->SetXTitle("p_{T #gamma}(GeV/c)"); + fListHistos->Add(hNTPC) ; + + TH2F * hJetTPCRatio = new TH2F + ("JetTPCRatio", "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}", + 240,0,120,200,0,10); + hJetTPCRatio->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}"); + hJetTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetTPCRatio) ; + + TH2F * hBkgTPCRatio = new TH2F + ("BkgTPCRatio","p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}", + 240,0,120,200,0,10); + hBkgTPCRatio->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}"); + hBkgTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgTPCRatio) ; + + TH2F * hJetTPCPt = new TH2F + ("JetTPCPt", "p_{T jet lead TPC} vs p_{T #gamma}",240,0,120,400,0,200); + hJetTPCPt->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}"); + hJetTPCPt->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetTPCPt) ; + + TH2F * hBkgTPCPt = new TH2F + ("BkgTPCPt", "p_{T bkg lead TPC} vs p_{T #gamma}",240,0,120,400,0,200); + hBkgTPCPt->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}"); + hBkgTPCPt->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgTPCPt) ; + + //JetDistributions + + TH2F * hJetTPCFragment = + new TH2F("JetTPCFragment","x = p_{T i charged}/p_{T #gamma}", + 240,0.,120.,1000,0.,1.2); + hJetTPCFragment->SetYTitle("x_{T}"); + hJetTPCFragment->SetXTitle("p_{T #gamma}"); + fListHistos->Add(hJetTPCFragment) ; + + TH2F * hBkgTPCFragment = new TH2F + ("BkgTPCFragment","x = p_{T i charged}/p_{T #gamma}", + 240,0.,120.,1000,0.,1.2); + hBkgTPCFragment->SetYTitle("x_{T}"); + hBkgTPCFragment->SetXTitle("p_{T #gamma}"); + fListHistos->Add(hBkgTPCFragment) ; + + + TH2F * hJetTPCPtDist = new TH2F("JetTPCPtDist", + "x = p_{T i charged}",240,0.,120.,400,0.,200.); + hJetTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetTPCPtDist) ; + + TH2F * hBkgTPCPtDist = new TH2F + ("BkgTPCPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); + hBkgTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgTPCPtDist) ; + + } + + + if(fAnyConeOrPt){ + //If we want to study the jet for different cones and pt. Old version + + TH2F * hJetRatios[5][5]; + TH2F * hJetTPCRatios[5][5]; + + TH2F * hJetPts[5][5]; + TH2F * hJetTPCPts[5][5]; + + TH2F * hBkgRatios[5][5]; + TH2F * hBkgTPCRatios[5][5]; + + TH2F * hBkgPts[5][5]; + TH2F * hBkgTPCPts[5][5]; + + TH2F * hNLeadings[5][5]; + TH2F * hNTPCLeadings[5][5]; + + TH1F * hNs[5][5]; + TH1F * hNTPCs[5][5]; + + TH1F * hNBkgs[5][5]; + TH1F * hNBkgTPCs[5][5]; + + TH2F * hJetFragments[5][5]; + TH2F * hBkgFragments[5][5]; + TH2F * hJetPtDists[5][5]; + TH2F * hBkgPtDists[5][5]; + + TH2F * hJetTPCFragments[5][5]; + TH2F * hBkgTPCFragments[5][5]; + TH2F * hJetTPCPtDists[5][5]; + TH2F * hBkgTPCPtDists[5][5]; + + + for(Int_t icone = 0; icone" +fNamePtThres[ipt]+" GeV/c", + 240,0,120,200,0,10); + hJetRatios[icone][ipt]-> + SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); + hJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetRatios[icone][ipt]) ; + + + hJetPts[icone][ipt] = new TH2F + ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" + +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", + 240,0,120,400,0,200); + hJetPts[icone][ipt]-> + SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); + hJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetPts[icone][ipt]) ; + + + //Bkg + hBkgRatios[icone][ipt] = new TH2F + ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" + +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", + 240,0,120,200,0,10); + hBkgRatios[icone][ipt]-> + SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}"); + hBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgRatios[icone][ipt]) ; + + + + hBkgPts[icone][ipt] = new TH2F + ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" + +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", + 240,0,120,400,0,200); + hBkgPts[icone][ipt]-> + SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); + hBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgPts[icone][ipt]) ; + + + if(fOnlyCharged){ + + hJetTPCRatios[icone][ipt] = new TH2F + ("JetTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone =" + +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", + 240,0,120,200,0,10); + hJetTPCRatios[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}"); + hJetTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetTPCRatios[icone][ipt]) ; + + hBkgTPCRatios[icone][ipt] = new TH2F + ("BkgTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone =" + +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", + 240,0,120,200,0,10); + hBkgTPCRatios[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}"); + hBkgTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgTPCRatios[icone][ipt]) ; + + hJetTPCPts[icone][ipt] = new TH2F + ("JetTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone =" + +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", + 240,0,120,400,0,200); + hJetTPCPts[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}"); + hJetTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetTPCPts[icone][ipt]) ; + + hBkgTPCPts[icone][ipt] = new TH2F + ("BkgTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone =" + +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", + 240,0,120,400,0,200); + hBkgTPCPts[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}"); + hBkgTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgTPCPts[icone][ipt]) ; + + } + + //Counts + hNBkgs[icone][ipt] = new TH1F + ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "bkg multiplicity cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",9000,0,9000); + hNBkgs[icone][ipt]->SetYTitle("counts"); + hNBkgs[icone][ipt]->SetXTitle("N"); + fListHistos->Add(hNBkgs[icone][ipt]) ; + + hNBkgTPCs[icone][ipt] = new TH1F + ("NBkgTPCCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "bkg multiplicity cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",9000,0,9000); + hNBkgTPCs[icone][ipt]->SetYTitle("counts"); + hNBkgTPCs[icone][ipt]->SetXTitle("N"); + fListHistos->Add(hNBkgTPCs[icone][ipt]) ; + + + hNLeadings[icone][ipt] = new TH2F + ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120); + hNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)"); + hNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)"); + fListHistos->Add(hNLeadings[icone][ipt]) ; + + hNs[icone][ipt] = new TH1F + ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "Number of neutral jets, cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",120,0,120); + hNs[icone][ipt]->SetYTitle("N"); + hNs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)"); + fListHistos->Add(hNs[icone][ipt]) ; + + //Fragmentation Function + hJetFragments[icone][ipt] = new TH2F + ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); + hJetFragments[icone][ipt]->SetYTitle("x_{T}"); + hJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}"); + fListHistos->Add(hJetFragments[icone][ipt]) ; + + hBkgFragments[icone][ipt] = new TH2F + ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); + hBkgFragments[icone][ipt]->SetYTitle("x_{T}"); + hBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}"); + fListHistos->Add(hBkgFragments[icone][ipt]) ; + + + //Jet particle distribution + + hJetPtDists[icone][ipt] = new TH2F + ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+ + " GeV/c",120,0.,120.,120,0.,120.); + hJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetPtDists[icone][ipt]) ; + + hBkgPtDists[icone][ipt] = new TH2F + ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+ + " GeV/c",120,0.,120.,120,0.,120.); + hBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgPtDists[icone][ipt]) ; + + + if(fOnlyCharged){ + hNTPCLeadings[icone][ipt] = new TH2F + ("NTPCLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "p_{T #gamma} vs p_{T charge} cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120); + hNTPCLeadings[icone][ipt]->SetYTitle("p_{T charge} (GeV/c)"); + hNTPCLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)"); + fListHistos->Add(hNTPCLeadings[icone][ipt]) ; + + hNTPCs[icone][ipt] = new TH1F + ("NTPCJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "Number of charged jets, cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",120,0,120); + hNTPCs[icone][ipt]->SetYTitle("N"); + hNTPCs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)"); + fListHistos->Add(hNTPCs[icone][ipt]) ; + + hJetTPCFragments[icone][ipt] = new TH2F + ("JetTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); + hJetTPCFragments[icone][ipt]->SetYTitle("x_{T}"); + hJetTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}"); + fListHistos->Add(hJetTPCFragments[icone][ipt]) ; + + hBkgTPCFragments[icone][ipt] = new TH2F + ("BkgTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); + hBkgTPCFragments[icone][ipt]->SetYTitle("x_{T}"); + hBkgTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}"); + fListHistos->Add(hBkgTPCFragments[icone][ipt]) ; + + hJetTPCPtDists[icone][ipt] = new TH2F + ("JetTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>" + +fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.); + hJetTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hJetTPCPtDists[icone][ipt]) ; + + hBkgTPCPtDists[icone][ipt] = new TH2F + ("BkgTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], + "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>" + + fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.); + hBkgTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); + fListHistos->Add(hBkgTPCPtDists[icone][ipt]) ; + } + }//ipt + } //icone + }//If we want to study any cone or pt threshold +} + + +//____________________________________________________________________________ +void AliPHOSGammaJet::MakeJet(TClonesArray * pl, + Double_t ptg, Double_t phig, + 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. ; + else ptcut = 0.5; + } + + TClonesArray * jetList = new TClonesArray("TParticle",1000); + TClonesArray * bkgList = new TClonesArray("TParticle",1000); + + TLorentzVector bkg(0,0,0,0); + TLorentzVector lv (0,0,0,0); + + Double_t ptjet = 0.0; + Double_t ptbkg = 0.0; + Int_t n0 = 0; + Int_t n1 = 0; + Bool_t b1 = kFALSE; + Bool_t b0 = kFALSE; + + TIter next(pl) ; + TParticle * particle = 0 ; + + while ( (particle = dynamic_cast(next())) ) { + + b0 = kFALSE; + b1 = kFALSE; + + //Particles in jet + SetJet(particle, b0, fCone, etal, phil) ; + + if(b0){ + new((*jetList)[n0++]) TParticle(*particle) ; + particle->Momentum(lv); + if(particle->Pt() > ptcut ){ + jet+=lv; + ptjet+=particle->Pt(); + } + } + + //Background around (phi_gamma-pi, eta_leading) + SetJet(particle, b1, fCone,etal, phig) ; + + if(b1) { + new((*bkgList)[n1++]) TParticle(*particle) ; + if(particle->Pt() > ptcut ){ + bkg+=lv; + ptbkg+=particle->Pt(); + } + } + } + + ptjet = jet.Pt(); + ptbkg = bkg.Pt(); + + if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all")) + Info("MakeJet","Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg); + + + //Fill histograms + + Double_t rat = ptjet/ptg ; + Double_t ratbkg = ptbkg/ptg ; + + dynamic_cast + (fListHistos->FindObject("Jet"+conf+"Ratio"))->Fill(ptg,rat); + dynamic_cast + (fListHistos->FindObject("Jet"+conf+"Pt")) ->Fill(ptg,ptjet); + dynamic_cast + (fListHistos->FindObject("Bkg"+conf+"Ratio"))->Fill(ptg,ratbkg); + dynamic_cast + (fListHistos->FindObject("Bkg"+conf+"Pt")) ->Fill(ptg,ptbkg); + + + if(IsJetSelected(ptg,ptjet,conf) || fSelect){ + if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all")) + Info("MakeJet","JetSelected"); + dynamic_cast(fListHistos->FindObject("N"+conf+"Jet"))-> + Fill(ptg); + dynamic_cast(fListHistos->FindObject("N"+conf+"Leading")) + ->Fill(ptg,ptl); + FillJetHistos(jetList, ptg, conf, "Jet"); + FillJetHistos(bkgList, ptg, conf, "Bkg"); + } + + jetList ->Delete(); + bkgList ->Delete(); +} + +//____________________________________________________________________________ +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); + + Double_t ptjet = 0.0; + Double_t ptbkg = 0.0; + + Int_t n1 = 0; + Int_t n0 = 0; + Bool_t b1 = kFALSE; + Bool_t b0 = kFALSE; + + //Create as many jets as cones and pt thresholds are defined + Double_t maxcut = fJetRatioMaxCut; + Double_t mincut = fJetRatioMinCut; + + if(conf == "TPC"){ + maxcut = fJetTPCRatioMaxCut; + mincut = fJetTPCRatioMinCut; + } + + Double_t ratjet = 0; + Double_t ratbkg = 0; + + for(Int_t icone = 0; icone(next())) ) { + b1 = kFALSE; + b0 = kFALSE; + + SetJet(particle, b0, fCones[icone],etal, phil) ; + SetJet(particle, b1, fCones[icone],etal, phig) ; + + if(b0){ + new((*jetList)[n0++]) TParticle(*particle) ; + if(particle->Pt() > fPtThres[ipt] ) + ptjet+=particle->Pt(); + } + if(b1) { + new((*bkgList)[n1++]) TParticle(*particle) ; + if(particle->Pt() > fPtThres[ipt] ) + ptbkg+=particle->Pt(); + } + + } + + //Fill histograms + if(ptjet > 0.) { + + 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); + } + + ratjet = ptjet /ptg; + ratbkg = ptbkg/ptg; + + dynamic_cast + (fListHistos->FindObject("Jet"+conf+"RatioCone"+cone+"Pt"+ptcut)) + ->Fill(ptg,ratjet); + dynamic_cast + (fListHistos->FindObject("Jet"+conf+"PtCone"+cone+"Pt"+ptcut)) + ->Fill(ptg,ptjet); + + dynamic_cast + (fListHistos->FindObject("Bkg"+conf+"RatioCone"+cone+"Pt"+ptcut)) + ->Fill(ptg,ratbkg); + + dynamic_cast + (fListHistos->FindObject("Bkg"+conf+"PtCone"+cone+"Pt"+ptcut)) + ->Fill(ptg,ptbkg); + + + //Select Jet + if((ratjet < maxcut) && (ratjet > mincut)){ + + dynamic_cast + (fListHistos->FindObject("N"+conf+"JetCone"+cone+"Pt"+ptcut))-> + Fill(ptg); + dynamic_cast + (fListHistos->FindObject("N"+conf+"LeadingCone"+cone+"Pt"+ptcut)) + ->Fill(ptg,ptl); + + FillJetHistosAnyConeOrPt + (jetList,ptg,conf,"Jet",fNameCones[icone],fNamePtThres[ipt]); + FillJetHistosAnyConeOrPt + (bkgList,ptg,conf,"Bkg",fNameCones[icone],fNamePtThres[ipt]); + + } + } //ptjet > 0 + jetList ->Delete(); + bkgList ->Delete(); + }//for pt threshold + }// for cone +} + +//____________________________________________________________________________ +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); + + // get the detected direction + TVector3 pos = particle.Vect(); + pos*=460./energy; + TVector3 modpos = MakePosition(energy, pos) ; + modpos *= modenergy / 460.; + + Float_t modtheta = modpos.Theta(); + Float_t modphi = modpos.Phi(); + + // Set the modified 4-momentum of the reconstructed particle + Float_t py = modenergy*TMath::Sin(modphi)*TMath::Sin(modtheta); + Float_t px = modenergy*TMath::Cos(modphi)*TMath::Sin(modtheta); + Float_t pz = modenergy*TMath::Cos(modtheta); + + particle.SetPxPyPzE(px,py,pz,modenergy); + +} + +//____________________________________________________________________________ +TVector3 AliPHOSGammaJet::MakePosition(const Double_t energy, const TVector3 pos) +{ + // Smears the impact position according to the energy dependent position resolution + // A gaussian position distribution is assumed + + TVector3 newpos ; + + Double_t sigma = SigmaP(energy) ; + Double_t x = fRan.Gaus( pos.X(), sigma ) ; + Double_t z = fRan.Gaus( pos.Z(), sigma ) ; + Double_t y = pos.Y() ; + + newpos.SetX(x) ; + newpos.SetY(y) ; + newpos.SetZ(z) ; + + // Info("MakePosition","Theta dif %f",pos.Theta()-newpos.Theta()); +// Info("MakePosition","Phi dif %f",pos.Phi()-newpos.Phi()); + return newpos ; +} + //____________________________________________________________________________ void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0, TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) { // Perform isotropic decay pi0 -> 2 photons - // p0 is pi0 4-momentum (input) + // p0 is pi0 4-momentum (inut) // p1 and p2 are photon 4-momenta (output) // cout<<"Boost vector"<(fOutputFile->Get(what)); + if(h){ + h->Draw(); + } + else if (what == "all") { + TCanvas * can = new TCanvas("GammaJet", "Gamma-Jet Study",10,40,1600,1200); + can->cd() ; + can->Range(0,0,22,20); + TPaveLabel *pl1 = new TPaveLabel(1,18,20,19.5,"Titre","br"); + pl1->SetFillColor(18); + pl1->SetTextFont(32); + pl1->SetTextColor(49); + pl1->Draw(); + Int_t index ; + TPad * pad = 0; + Float_t begx = -0.29, begy = 0., endx = 0., endy = 0.30 ; + for (index = 0 ; index < fListHistos->GetEntries() ; index++) { + TString name("pad") ; + name += index ; + begx += 0.30 ; + endx += 0.30 ; + if (begx >= 1.0 || endx >= 1.0) { + begx = 0.01 ; + endx = 0.30 ; + begy += 0.30 ; + endy += 0.30 ; + } + printf("%f %f %f %f \n", begx, begy, endx, endy) ; + pad = new TPad(name,"This is a pad",begx,begy,endx,endy,33); + pad->Draw(); + pad->cd() ; + fListHistos->At(index)->Draw(option) ; + pad->Modified() ; + pad->Update() ; + can->cd() ; + } + + } + else{ + Info("Draw", "Histogram %s does not exist or unknown option", what.Data()) ; + fOutputFile->ls(); + } +} + +//____________________________________________________________________________ +void AliPHOSGammaJet::Print(const Option_t * opt) const +{ + + //Print some relevant parameters set for the analysis + if(! opt) + return; + + Info("Print", "%s %s", GetName(), GetTitle() ) ; + + printf("Eta cut : %f\n", fEtaCut) ; + printf("D phi max cut : %f\n", fPhiMaxCut) ; + printf("D phi min cut : %f\n", fPhiMinCut) ; + printf("Leading Ratio max cut : %f\n", fRatioMaxCut) ; + printf("Leading Ratio min cut : %f\n", fRatioMinCut) ; + printf("Jet Ratio max cut : %f\n", fJetRatioMaxCut) ; + printf("Jet Ratio min cut : %f\n", fJetRatioMinCut) ; + printf("Jet TPC Ratio max cut : %f\n", fJetTPCRatioMaxCut) ; + printf("Jet TPC Ratio min cut : %f\n", fJetTPCRatioMinCut) ; + printf("Fast recons : %d\n", fOptFast); + printf("Inv Mass max cut : %f\n", fInvMassMaxCut) ; + printf("Inv Mass min cut : %f\n", fInvMassMinCut) ; + +} + +//__________________________________________________________________________ +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) ) ; - Int_t iPrimary=-1, id_motherg = -1; + // 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 - TIter next(&particleList); - while ( (particle = (TParticle*)next()) ) { - iPrimary++; - Int_t ksCode = particle->GetStatusCode(); - Int_t iMother= particle->GetMother(0); + AliInfo( Form(" %d events read", cESDTree->GetEntriesFast()) ) ; + rv = cESDTree ; - if (ksCode == 21 && iMother == -1) - if(particle->GetPdgCode()==22){ - id_motherg = iPrimary; - //cout<<"idmother "<GetMother(0)){ - particle->Momentum(gamma); - id = iPrimary; - break; - } - }// kscode == 1 - }// while + } // 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 -//____________________________________________________________________________ -void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id) + 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) { - // Gets the leading particle from the list of charged particles - TParticle *particle = 0x0; - - Int_t iPrimary=-1; - Double_t ptmax = 0, pti = 0; - TIter next(&particleList); - while ( (particle = (TParticle*)next()) ) { - iPrimary++; - Int_t ksCode = particle->GetStatusCode(); - - if((ksCode == 1)&&(id != iPrimary) - &&(particle->GetPDG(0)->Charge()!=0)){ - pti = particle->Pt(); - if(pti> ptmax){ - ptmax = pti; - particle->Momentum(charge); - }//ptmax - }// kscode == 1 - }// while + + //Check if the particle is inside the cone defined by the leading particle + b = kFALSE; + + if(phi > TMath::TwoPi()) + phi-=TMath::TwoPi(); + if(phi < 0.) + phi+=TMath::TwoPi(); + + Double_t rad = 10000 + cone; + + if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone)) + rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+ + TMath::Power(part->Phi()-phi,2)); + else{ + if(part->Phi()-phi > TMath::TwoPi() - cone) + rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+ + TMath::Power((part->Phi()-TMath::TwoPi())-phi,2)); + if(part->Phi()-phi < -(TMath::TwoPi() - cone)) + rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+ + TMath::Power((part->Phi()+TMath::TwoPi())-phi,2)); + } + + if(rad < cone ) + b = kTRUE; + } //____________________________________________________________________________ -void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0) +Double_t AliPHOSGammaJet::SigmaE(Double_t energy) { - // Gets the leading pi0 from the list of particles - TParticle *particle = 0x0; - - Int_t iPrimary=-1; - Double_t ptmax = 0, pti = 0; - TIter next(&particleList); - while ( (particle = (TParticle*)next()) ) { - iPrimary++; - Int_t ksCode = particle->GetStatusCode(); - - if((ksCode == 2)) - { - pti = particle->Pt(); - if(pti> ptmax){ - ptmax = pti; - particle->Momentum(pi0); - }//ptmax - }// kscode == 1 - }// while + // Calculates the energy dependent energy resolution + + Double_t rv = -1 ; + + rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2) + + TMath::Power(fResPara2/TMath::Sqrt(energy), 2) + + TMath::Power(fResPara3, 2) ) ; + + return rv * energy ; } //____________________________________________________________________________ -// void AliPHOSGammaJet::GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id, -// Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2, -// Double_t & invmasscut1,Double_t & invmasscut2) -// { -// TParticle *particle = 0x0; - -// Int_t iPrimary=-1; -// Double_t El = 0, E12 = 0; -// TLorentzVector gamma_i,gamma_j; -// TIter next(&particleList); -// while ( (particle = (TParticle*)next()) ) { -// iPrimary++; -// Int_t ksCode = particle->GetStatusCode(); -// Int_t ksPdg = particle->GetPdgCode(); -// if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){ -// particle->Momentum(gamma_i); -// Int_t jPrimary=-1; -// TIter next2(&particleList); -// while ( (particle = (TParticle*)next2()) ) { -// jPrimary++; -// if(jPrimary>iPrimary){ -// ksCode = particle->GetStatusCode(); -// ksPdg = particle->GetPdgCode(); -// if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){ -// particle->Momentum(gamma_j); -// if(gamma_j.Angle(gamma_i.Vect())Fill(Eg,invmass); -// if((invmass>invmasscut1) && (invmassEl && (E12/Eg>ratiocut1) && (E12/EgEl && (E12/Eg>0.2 && E12/Eg<1.) -// }//(invmass>0.125) && (invmass<0.145) -// }//gamma_j.Angle(gamma_i.Vect())<0.04 -// }//(ksCode == 1) -// } -// }//while -// // cout<<"jPrimary "<