// Class for the photon identification.
// Clusters from calorimeters are identified as photons
// and kept in the AOD. Few histograms produced.
+// Produces input for other analysis classes like AliAnaPi0,
+// AliAnaParticleHadronCorrelation ...
//
// -- Author: Gustavo Conesa (LNF-INFN)
//////////////////////////////////////////////////////////////////////////////
// --- ROOT system ---
#include <TH2F.h>
+#include <TH3D.h>
#include <TClonesArray.h>
#include <TObjString.h>
//#include <Riostream.h>
+#include "TParticle.h"
+#include "TDatabasePDG.h"
// --- Analysis system ---
#include "AliAnaPhoton.h"
#include "AliCaloTrackReader.h"
+#include "AliStack.h"
#include "AliCaloPID.h"
#include "AliMCAnalysisUtils.h"
-#include "AliFidutialCut.h"
-#include "AliAODCaloCluster.h"
+#include "AliFiducialCut.h"
+#include "AliVCluster.h"
+#include "AliAODMCParticle.h"
+#include "AliMixedEvent.h"
+
ClassImp(AliAnaPhoton)
//____________________________________________________________________________
- AliAnaPhoton::AliAnaPhoton() :
- AliAnaPartCorrBaseClass(), fCalorimeter(""),
- fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
- fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
- //MC
- fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
- fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
- fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
- fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
- fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
- fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
- fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
+AliAnaPhoton::AliAnaPhoton() :
+ AliAnaPartCorrBaseClass(), fCalorimeter(""),
+ fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
+ fRejectTrackMatch(0), fTimeCutMin(-1), fTimeCutMax(999999),
+ fNCellsCut(0), fFillSSHistograms(kFALSE),
+ fNOriginHistograms(8), fNPrimaryHistograms(4),
+ fCheckConversion(kFALSE), fRemoveConvertedPair(kFALSE),
+ fAddConvertedPairsToAOD(kFALSE),
+ fMassCut(0), fConvAsymCut(1.), fConvDEtaCut(2.),
+ fConvDPhiMinCut(-1.), fConvDPhiMaxCut(7.),
+
+ // Histograms
+ fhNCellsE(0), fhMaxCellDiffClusterE(0),
+ fhEPhoton(0), fhPtPhoton(0),
+ fhPhiPhoton(0), fhEtaPhoton(0),
+ fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
+
+ // Conversion histograms
+ fhPtPhotonConv(0), fhEtaPhiPhotonConv(0), fhEtaPhi05PhotonConv(0),
+ fhConvDeltaEta(0), fhConvDeltaPhi(0), fhConvDeltaEtaPhi(0),
+ fhConvAsym(0), fhConvPt(0),
+ fhConvDistEta(0), fhConvDistEn(0), fhConvDistMass(0),
+ fhConvDistEtaCutEta(0), fhConvDistEnCutEta(0), fhConvDistMassCutEta(0),
+ fhConvDistEtaCutMass(0), fhConvDistEnCutMass(0),
+ fhConvDistEtaCutAsy(0), fhConvDistEnCutAsy(0),
+
+ //Shower shape histograms
+ fhDispE(0), fhLam0E(0), fhLam1E(0),
+ fhdDispE(0), fhdLam0E(0), fhdLam1E(0),
+ fhDispETRD(0), fhLam0ETRD(0), fhLam1ETRD(0),
+ fhdDispETRD(0), fhdLam0ETRD(0), fhdLam1ETRD(0),
+
+ fhNCellsLam0LowE(0), fhNCellsLam1LowE(0), fhNCellsDispLowE(0),
+ fhNCellsLam0HighE(0), fhNCellsLam1HighE(0), fhNCellsDispHighE(0),
+ fhNCellsdLam0LowE(0), fhNCellsdLam1LowE(0), fhNCellsdDispLowE(0),
+ fhNCellsdLam0HighE(0), fhNCellsdLam1HighE(0), fhNCellsdDispHighE(0),
+
+ fhEtaLam0LowE(0), fhPhiLam0LowE(0),
+ fhEtaLam0HighE(0), fhPhiLam0HighE(0),
+ fhLam0DispLowE(0), fhLam0DispHighE(0),
+ fhLam1Lam0LowE(0), fhLam1Lam0HighE(0),
+ fhDispLam1LowE(0), fhDispLam1HighE(0),
+ fhEtadLam0LowE(0), fhPhidLam0LowE(0),
+ fhEtadLam0HighE(0), fhPhidLam0HighE(0),
+ fhdLam0dDispLowE(0), fhdLam0dDispHighE(0),
+ fhdLam1dLam0LowE(0), fhdLam1dLam0HighE(0),
+ fhdDispdLam1LowE(0), fhdDispdLam1HighE(0),
+
+ //MC histograms
+ fhDeltaE(0), fhDeltaPt(0),
+ fhRatioE(0), fhRatioPt(0),
+ fh2E(0), fh2Pt(0),
+
+ // Conversion MC histograms
+ fhPtConversionTagged(0), fhPtAntiNeutronTagged(0),
+ fhPtAntiProtonTagged(0), fhPtUnknownTagged(0),
+ fhEtaPhiConversion(0), fhEtaPhi05Conversion(0),
+
+ fhConvDeltaEtaMCConversion(0), fhConvDeltaPhiMCConversion(0), fhConvDeltaEtaPhiMCConversion(0),
+ fhConvAsymMCConversion(0), fhConvPtMCConversion(0),
+ fhConvDispersionMCConversion(0), fhConvM02MCConversion(0),
+
+ fhConvDeltaEtaMCAntiNeutron(0), fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0),
+ fhConvAsymMCAntiNeutron(0), fhConvPtMCAntiNeutron(0),
+ fhConvDispersionMCAntiNeutron(0), fhConvM02MCAntiNeutron(0),
+ fhConvDeltaEtaMCAntiProton(0), fhConvDeltaPhiMCAntiProton(0), fhConvDeltaEtaPhiMCAntiProton(0),
+ fhConvAsymMCAntiProton(0), fhConvPtMCAntiProton(0),
+ fhConvDispersionMCAntiProton(0), fhConvM02MCAntiProton(0),
+ fhConvDeltaEtaMCString(0), fhConvDeltaPhiMCString(0), fhConvDeltaEtaPhiMCString(0),
+ fhConvAsymMCString(0), fhConvPtMCString(0),
+ fhConvDispersionMCString(0), fhConvM02MCString(0),
+ fhConvDistMCConversion(0), fhConvDistMCConversionCuts(0),
+
+ // Photon SS MC histograms
+ fhMCPhotonELambda0NoOverlap(0), fhMCPhotonELambda0TwoOverlap(0), fhMCPhotonELambda0NOverlap(0),
+ fhMCPhotonEdLambda0NoOverlap(0), fhMCPhotonEdLambda0TwoOverlap(0), fhMCPhotonEdLambda0NOverlap(0),
+
+ //Embedding
+ fhEmbeddedSignalFractionEnergy(0),
+ fhEmbedPhotonELambda0FullSignal(0), fhEmbedPhotonEdLambda0FullSignal(0),
+ fhEmbedPhotonELambda0MostlySignal(0), fhEmbedPhotonEdLambda0MostlySignal(0),
+ fhEmbedPhotonELambda0MostlyBkg(0), fhEmbedPhotonEdLambda0MostlyBkg(0),
+ fhEmbedPhotonELambda0FullBkg(0), fhEmbedPhotonEdLambda0FullBkg(0),
+ fhEmbedPi0ELambda0FullSignal(0), fhEmbedPi0EdLambda0FullSignal(0),
+ fhEmbedPi0ELambda0MostlySignal(0), fhEmbedPi0EdLambda0MostlySignal(0),
+ fhEmbedPi0ELambda0MostlyBkg(0), fhEmbedPi0EdLambda0MostlyBkg(0),
+ fhEmbedPi0ELambda0FullBkg(0), fhEmbedPi0EdLambda0FullBkg(0)
{
//default ctor
+ for(Int_t i = 0; i < 14; i++){
+ fhPtMC [i] = 0;
+ fhMCE [i] = 0;
+ fhPhiMC[i] = 0;
+ fhEtaMC[i] = 0;
+ }
+
+ for(Int_t i = 0; i < 7; i++){
+ fhPtPrimMC [i] = 0;
+ fhEPrimMC [i] = 0;
+ fhPhiPrimMC[i] = 0;
+ fhYPrimMC [i] = 0;
+
+ fhPtPrimMCAcc [i] = 0;
+ fhEPrimMCAcc [i] = 0;
+ fhPhiPrimMCAcc[i] = 0;
+ fhYPrimMCAcc [i] = 0;
+ }
+
+ for(Int_t i = 0; i < 6; i++){
+ fhMCELambda0 [i] = 0;
+ fhMCELambda1 [i] = 0;
+ fhMCEDispersion [i] = 0;
+ fhMCEdLambda0 [i] = 0;
+ fhMCEdLambda1 [i] = 0;
+ fhMCEdDispersion[i] = 0;
+ fhMCNCellsE [i] = 0;
+ fhMCMaxCellDiffClusterE[i] = 0;
+ fhMCLambda0vsClusterMaxCellDiffE0[i] = 0;
+ fhMCLambda0vsClusterMaxCellDiffE2[i] = 0;
+ fhMCLambda0vsClusterMaxCellDiffE6[i] = 0;
+ fhMCNCellsvsClusterMaxCellDiffE0 [i] = 0;
+ fhMCNCellsvsClusterMaxCellDiffE2 [i] = 0;
+ fhMCNCellsvsClusterMaxCellDiffE6 [i] = 0;
+ }
+
//Initialize parameters
InitParameters();
+}//____________________________________________________________________________
+AliAnaPhoton::~AliAnaPhoton()
+{
+ //dtor
+
}
-//____________________________________________________________________________
-AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) :
- AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),
- fMinDist(g.fMinDist),fMinDist2(g.fMinDist2), fMinDist3(g.fMinDist3),
- fRejectTrackMatch(g.fRejectTrackMatch),
- fhPtPhoton(g.fhPtPhoton),fhPhiPhoton(g.fhPhiPhoton),fhEtaPhoton(g.fhEtaPhoton),
- //MC
- fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
- fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
- fhPtISR(g.fhPtISR),fhPhiISR(g.fhPhiISR),fhEtaISR(g.fhEtaISR),
- fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
- fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
- fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
- fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
-
+//__________________________________________________________________
+Bool_t AliAnaPhoton::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
{
- // cpy ctor
+ //Select clusters if they pass different cuts
+ if(GetDebug() > 2)
+ printf("AliAnaPhoton::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
+ GetReader()->GetEventNumber(),
+ mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
+
+ //.......................................
+ //If too small or big energy, skip it
+ if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
+ if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
+
+ //.......................................
+ // TOF cut, BE CAREFUL WITH THIS CUT
+ Double_t tof = calo->GetTOF()*1e9;
+ if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
+ if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
+
+ //.......................................
+ if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
+ if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
+
+ //.......................................
+ //Check acceptance selection
+ if(IsFiducialCutOn()){
+ Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
+ if(! in ) return kFALSE ;
+ }
+ if(GetDebug() > 2) printf("Fiducial cut passed \n");
+ //.......................................
+ //Skip matched clusters with tracks
+ if(fRejectTrackMatch){
+ if(IsTrackMatched(calo)) {
+ if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
+ return kFALSE ;
+ }
+ else
+ if(GetDebug() > 2) printf(" Track-matching cut passed \n");
+ }// reject matched clusters
+
+ //.......................................
+ //Check Distance to Bad channel, set bit.
+ Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
+ if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
+ if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
+ return kFALSE ;
+ }
+ else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
+ //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
+
+ if(GetDebug() > 0)
+ printf("AliAnaPhoton::ClusterSelected() Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
+ GetReader()->GetEventNumber(),
+ mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
+
+ //All checks passed, cluster selected
+ return kTRUE;
+
}
-//_________________________________________________________________________
-AliAnaPhoton & AliAnaPhoton::operator = (const AliAnaPhoton & g)
-{
- // assignment operator
-
- if(&g == this) return *this;
-
- fCalorimeter = g.fCalorimeter ;
- fMinDist = g.fMinDist;
- fMinDist2 = g.fMinDist2;
- fMinDist3 = g.fMinDist3;
- fRejectTrackMatch = g.fRejectTrackMatch;
-
- fhPtPhoton = g.fhPtPhoton ;
- fhPhiPhoton = g.fhPhiPhoton ;
- fhEtaPhoton = g.fhEtaPhoton ;
-
- fhPtPrompt = g.fhPtPrompt;
- fhPhiPrompt = g.fhPhiPrompt;
- fhEtaPrompt = g.fhEtaPrompt;
- fhPtFragmentation = g.fhPtFragmentation;
- fhPhiFragmentation = g.fhPhiFragmentation;
- fhEtaFragmentation = g.fhEtaFragmentation;
- fhPtISR = g.fhPtISR;
- fhPhiISR = g.fhPhiISR;
- fhEtaISR = g.fhEtaISR;
- fhPtPi0Decay = g.fhPtPi0Decay;
- fhPhiPi0Decay = g.fhPhiPi0Decay;
- fhEtaPi0Decay = g.fhEtaPi0Decay;
- fhPtOtherDecay = g.fhPtOtherDecay;
- fhPhiOtherDecay = g.fhPhiOtherDecay;
- fhEtaOtherDecay = g.fhEtaOtherDecay;
- fhPtConversion = g. fhPtConversion;
- fhPhiConversion = g.fhPhiConversion;
- fhEtaConversion = g.fhEtaConversion;
- fhPtUnknown = g.fhPtUnknown;
- fhPhiUnknown = g.fhPhiUnknown;
- fhEtaUnknown = g.fhEtaUnknown;
-
- return *this;
+//_____________________________________________________________
+void AliAnaPhoton::FillAcceptanceHistograms(){
+ //Fill acceptance histograms if MC data is available
+ if(GetReader()->ReadStack()){
+ AliStack * stack = GetMCStack();
+ if(stack){
+ for(Int_t i=0 ; i<stack->GetNtrack(); i++){
+ TParticle * prim = stack->Particle(i) ;
+ Int_t pdg = prim->GetPdgCode();
+ //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
+ // prim->GetName(), prim->GetPdgCode());
+
+ if(pdg == 22){
+
+ // Get tag of this particle photon from fragmentation, decay, prompt ...
+ Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
+ if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
+ //A conversion photon from a hadron, skip this kind of photon
+ // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!: tag %d, conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d\n",tag,
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
+ // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
+
+ return;
+ }
+
+ //Get photon kinematics
+ if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
+
+ Double_t photonY = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
+ Double_t photonE = prim->Energy() ;
+ Double_t photonPt = prim->Pt() ;
+ Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
+ if(photonPhi < 0) photonPhi+=TMath::TwoPi();
+ Double_t photonEta = prim->Eta() ;
+
+
+ //Check if photons hit the Calorimeter
+ TLorentzVector lv;
+ prim->Momentum(lv);
+ Bool_t inacceptance = kFALSE;
+ if (fCalorimeter == "PHOS"){
+ if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
+ Int_t mod ;
+ Double_t x,z ;
+ if(GetPHOSGeometry()->ImpactOnEmc(prim,mod,z,x))
+ inacceptance = kTRUE;
+ if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ else{
+ if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
+ inacceptance = kTRUE ;
+ if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ }
+ else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
+ if(GetEMCALGeometry()){
+
+ Int_t absID=0;
+
+ GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
+
+ if( absID >= 0)
+ inacceptance = kTRUE;
+
+ // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
+ // inacceptance = kTRUE;
+ if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ else{
+ if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
+ inacceptance = kTRUE ;
+ if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ } //In EMCAL
+
+ //Fill histograms
+
+ fhYPrimMC[mcPPhoton]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0)
+ {
+ fhEPrimMC [mcPPhoton]->Fill(photonE ) ;
+ fhPtPrimMC [mcPPhoton]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPPhoton]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPPhoton] ->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPPhoton] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPPhoton] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPPhoton]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPPhoton] ->Fill(photonE , photonY) ;
+ }//Accepted
+
+ //Origin of photon
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[mcPPrompt])
+ {
+ fhYPrimMC[mcPPrompt]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPPrompt]->Fill(photonE ) ;
+ fhPtPrimMC [mcPPrompt]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPPrompt]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPPrompt] ->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPPrompt] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPPrompt] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPPrompt]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPPrompt] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[mcPFragmentation])
+ {
+ fhYPrimMC[mcPFragmentation]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPFragmentation]->Fill(photonE ) ;
+ fhPtPrimMC [mcPFragmentation]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPFragmentation]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPFragmentation] ->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPFragmentation] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPFragmentation] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPFragmentation]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPFragmentation] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[mcPISR])
+ {
+ fhYPrimMC[mcPISR]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPISR]->Fill(photonE ) ;
+ fhPtPrimMC [mcPISR]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPISR]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPISR]->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPISR] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPISR] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPISR]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPISR] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[mcPPi0Decay])
+ {
+ fhYPrimMC[mcPPi0Decay]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPPi0Decay]->Fill(photonE ) ;
+ fhPtPrimMC [mcPPi0Decay]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPPi0Decay]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPPi0Decay] ->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPPi0Decay] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPPi0Decay] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPPi0Decay]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPPi0Decay] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
+ GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) && fhEPrimMC[mcPOtherDecay])
+ {
+ fhYPrimMC[mcPOtherDecay]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPOtherDecay]->Fill(photonE ) ;
+ fhPtPrimMC [mcPOtherDecay]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPOtherDecay]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPOtherDecay] ->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPOtherDecay] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPOtherDecay] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPOtherDecay]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPOtherDecay] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(fhEPrimMC[mcPOther])
+ {
+ fhYPrimMC[mcPOther]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPOther]->Fill(photonE ) ;
+ fhPtPrimMC [mcPOther]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPOther]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPOther] ->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPOther] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPOther] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPOther]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPOther] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }//Other origin
+ }// Primary photon
+ }//loop on primaries
+ }//stack exists and data is MC
+ }//read stack
+ else if(GetReader()->ReadAODMCParticles()){
+ TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
+ if(mcparticles){
+ Int_t nprim = mcparticles->GetEntriesFast();
+
+ for(Int_t i=0; i < nprim; i++)
+ {
+ AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
+
+ Int_t pdg = prim->GetPdgCode();
+
+ if(pdg == 22){
+
+ // Get tag of this particle photon from fragmentation, decay, prompt ...
+ Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
+ if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
+ //A conversion photon from a hadron, skip this kind of photon
+// printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!: tag %d, conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d\n",tag,
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
+// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
+
+ return;
+ }
+
+ //Get photon kinematics
+ if(prim->E() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
+
+ Double_t photonY = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
+ Double_t photonE = prim->E() ;
+ Double_t photonPt = prim->Pt() ;
+ Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
+ if(photonPhi < 0) photonPhi+=TMath::TwoPi();
+ Double_t photonEta = prim->Eta() ;
+
+ //Check if photons hit the Calorimeter
+ TLorentzVector lv;
+ lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
+ Bool_t inacceptance = kFALSE;
+ if (fCalorimeter == "PHOS"){
+ if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
+ Int_t mod ;
+ Double_t x,z ;
+ Double_t vtx[]={prim->Xv(),prim->Yv(),prim->Zv()};
+ if(GetPHOSGeometry()->ImpactOnEmc(vtx, prim->Theta(),prim->Phi(),mod,z,x))
+ inacceptance = kTRUE;
+ if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ else{
+ if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
+ inacceptance = kTRUE ;
+ if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ }
+ else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
+ if(GetEMCALGeometry()){
+
+ Int_t absID=0;
+
+ GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
+
+ if( absID >= 0)
+ inacceptance = kTRUE;
+
+ if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ else{
+ if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
+ inacceptance = kTRUE ;
+ if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+ }
+ } //In EMCAL
+
+ //Fill histograms
+
+ fhYPrimMC[mcPPhoton]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0)
+ {
+ fhEPrimMC [mcPPhoton]->Fill(photonE ) ;
+ fhPtPrimMC [mcPPhoton]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPPhoton]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPPhoton]->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPPhoton] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPPhoton] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPPhoton]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPPhoton] ->Fill(photonE , photonY) ;
+ }//Accepted
+
+
+ //Origin of photon
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[mcPPrompt])
+ {
+ fhYPrimMC[mcPPrompt]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPPrompt]->Fill(photonE ) ;
+ fhPtPrimMC [mcPPrompt]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPPrompt]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPPrompt]->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPPrompt] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPPrompt] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPPrompt]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPPrompt] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[mcPFragmentation] )
+ {
+ fhYPrimMC[mcPFragmentation]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPFragmentation]->Fill(photonE ) ;
+ fhPtPrimMC [mcPFragmentation]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPFragmentation]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPFragmentation]->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPFragmentation] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPFragmentation] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPFragmentation]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPFragmentation] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[mcPISR])
+ {
+ fhYPrimMC[mcPISR]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPISR]->Fill(photonE ) ;
+ fhPtPrimMC [mcPISR]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPISR]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPISR]->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPISR] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPISR] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPISR]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPISR] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[mcPPi0Decay])
+ {
+ fhYPrimMC[mcPPi0Decay]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPPi0Decay]->Fill(photonE ) ;
+ fhPtPrimMC [mcPPi0Decay]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPPi0Decay]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPPi0Decay]->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPPi0Decay] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPPi0Decay] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPPi0Decay]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPPi0Decay] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
+ GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) && fhEPrimMC[mcPOtherDecay])
+ {
+ fhYPrimMC[mcPOtherDecay]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPOtherDecay]->Fill(photonE ) ;
+ fhPtPrimMC [mcPOtherDecay]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPOtherDecay]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPOtherDecay]->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPOtherDecay] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPOtherDecay] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPOtherDecay]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPOtherDecay] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }
+ else if(fhEPrimMC[mcPOther])
+ {
+ fhYPrimMC[mcPOther]->Fill(photonPt, photonY) ;
+ if(TMath::Abs(photonY) < 1.0){
+ fhEPrimMC [mcPOther]->Fill(photonE ) ;
+ fhPtPrimMC [mcPOther]->Fill(photonPt) ;
+ fhPhiPrimMC[mcPOther]->Fill(photonE , photonPhi) ;
+ fhYPrimMC[mcPOther]->Fill(photonE , photonEta) ;
+ }
+ if(inacceptance){
+ fhEPrimMCAcc[mcPOther] ->Fill(photonE ) ;
+ fhPtPrimMCAcc[mcPOther] ->Fill(photonPt) ;
+ fhPhiPrimMCAcc[mcPOther]->Fill(photonE , photonPhi) ;
+ fhYPrimMCAcc[mcPOther] ->Fill(photonE , photonY) ;
+ }//Accepted
+ }//Other origin
+ }// Primary photon
+ }//loop on primaries
+
+ }//mc array exists and data is MC
+ } // read AOD MC
}
-//____________________________________________________________________________
-AliAnaPhoton::~AliAnaPhoton()
-{
- //dtor
+//__________________________________________________________________
+void AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag){
+
+ //Fill cluster Shower Shape histograms
+
+ if(!fFillSSHistograms || GetMixedEvent()) return;
+
+ Float_t energy = cluster->E();
+ Int_t ncells = cluster->GetNCells();
+ Int_t ncells2 = ncells*ncells;
+ Float_t lambda0 = cluster->GetM02();
+ Float_t lambda1 = cluster->GetM20();
+ Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
+
+ TLorentzVector mom;
+ if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
+ cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
+ else{
+ Double_t vertex[]={0,0,0};
+ cluster->GetMomentum(mom,vertex) ;
+ }
+
+ Float_t eta = mom.Eta();
+ Float_t phi = mom.Phi();
+ if(phi < 0) phi+=TMath::TwoPi();
+
+ fhLam0E ->Fill(energy,lambda0);
+ fhLam1E ->Fill(energy,lambda1);
+ fhDispE ->Fill(energy,disp);
+ fhdLam0E->Fill(energy,lambda0/ncells2);
+ fhdLam1E->Fill(energy,lambda1/ncells2);
+ fhdDispE->Fill(energy,disp/ncells2);
+
+ if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5){
+ fhLam0ETRD->Fill(energy,lambda0);
+ fhLam1ETRD->Fill(energy,lambda1);
+ fhDispETRD->Fill(energy,disp);
+ fhdLam0ETRD->Fill(energy,lambda0/ncells2);
+ fhdLam1ETRD->Fill(energy,lambda1/ncells2);
+ fhdDispETRD->Fill(energy,disp/ncells2);
+ }
+
+ if(energy < 2){
+ fhNCellsLam0LowE ->Fill(ncells,lambda0);
+ fhNCellsLam1LowE ->Fill(ncells,lambda1);
+ fhNCellsDispLowE ->Fill(ncells,disp);
+ fhNCellsdLam0LowE->Fill(ncells,lambda0/ncells2);
+ fhNCellsdLam1LowE->Fill(ncells,lambda1/ncells2);
+ fhNCellsdDispLowE->Fill(ncells,disp/ncells2);
+
+ fhLam1Lam0LowE ->Fill(lambda1,lambda0);
+ fhLam0DispLowE ->Fill(lambda0,disp);
+ fhDispLam1LowE ->Fill(disp,lambda1);
+ fhEtaLam0LowE ->Fill(eta,lambda0);
+ fhPhiLam0LowE ->Fill(phi,lambda0);
+
+ fhdLam1dLam0LowE->Fill(lambda1/ncells2,lambda0/ncells2);
+ fhdLam0dDispLowE->Fill(lambda0/ncells2,disp/ncells2);
+ fhdDispdLam1LowE->Fill(disp/ncells2,lambda1/ncells2);
+ fhEtadLam0LowE ->Fill(eta,lambda0/ncells2);
+ fhPhidLam0LowE ->Fill(phi,lambda0/ncells2);
+ }
+ else {
+ fhNCellsLam0HighE ->Fill(ncells,lambda0);
+ fhNCellsLam1HighE ->Fill(ncells,lambda1);
+ fhNCellsDispHighE ->Fill(ncells,disp);
+ fhNCellsdLam0HighE->Fill(ncells,lambda0/ncells2);
+ fhNCellsdLam1HighE->Fill(ncells,lambda1/ncells2);
+ fhNCellsdDispHighE->Fill(ncells,disp/ncells2);
+
+ fhLam1Lam0HighE ->Fill(lambda1,lambda0);
+ fhLam0DispHighE ->Fill(lambda0,disp);
+ fhDispLam1HighE ->Fill(disp,lambda1);
+ fhEtaLam0HighE ->Fill(eta, lambda0);
+ fhPhiLam0HighE ->Fill(phi, lambda0);
+
+ fhdLam1dLam0HighE->Fill(lambda1/ncells2,lambda0/ncells2);
+ fhdLam0dDispHighE->Fill(lambda0/ncells2,disp/ncells2);
+ fhdDispdLam1HighE->Fill(disp/ncells2,lambda1/ncells2);
+ fhEtadLam0HighE ->Fill(eta, lambda0/ncells2);
+ fhPhidLam0HighE ->Fill(phi, lambda0/ncells2);
+ }
+
+ if(IsDataMC()){
+
+ AliVCaloCells* cells = 0;
+ if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
+ else cells = GetPHOSCells();
+
+ //Fill histograms to check shape of embedded clusters
+ Float_t fraction = 0;
+ if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
+
+ Float_t clusterE = 0; // recalculate in case corrections applied.
+ Float_t cellE = 0;
+ for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
+ cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
+ clusterE+=cellE;
+ fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
+ }
+
+ //Fraction of total energy due to the embedded signal
+ fraction/=clusterE;
+
+ if(GetDebug() > 1 ) printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
+
+ fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
+
+ } // embedded fraction
+
+ // Get the fraction of the cluster energy that carries the cell with highest energy
+ Int_t absID =-1 ;
+ Float_t maxCellFraction = 0.;
+
+ absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
+
+ // Check the origin and fill histograms
+ if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
+ !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
+ !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
+ !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)){
+ fhMCELambda0[mcssPhoton] ->Fill(energy, lambda0);
+ fhMCEdLambda0[mcssPhoton] ->Fill(energy, lambda0/ncells2);
+ fhMCELambda1[mcssPhoton] ->Fill(energy, lambda1);
+ fhMCEdLambda1[mcssPhoton] ->Fill(energy, lambda1/ncells2);
+ fhMCEDispersion[mcssPhoton] ->Fill(energy, disp);
+ fhMCEdDispersion[mcssPhoton]->Fill(energy, disp/ncells2);
+ fhMCNCellsE[mcssPhoton] ->Fill(energy, ncells);
+ fhMCMaxCellDiffClusterE[mcssPhoton]->Fill(energy,maxCellFraction);
+
+ if (energy < 2.){
+ fhMCLambda0vsClusterMaxCellDiffE0[mcssPhoton]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE0[mcssPhoton] ->Fill(ncells, maxCellFraction);
+ }
+ else if(energy < 6.){
+ fhMCLambda0vsClusterMaxCellDiffE2[mcssPhoton]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE2[mcssPhoton] ->Fill(ncells, maxCellFraction);
+ }
+ else{
+ fhMCLambda0vsClusterMaxCellDiffE6[mcssPhoton]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE6[mcssPhoton] ->Fill(ncells, maxCellFraction);
+ }
+
+ if(!GetReader()->IsEmbeddedClusterSelectionOn()){
+ //Check particle overlaps in cluster
+
+ //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
+ Int_t ancPDG = 0, ancStatus = -1;
+ TLorentzVector momentum; TVector3 prodVertex;
+ Int_t ancLabel = 0;
+ Int_t noverlaps = 1;
+ for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
+ ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
+ if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
+ }
+
+ if(noverlaps == 1){
+ fhMCPhotonELambda0NoOverlap ->Fill(energy, lambda0);
+ fhMCPhotonEdLambda0NoOverlap ->Fill(energy, lambda0/ncells2);
+ }
+ else if(noverlaps == 2){
+ fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
+ fhMCPhotonEdLambda0TwoOverlap->Fill(energy, lambda0/ncells2);
+ }
+ else if(noverlaps > 2){
+ fhMCPhotonELambda0NOverlap ->Fill(energy, lambda0);
+ fhMCPhotonEdLambda0NOverlap ->Fill(energy, lambda0/ncells2);
+ }
+ else {
+ printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
+ }
+ }//No embedding
+
+ //Fill histograms to check shape of embedded clusters
+ if(GetReader()->IsEmbeddedClusterSelectionOn()){
+
+ if (fraction > 0.9)
+ {
+ fhEmbedPhotonELambda0FullSignal ->Fill(energy, lambda0);
+ fhEmbedPhotonEdLambda0FullSignal ->Fill(energy, lambda0/ncells2);
+ }
+ else if(fraction > 0.5)
+ {
+ fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
+ fhEmbedPhotonEdLambda0MostlySignal->Fill(energy, lambda0/ncells2);
+ }
+ else if(fraction > 0.1)
+ {
+ fhEmbedPhotonELambda0MostlyBkg ->Fill(energy, lambda0);
+ fhEmbedPhotonEdLambda0MostlyBkg ->Fill(energy, lambda0/ncells2);
+ }
+ else
+ {
+ fhEmbedPhotonELambda0FullBkg ->Fill(energy, lambda0);
+ fhEmbedPhotonEdLambda0FullBkg ->Fill(energy, lambda0/ncells2);
+ }
+ } // embedded
+
+ }//photon no conversion
+ else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)){
+ fhMCELambda0[mcssElectron] ->Fill(energy, lambda0);
+ fhMCEdLambda0[mcssElectron] ->Fill(energy, lambda0/ncells2);
+ fhMCELambda1[mcssElectron] ->Fill(energy, lambda1);
+ fhMCEdLambda1[mcssElectron] ->Fill(energy, lambda1/ncells2);
+ fhMCEDispersion[mcssElectron] ->Fill(energy, disp);
+ fhMCEdDispersion[mcssElectron]->Fill(energy, disp/ncells2);
+ fhMCNCellsE[mcssElectron] ->Fill(energy, ncells);
+ fhMCMaxCellDiffClusterE[mcssElectron]->Fill(energy,maxCellFraction);
+
+ if (energy < 2.){
+ fhMCLambda0vsClusterMaxCellDiffE0[mcssElectron]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE0[mcssElectron] ->Fill(ncells, maxCellFraction);
+ }
+ else if(energy < 6.){
+ fhMCLambda0vsClusterMaxCellDiffE2[mcssElectron]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE2[mcssElectron] ->Fill(ncells, maxCellFraction);
+ }
+ else{
+ fhMCLambda0vsClusterMaxCellDiffE6[mcssElectron]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE6[mcssElectron] ->Fill(ncells, maxCellFraction);
+ }
+ }//electron
+ else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
+ GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) ){
+ fhMCELambda0[mcssConversion] ->Fill(energy, lambda0);
+ fhMCEdLambda0[mcssConversion] ->Fill(energy, lambda0/ncells2);
+ fhMCELambda1[mcssConversion] ->Fill(energy, lambda1);
+ fhMCEdLambda1[mcssConversion] ->Fill(energy, lambda1/ncells2);
+ fhMCEDispersion[mcssConversion] ->Fill(energy, disp);
+ fhMCEdDispersion[mcssConversion]->Fill(energy, disp/ncells2);
+ fhMCNCellsE[mcssConversion] ->Fill(energy, ncells);
+ fhMCMaxCellDiffClusterE[mcssConversion]->Fill(energy,maxCellFraction);
+
+ if (energy < 2.){
+ fhMCLambda0vsClusterMaxCellDiffE0[mcssConversion]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE0[mcssConversion] ->Fill(ncells, maxCellFraction);
+ }
+ else if(energy < 6.){
+ fhMCLambda0vsClusterMaxCellDiffE2[mcssConversion]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE2[mcssConversion] ->Fill(ncells, maxCellFraction);
+ }
+ else{
+ fhMCLambda0vsClusterMaxCellDiffE6[mcssConversion]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE6[mcssConversion] ->Fill(ncells, maxCellFraction);
+ }
+
+ }//conversion photon
+ else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) ){
+ fhMCELambda0[mcssPi0] ->Fill(energy, lambda0);
+ fhMCEdLambda0[mcssPi0] ->Fill(energy, lambda0/ncells2);
+ fhMCELambda1[mcssPi0] ->Fill(energy, lambda1);
+ fhMCEdLambda1[mcssPi0] ->Fill(energy, lambda1/ncells2);
+ fhMCEDispersion[mcssPi0] ->Fill(energy, disp);
+ fhMCEdDispersion[mcssPi0]->Fill(energy, disp/ncells2);
+ fhMCNCellsE[mcssPi0] ->Fill(energy, ncells);
+ fhMCMaxCellDiffClusterE[mcssPi0]->Fill(energy,maxCellFraction);
+
+ if (energy < 2.){
+ fhMCLambda0vsClusterMaxCellDiffE0[mcssPi0]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE0[mcssPi0] ->Fill(ncells, maxCellFraction);
+ }
+ else if(energy < 6.){
+ fhMCLambda0vsClusterMaxCellDiffE2[mcssPi0]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE2[mcssPi0] ->Fill(ncells, maxCellFraction);
+ }
+ else{
+ fhMCLambda0vsClusterMaxCellDiffE6[mcssPi0]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE6[mcssPi0] ->Fill(ncells, maxCellFraction);
+ }
+
+ //Fill histograms to check shape of embedded clusters
+ if(GetReader()->IsEmbeddedClusterSelectionOn()){
+
+ if (fraction > 0.9)
+ {
+ fhEmbedPi0ELambda0FullSignal ->Fill(energy, lambda0);
+ fhEmbedPi0EdLambda0FullSignal ->Fill(energy, lambda0/ncells2);
+ }
+ else if(fraction > 0.5)
+ {
+ fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
+ fhEmbedPi0EdLambda0MostlySignal->Fill(energy, lambda0/ncells2);
+ }
+ else if(fraction > 0.1)
+ {
+ fhEmbedPi0ELambda0MostlyBkg ->Fill(energy, lambda0);
+ fhEmbedPi0EdLambda0MostlyBkg ->Fill(energy, lambda0/ncells2);
+ }
+ else
+ {
+ fhEmbedPi0ELambda0FullBkg ->Fill(energy, lambda0);
+ fhEmbedPi0EdLambda0FullBkg ->Fill(energy, lambda0/ncells2);
+ }
+ } // embedded
+
+ }//pi0
+ else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) ){
+ fhMCELambda0[mcssEta] ->Fill(energy, lambda0);
+ fhMCEdLambda0[mcssEta] ->Fill(energy, lambda0/ncells2);
+ fhMCELambda1[mcssEta] ->Fill(energy, lambda1);
+ fhMCEdLambda1[mcssEta] ->Fill(energy, lambda1/ncells2);
+ fhMCEDispersion[mcssEta] ->Fill(energy, disp);
+ fhMCEdDispersion[mcssEta]->Fill(energy, disp/ncells2);
+ fhMCNCellsE[mcssEta] ->Fill(energy, ncells);
+ fhMCMaxCellDiffClusterE[mcssEta]->Fill(energy,maxCellFraction);
+
+ if (energy < 2.){
+ fhMCLambda0vsClusterMaxCellDiffE0[mcssEta]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE0[mcssEta] ->Fill(ncells, maxCellFraction);
+ }
+ else if(energy < 6.){
+ fhMCLambda0vsClusterMaxCellDiffE2[mcssEta]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE2[mcssEta] ->Fill(ncells, maxCellFraction);
+ }
+ else{
+ fhMCLambda0vsClusterMaxCellDiffE6[mcssEta]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE6[mcssEta] ->Fill(ncells, maxCellFraction);
+ }
+
+ }//eta
+ else {
+ fhMCELambda0[mcssOther] ->Fill(energy, lambda0);
+ fhMCEdLambda0[mcssOther] ->Fill(energy, lambda0/ncells2);
+ fhMCELambda1[mcssOther] ->Fill(energy, lambda1);
+ fhMCEdLambda1[mcssOther] ->Fill(energy, lambda1/ncells2);
+ fhMCEDispersion[mcssOther] ->Fill(energy, disp);
+ fhMCEdDispersion[mcssOther]->Fill(energy, disp/ncells2);
+ fhMCNCellsE[mcssOther] ->Fill(energy, ncells);
+ fhMCMaxCellDiffClusterE[mcssOther]->Fill(energy,maxCellFraction);
+ if (energy < 2.){
+ fhMCLambda0vsClusterMaxCellDiffE0[mcssOther]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE0[mcssOther] ->Fill(ncells, maxCellFraction);
+ }
+ else if(energy < 6.){
+ fhMCLambda0vsClusterMaxCellDiffE2[mcssOther]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE2[mcssOther] ->Fill(ncells, maxCellFraction);
+ }
+ else{
+ fhMCLambda0vsClusterMaxCellDiffE6[mcssOther]->Fill(lambda0, maxCellFraction);
+ fhMCNCellsvsClusterMaxCellDiffE6[mcssOther] ->Fill(ncells, maxCellFraction);
+ }
+
+ }//other particles
+
+ }//MC data
+
}
+//________________________________________________________________________
+TObjString * AliAnaPhoton::GetAnalysisCuts()
+{
+ //Save parameters used for analysis
+ TString parList ; //this will be list of parameters used for this analysis.
+ const Int_t buffersize = 255;
+ char onePar[buffersize] ;
+
+ snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f)\n",
+ fConvAsymCut, fConvDEtaCut, fConvDPhiMinCut, fConvDPhiMaxCut) ;
+ parList+=onePar ;
+
+ //Get parameters set in base class.
+ parList += GetBaseParametersList() ;
+
+ //Get parameters set in PID class.
+ parList += GetCaloPID()->GetPIDParametersList() ;
+
+ //Get parameters set in FiducialCut class (not available yet)
+ //parlist += GetFidCut()->GetFidCutParametersList()
+
+ return new TObjString(parList) ;
+}
//________________________________________________________________________
TList * AliAnaPhoton::GetCreateOutputObjects()
// store them in outputContainer
TList * outputContainer = new TList() ;
outputContainer->SetName("PhotonHistos") ;
+
+ Int_t nptbins = GetHistoPtBins(); Float_t ptmax = GetHistoPtMax(); Float_t ptmin = GetHistoPtMin();
+ Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin();
+ Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();
+ Int_t ssbins = GetHistoShowerShapeBins(); Float_t ssmax = GetHistoShowerShapeMax(); Float_t ssmin = GetHistoShowerShapeMin();
+ Int_t nbins = GetHistoNClusterCellBins(); Int_t nmax = GetHistoNClusterCellMax(); Int_t nmin = GetHistoNClusterCellMin();
+
+ fhNCellsE = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
+ fhNCellsE->SetXTitle("E (GeV)");
+ fhNCellsE->SetYTitle("# of cells in cluster");
+ outputContainer->Add(fhNCellsE);
- Int_t nptbins = GetHistoNPtBins();
- Int_t nphibins = GetHistoNPhiBins();
- Int_t netabins = GetHistoNEtaBins();
- Float_t ptmax = GetHistoPtMax();
- Float_t phimax = GetHistoPhiMax();
- Float_t etamax = GetHistoEtaMax();
- Float_t ptmin = GetHistoPtMin();
- Float_t phimin = GetHistoPhiMin();
- Float_t etamin = GetHistoEtaMin();
-
- //Histograms of highest Photon identified in Event
- fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
+ fhMaxCellDiffClusterE = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
+ nptbins,ptmin,ptmax, 500,0,1.);
+ fhMaxCellDiffClusterE->SetXTitle("E_{cluster} (GeV) ");
+ fhMaxCellDiffClusterE->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
+ outputContainer->Add(fhMaxCellDiffClusterE);
+
+ fhEPhoton = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
+ fhEPhoton->SetYTitle("N");
+ fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
+ outputContainer->Add(fhEPhoton) ;
+
+ fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax);
fhPtPhoton->SetYTitle("N");
fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
outputContainer->Add(fhPtPhoton) ;
fhPhiPhoton = new TH2F
- ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiPhoton->SetYTitle("#phi");
+ ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
+ fhPhiPhoton->SetYTitle("#phi (rad)");
fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
outputContainer->Add(fhPhiPhoton) ;
fhEtaPhoton = new TH2F
- ("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
+ ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
fhEtaPhoton->SetYTitle("#eta");
fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
outputContainer->Add(fhEtaPhoton) ;
- if(IsDataMC()){
-
- fhPtPrompt = new TH1F("hPtMCPrompt","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtPrompt->SetYTitle("N");
- fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtPrompt) ;
-
- fhPhiPrompt = new TH2F
- ("hPhiMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiPrompt->SetYTitle("#phi");
- fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiPrompt) ;
-
- fhEtaPrompt = new TH2F
- ("hEtaMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaPrompt->SetYTitle("#eta");
- fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaPrompt) ;
-
- fhPtFragmentation = new TH1F("hPtMCFragmentation","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtFragmentation->SetYTitle("N");
- fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtFragmentation) ;
-
- fhPhiFragmentation = new TH2F
- ("hPhiMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiFragmentation->SetYTitle("#phi");
- fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiFragmentation) ;
-
- fhEtaFragmentation = new TH2F
- ("hEtaMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaFragmentation->SetYTitle("#eta");
- fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaFragmentation) ;
-
- fhPtISR = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtISR->SetYTitle("N");
- fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtISR) ;
-
- fhPhiISR = new TH2F
- ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiISR->SetYTitle("#phi");
- fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiISR) ;
-
- fhEtaISR = new TH2F
- ("hEtaMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaISR->SetYTitle("#eta");
- fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaISR) ;
-
- fhPtPi0Decay = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtPi0Decay->SetYTitle("N");
- fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtPi0Decay) ;
-
- fhPhiPi0Decay = new TH2F
- ("hPhiMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiPi0Decay->SetYTitle("#phi");
- fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiPi0Decay) ;
-
- fhEtaPi0Decay = new TH2F
- ("hEtaMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaPi0Decay->SetYTitle("#eta");
- fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaPi0Decay) ;
-
- fhPtOtherDecay = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtOtherDecay->SetYTitle("N");
- fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtOtherDecay) ;
-
- fhPhiOtherDecay = new TH2F
- ("hPhiMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiOtherDecay->SetYTitle("#phi");
- fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiOtherDecay) ;
-
- fhEtaOtherDecay = new TH2F
- ("hEtaMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaOtherDecay->SetYTitle("#eta");
- fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaOtherDecay) ;
-
- fhPtConversion = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtConversion->SetYTitle("N");
- fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtConversion) ;
-
- fhPhiConversion = new TH2F
- ("hPhiMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiConversion->SetYTitle("#phi");
- fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiConversion) ;
-
- fhEtaConversion = new TH2F
- ("hEtaMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaConversion->SetYTitle("#eta");
- fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaConversion) ;
-
- fhPtUnknown = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtUnknown->SetYTitle("N");
- fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtUnknown) ;
-
- fhPhiUnknown = new TH2F
- ("hPhiMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiUnknown->SetYTitle("#phi");
- fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiUnknown) ;
-
- fhEtaUnknown = new TH2F
- ("hEtaMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaUnknown->SetYTitle("#eta");
- fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaUnknown) ;
-
- }//Histos with MC
-
- //Save parameters used for analysis
- TString parList ; //this will be list of parameters used for this analysis.
- char onePar[255] ;
+ fhEtaPhiPhoton = new TH2F
+ ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
+ fhEtaPhiPhoton->SetYTitle("#phi (rad)");
+ fhEtaPhiPhoton->SetXTitle("#eta");
+ outputContainer->Add(fhEtaPhiPhoton) ;
+ if(GetMinPt() < 0.5){
+ fhEtaPhi05Photon = new TH2F
+ ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
+ fhEtaPhi05Photon->SetYTitle("#phi (rad)");
+ fhEtaPhi05Photon->SetXTitle("#eta");
+ outputContainer->Add(fhEtaPhi05Photon) ;
+ }
- sprintf(onePar,"--- AliAnaPhoton ---\n") ;
- parList+=onePar ;
- sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
- parList+=onePar ;
- sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
- parList+=onePar ;
- sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
- parList+=onePar ;
- sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
- parList+=onePar ;
- sprintf(onePar,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
- parList+=onePar ;
+ //Conversion
+ if(fCheckConversion){
+
+ fhEtaPhiConversion = new TH2F
+ ("hEtaPhiConversion","#eta vs #phi for converted clusters",netabins,etamin,etamax,nphibins,phimin,phimax);
+ fhEtaPhiConversion->SetYTitle("#phi (rad)");
+ fhEtaPhiConversion->SetXTitle("#eta");
+ outputContainer->Add(fhEtaPhiConversion) ;
+ if(GetMinPt() < 0.5){
+ fhEtaPhi05Conversion = new TH2F
+ ("hEtaPhi05Conversion","#eta vs #phi, E > 0.5, for converted clusters",netabins,etamin,etamax,nphibins,phimin,phimax);
+ fhEtaPhi05Conversion->SetYTitle("#phi (rad)");
+ fhEtaPhi05Conversion->SetXTitle("#eta");
+ outputContainer->Add(fhEtaPhi05Conversion) ;
+ }
+
+ fhPtPhotonConv = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax);
+ fhPtPhotonConv->SetYTitle("N");
+ fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
+ outputContainer->Add(fhPtPhotonConv) ;
+
+ fhEtaPhiPhotonConv = new TH2F
+ ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
+ fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
+ fhEtaPhiPhotonConv->SetXTitle("#eta");
+ outputContainer->Add(fhEtaPhiPhotonConv) ;
+ if(GetMinPt() < 0.5){
+ fhEtaPhi05PhotonConv = new TH2F
+ ("hEtaPhi05PhotonConv","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
+ fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
+ fhEtaPhi05PhotonConv->SetXTitle("#eta");
+ outputContainer->Add(fhEtaPhi05PhotonConv) ;
+ }
+
+ fhConvDeltaEta = new TH2F
+ ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5);
+ fhConvDeltaEta->SetYTitle("#Delta #eta");
+ fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaEta) ;
+
+ fhConvDeltaPhi = new TH2F
+ ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5);
+ fhConvDeltaPhi->SetYTitle("#Delta #phi");
+ fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaPhi) ;
+
+ fhConvDeltaEtaPhi = new TH2F
+ ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5);
+ fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
+ fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
+ outputContainer->Add(fhConvDeltaEtaPhi) ;
+
+ fhConvAsym = new TH2F
+ ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1);
+ fhConvAsym->SetYTitle("Asymmetry");
+ fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvAsym) ;
+
+ fhConvPt = new TH2F
+ ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.);
+ fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
+ fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvPt) ;
+
+ fhConvDistEta = new TH2F
+ ("hConvDistEta","distance to conversion vertex",100,-0.7,0.7,100,0.,5.);
+ fhConvDistEta->SetXTitle("#eta");
+ fhConvDistEta->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistEta) ;
+
+ fhConvDistEn = new TH2F
+ ("hConvDistEn","distance to conversion vertex",nptbins,ptmin,ptmax,100,0.,5.);
+ fhConvDistEn->SetXTitle("E (GeV)");
+ fhConvDistEn->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistEn) ;
+
+ fhConvDistMass = new TH2F
+ ("hConvDistMass","distance to conversion vertex",100,0,fMassCut,100,0.,5.);
+ fhConvDistMass->SetXTitle("m (GeV/c^2)");
+ fhConvDistMass->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistMass) ;
+
+ fhConvDistEtaCutEta = new TH2F
+ ("hConvDistEtaCutEta","distance to conversion vertex, dEta < 0.05",100,-0.7,0.7,100,0.,5.);
+ fhConvDistEtaCutEta->SetXTitle("#eta");
+ fhConvDistEtaCutEta->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistEtaCutEta) ;
+
+ fhConvDistEnCutEta = new TH2F
+ ("hConvDistEnCutEta","distance to conversion vertex, dEta < 0.05",nptbins,ptmin,ptmax,100,0.,5.);
+ fhConvDistEnCutEta->SetXTitle("E (GeV)");
+ fhConvDistEnCutEta->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistEnCutEta) ;
+
+ fhConvDistMassCutEta = new TH2F
+ ("hConvDistMassCutEta","distance to conversion vertex, dEta < 0.05",100,0,fMassCut,100,0.,5.);
+ fhConvDistMassCutEta->SetXTitle("m (GeV/c^2)");
+ fhConvDistMassCutEta->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistMassCutEta) ;
+
+ fhConvDistEtaCutMass = new TH2F
+ ("hConvDistEtaCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",100,-0.7,0.7,100,0.,5.);
+ fhConvDistEtaCutMass->SetXTitle("#eta");
+ fhConvDistEtaCutMass->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistEtaCutMass) ;
+
+ fhConvDistEnCutMass = new TH2F
+ ("hConvDistEnCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",nptbins,ptmin,ptmax,100,0.,5.);
+ fhConvDistEnCutMass->SetXTitle("E (GeV)");
+ fhConvDistEnCutMass->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistEnCutMass) ;
+
+ fhConvDistEtaCutAsy = new TH2F
+ ("hConvDistEtaCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",100,-0.7,0.7,100,0.,5.);
+ fhConvDistEtaCutAsy->SetXTitle("#eta");
+ fhConvDistEtaCutAsy->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistEtaCutAsy) ;
+
+ fhConvDistEnCutAsy = new TH2F
+ ("hConvDistEnCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",nptbins,ptmin,ptmax,100,0.,5.);
+ fhConvDistEnCutAsy->SetXTitle("E (GeV)");
+ fhConvDistEnCutAsy->SetYTitle(" distance (m)");
+ outputContainer->Add(fhConvDistEnCutAsy) ;
+
+ } // check conversion
- //Get parameters set in base class.
- parList += GetBaseParametersList() ;
- //Get parameters set in PID class.
- parList += GetCaloPID()->GetPIDParametersList() ;
+ //Shower shape
+ if(fFillSSHistograms){
+
+ fhLam0E = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhLam0E->SetYTitle("#lambda_{0}^{2}");
+ fhLam0E->SetXTitle("E (GeV)");
+ outputContainer->Add(fhLam0E);
+
+ fhLam1E = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhLam1E->SetYTitle("#lambda_{1}^{2}");
+ fhLam1E->SetXTitle("E (GeV)");
+ outputContainer->Add(fhLam1E);
+
+ fhDispE = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhDispE->SetYTitle("D^{2}");
+ fhDispE->SetXTitle("E (GeV) ");
+ outputContainer->Add(fhDispE);
+
+ fhdLam0E = new TH2F ("hdLam0E","#lambda_{0}^{2}/N_{cells}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhdLam0E->SetYTitle("d#lambda_{0}^{2}");
+ fhdLam0E->SetXTitle("E (GeV)");
+ outputContainer->Add(fhdLam0E);
+
+ fhdLam1E = new TH2F ("hdLam1E","#lambda_{1}^{2}/N_{cells}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhdLam1E->SetYTitle("d#lambda_{1}^{2}");
+ fhdLam1E->SetXTitle("E (GeV)");
+ outputContainer->Add(fhdLam1E);
+
+ fhdDispE = new TH2F ("hdDispE"," dispersion^{2}/N_{cells}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhdDispE->SetYTitle("dD^{2}");
+ fhdDispE->SetXTitle("E (GeV) ");
+ outputContainer->Add(fhdDispE);
+
+
+ if(fCalorimeter == "EMCAL"){
+ fhLam0ETRD = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
+ fhLam0ETRD->SetXTitle("E (GeV)");
+ outputContainer->Add(fhLam0ETRD);
+
+ fhLam1ETRD = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
+ fhLam1ETRD->SetXTitle("E (GeV)");
+ outputContainer->Add(fhLam1ETRD);
+
+ fhDispETRD = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhDispETRD->SetYTitle("Dispersion^{2}");
+ fhDispETRD->SetXTitle("E (GeV) ");
+ outputContainer->Add(fhDispETRD);
+
+ fhdLam0ETRD = new TH2F ("hdLam0ETRD","#lambda_{0}^{2}/N_{cells}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhdLam0ETRD->SetYTitle("d#lambda_{0}^{2}");
+ fhdLam0ETRD->SetXTitle("E (GeV)");
+ outputContainer->Add(fhdLam0ETRD);
+
+ fhdLam1ETRD = new TH2F ("hdLam1ETRD","#lambda_{1}^{2}/N_{cells}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhdLam1ETRD->SetYTitle("d#lambda_{1}^{2}");
+ fhdLam1ETRD->SetXTitle("E (GeV)");
+ outputContainer->Add(fhdLam1ETRD);
+
+ fhdDispETRD = new TH2F ("hdDispETRD"," dispersion^{2}/N_{cells}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhdDispETRD->SetYTitle("dD^{2}");
+ fhdDispETRD->SetXTitle("E (GeV) ");
+ outputContainer->Add(fhdDispETRD);
+
+ }
+
+ fhNCellsLam0LowE = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
+ fhNCellsLam0LowE->SetXTitle("N_{cells}");
+ fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
+ outputContainer->Add(fhNCellsLam0LowE);
+
+ fhNCellsLam0HighE = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
+ fhNCellsLam0HighE->SetXTitle("N_{cells}");
+ fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
+ outputContainer->Add(fhNCellsLam0HighE);
+
+ fhNCellsLam1LowE = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
+ fhNCellsLam1LowE->SetXTitle("N_{cells}");
+ fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
+ outputContainer->Add(fhNCellsLam1LowE);
+
+ fhNCellsLam1HighE = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
+ fhNCellsLam1HighE->SetXTitle("N_{cells}");
+ fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
+ outputContainer->Add(fhNCellsLam1HighE);
+
+ fhNCellsDispLowE = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
+ fhNCellsDispLowE->SetXTitle("N_{cells}");
+ fhNCellsDispLowE->SetYTitle("D^{2}");
+ outputContainer->Add(fhNCellsDispLowE);
+
+ fhNCellsDispHighE = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
+ fhNCellsDispHighE->SetXTitle("N_{cells}");
+ fhNCellsDispHighE->SetYTitle("D^{2}");
+ outputContainer->Add(fhNCellsDispHighE);
+
+
+ fhNCellsdLam0LowE = new TH2F ("hNCellsdLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
+ fhNCellsdLam0LowE->SetXTitle("N_{cells}");
+ fhNCellsdLam0LowE->SetYTitle("#lambda_{0}^{2}");
+ outputContainer->Add(fhNCellsdLam0LowE);
+
+ fhNCellsdLam0HighE = new TH2F ("hNCellsdLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}/N_{cells}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
+ fhNCellsdLam0HighE->SetXTitle("N_{cells}");
+ fhNCellsdLam0HighE->SetYTitle("#lambda_{0}^{2}");
+ outputContainer->Add(fhNCellsdLam0HighE);
+
+ fhNCellsdLam1LowE = new TH2F ("hNCellsdLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
+ fhNCellsdLam1LowE->SetXTitle("N_{cells}");
+ fhNCellsdLam1LowE->SetYTitle("#lambda_{0}^{2}");
+ outputContainer->Add(fhNCellsdLam1LowE);
+
+ fhNCellsdLam1HighE = new TH2F ("hNCellsdLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}/N_{cells}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
+ fhNCellsdLam1HighE->SetXTitle("N_{cells}");
+ fhNCellsdLam1HighE->SetYTitle("#lambda_{0}^{2}");
+ outputContainer->Add(fhNCellsdLam1HighE);
+
+ fhNCellsdDispLowE = new TH2F ("hNCellsdDispLowE","N_{cells} in cluster vs dispersion^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
+ fhNCellsdDispLowE->SetXTitle("N_{cells}");
+ fhNCellsdDispLowE->SetYTitle("D^{2}");
+ outputContainer->Add(fhNCellsdDispLowE);
+
+ fhNCellsdDispHighE = new TH2F ("hNCellsdDispHighE","N_{cells} in cluster vs dispersion^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
+ fhNCellsdDispHighE->SetXTitle("N_{cells}");
+ fhNCellsdDispHighE->SetYTitle("D^{2}");
+ outputContainer->Add(fhNCellsdDispHighE);
+
+
+ fhEtaLam0LowE = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
+ fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
+ fhEtaLam0LowE->SetXTitle("#eta");
+ outputContainer->Add(fhEtaLam0LowE);
+
+ fhPhiLam0LowE = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
+ fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
+ fhPhiLam0LowE->SetXTitle("#phi");
+ outputContainer->Add(fhPhiLam0LowE);
+
+ fhEtaLam0HighE = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, E > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
+ fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
+ fhEtaLam0HighE->SetXTitle("#eta");
+ outputContainer->Add(fhEtaLam0HighE);
+
+ fhPhiLam0HighE = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
+ fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
+ fhPhiLam0HighE->SetXTitle("#phi");
+ outputContainer->Add(fhPhiLam0HighE);
+
+ fhLam1Lam0LowE = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
+ fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
+ fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
+ outputContainer->Add(fhLam1Lam0LowE);
+
+ fhLam1Lam0HighE = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
+ fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
+ fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
+ outputContainer->Add(fhLam1Lam0HighE);
+
+ fhLam0DispLowE = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
+ fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
+ fhLam0DispLowE->SetYTitle("D^{2}");
+ outputContainer->Add(fhLam0DispLowE);
+
+ fhLam0DispHighE = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
+ fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
+ fhLam0DispHighE->SetYTitle("D^{2}");
+ outputContainer->Add(fhLam0DispHighE);
+
+ fhDispLam1LowE = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
+ fhDispLam1LowE->SetXTitle("D^{2}");
+ fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
+ outputContainer->Add(fhDispLam1LowE);
+
+ fhDispLam1HighE = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of E > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
+ fhDispLam1HighE->SetXTitle("D^{2}");
+ fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
+ outputContainer->Add(fhDispLam1HighE);
+
+
+
+ fhEtadLam0LowE = new TH2F ("hEtadLam0LowE","#eta vs #lambda_{0}^{2}/N_{cells}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax/50);
+ fhEtadLam0LowE->SetYTitle("d#lambda_{0}^{2}");
+ fhEtadLam0LowE->SetXTitle("#eta");
+ outputContainer->Add(fhEtadLam0LowE);
+
+ fhPhidLam0LowE = new TH2F ("hPhidLam0LowE","#phi vs #lambda_{0}^{2}/N_{cells}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax/50);
+ fhPhidLam0LowE->SetYTitle("d#lambda_{0}^{2}");
+ fhPhidLam0LowE->SetXTitle("#phi");
+ outputContainer->Add(fhPhidLam0LowE);
+
+ fhEtadLam0HighE = new TH2F ("hEtadLam0HighE","#eta vs #lambda_{0}^{2}/N_{cells}^{2}", netabins,etamin,etamax, ssbins,ssmin,ssmax/50);
+ fhEtadLam0HighE->SetYTitle("d#lambda_{0}^{2}");
+ fhEtadLam0HighE->SetXTitle("#eta");
+ outputContainer->Add(fhEtadLam0HighE);
+
+ fhPhidLam0HighE = new TH2F ("hPhidLam0HighE","#phi vs #lambda_{0}^{2}/N_{cells}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax/50);
+ fhPhidLam0HighE->SetYTitle("d#lambda_{0}^{2}");
+ fhPhidLam0HighE->SetXTitle("#phi");
+ outputContainer->Add(fhPhidLam0HighE);
+
+ fhdLam1dLam0LowE = new TH2F ("hdLam1dLam0LowE","#lambda_{0}^{2}/N_{cells}^{2} vs #lambda_{1}^{2}/N_{cells}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
+ fhdLam1dLam0LowE->SetYTitle("d#lambda_{0}^{2}");
+ fhdLam1dLam0LowE->SetXTitle("d#lambda_{1}^{2}");
+ outputContainer->Add(fhdLam1dLam0LowE);
+
+ fhdLam1dLam0HighE = new TH2F ("hdLam1dLam0HighE","#lambda_{0}^{2}/N_{cells}^{2}vs #lambda_{1}^{2}/N_{cells}^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
+ fhdLam1dLam0HighE->SetYTitle("d#lambda_{0}^{2}");
+ fhdLam1dLam0HighE->SetXTitle("d#lambda_{1}^{2}");
+ outputContainer->Add(fhdLam1dLam0HighE);
+
+ fhdLam0dDispLowE = new TH2F ("hdLam0dDispLowE","#lambda_{0}^{2}/N_{cells}^{2} vs dispersion^{2}/N_{cells}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
+ fhdLam0dDispLowE->SetXTitle("d#lambda_{0}^{2}");
+ fhdLam0dDispLowE->SetYTitle("dD^{2}");
+ outputContainer->Add(fhdLam0dDispLowE);
+
+ fhdLam0dDispHighE = new TH2F ("hdLam0dDispHighE","#lambda_{0}^{2}/N_{cells}^{2} vs dispersion^{2}/N_{cells}^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
+ fhdLam0dDispHighE->SetXTitle("d#lambda_{0}^{2}");
+ fhdLam0dDispHighE->SetYTitle("dD^{2}");
+ outputContainer->Add(fhdLam0dDispHighE);
+
+ fhdDispdLam1LowE = new TH2F ("hdDispddLam1LowE","Dispersion^{2}/N_{cells}^{2} vs #lambda_{1}^{2}/N_{cells}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
+ fhdDispdLam1LowE->SetXTitle("dD^{2}");
+ fhdDispdLam1LowE->SetYTitle("d#lambda_{1}^{2}");
+ outputContainer->Add(fhdDispdLam1LowE);
+
+ fhdDispdLam1HighE = new TH2F ("hdDispdLam1HighE","Dispersion^{2}/N_{cells}^{2} vs #lambda_{1^{2}}/N_{cells}^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
+ fhdDispdLam1HighE->SetXTitle("dD^{2}");
+ fhdDispdLam1HighE->SetYTitle("d#lambda_{1}^{2}");
+ outputContainer->Add(fhdDispdLam1HighE);
+
+
+ } // Shower shape
- //Get parameters set in FidutialCut class (not available yet)
- //parlist += GetFidCut()->GetFidCutParametersList()
- TObjString *oString= new TObjString(parList) ;
- outputContainer->Add(oString);
+ if(IsDataMC()){
+ fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
+ fhDeltaE->SetXTitle("#Delta E (GeV)");
+ outputContainer->Add(fhDeltaE);
+
+ fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
+ fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
+ outputContainer->Add(fhDeltaPt);
+
+ fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
+ fhRatioE->SetXTitle("E_{reco}/E_{gen}");
+ outputContainer->Add(fhRatioE);
+
+ fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
+ fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
+ outputContainer->Add(fhRatioPt);
+
+ fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+ fh2E->SetXTitle("E_{rec} (GeV)");
+ fh2E->SetYTitle("E_{gen} (GeV)");
+ outputContainer->Add(fh2E);
+
+ fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+ fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
+ fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
+ outputContainer->Add(fh2Pt);
+
+ TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
+ "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P",
+ "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","String" } ;
+
+ TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
+ "Conversion", "Hadron", "AntiNeutron","AntiProton",
+ "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
+
+ for(Int_t i = 0; i < fNOriginHistograms; i++){
+
+ fhMCE[i] = new TH1F(Form("hE_MC%s",pname[i].Data()),
+ Form("cluster from %s : E ",ptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhMCE[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCE[i]) ;
+
+ fhPtMC[i] = new TH1F(Form("hPt_MC%s",pname[i].Data()),
+ Form("cluster from %s : p_{T} ",ptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhPtMC[i]->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhPtMC[i]) ;
+
+ fhEtaMC[i] = new TH2F(Form("hEta_MC%s",pname[i].Data()),
+ Form("cluster from %s : #eta ",ptype[i].Data()),
+ nptbins,ptmin,ptmax,netabins,etamin,etamax);
+ fhEtaMC[i]->SetYTitle("#eta");
+ fhEtaMC[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEtaMC[i]) ;
+
+ fhPhiMC[i] = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
+ Form("cluster from %s : #phi ",ptype[i].Data()),
+ nptbins,ptmin,ptmax,nphibins,phimin,phimax);
+ fhPhiMC[i]->SetYTitle("#phi (rad)");
+ fhPhiMC[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhPhiMC[i]) ;
+
+ }
+
+ TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}","hadron?",
+ "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ;
+
+ TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Hadron",
+ "PhotonPrompt","PhotonFragmentation","PhotonISR"} ;
+
+ for(Int_t i = 0; i < fNPrimaryHistograms; i++){
+ fhEPrimMC[i] = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
+ Form("primary photon %s : E ",pptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhEPrimMC[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEPrimMC[i]) ;
+
+ fhPtPrimMC[i] = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
+ Form("primary photon %s : p_{T} ",pptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhPtPrimMC[i]->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhPtPrimMC[i]) ;
+
+ fhYPrimMC[i] = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
+ Form("primary photon %s : Rapidity ",pptype[i].Data()),
+ nptbins,ptmin,ptmax,800,-8,8);
+ fhYPrimMC[i]->SetYTitle("Rapidity");
+ fhYPrimMC[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhYPrimMC[i]) ;
+
+ fhPhiPrimMC[i] = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
+ Form("primary photon %s : #phi ",pptype[i].Data()),
+ nptbins,ptmin,ptmax,nphibins,phimin,phimax);
+ fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
+ fhPhiPrimMC[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhPhiPrimMC[i]) ;
+
+
+ fhEPrimMCAcc[i] = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
+ Form("primary photon %s in acceptance: E ",pptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhEPrimMCAcc[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEPrimMCAcc[i]) ;
+
+ fhPtPrimMCAcc[i] = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
+ Form("primary photon %s in acceptance: p_{T} ",pptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhPtPrimMCAcc[i]->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhPtPrimMCAcc[i]) ;
+
+ fhYPrimMCAcc[i] = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
+ Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
+ nptbins,ptmin,ptmax,100,-1,1);
+ fhYPrimMCAcc[i]->SetYTitle("Rapidity");
+ fhYPrimMCAcc[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhYPrimMCAcc[i]) ;
+
+ fhPhiPrimMCAcc[i] = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
+ Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
+ nptbins,ptmin,ptmax,nphibins,phimin,phimax);
+ fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
+ fhPhiPrimMCAcc[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhPhiPrimMCAcc[i]) ;
+
+ }
+
+ if(fCheckConversion){
+ fhPtConversionTagged = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
+ fhPtConversionTagged->SetYTitle("N");
+ fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
+ outputContainer->Add(fhPtConversionTagged) ;
+
+
+ fhPtAntiNeutronTagged = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
+ fhPtAntiNeutronTagged->SetYTitle("N");
+ fhPtAntiNeutronTagged->SetXTitle("p_{T #gamma}(GeV/c)");
+ outputContainer->Add(fhPtAntiNeutronTagged) ;
+
+ fhPtAntiProtonTagged = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
+ fhPtAntiProtonTagged->SetYTitle("N");
+ fhPtAntiProtonTagged->SetXTitle("p_{T #gamma}(GeV/c)");
+ outputContainer->Add(fhPtAntiProtonTagged) ;
+
+ fhPtUnknownTagged = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
+ fhPtUnknownTagged->SetYTitle("N");
+ fhPtUnknownTagged->SetXTitle("p_{T #gamma}(GeV/c)");
+ outputContainer->Add(fhPtUnknownTagged) ;
+
+ fhConvDeltaEtaMCConversion = new TH2F
+ ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5);
+ fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
+ fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaEtaMCConversion) ;
+
+ fhConvDeltaPhiMCConversion = new TH2F
+ ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5);
+ fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
+ fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaPhiMCConversion) ;
+
+ fhConvDeltaEtaPhiMCConversion = new TH2F
+ ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5);
+ fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
+ fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
+ outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
+
+ fhConvAsymMCConversion = new TH2F
+ ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1);
+ fhConvAsymMCConversion->SetYTitle("Asymmetry");
+ fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvAsymMCConversion) ;
+
+ fhConvPtMCConversion = new TH2F
+ ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.);
+ fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
+ fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvPtMCConversion) ;
+
+ fhConvDispersionMCConversion = new TH2F
+ ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.);
+ fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
+ fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
+ outputContainer->Add(fhConvDispersionMCConversion) ;
+
+ fhConvM02MCConversion = new TH2F
+ ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
+ fhConvM02MCConversion->SetYTitle("M02 cluster 1");
+ fhConvM02MCConversion->SetXTitle("M02 cluster 2");
+ outputContainer->Add(fhConvM02MCConversion) ;
+
+ fhConvDeltaEtaMCAntiNeutron = new TH2F
+ ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5);
+ fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
+ fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
+
+ fhConvDeltaPhiMCAntiNeutron = new TH2F
+ ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5);
+ fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
+ fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
+
+ fhConvDeltaEtaPhiMCAntiNeutron = new TH2F
+ ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
+ fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
+ fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
+ outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;
+
+ fhConvAsymMCAntiNeutron = new TH2F
+ ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1);
+ fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
+ fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvAsymMCAntiNeutron) ;
+
+ fhConvPtMCAntiNeutron = new TH2F
+ ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.);
+ fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
+ fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvPtMCAntiNeutron) ;
+
+ fhConvDispersionMCAntiNeutron = new TH2F
+ ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.);
+ fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
+ fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
+ outputContainer->Add(fhConvDispersionMCAntiNeutron) ;
+
+ fhConvM02MCAntiNeutron = new TH2F
+ ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
+ fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
+ fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
+ outputContainer->Add(fhConvM02MCAntiNeutron) ;
+
+ fhConvDeltaEtaMCAntiProton = new TH2F
+ ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5);
+ fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
+ fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
+
+ fhConvDeltaPhiMCAntiProton = new TH2F
+ ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5);
+ fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
+ fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
+
+ fhConvDeltaEtaPhiMCAntiProton = new TH2F
+ ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
+ fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
+ fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
+ outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;
+
+ fhConvAsymMCAntiProton = new TH2F
+ ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1);
+ fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
+ fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvAsymMCAntiProton) ;
+
+ fhConvPtMCAntiProton = new TH2F
+ ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.);
+ fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
+ fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvPtMCAntiProton) ;
+
+ fhConvDispersionMCAntiProton = new TH2F
+ ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.);
+ fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
+ fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
+ outputContainer->Add(fhConvDispersionMCAntiProton) ;
+
+ fhConvM02MCAntiProton = new TH2F
+ ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
+ fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
+ fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
+ outputContainer->Add(fhConvM02MCAntiProton) ;
+
+ fhConvDeltaEtaMCString = new TH2F
+ ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5);
+ fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
+ fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaEtaMCString) ;
+
+ fhConvDeltaPhiMCString = new TH2F
+ ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5);
+ fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
+ fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvDeltaPhiMCString) ;
+
+ fhConvDeltaEtaPhiMCString = new TH2F
+ ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5);
+ fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
+ fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
+ outputContainer->Add(fhConvDeltaEtaPhiMCString) ;
+
+ fhConvAsymMCString = new TH2F
+ ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1);
+ fhConvAsymMCString->SetYTitle("Asymmetry");
+ fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvAsymMCString) ;
+
+ fhConvPtMCString = new TH2F
+ ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.);
+ fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
+ fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
+ outputContainer->Add(fhConvPtMCString) ;
+
+ fhConvDispersionMCString = new TH2F
+ ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
+ fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
+ fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
+ outputContainer->Add(fhConvDispersionMCString) ;
+
+ fhConvM02MCString = new TH2F
+ ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
+ fhConvM02MCString->SetYTitle("M02 cluster 1");
+ fhConvM02MCString->SetXTitle("M02 cluster 2");
+ outputContainer->Add(fhConvM02MCString) ;
+
+ fhConvDistMCConversion = new TH2F
+ ("hConvDistMCConversion","calculated conversion distance vs real vertes for MC conversion",100,0.,5.,100,0.,5.);
+ fhConvDistMCConversion->SetYTitle("distance");
+ fhConvDistMCConversion->SetXTitle("vertex R");
+ outputContainer->Add(fhConvDistMCConversion) ;
+
+ fhConvDistMCConversionCuts = new TH2F
+ ("hConvDistMCConversionCuts","calculated conversion distance vs real vertes for MC conversion, deta < 0.05, m < 10 MeV, asym < 0.1",100,0.,5.,100,0.,5.);
+ fhConvDistMCConversionCuts->SetYTitle("distance");
+ fhConvDistMCConversionCuts->SetXTitle("vertex R");
+ outputContainer->Add(fhConvDistMCConversionCuts) ;
+ }
+
+ if(fFillSSHistograms){
+
+ TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
+
+ TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
+
+ for(Int_t i = 0; i < 6; i++){
+
+ fhMCELambda0[i] = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
+ fhMCELambda0[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCELambda0[i]) ;
+
+ fhMCEdLambda0[i] = new TH2F(Form("hEdLambda0_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : E vs #lambda_{0}^{2}/N_{cells}^{2}",ptypess[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhMCEdLambda0[i]->SetYTitle("d#lambda_{0}^{2}");
+ fhMCEdLambda0[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCEdLambda0[i]) ;
+
+ fhMCELambda1[i] = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
+ fhMCELambda1[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCELambda1[i]) ;
+
+ fhMCEdLambda1[i] = new TH2F(Form("hEdLambda1_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : E vs d#lambda_{1}^{2}/N_{cells}^{2}",ptypess[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhMCEdLambda1[i]->SetYTitle("d#lambda_{1}^{2}");
+ fhMCEdLambda1[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCEdLambda1[i]) ;
+
+ fhMCEDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhMCEDispersion[i]->SetYTitle("D^{2}");
+ fhMCEDispersion[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCEDispersion[i]) ;
+
+ fhMCEdDispersion[i] = new TH2F(Form("hEdDispersion_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : E vs dispersion^{2}/N_{cells}^{2}",ptypess[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhMCEdDispersion[i]->SetYTitle("dD^{2}");
+ fhMCEdDispersion[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCEdDispersion[i]) ;
+
+ fhMCNCellsE[i] = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
+ Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
+ nptbins,ptmin,ptmax, nbins,nmin,nmax);
+ fhMCNCellsE[i]->SetXTitle("E (GeV)");
+ fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
+ outputContainer->Add(fhMCNCellsE[i]);
+
+ fhMCMaxCellDiffClusterE[i] = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
+ Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
+ nptbins,ptmin,ptmax, 500,0,1.);
+ fhMCMaxCellDiffClusterE[i]->SetXTitle("E_{cluster} (GeV) ");
+ fhMCMaxCellDiffClusterE[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
+ outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
+
+ fhMCLambda0vsClusterMaxCellDiffE0[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
+ ssbins,ssmin,ssmax,500,0,1.);
+ fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
+ fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
+ outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
+
+ fhMCLambda0vsClusterMaxCellDiffE2[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
+ ssbins,ssmin,ssmax,500,0,1.);
+ fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
+ fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
+ outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
+
+ fhMCLambda0vsClusterMaxCellDiffE6[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
+ ssbins,ssmin,ssmax,500,0,1.);
+ fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
+ fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
+ outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
+
+ fhMCNCellsvsClusterMaxCellDiffE0[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
+ nbins/5,nmin,nmax/5,500,0,1.);
+ fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
+ fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
+ outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
+
+ fhMCNCellsvsClusterMaxCellDiffE2[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
+ nbins/5,nmin,nmax/5,500,0,1.);
+ fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
+ fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
+ outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
+
+ fhMCNCellsvsClusterMaxCellDiffE6[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
+ Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
+ nbins/5,nmin,nmax/5,500,0,1.);
+ fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
+ fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("E (GeV)");
+ outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
+
+ }// loop
+
+ if(!GetReader()->IsEmbeddedClusterSelectionOn())
+ {
+ fhMCPhotonELambda0NoOverlap = new TH2F("hELambda0_MCPhoton_NoOverlap",
+ "cluster from Photon : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
+ fhMCPhotonELambda0NoOverlap->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
+
+ fhMCPhotonEdLambda0NoOverlap = new TH2F("hEdLambda0_MCPhoton_NoOverlap",
+ "cluster from Photon : E vs #lambda_{0}^{2}/N_{cells}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhMCPhotonEdLambda0NoOverlap->SetYTitle("d#lambda_{0}^{2}");
+ fhMCPhotonEdLambda0NoOverlap->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPhotonEdLambda0NoOverlap) ;
+
+ fhMCPhotonELambda0TwoOverlap = new TH2F("hELambda0_MCPhoton_TwoOverlap",
+ "cluster from Photon : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
+ fhMCPhotonELambda0TwoOverlap->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
+
+ fhMCPhotonEdLambda0TwoOverlap = new TH2F("hEdLambda0_MCPhoton_TwoOverlap",
+ "cluster from Photon : E vs #lambda_{0}^{2}/N_{cells}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhMCPhotonEdLambda0TwoOverlap->SetYTitle("d#lambda_{0}^{2}");
+ fhMCPhotonEdLambda0TwoOverlap->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPhotonEdLambda0TwoOverlap) ;
+
+ fhMCPhotonELambda0NOverlap = new TH2F("hELambda0_MCPhoton_NOverlap",
+ "cluster from Photon : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
+ fhMCPhotonELambda0NOverlap->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
+
+ fhMCPhotonEdLambda0NOverlap = new TH2F("hEdLambda0_MCPhoton_NOverlap",
+ "cluster from Photon : E vs #lambda_{0}^{2}/N_{cells}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
+ fhMCPhotonEdLambda0NOverlap->SetYTitle("d#lambda_{0}^{2}");
+ fhMCPhotonEdLambda0NOverlap->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCPhotonEdLambda0NOverlap) ;
+ } //No embedding
+
+ //Fill histograms to check shape of embedded clusters
+ if(GetReader()->IsEmbeddedClusterSelectionOn())
+ {
+
+ fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
+ "Energy Fraction of embedded signal versus cluster energy",
+ nptbins,ptmin,ptmax,100,0.,1.);
+ fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
+ fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
+
+ fhEmbedPhotonELambda0FullSignal = new TH2F("hELambda0_EmbedPhoton_FullSignal",
+ "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPhotonELambda0FullSignal->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
+
+ fhEmbedPhotonEdLambda0FullSignal = new TH2F("hEdLambda0_EmbedPhoton_FullSignal",
+ "cluster from Photon with more than 90% energy in cluster: E vs d#lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPhotonEdLambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPhotonEdLambda0FullSignal->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPhotonEdLambda0FullSignal) ;
+
+
+ fhEmbedPhotonELambda0MostlySignal = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
+ "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPhotonELambda0MostlySignal->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
+
+ fhEmbedPhotonEdLambda0MostlySignal = new TH2F("hEdLambda0_EmbedPhoton_MostlySignal",
+ "cluster from Photon with 50% to 90% energy in cluster: E vs d#lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPhotonEdLambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPhotonEdLambda0MostlySignal->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPhotonEdLambda0MostlySignal) ;
+
+
+ fhEmbedPhotonELambda0MostlyBkg = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
+ "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPhotonELambda0MostlyBkg->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
+
+ fhEmbedPhotonEdLambda0MostlyBkg = new TH2F("hEdLambda0_EmbedPhoton_MostlyBkg",
+ "cluster from Photon with 10% to 50% energy in cluster: E vs d#lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPhotonEdLambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPhotonEdLambda0MostlyBkg->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPhotonEdLambda0MostlyBkg) ;
+
+
+ fhEmbedPhotonELambda0FullBkg = new TH2F("hELambda0_EmbedPhoton_FullBkg",
+ "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPhotonELambda0FullBkg->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
+
+ fhEmbedPhotonEdLambda0FullBkg = new TH2F("hEdLambda0_EmbedPhoton_FullBkg",
+ "cluster from Photon with 0% to 10% energy in cluster: E vs d#lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPhotonEdLambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPhotonEdLambda0FullBkg->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPhotonEdLambda0FullBkg) ;
+
+
+ fhEmbedPi0ELambda0FullSignal = new TH2F("hELambda0_EmbedPi0_FullSignal",
+ "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPi0ELambda0FullSignal->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
+
+ fhEmbedPi0EdLambda0FullSignal = new TH2F("hEdLambda0_EmbedPi0_FullSignal",
+ "cluster from Pi0 with more than 90% energy in cluster: E vs d#lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPi0EdLambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPi0EdLambda0FullSignal->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPi0EdLambda0FullSignal) ;
+
+
+ fhEmbedPi0ELambda0MostlySignal = new TH2F("hELambda0_EmbedPi0_MostlySignal",
+ "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPi0ELambda0MostlySignal->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
+
+ fhEmbedPi0EdLambda0MostlySignal = new TH2F("hEdLambda0_EmbedPi0_MostlySignal",
+ "cluster from Pi0 with 50% to 90% energy in cluster: E vs d#lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPi0EdLambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPi0EdLambda0MostlySignal->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPi0EdLambda0MostlySignal) ;
+
+
+ fhEmbedPi0ELambda0MostlyBkg = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
+ "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPi0ELambda0MostlyBkg->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
+
+ fhEmbedPi0EdLambda0MostlyBkg = new TH2F("hEdLambda0_EmbedPi0_MostlyBkg",
+ "cluster from Pi0 with 10% to 50% energy in cluster: E vs d#lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPi0EdLambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPi0EdLambda0MostlyBkg->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPi0EdLambda0MostlyBkg) ;
+
+
+ fhEmbedPi0ELambda0FullBkg = new TH2F("hELambda0_EmbedPi0_FullBkg",
+ "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPi0ELambda0FullBkg->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
+
+ fhEmbedPi0EdLambda0FullBkg = new TH2F("hEdLambda0_EmbedPi0_FullBkg",
+ "cluster from Pi0 with 0% to 10% energy in cluster: E vs d#lambda_{0}^{2}",
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEmbedPi0EdLambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
+ fhEmbedPi0EdLambda0FullBkg->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEmbedPi0EdLambda0FullBkg) ;
+
+ }// embedded histograms
+
+
+ }// Fill SS MC histograms
+
+ }//Histos with MC
+
+ //Store calo PID histograms
+ if(fRejectTrackMatch){
+ TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
+ for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
+ outputContainer->Add(caloPIDHistos->At(i)) ;
+ }
+ delete caloPIDHistos;
+ }
return outputContainer ;
//Init
//Do some checks
- if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){
- printf("AliAnaPhoton::Init() - !!ABORT: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
+ if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
+ printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
abort();
}
- else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
- printf("AliAnaPhoton::Init() - !!ABORT: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
+ else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
+ printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
abort();
}
}
-
//____________________________________________________________________________
void AliAnaPhoton::InitParameters()
{
//Initialize the parameters of the analysis.
- SetOutputAODClassName("AliAODPWG4Particle");
- SetOutputAODName("PWG4Particle");
-
AddToHistogramsName("AnaPhoton_");
-
- fCalorimeter = "PHOS" ;
- fMinDist = 2.;
- fMinDist2 = 4.;
- fMinDist3 = 5.;
- fRejectTrackMatch = kTRUE ;
+
+ fCalorimeter = "EMCAL" ;
+ fMinDist = 2.;
+ fMinDist2 = 4.;
+ fMinDist3 = 5.;
+ fMassCut = 0.03; //30 MeV
+
+ fTimeCutMin = -1;
+ fTimeCutMax = 9999999;
+ fNCellsCut = 0;
+
+ fRejectTrackMatch = kTRUE ;
+ fCheckConversion = kFALSE;
+ fRemoveConvertedPair = kFALSE;
+ fAddConvertedPairsToAOD = kFALSE;
+
}
//__________________________________________________________________
void AliAnaPhoton::MakeAnalysisFillAOD()
{
- //Do analysis and fill aods
- //Search for photons in fCalorimeter
- TRefArray * pl = new TRefArray();
+ //Do photon analysis and fill aods
+
+ //Get the vertex
+ Double_t v[3] = {0,0,0}; //vertex ;
+ GetReader()->GetVertex(v);
- //Get vertex for photon momentum calculation
- Double_t vertex[]={0,0,0} ; //vertex ;
- if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
//Select the Calorimeter of the photon
+ TObjArray * pl = 0x0;
if(fCalorimeter == "PHOS")
- pl = GetAODPHOS();
+ pl = GetPHOSClusters();
else if (fCalorimeter == "EMCAL")
- pl = GetAODEMCAL();
-
- //Fill AODCaloClusters and AODParticle with PHOS aods
- TLorentzVector mom ;
-
- for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
- AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));
- //Cluster selection, not charged, with photon id and in fidutial cut
- //Get Momentum vector,
- calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
- //If too small or big pt, skip it
- if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
- //Check acceptance selection
- if(IsFidutialCutOn()){
- Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;
- if(! in ) continue ;
+ pl = GetEMCALClusters();
+
+ if(!pl) {
+ Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
+ return;
+ }
+
+ //Init arrays, variables, get number of clusters
+ TLorentzVector mom, mom2 ;
+ Int_t nCaloClusters = pl->GetEntriesFast();
+ //List to be used in conversion analysis, to tag the cluster as candidate for conversion
+ Bool_t * indexConverted = 0x0;
+ if(fCheckConversion){
+ indexConverted = new Bool_t[nCaloClusters];
+ for (Int_t i = 0; i < nCaloClusters; i++)
+ indexConverted[i] = kFALSE;
+ }
+
+ if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
+
+ //----------------------------------------------------
+ // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
+ //----------------------------------------------------
+ // Loop on clusters
+ for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
+
+ AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
+ //printf("calo %d, %f\n",icalo,calo->E());
+
+ //Get the index where the cluster comes, to retrieve the corresponding vertex
+ Int_t evtIndex = 0 ;
+ if (GetMixedEvent()) {
+ evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
+ //Get the vertex and check it is not too large in z
+ if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
}
+ //Cluster selection, not charged, with photon id and in fiducial cut
+ if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
+ calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
+ else{
+ Double_t vertex[]={0,0,0};
+ calo->GetMomentum(mom,vertex) ;
+ }
+
+ //--------------------------------------
+ // Cluster selection
+ //--------------------------------------
+ if(!ClusterSelected(calo,mom)) continue;
+
+ //----------------------------
//Create AOD for analysis
+ //----------------------------
AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
- aodph.SetLabel(calo->GetLabel(0));
- //printf("Index %d, Id %d\n",icalo, calo->GetID());
- //Set the indeces of the original caloclusters
+
+ //...............................................
+ //Set the indeces of the original caloclusters (MC, ID), and calorimeter
+ Int_t label = calo->GetLabel();
+ aodph.SetLabel(label);
aodph.SetCaloLabel(calo->GetID(),-1);
aodph.SetDetector(fCalorimeter);
- if(GetDebug() > 1)
- printf("AliAnaPhoton::MakeAnalysisFillAOD() - Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());
+ //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
- //Check Distance to Bad channel, set bit.
- Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
- if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
- if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
- continue ;
+ //...............................................
+ //Set bad channel distance bit
+ Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
+ if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
+ else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
+ else aodph.SetDistToBad(0) ;
+ //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
- if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
+ //--------------------------------------------------------------------------------------
+ //Play with the MC stack if available
+ //--------------------------------------------------------------------------------------
- if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
- else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
- else aodph.SetDistToBad(0) ;
+ //Check origin of the candidates
+ if(IsDataMC()){
+ aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
+
+ if(GetDebug() > 0)
+ printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
+ }//Work with stack also
+
+ //--------------------------------------------------------------------------------------
+ //Fill some shower shape histograms before PID is applied
+ //--------------------------------------------------------------------------------------
+
+ FillShowerShapeHistograms(calo,aodph.GetTag());
- //Check PID
+ //-------------------------------------
//PID selection or bit setting
+ //-------------------------------------
+ // MC
if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
//Get most probable PID, check PID weights (in MC this option is mandatory)
- aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
- if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
+ aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
+ if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
//If primary is not photon, skip it.
- if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
- }
+ if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
+ }
+ //...............................................
+ // Data, PID check on
else if(IsCaloPIDOn()){
- //Skip matched clusters with tracks
- if(fRejectTrackMatch && calo->GetNTracksMatched() > 0) continue ;
-
//Get most probable PID, 2 options check PID weights
- //or redo PID, recommended option for EMCal.
+ //or redo PID, recommended option for MCEal.
if(!IsCaloPIDRecalculationOn())
- aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
+ aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
else
- aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
+ aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
- if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
+ if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
//If cluster does not pass pid, not photon, skip it.
- if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
+ if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
}
+ //...............................................
+ // Data, PID check off
else{
//Set PID bits for later selection (AliAnaPi0 for example)
//GetPDG already called in SetPIDBits.
- GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
+ GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
}
- if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
+ if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
- //Play with the MC stack if available
- //Check origin of the candidates
- if(IsDataMC()){
- aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
- if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate %d\n",aodph.GetTag());
- }//Work with stack also
+ //--------------------------------------------------------------------------------------
+ // Conversions pairs analysis
+ // Check if cluster comes from a conversion in the material in front of the calorimeter
+ // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
+ //--------------------------------------------------------------------------------------
+
+ // Do analysis only if there are more than one cluster
+ if( nCaloClusters > 1 && fCheckConversion){
+ Bool_t bConverted = kFALSE;
+ Int_t id2 = -1;
+
+ //Check if set previously as converted couple, if so skip its use.
+ if (indexConverted[icalo]) continue;
+
+ // Second cluster loop
+ for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
+ //Check if set previously as converted couple, if so skip its use.
+ if (indexConverted[jcalo]) continue;
+ //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
+ AliVCluster * calo2 = (AliVCluster*) (pl->At(jcalo)); //Get cluster kinematics
+
+
+ //Mixed event, get index of event
+ Int_t evtIndex2 = 0 ;
+ if (GetMixedEvent()) {
+ evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
+
+ }
+
+ //Get kinematics of second cluster
+ if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
+ calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
+ else{
+ Double_t vertex[]={0,0,0};
+ calo2->GetMomentum(mom2,vertex) ;
+ }
+
+ //--------------------------------------
+ // Cluster selection
+ //--------------------------------------
+
+ if(!ClusterSelected(calo2,mom2)) continue;
+
+ //................................................
+ // Get TOF of each cluster in pair, calculate difference if small,
+ // take this pair. Way to reject clusters from hadrons (or pileup?)
+ Double_t t12diff = calo2->GetTOF()-calo->GetTOF()*1e9;
+ if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
+
+ //................................................
+ //Get mass of pair, if small, take this pair.
+ Float_t pairM = (mom+mom2).M();
+ //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
+ if(pairM < fMassCut){
+ aodph.SetTagged(kFALSE);
+ id2 = calo2->GetID();
+ indexConverted[icalo]=kTRUE;
+ indexConverted[jcalo]=kTRUE;
+ Float_t asymmetry = TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E());
+ Float_t dPhi = mom.Phi()-mom2.Phi();
+ Float_t dEta = mom.Eta()-mom2.Eta();
+
+ //...............................................
+ //Fill few histograms with kinematics of the pair
+ //FIXME, move all this to MakeAnalysisFillHistograms ...
+
+ fhConvDeltaEta ->Fill( pairM, dPhi );
+ fhConvDeltaPhi ->Fill( pairM, dEta );
+ fhConvAsym ->Fill( pairM, asymmetry );
+ fhConvDeltaEtaPhi->Fill( dEta , dPhi );
+ fhConvPt ->Fill( pairM, (mom+mom2).Pt());
+
+ //Estimate conversion distance, T. Awes, M. Ivanov
+ //Under the assumption that the pair has zero mass, and that each electron
+ //of the pair has the same momentum, they will each have the same bend radius
+ //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
+ //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
+ //R = E/1.5 [cm]. Under these assumptions, the distance from the conversion
+ //point to the MCEal can be related to the separation distance, L=2y, on the MCEal
+ //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
+ //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
+ //the clusters.
+ Float_t pos1[3];
+ calo->GetPosition(pos1);
+ Float_t pos2[3];
+ calo2->GetPosition(pos2);
+ Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
+ (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
+ (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
+
+ Float_t convDist = TMath::Sqrt(mom.E() *clustDist*0.01/0.15);
+ Float_t convDist2 = TMath::Sqrt(mom2.E()*clustDist*0.01/0.15);
+ //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,mom.E(),convDist,mom2.E(),convDist2);
+ if(GetDebug() > 2)
+ printf("AliAnaPhoton::MakeAnalysisFillAOD(): Pair with mass %2.3f < %2.3f, %1.2f < dPhi %2.2f < %2.2f, dEta %f < %2.2f, asymmetry %2.2f< %2.2f; \n cluster1 id %d, e %2.3f SM %d, eta %2.3f, phi %2.3f ; \n cluster2 id %d, e %2.3f, SM %d,eta %2.3f, phi %2.3f\n",
+ pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
+ calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
+ id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
+
+ fhConvDistEta ->Fill(mom .Eta(), convDist );
+ fhConvDistEta ->Fill(mom2.Eta(), convDist2);
+ fhConvDistEn ->Fill(mom .E(), convDist );
+ fhConvDistEn ->Fill(mom2.E(), convDist2);
+ fhConvDistMass->Fill((mom+mom2).M(), convDist );
+ //dEta cut
+ if(dEta<0.05){
+ fhConvDistEtaCutEta ->Fill(mom .Eta(), convDist );
+ fhConvDistEtaCutEta ->Fill(mom2.Eta(), convDist2);
+ fhConvDistEnCutEta ->Fill(mom .E(), convDist );
+ fhConvDistEnCutEta ->Fill(mom2.E(), convDist2);
+ fhConvDistMassCutEta->Fill((mom+mom2).M(), convDist );
+ //mass cut
+ if(pairM<0.01){//10 MeV
+ fhConvDistEtaCutMass ->Fill(mom .Eta(), convDist );
+ fhConvDistEtaCutMass ->Fill(mom2.Eta(), convDist2);
+ fhConvDistEnCutMass ->Fill(mom .E(), convDist );
+ fhConvDistEnCutMass ->Fill(mom2.E(), convDist2);
+ // asymmetry cut
+ if(asymmetry<0.1){
+ fhConvDistEtaCutAsy ->Fill(mom .Eta(), convDist );
+ fhConvDistEtaCutAsy ->Fill(mom2.Eta(), convDist2);
+ fhConvDistEnCutAsy ->Fill(mom .E(), convDist );
+ fhConvDistEnCutAsy ->Fill(mom2.E(), convDist2);
+ }//asymmetry cut
+ }//mass cut
+ }//dEta cut
+
+ //...............................................
+ //Select pairs in a eta-phi window
+ if(TMath::Abs(dEta) < fConvDEtaCut &&
+ TMath::Abs(dPhi) < fConvDPhiMaxCut &&
+ TMath::Abs(dPhi) > fConvDPhiMinCut &&
+ asymmetry < fConvAsymCut ){
+ bConverted = kTRUE;
+ }
+ //printf("Accepted? %d\n",bConverted);
+ //...........................................
+ //Fill more histograms, simulated data
+ //FIXME, move all this to MakeAnalysisFillHistograms ...
+ if(IsDataMC()){
+
+ //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
+ Int_t ancPDG = 0;
+ Int_t ancStatus = 0;
+ TLorentzVector momentum;
+ TVector3 prodVertex;
+ Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(),
+ GetReader(), ancPDG, ancStatus, momentum, prodVertex);
+
+ // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
+ // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
+
+ Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabels(),calo2->GetNLabels(),GetReader(), 0);
+ if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCConversion)){
+ if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
+ fhConvDeltaEtaMCConversion ->Fill( pairM, dEta );
+ fhConvDeltaPhiMCConversion ->Fill( pairM, dPhi );
+ fhConvAsymMCConversion ->Fill( pairM, asymmetry );
+ fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi );
+ fhConvPtMCConversion ->Fill( pairM, (mom+mom2).Pt());
+ fhConvDispersionMCConversion ->Fill( calo->GetDispersion(), calo2->GetDispersion());
+ fhConvM02MCConversion ->Fill( calo->GetM02(), calo2->GetM02());
+ fhConvDistMCConversion ->Fill( convDist , prodVertex.Mag() );
+ fhConvDistMCConversion ->Fill( convDist2, prodVertex.Mag() );
+
+ if(dEta<0.05 && pairM<0.01 && asymmetry<0.1){
+ fhConvDistMCConversionCuts->Fill( convDist , prodVertex.Mag() );
+ fhConvDistMCConversionCuts->Fill( convDist2, prodVertex.Mag() );
+ }
+
+ }
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiNeutron)){
+ if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
+ fhConvDeltaEtaMCAntiNeutron ->Fill( pairM, dEta );
+ fhConvDeltaPhiMCAntiNeutron ->Fill( pairM, dPhi );
+ fhConvAsymMCAntiNeutron ->Fill( pairM, asymmetry );
+ fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi );
+ fhConvPtMCAntiNeutron ->Fill( pairM, (mom+mom2).Pt());
+ fhConvDispersionMCAntiNeutron ->Fill( calo->GetDispersion(), calo2->GetDispersion());
+ fhConvM02MCAntiNeutron ->Fill( calo->GetM02(), calo2->GetM02());
+ }
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiProton)){
+ if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
+ fhConvDeltaEtaMCAntiProton ->Fill( pairM, dEta );
+ fhConvDeltaPhiMCAntiProton ->Fill( pairM, dPhi );
+ fhConvAsymMCAntiProton ->Fill( pairM, asymmetry );
+ fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi );
+ fhConvPtMCAntiProton ->Fill( pairM, (mom+mom2).Pt());
+ fhConvDispersionMCAntiProton ->Fill( calo->GetDispersion(), calo2->GetDispersion());
+ fhConvM02MCAntiProton ->Fill( calo->GetM02(), calo2->GetM02());
+ }
+ }
+
+ //Pairs coming from fragmenting pairs.
+ if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
+ fhConvDeltaEtaMCString ->Fill( pairM, dPhi);
+ fhConvDeltaPhiMCString ->Fill( pairM, dPhi);
+ fhConvAsymMCString ->Fill( pairM, TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
+ fhConvDeltaEtaPhiMCString ->Fill( dPhi, dPhi );
+ fhConvPtMCString ->Fill( pairM, (mom+mom2).Pt());
+ fhConvDispersionMCString ->Fill( calo->GetDispersion(), calo2->GetDispersion());
+ fhConvM02MCString ->Fill( calo->GetM02(), calo2->GetM02());
+ }
+
+ }// Data MC
+
+ break;
+ }
+
+ }//Mass loop
+
+ //..........................................................................................................
+ //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
+ if(bConverted){
+ //Add to AOD
+ if(fAddConvertedPairsToAOD){
+ //Create AOD of pair analysis
+ TLorentzVector mpair = mom+mom2;
+ AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
+ aodpair.SetLabel(aodph.GetLabel());
+ //aodpair.SetInputFileIndex(input);
+
+ //printf("Index %d, Id %d\n",icalo, calo->GetID());
+ //Set the indeces of the original caloclusters
+ aodpair.SetCaloLabel(calo->GetID(),id2);
+ aodpair.SetDetector(fCalorimeter);
+ aodpair.SetIdentifiedParticleType(aodph.GetIdentifiedParticleType());
+ aodpair.SetTag(aodph.GetTag());
+ aodpair.SetTagged(kTRUE);
+ //Add AOD with pair object to aod branch
+ AddAODParticle(aodpair);
+ //printf("\t \t both added pair\n");
+ }
+
+ //Do not add the current calocluster
+ if(fRemoveConvertedPair) continue;
+ else {
+ //printf("TAGGED\n");
+ //Tag this cluster as likely conversion
+ aodph.SetTagged(kTRUE);
+ }
+ }//converted pair
+ }//check conversion
+ //printf("\t \t added single cluster %d\n",icalo);
+
+ //FIXME, this to MakeAnalysisFillHistograms ...
+ Int_t absID = 0;
+ Float_t maxCellFraction = 0;
+ AliVCaloCells* cells = 0;
+
+ if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
+ else cells = GetPHOSCells();
+
+ absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
+
+ fhMaxCellDiffClusterE->Fill(aodph.E(),maxCellFraction);
+ fhNCellsE ->Fill(aodph.E(),calo->GetNCells());
+
//Add AOD with photon object to aod branch
AddAODParticle(aodph);
}//loop
- if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs \n");
+ delete [] indexConverted;
+
+ if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
}
//__________________________________________________________________
void AliAnaPhoton::MakeAnalysisFillHistograms()
{
- //Do analysis and fill histograms
+ //Fill histograms
+
+ //-------------------------------------------------------------------
+ // Access MC information in stack if requested, check that it exists.
+ AliStack * stack = 0x0;
+ TParticle * primary = 0x0;
+ TClonesArray * mcparticles = 0x0;
+ AliAODMCParticle * aodprimary = 0x0;
+
+ if(IsDataMC()){
+
+ if(GetReader()->ReadStack()){
+ stack = GetMCStack() ;
+ if(!stack) {
+ printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
+ abort();
+ }
+
+ }
+ else if(GetReader()->ReadAODMCParticles()){
+
+ //Get the list of MC particles
+ mcparticles = GetReader()->GetAODMCParticles(0);
+ if(!mcparticles && GetDebug() > 0) {
+ printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
+ }
+ }
+ }// is data and MC
+
- //Get vertex for photon momentum calculation
- //Double_t v[]={0,0,0} ; //vertex ;
- //if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC)
- //GetReader()->GetVertex(v);
+ // Get vertex
+ Double_t v[3] = {0,0,0}; //vertex ;
+ GetReader()->GetVertex(v);
+ //fhVertex->Fill(v[0],v[1],v[2]);
+ if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
+ //----------------------------------
//Loop on stored AOD photons
Int_t naod = GetOutputAODBranch()->GetEntriesFast();
if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
for(Int_t iaod = 0; iaod < naod ; iaod++){
AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
- Int_t pdg = ph->GetPdg();
+ Int_t pdg = ph->GetIdentifiedParticleType();
if(GetDebug() > 3)
- printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
+ printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
//If PID used, fill histos with photons in Calorimeter fCalorimeter
- if(IsCaloPIDOn())
- if(pdg != AliCaloPID::kPhoton) continue;
+ if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
if(ph->GetDetector() != fCalorimeter) continue;
- if(GetDebug() > 1)
+ if(GetDebug() > 2)
printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
+ //................................
//Fill photon histograms
- Float_t ptcluster = ph->Pt();
+ Float_t ptcluster = ph->Pt();
Float_t phicluster = ph->Phi();
Float_t etacluster = ph->Eta();
+ Float_t ecluster = ph->E();
- fhPtPhoton ->Fill(ptcluster);
+ fhEPhoton ->Fill(ecluster);
+ fhPtPhoton ->Fill(ptcluster);
fhPhiPhoton ->Fill(ptcluster,phicluster);
- fhEtaPhoton ->Fill(ptcluster,etacluster);
+ fhEtaPhoton ->Fill(ptcluster,etacluster);
+ if (ecluster > 0.5) fhEtaPhiPhoton ->Fill(etacluster, phicluster);
+ else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
+
+ if(fCheckConversion &&ph->IsTagged()){
+ fhPtPhotonConv->Fill(ptcluster);
+ if(ecluster > 0.5) fhEtaPhiPhotonConv ->Fill(etacluster, phicluster);
+ else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
+ }
+ //.......................................
+ //Play with the MC data if available
if(IsDataMC()){
+
+ FillAcceptanceHistograms();
+
Int_t tag =ph->GetTag();
- if(tag == AliMCAnalysisUtils::kMCPrompt){
- fhPtPrompt ->Fill(ptcluster);
- fhPhiPrompt ->Fill(ptcluster,phicluster);
- fhEtaPrompt ->Fill(ptcluster,etacluster);
+
+ if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[mcPhoton])
+ {
+ fhMCE [mcPhoton] ->Fill(ecluster);
+ fhPtMC [mcPhoton] ->Fill(ptcluster);
+ fhPhiMC[mcPhoton] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcPhoton] ->Fill(ecluster,etacluster);
+
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[mcConversion])
+ {
+ fhMCE [mcConversion] ->Fill(ecluster);
+ fhPtMC [mcConversion] ->Fill(ptcluster);
+ fhPhiMC[mcConversion] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcConversion] ->Fill(ecluster,etacluster);
+
+ if(fCheckConversion){
+ if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
+ if(ptcluster > 0.5)fhEtaPhiConversion ->Fill(etacluster,phicluster);
+ else fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
+ }
+ }
+
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhMCE[mcPrompt]){
+ fhMCE [mcPrompt] ->Fill(ecluster);
+ fhPtMC [mcPrompt] ->Fill(ptcluster);
+ fhPhiMC[mcPrompt] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcPrompt] ->Fill(ecluster,etacluster);
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)&& fhMCE[mcFragmentation])
+ {
+ fhMCE [mcFragmentation] ->Fill(ecluster);
+ fhPtMC [mcFragmentation] ->Fill(ptcluster);
+ fhPhiMC[mcFragmentation] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcFragmentation] ->Fill(ecluster,etacluster);
+
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[mcISR])
+ {
+ fhMCE [mcISR] ->Fill(ecluster);
+ fhPtMC [mcISR] ->Fill(ptcluster);
+ fhPhiMC[mcISR] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcISR] ->Fill(ecluster,etacluster);
+ }
+ else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
+ !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[mcPi0Decay])
+ {
+ fhMCE [mcPi0Decay] ->Fill(ecluster);
+ fhPtMC [mcPi0Decay] ->Fill(ptcluster);
+ fhPhiMC[mcPi0Decay] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcPi0Decay] ->Fill(ecluster,etacluster);
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
+ GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) && fhMCE[mcOtherDecay])
+ {
+ fhMCE [mcOtherDecay] ->Fill(ecluster);
+ fhPtMC [mcOtherDecay] ->Fill(ptcluster);
+ fhPhiMC[mcOtherDecay] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcOtherDecay] ->Fill(ecluster,etacluster);
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [mcPi0])
+ {
+ fhMCE [mcPi0] ->Fill(ecluster);
+ fhPtMC [mcPi0] ->Fill(ptcluster);
+ fhPhiMC[mcPi0] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcPi0] ->Fill(ecluster,etacluster);
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[mcEta])
+ {
+ fhMCE [mcEta] ->Fill(ecluster);
+ fhPtMC [mcEta] ->Fill(ptcluster);
+ fhPhiMC[mcEta] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcEta] ->Fill(ecluster,etacluster);
+ }
}
- else if(tag==AliMCAnalysisUtils::kMCFragmentation)
- {
- fhPtFragmentation ->Fill(ptcluster);
- fhPhiFragmentation ->Fill(ptcluster,phicluster);
- fhEtaFragmentation ->Fill(ptcluster,etacluster);
- }
- else if(tag==AliMCAnalysisUtils::kMCISR)
- {
- fhPtISR ->Fill(ptcluster);
- fhPhiISR ->Fill(ptcluster,phicluster);
- fhEtaISR ->Fill(ptcluster,etacluster);
- }
- else if(tag==AliMCAnalysisUtils::kMCPi0Decay)
- {
- fhPtPi0Decay ->Fill(ptcluster);
- fhPhiPi0Decay ->Fill(ptcluster,phicluster);
- fhEtaPi0Decay ->Fill(ptcluster,etacluster);
- }
- else if(tag==AliMCAnalysisUtils::kMCEtaDecay || tag==AliMCAnalysisUtils::kMCOtherDecay)
- {
- fhPtOtherDecay ->Fill(ptcluster);
- fhPhiOtherDecay ->Fill(ptcluster,phicluster);
- fhEtaOtherDecay ->Fill(ptcluster,etacluster);
- }
- else if(tag==AliMCAnalysisUtils::kMCConversion)
- {
- fhPtConversion ->Fill(ptcluster);
- fhPhiConversion ->Fill(ptcluster,phicluster);
- fhEtaConversion ->Fill(ptcluster,etacluster);
-
- }
- else{
- fhPtUnknown ->Fill(ptcluster);
- fhPhiUnknown ->Fill(ptcluster,phicluster);
- fhEtaUnknown ->Fill(ptcluster,etacluster);
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[mcAntiNeutron])
+ {
+ fhMCE [mcAntiNeutron] ->Fill(ecluster);
+ fhPtMC [mcAntiNeutron] ->Fill(ptcluster);
+ fhPhiMC[mcAntiNeutron] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcAntiNeutron] ->Fill(ecluster,etacluster);
+ if(ph->IsTagged() && fCheckConversion) fhPtAntiNeutronTagged ->Fill(ptcluster);
+
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[mcAntiProton])
+ {
+ fhMCE [mcAntiProton] ->Fill(ecluster);
+ fhPtMC [mcAntiProton] ->Fill(ptcluster);
+ fhPhiMC[mcAntiProton] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcAntiProton] ->Fill(ecluster,etacluster);
+ if(ph->IsTagged() && fCheckConversion) fhPtAntiProtonTagged ->Fill(ptcluster);
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[mcElectron])
+ {
+ fhMCE [mcElectron] ->Fill(ecluster);
+ fhPtMC [mcElectron] ->Fill(ptcluster);
+ fhPhiMC[mcElectron] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcElectron] ->Fill(ecluster,etacluster);
+ }
+ else if( fhMCE[mcOther]){
+ fhMCE [mcOther] ->Fill(ecluster);
+ fhPtMC [mcOther] ->Fill(ptcluster);
+ fhPhiMC[mcOther] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcOther] ->Fill(ecluster,etacluster);
+ if(ph->IsTagged() && fCheckConversion) fhPtUnknownTagged ->Fill(ptcluster);
+
+
+ // printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
+ // ph->GetLabel(),ph->Pt());
+ // for(Int_t i = 0; i < 20; i++) {
+ // if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
+ // }
+ // printf("\n");
+
}
+
+ //....................................................................
+ // Access MC information in stack if requested, check that it exists.
+ Int_t label =ph->GetLabel();
+ if(label < 0) {
+ if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
+ continue;
+ }
+
+ Float_t eprim = 0;
+ Float_t ptprim = 0;
+ if(GetReader()->ReadStack()){
+
+ if(label >= stack->GetNtrack()) {
+ if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
+ continue ;
+ }
+
+ primary = stack->Particle(label);
+ if(!primary){
+ printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
+ continue;
+ }
+ eprim = primary->Energy();
+ ptprim = primary->Pt();
+
+ }
+ else if(GetReader()->ReadAODMCParticles()){
+ //Check which is the input
+ if(ph->GetInputFileIndex() == 0){
+ if(!mcparticles) continue;
+ if(label >= mcparticles->GetEntriesFast()) {
+ if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
+ label, mcparticles->GetEntriesFast());
+ continue ;
+ }
+ //Get the particle
+ aodprimary = (AliAODMCParticle*) mcparticles->At(label);
+
+ }
+
+ if(!aodprimary){
+ printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
+ continue;
+ }
+
+ eprim = aodprimary->E();
+ ptprim = aodprimary->Pt();
+
+ }
+
+ fh2E ->Fill(ecluster, eprim);
+ fh2Pt ->Fill(ptcluster, ptprim);
+ fhDeltaE ->Fill(eprim-ecluster);
+ fhDeltaPt->Fill(ptprim-ptcluster);
+ if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
+ if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
+
}//Histograms with MC
}// aod loop
printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
-
+ printf("Check Pair Conversion = %d\n",fCheckConversion);
+ printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
+ printf("Conversion pair mass cut = %f\n",fMassCut);
+ printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
+ fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
+ printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
+ printf("Number of cells in cluster is > %d \n", fNCellsCut);
printf(" \n") ;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff