// --- ROOT system ---
-#include "TRandom.h"
+#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 "AliPHOSGammaJet.h"
#include "AliPHOSGetter.h"
+#include "AliPHOSGeometry.h"
+#include "TF1.h"
+#include "Riostream.h"
+#include "../PYTHIA6/AliGenPythia.h"
ClassImp(AliPHOSGammaJet)
//____________________________________________________________________________
-AliPHOSGammaJet::AliPHOSGammaJet() {
+AliPHOSGammaJet::AliPHOSGammaJet() : TTask()
+{
// ctor
+ fAngleMaxParam.Set(4) ;
+ fAngleMaxParam.Reset(0.);
+ fAnyConeOrPt = kFALSE ;
+ fEtaCut = 0. ;
+ fFastRec = 0 ;
+ fInvMassMaxCut = 0. ;
+ fInvMassMinCut = 0. ;
+ fJetRatioMaxCut = 0. ;
+ fJetRatioMinCut = 0. ;
+ fJetTPCRatioMaxCut = 0. ;
+ fJetTPCRatioMinCut = 0. ;
+ fMinDistance = 0. ;
+ fNEvent = 0 ;
+ fNCone = 0 ;
+ fNPt = 0 ;
+ fOnlyCharged = kFALSE ;
+ fOptFast = kFALSE ;
+ fPhiMaxCut = 0. ;
+ fPhiMinCut = 0. ;
+ fPtCut = 0. ;
+ fNeutralPtCut = 0. ;
+ fChargedPtCut = 0. ;
+ fRatioMaxCut = 0. ;
+ fRatioMinCut = 0. ;
+ fCone = 0 ;
+ fPtThreshold = 0 ;
+ fTPCCutsLikeEMCAL = kFALSE ;
+ fSelect = kFALSE ;
+ 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 ;
+ }
+ }
+ }
+
+ fOption = "" ;
+ 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<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
+ fListHistos->Add(h) ;
+ }
+ List() ;
}
//____________________________________________________________________________
-AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
+AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) :
+ TTask("GammaJet","Analysis of gamma-jet correlations")
+{
// ctor
- AliPHOSGetter::Instance(inputfilename) ;
+ 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;
+ fEtaCut = gj.fEtaCut ;
+ fInvMassMaxCut = gj.fInvMassMaxCut ;
+ fInvMassMinCut = gj.fInvMassMinCut ;
+ fFastRec = gj.fFastRec ;
+ fOption = gj.fOption ;
+ 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 ;
+
+ 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,
+ const AliPHOSGeometry * geom )
{
- // does the job
-
- if (! opt)
- return ;
- AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
- Int_t nEvent1 = gime->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--------------------
+ // List of particles copied to a root file.
+
+// Char_t sb[] = "bgrd/";
+// // cout<<sb<<endl;
+// Char_t si[10];
+// Char_t sf[] = "/list.root";
+// // cout<<sf<<endl;
+// sprintf(si,"%d",iEvent);
+// strcat(sb,si) ;
+// // cout<<si<<endl;
+// strcat(sb,sf) ;
+// // cout<<si<<endl;
+// TFile * f = TFile::Open(sb) ;
+// //cout<<f->GetName()<<endl;
+
+ Char_t fi[100];
+ sprintf(fi,"bgrd/%d/list.root",iEvent);
+ TFile * f = TFile::Open(fi) ;
+ cout<<f->GetName()<<endl;
+
+ TParticle *particle = new TParticle();
+ TTree *T = (TTree*) f->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 "<<T->GetEntries()<<endl;
+ 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())<fEtaCut){
+ new((*particleList)[index]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))
+ ->SetStatusCode(0) ;
+ index++ ;
+
+ new((*plCh)[indexCh]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
+ indexCh++ ;
+
+ if(strstr(fOption,"deb all")||strstr(fOption,"deb")){
+ dynamic_cast<TH1F*>(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(fOption,"deb all")|| strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(particle->Pt());
+
+ new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
+ indexNePHOS++ ;
+ }
+
+ if((particle->Phi()>fPhiEMCALCut[0]) &&
+ (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
+ {//Is in EMCAL
+ if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(particle->Pt());
+
+ new((*particleList)[index]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))
+ ->SetStatusCode(0) ;
+ index++ ;
+ new((*plNe)[indexNe]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(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())<fEtaCut)){
+
+ new((*particleList)[index]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))
+ ->SetStatusCode(0) ;
+ index++ ;
+
+ new((*plCh)[indexCh]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
+ indexCh++ ;
+ }
+ else if(charge == 0){
+ geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
+ if((particle->GetPdgCode() != 111) && (particle->Pt() > fNeutralPtCut) &&
+ (TMath::Abs(particle->Eta())<fEtaCut) ){
+ TLorentzVector part(particle->Px(),particle->Py(),
+ particle->Pz(),particle->Energy());
+ MakePhoton(part);
+ if(part.Pt() > fNeutralPtCut){
+
+ if(particle->Phi()>fPhiEMCALCut[0] &&
+ particle->Phi()<fPhiEMCALCut[1] && m == 0)
+ {
+ new((*particleList)[index]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))
+ ->SetStatusCode(0) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ index++ ;
+
+
+ new((*plNe)[indexNe]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plNe->At(indexNe)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
+ indexNe++ ;
+ }
+ if(m != 0)
+ {
+ new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
+ indexNePHOS++ ;
+ }
+ }
+ }
+ if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
+ (TMath::Abs(particle->Eta())<fEtaCut+1))
+ {
+
+ pPi0.SetPxPyPzE(particle->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(fOption,"deb all")|| strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(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()<fPhiEMCALCut[1] && m == 0){
+
+ new((*particleList)[index]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ index++ ;
+
+ new((*plNe)[indexNe]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
+ (dynamic_cast<TParticle*>(plNe->At(indexNe)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ indexNe++ ;
+ }
+ if(m != 0){
+ if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(particle->Pt());
+ new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
+ (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ indexNePHOS++;
+ }//In PHOS
+ }
+ }// if cell<distance
+ else {
+
+ dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
+ ->Fill(pPi0.E(),angle);
+
+ p1 = kFALSE;
+ if(pGamma1.Pt() > 0. && TMath::Abs(pGamma1.Eta())<fEtaCut){
+
+ MakePhoton(pGamma1);
+
+ if(pGamma1.Pt() > 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()<fPhiEMCALCut[1]
+ && m == 0){
+ if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(photon1->Pt());
+ new((*particleList)[index]) TParticle(*photon1) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
+ index++ ;
+
+ new((*plNe)[indexNe]) TParticle(*photon1) ;
+ (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
+ indexNe++ ;
+ p1 = kTRUE;
+ }
+ if(m != 0){
+ if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(photon1->Pt());
+ new((*plNePHOS)[indexNePHOS]) TParticle(*photon1) ;
+ (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
+ indexNePHOS++;
+ p1 = kTRUE;
+ }//Is in PHOS
+ }
+ }
+ if(pGamma2.Pt() > 0. && TMath::Abs(pGamma2.Eta())<fEtaCut){
+
+ MakePhoton(pGamma2);
+ if(pGamma2.Pt() > 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()<fPhiEMCALCut[1] && m == 0){
+ if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(photon2->Pt());
+ new((*particleList)[index]) TParticle(*photon2) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
+ index++ ;
+
+ new((*plNe)[indexNe]) TParticle(*photon2) ;
+ (dynamic_cast<TParticle*>(plNe->At(indexNe))) ->SetStatusCode(0) ;
+ indexNe++ ;
+ }
+ if(m != 0){
+ if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(photon2->Pt());
+ new((*plNePHOS)[indexNePHOS]) TParticle(*photon2) ;
+ (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
+ indexNePHOS++;
+ }
+ if(p1){
+ // e = (pGamma1+pGamma2).E();
+ // if(IsAngleInWindow(angle,e))
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("AnglePairAccepted"))->
+ Fill(pPi0.E(),angle);
+ }
+ }
+ }//photon2 in acceptance
+ }//if angle > mindist
+ }//if pi0
}
- //Keep Stable particles within eta range
- Float_t etacut = 0. ;
- if((particle->GetStatusCode() == 1)&&
- (TMath::Abs(particle->Eta())<etacut)){
- // Keep particles different from Pi0
- if(particle->GetPdgCode() != 111){
- particleList.Add(particle);
- n++;
- }
- //Decay Pi0 and keep it with different status name
- //Keep decay photons
- if(particle->GetPdgCode() == 111) {
- //cout<<"Pi0 "<<endl;
- n += 3 ;
- particle->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
}//for (iParticle<nParticle)
- TLorentzVector gamma,charge,pi0, gammapair;
- Int_t idg = -1;
- GetGammaJet(particleList,gamma, idg);
- GetLeadingCharge(particleList,charge, idg);
- GetLeadingPi0(particleList,pi0);
- Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), idg) ;
- Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
- Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
- // GetLeadingGammaPair(particleList, gammapair, idg,
- // 0.04,0.2, 1.0,0.13,0.14);
- Info("Pi0Decay", "Pair: %f", gammapair.Energy()) ;
+ }
+
+ //Info("AddHIJINGToList","End HIJING");
+}
+
+//____________________________________________________________________________
+Double_t AliPHOSGammaJet::CalculateJetRatioLimit(const Double_t ptg,
+ const Double_t *par,
+ 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 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);
+ return rat ;
+}
+
+//____________________________________________________________________________
+void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
+ TClonesArray * particleList,
+ TClonesArray * plCh,
+ TClonesArray * plNe,
+ TClonesArray * plNePHOS,
+ const AliPHOSGeometry * geom )
+{
+ //Info("CreateParticleList","Inside");
+ AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
+ 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())<fEtaCut){
+ //Fill lists
+
+ charge = TDatabasePDG::Instance()
+ ->GetParticle(particle->GetPdgCode())->Charge();
+ if((charge != 0) && (particle->Pt() > fChargedPtCut)){
+
+ if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(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(fOption,"deb all")|| strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(particle->Pt());
+
+ new((*plNePHOS)[indexNePHOS++]) TParticle(*particle) ;
+ }
+ if((particle->Phi()>fPhiEMCALCut[0]) &&
+ (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
+ {//Is in EMCAL
+ if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(particle->Pt());
+ new((*plNe)[indexNe++]) TParticle(*particle) ;
+ new((*particleList)[index++]) TParticle(*particle) ;
+ }
+ }
+ }
+ }//final particle, etacut
+ }//for (iParticle<nParticle)
+ }// No Op
+ else
+ {
+ Double_t mass = TDatabasePDG::Instance()->GetParticle(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())<fEtaCut)){
+ new((*plCh)[indexCh++]) TParticle(*particle) ;
+ new((*particleList)[index++]) TParticle(*particle) ;
+ }
+ else if(charge == 0) {
+ geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
+ if((particle->GetPdgCode() != 111) && particle->Pt() > 0 &&
+ (TMath::Abs(particle->Eta())<fEtaCut))
+ {
+
+ TLorentzVector part(particle->Px(),particle->Py(),
+ particle->Pz(),particle->Energy());
+
+ MakePhoton(part);
+
+ if(part.Pt() > fNeutralPtCut){
+ if(particle->Phi()>fPhiEMCALCut[0] &&
+ particle->Phi()<fPhiEMCALCut[1] && m == 0)
+ {
+ new((*particleList)[index]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ index++ ;
+
+ new((*plNe)[indexNe]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plNe->At(indexNe)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ indexNe++ ;
+ }
+ if(m != 0)
+ {
+ new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(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())<fEtaCut+1))
+ {
+
+ pPi0.SetPxPyPzE(particle->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()<fPhiEMCALCut[1] && m == 0){
+ if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(particle->Pt());
+
+ new((*plNe)[indexNe]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plNe->At(indexNe)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ new((*particleList)[index]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(particleList->At(index)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ index++;
+ indexNe++;
+ }//InEMCAL
+ if(m != 0){
+ if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(particle->Pt());
+ new((*plNePHOS)[indexNePHOS]) TParticle(*particle) ;
+ (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
+ ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
+ indexNePHOS++;
+ }//In PHOS
+ }//Small Pt
+ }// if cell<distance
+ else {
+
+ dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
+ ->Fill(pPi0.E(),angle);
+
+ Bool_t p1 = kFALSE;
+
+ if(pGamma1.Pt() > 0 && TMath::Abs(pGamma1.Eta())<fEtaCut){
+ MakePhoton(pGamma1);
+
+ if(pGamma1.Pt() > 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()<fPhiEMCALCut[1]
+ && m == 0){
+ if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(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(fOption,"deb all") || strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
+ Fill(photon1->Pt());
+ new((*plNePHOS)[indexNePHOS++]) TParticle(*photon1) ;
+ p1 = kTRUE;
+ }
+ }
+ }
+ if(pGamma2.Pt() > 0 && TMath::Abs(pGamma2.Eta())<fEtaCut){
+ MakePhoton(pGamma2);
+
+ if(pGamma2.Pt() > 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()<fPhiEMCALCut[1] && m == 0){
+ if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+ dynamic_cast<TH1F*>(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"))
+ dynamic_cast<TH1F*>(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<TH2F*>
+ (fListHistos->FindObject("AnglePairAccepted"))->
+ Fill(pPi0.E(),angle);
+ }
+ }
+ }//photon2 in acceptance
+ }//if angle > mindist
+ }//if pi0
+ }//If neutral
+ }//final particle, etacut
+ }//for (iParticle<nParticle)
+ }//OptFast
+ //gime->Delete() ;
+}
+
+
+
+//____________________________________________________________________________
+void AliPHOSGammaJet::Exec(Option_t *option)
+{
+ // does the job
+ fOption = option;
+ MakeHistos() ;
+
+ AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
+ const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
+
+
+// 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(fOption,"deb")||strstr(fOption,"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);
+
+ CreateParticleList(iEvent, particleList, plCh,plNe,plNePHOS, geom );
+
+// TLorentzVector pyjet(0,0,0,0);
+
+// 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());
+// }
+
+// }
+
+ if(fHIJING)
+ AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS, geom);
+
+
+ Bool_t IsInPHOS = kFALSE ;
+ GetGammaJet(plNePHOS, ptg, phig, etag, IsInPHOS) ;
+
+ if(IsInPHOS){
+
+ //Info("Exec"," In PHOS") ;
+ dynamic_cast<TH1F*>(fListHistos->FindObject("NGamma"))->Fill(ptg);
+ dynamic_cast<TH2F*>(fListHistos->FindObject("PhiGamma"))
+ ->Fill(ptg,phig);
+ dynamic_cast<TH2F*>(fListHistos->FindObject("EtaGamma"))
+ ->Fill(ptg,etag);
+ if(strstr(fOption,"deb")||strstr(fOption,"deb all"))
+ Info("Exec", "Gamma: pt %f, phi %f, eta %f", ptg,
+ phig,etag) ;
+
+// cout<<"n charged "<<plCh->GetEntries()<<endl;
+// cout<<"n neutral "<<plNe->GetEntries()<<endl;
+// cout<<"n All "<<particleList->GetEntries()<<endl;
+
+ GetLeadingCharge(plCh, ptg, phig, ptch, etach, phich) ;
+ GetLeadingPi0 (plNe, ptg, phig, ptpi, etapi, phipi) ;
+
+// cout<<"n2 charged "<<plCh->GetEntries()<<endl;
+// cout<<"n2 neutral "<<plNe->GetEntries()<<endl;
+// cout<<"n2 All "<<particleList->GetEntries()<<endl;
+
+
+ //TPC+EMCAL
+
+ //Is the leading cone inside EMCAL?
+ Bool_t insidech = kFALSE ;
+ if((phich - fCone) > 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<TH2F*>(fListHistos->FindObject("ChargeRatio"))
+ ->Fill(ptg,ptch/ptg);
+ dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiCharge"))
+ ->Fill(ptg,phig-phich);
+ dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaCharge"))
+ ->Fill(ptg,etag-etach);
+ if(strstr(fOption,"deb"))
+ Info("Exec"," Charged Leading") ;
+ }
+ if((ptpi > ptch) && insidepi){
+ phil = phipi ;
+ etal = etapi ;
+ ptl = ptpi ;
+
+ dynamic_cast<TH2F*>(fListHistos->FindObject("Pi0Ratio"))
+ ->Fill(ptg,ptpi/ptg);
+ dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiPi0"))
+ ->Fill(ptg,phig-phipi);
+ dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaPi0"))
+ ->Fill(ptg,etag-etapi);
+
+ if(ptpi > 0. && strstr(fOption,"deb"))
+ Info("Exec"," Pi0 Leading") ;
+ }
+
+ if(strstr(fOption,"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")){
+// 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<TH2F*>(fListHistos->FindObject("DeltaPhiJet"))
+// ->Fill(ptg,pyjet.Phi()-jet.Phi());
+// dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaJet"))
+// ->Fill(ptg,pyjet.Eta()-jet.Eta());
+// dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtJet"))
+// ->Fill(ptg,pyjet.Pt()-jet.Pt());
+ }
+ else
+ MakeJetAnyConeOrPt(particleList, ptg, phig, ptl, phil, etal, "");
+ }
+
+ //TPC
+ if(fOnlyCharged && ptch > 0.)
+ {
+ if(strstr(fOption,"deb"))
+ Info("Exec","Leading TPC pt %f, phi %f",ptch,phich);
+
+ dynamic_cast<TH2F*>(fListHistos->FindObject("TPCRatio"))
+ ->Fill(ptg,ptch/ptg);
+ dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPC"))
+ ->Fill(ptg,phig-phich);
+ dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPC"))
+ ->Fill(ptg,etag-etach);
+
+ if(!fAnyConeOrPt){
+
+ MakeJet(plCh, ptg, phig, ptch, phich, etach, "TPC",jettpc);
+
+ if(strstr(fOption,"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<TH2F*>(fListHistos->FindObject("DeltaPhiTPCJet"))
+// ->Fill(ptg,pyjet.Phi()-jettpc.Phi());
+// dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPCJet"))
+// ->Fill(ptg,pyjet.Eta()-jettpc.Eta());
+// dynamic_cast<TH2F*>(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)
+{
+
+ TParticle * particle = 0 ;
+ Int_t ipr = -1 ;
+ Float_t charge = 0;
+
+ TIter next(pl) ;
+ while ( (particle = dynamic_cast<TParticle*>(next())) ) {
+ ipr++ ;
+ Double_t pt = particle->Pt();
+
+ charge = TDatabasePDG::Instance()
+ ->GetParticle(particle->GetPdgCode())->Charge();
+ if(charge != 0){//Only jet Charged particles
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject(type+conf+"Fragment"))
+ ->Fill(ptg,pt/ptg);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject(type+conf+"PtDist"))
+ ->Fill(ptg,pt);
+ }
+ }
+ if(type == "Bkg"){
+ dynamic_cast<TH1F*>
+ (fListHistos->FindObject("NBkg"+conf))->Fill(ipr);
+ }
+}
+//____________________________________________________________________________
+void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg,
+ TString conf, TString type,
+ TString cone, TString ptcut)
+{
+
+ TParticle *particle = 0;
+ Int_t ipr=-1;
+ Float_t charge = 0;
+
+ TIter next(pl) ;
+ while ( (particle = dynamic_cast<TParticle*>(next())) ) {
+ ipr++;
+ Double_t pt = particle->Pt();
+ charge = TDatabasePDG::Instance()
+ ->GetParticle(particle->GetPdgCode())->Charge();
+ if(charge != 0){//Only jet Charged particles
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject(type+conf+"FragmentCone"+cone+"Pt"+ptcut))
+ ->Fill(ptg,pt/ptg);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject(type+conf+"PtDistCone"+cone+"Pt"+ptcut))
+ ->Fill(ptg,pt);
+ }
+ }//while
+
+ if(type == "Bkg"){
+ dynamic_cast<TH1F*>
+ (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)
+
+{
+ pt = -10.;
+ eta = -10.;
+ phi = -10.;
+
+ for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
+ TParticle * particle = dynamic_cast<TParticle *>(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;
+ }
+ }
+ }
+}
+
+//____________________________________________________________________________
+void AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl,
+ Double_t ptg, Double_t phig,
+ Double_t &pt, Double_t &eta, Double_t &phi)
+{
+
+ pt = -100.;
+ eta = -100;
+ phi = -100;
+
+ for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
+
+ TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
+
+ Double_t ptl = particle->Pt();
+ Double_t rat = ptl/ptg ;
+ Double_t phil = particle->Phi() ;
+
+ if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
+ (rat > 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)
+{
+ 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<TH2F*>
+ (fListHistos->FindObject("AnglePairNoCut"))->
+ Fill(e,0.1);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("InvMassPairNoCut"))->
+ Fill(ptg,1.35);
+
+ if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
+ (rat > fRatioMinCut) && (rat < fRatioMaxCut)) {
+
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("AnglePairLeadingCut"))->
+ Fill(e,0.1);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("InvMassPairLeadingCut"))->
+ Fill(ptg,1.35);
+
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("AnglePairAngleCut"))->
+ Fill(e,0.15);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("InvMassPairAngleCut"))->
+ Fill(ptg,1.36);
+
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("InvMassPairAllCut"))->
+ Fill(ptg,0.27);
+ dynamic_cast<TH2F*>
+ (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<TH2F*>
+ (fListHistos->FindObject("InvMassPairLeading"))->
+ Fill(ptg,1.35);
+ dynamic_cast<TH2F*>
+ (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++;
+ // if(particle->GetStatusCode() == 1){
+
+ 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)<fPhiMaxCut)){
+ eta = particle->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<TH2F*>
+ (fListHistos->FindObject("AnglePairNoCut"))->
+ Fill(e,angle);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("InvMassPairNoCut"))->
+ Fill(ptg,invmass);
+
+ if((rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
+ ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
+
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("AnglePairLeadingCut"))->
+ Fill(e,angle);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("InvMassPairLeadingCut"))->
+ Fill(ptg,invmass);
+
+ if(IsAngleInWindow(angle,e)){
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("AnglePairAngleCut"))->
+ Fill(e,angle);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("InvMassPairAngleCut"))->
+ Fill(ptg,invmass);
+
+ //Info("GetLeadingPi0","InvMass %f", invmass);
+ if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("InvMassPairAllCut"))->
+ Fill(ptg,invmass);
+ dynamic_cast<TH2F*>
+ (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<jprimary
+ }//while
+ }// if pdg = 22
+ // }
+ }//while
+
+ if(ef > 0.){
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("InvMassPairLeading"))->
+ Fill(ptg,invmassf);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("AnglePairLeading"))->
+ Fill(ef,anglef);
+ }
+ }//fOptFast
+ // printf("GetLeadingPi0: %f %f %f \n", pt, eta, phi) ;
+}
+
+
+//____________________________________________________________________________
+void AliPHOSGammaJet::InitParameters()
+{
+ 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" ;
+ fOption = "";
+ fHIJINGFileName = "galice.root" ;
+ fHIJING = kFALSE ;
+ fMinDistance = 3.6 ;
+ fEtaCut = 0.7 ;
+ fInvMassMaxCut = 0.15 ;
+ fInvMassMinCut = 0.12 ;
+ fOnlyCharged = kFALSE ;
+ fOptFast = kFALSE ;
+ 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 ;
+
+ //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){
+ 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<max)&&(angle>=min))
+ result = kTRUE;
+
+ return result;
+}
+
+//__________________________________________________________________________-
+Bool_t AliPHOSGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj,
+ const TString type ){
+
+ 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<TH2F*>(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<fNCone; icone++){
+ for(Int_t ipt = 0; ipt<fNPt;ipt++){
+ //Jet Pt / Gamma Pt
+
+ //Jet
+
+ hJetRatios[icone][ipt] = new TH2F
+ ("JetRatioCone"+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,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)
+{
+
+ 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<TParticle*>(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(fOption,"deb") || strstr(fOption,"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<TH2F*>
+ (fListHistos->FindObject("Jet"+conf+"Ratio"))->Fill(ptg,rat);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("Jet"+conf+"Pt")) ->Fill(ptg,ptjet);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("Bkg"+conf+"Ratio"))->Fill(ptg,ratbkg);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("Bkg"+conf+"Pt")) ->Fill(ptg,ptbkg);
+
+
+ if(IsJetSelected(ptg,ptjet,conf) || fSelect){
+ if(strstr(fOption,"deb") || strstr(fOption,"deb all"))
+ Info("MakeJet","JetSelected");
+ dynamic_cast<TH1F*>(fListHistos->FindObject("N"+conf+"Jet"))->
+ Fill(ptg);
+ dynamic_cast<TH2F*>(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){
+
+ 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<fNCone; icone++) {
+ for(Int_t ipt = 0; ipt<fNPt;ipt++) {
+
+ TString cone = fNameCones[icone] ;
+ TString ptcut = fNamePtThres[ipt] ;
+
+ TIter next(pl) ;
+ TParticle * particle = 0 ;
+
+ ptjet = 0 ;
+ ptbkg = 0 ;
+
+ while ( (particle = dynamic_cast<TParticle*>(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(fOption,"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<TH2F*>
+ (fListHistos->FindObject("Jet"+conf+"RatioCone"+cone+"Pt"+ptcut))
+ ->Fill(ptg,ratjet);
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("Jet"+conf+"PtCone"+cone+"Pt"+ptcut))
+ ->Fill(ptg,ptjet);
+
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("Bkg"+conf+"RatioCone"+cone+"Pt"+ptcut))
+ ->Fill(ptg,ratbkg);
+
+ dynamic_cast<TH2F*>
+ (fListHistos->FindObject("Bkg"+conf+"PtCone"+cone+"Pt"+ptcut))
+ ->Fill(ptg,ptbkg);
+
+
+ //Select Jet
+ if((ratjet < maxcut) && (ratjet > mincut)){
+
+ dynamic_cast<TH1F*>
+ (fListHistos->FindObject("N"+conf+"JetCone"+cone+"Pt"+ptcut))->
+ Fill(ptg);
+ dynamic_cast<TH2F*>
+ (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)
+{
+
+ 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"<<endl;
TVector3 b = p0.BoostVector();
angle = p1.Angle(p2.Vect());
//cout<<angle<<endl;
}
+
+//____________________________________________________________________________
+void AliPHOSGammaJet::Plot(TString what, Option_t * option) const
+{
+ TH2F * h = dynamic_cast<TH2F*>(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::GetGammaJet(TList &particleList, TLorentzVector &gamma, Int_t & id)
+void AliPHOSGammaJet::Print(char * opt)
+{
+ if(! opt)
+ return;
+
+ // Print
+ 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) ;
+
+}
+
+//___________________________________________________________________
+void AliPHOSGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone,
+ Double_t eta, Double_t phi)
{
- // Get the lists of jet particles and gamma
- TParticle *particle = 0x0;
- Int_t iPrimary=-1, id_motherg = -1;
+ 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;
- TIter next(&particleList);
- while ( (particle = (TParticle*)next()) ) {
- iPrimary++;
- Int_t ksCode = particle->GetStatusCode();
- Int_t iMother= particle->GetMother(0);
-
- if (ksCode == 21 && iMother == -1)
- if(particle->GetPdgCode()==22){
- id_motherg = iPrimary;
- //cout<<"idmother "<<id_motherg<<endl;
- }
- if(ksCode == 1){
-
- if(id_motherg == particle->GetMother(0)){
- particle->Momentum(gamma);
- id = iPrimary;
- break;
- }
- }// kscode == 1
- }// while
}
-
+
//____________________________________________________________________________
-void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id)
+Double_t AliPHOSGammaJet::SigmaE(Double_t energy)
{
- // 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
+ // 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::GetLeadingPi0(TList &particleList, TLorentzVector &pi0)
+Double_t AliPHOSGammaJet::SigmaP(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 position resolution
+
+ Double_t sigma = TMath::Sqrt(TMath::Power(fPosParaA,2) +
+ TMath::Power(fPosParaB,2) / energy) ;
+
+
+ return sigma ; // in cm
}
-//____________________________________________________________________________
-// 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())<thetacut){
-// Float_t invmass = (gamma_i+gamma_j).M();
-// h_invmass->Fill(Eg,invmass);
-// if((invmass>invmasscut1) && (invmass<invmasscut2)){
-// E12 = (gamma_i+gamma_j).Energy();
-// if(E12>El && (E12/Eg>ratiocut1) && (E12/Eg<ratiocut2)){
-// //cout<<E12<<" "<<E12/Eg<<endl;
-// El = E12;
-// id_i = iPrimary;
-// id_j = jPrimary;
-// gammapair = gamma_i+gamma_j;
-// }//E12>El && (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 "<<jPrimary<<endl;
-// }// if kscode 1
-// }//while
-// // cout<<"iPrimary "<<iPrimary<<endl;
-// }
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff