ClassImp(AliAnaPhoton)
//____________________________________________________________________________
- AliAnaPhoton::AliAnaPhoton() :
- AliAnaPartCorrBaseClass(), fCalorimeter(""),
- fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
- fTimeCutMin(-1), fTimeCutMax(9999999), fNCellsCut(0),
- fCheckConversion(kFALSE), fRemoveConvertedPair(kFALSE), fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
- fConvAsymCut(1.), fConvDEtaCut(2.),fConvDPhiMinCut(-1.), fConvDPhiMaxCut(7.),
- //fhVertex(0),
- fhNtraNclu(0), fhNCellsPt(0),
- fhEPhoton(0), fhPtPhoton(0), fhPhiPhoton(0), fhEtaPhoton(0), fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
- fhPtPhotonConv(0), fhEtaPhiPhotonConv(0),fhEtaPhi05PhotonConv(0),
- fhConvDeltaEta(0), fhConvDeltaPhi(0), fhConvDeltaEtaPhi(0), fhConvAsym(0), fhConvPt(0),
- //MC
- fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
- fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),
- 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),fhEtaPhiConversion(0),fhEtaPhi05Conversion(0),
- fhPtAntiNeutron(0), fhPhiAntiNeutron(0), fhEtaAntiNeutron(0),
- fhPtAntiProton(0), fhPhiAntiProton(0), fhEtaAntiProton(0),
- fhPtUnknown(0), fhPhiUnknown(0), fhEtaUnknown(0),
- fhPtConversionTagged(0), fhPtAntiNeutronTagged(0), fhPtAntiProtonTagged(0), fhPtUnknownTagged(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)
+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();
}
+//__________________________________________________________________
+Bool_t AliAnaPhoton::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
+{
+ //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;
+
+}
+
+//_____________________________________________________________
+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
+}
+
+//__________________________________________________________________
+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()
{
return new TObjString(parList) ;
}
-
//________________________________________________________________________
TList * AliAnaPhoton::GetCreateOutputObjects()
{
TList * outputContainer = new TList() ;
outputContainer->SetName("PhotonHistos") ;
- Int_t nptbins = GetHistoPtBins();
- Int_t nphibins = GetHistoPhiBins();
- Int_t netabins = GetHistoEtaBins();
- 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
-// fhVertex = new TH3D ("Vertex","vertex position", 20,-10.,10., 20,-10.,10., 80,-40.,40.);
-// fhVertex->SetXTitle("X");
-// fhVertex->SetYTitle("Y");
-// fhVertex->SetZTitle("Z");
-// outputContainer->Add(fhVertex);
-
- fhNtraNclu = new TH2F ("hNtracksNcluster","# of tracks vs # of clusters", 500,0,500, 500,0,500);
- fhNtraNclu->SetXTitle("# of tracks");
- fhNtraNclu->SetYTitle("# of clusters");
- outputContainer->Add(fhNtraNclu);
-
- fhNCellsPt = new TH2F ("hNCellsPt","# of cells in cluster vs E of clusters", nptbins,ptmin, ptmax, 100,0,100);
- fhNCellsPt->SetXTitle("p_{T} (GeV/c)");
- fhNCellsPt->SetYTitle("# of cells in cluster");
- outputContainer->Add(fhNCellsPt);
+ 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);
+
+ 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");
//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)");
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
+
+
+ //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
+
if(IsDataMC()){
fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
outputContainer->Add(fh2Pt);
- fhPtMCPhoton = new TH1F("hPtMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtMCPhoton->SetYTitle("N");
- fhPtMCPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtMCPhoton) ;
-
- fhPhiMCPhoton = new TH2F
- ("hPhiMCPhoton","#phi_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiMCPhoton->SetYTitle("#phi");
- fhPhiMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiMCPhoton) ;
-
- fhEtaMCPhoton = new TH2F
- ("hEtaMCPhoton","#eta_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaMCPhoton->SetYTitle("#eta");
- fhEtaMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaMCPhoton) ;
-
- fhPtPrompt = new TH1F("hPtMCPrompt","Number of prompt #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}, prompt #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiPrompt->SetYTitle("#phi");
- fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiPrompt) ;
-
- fhEtaPrompt = new TH2F
- ("hEtaMCPrompt","#eta_{#gamma}, prompt #gamma in MC",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 fragmentation #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}, fragmentation #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiFragmentation->SetYTitle("#phi");
- fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiFragmentation) ;
-
- fhEtaFragmentation = new TH2F
- ("hEtaMCFragmentation","#eta_{#gamma}, fragmentation #gamma in MC",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","#eta_{#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("hPtMCPi0Decay","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}, #pi^{0} decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiPi0Decay->SetYTitle("#phi");
- fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiPi0Decay) ;
-
- fhEtaPi0Decay = new TH2F
- ("hEtaMCPi0Decay","#eta_{#gamma}, #pi^{0} #gamma in MC",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}, other decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiOtherDecay->SetYTitle("#phi");
- fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiOtherDecay) ;
-
- fhEtaOtherDecay = new TH2F
- ("hEtaMCOtherDecay","#eta_{#gamma}, other decay #gamma in MC",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}, conversion #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiConversion->SetYTitle("#phi");
- fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiConversion) ;
-
- fhEtaConversion = new TH2F
- ("hEtaMCConversion","#eta_{#gamma}, conversion #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaConversion->SetYTitle("#eta");
- fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaConversion) ;
+ 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" } ;
- fhEtaPhiConversion = new TH2F
- ("hEtaPhiConversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
- fhEtaPhiConversion->SetYTitle("#phi (rad)");
- fhEtaPhiConversion->SetXTitle("#eta");
- outputContainer->Add(fhEtaPhiConversion) ;
+ TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
+ "Conversion", "Hadron", "AntiNeutron","AntiProton",
+ "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
- fhEtaPhi05Conversion = new TH2F
- ("hEtaPhi05Conversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
- fhEtaPhi05Conversion->SetYTitle("#phi (rad)");
- fhEtaPhi05Conversion->SetXTitle("#eta");
- outputContainer->Add(fhEtaPhi05Conversion) ;
-
- fhPtAntiNeutron = new TH1F("hPtMCAntiNeutron","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtAntiNeutron->SetYTitle("N");
- fhPtAntiNeutron->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtAntiNeutron) ;
-
- fhPhiAntiNeutron = new TH2F
- ("hPhiMCAntiNeutron","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiAntiNeutron->SetYTitle("#phi");
- fhPhiAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiAntiNeutron) ;
-
- fhEtaAntiNeutron = new TH2F
- ("hEtaMCAntiNeutron","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaAntiNeutron->SetYTitle("#eta");
- fhEtaAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaAntiNeutron) ;
-
- fhPtAntiProton = new TH1F("hPtMCAntiProton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
- fhPtAntiProton->SetYTitle("N");
- fhPtAntiProton->SetXTitle("p_{T #gamma}(GeV/c)");
- outputContainer->Add(fhPtAntiProton) ;
-
- fhPhiAntiProton = new TH2F
- ("hPhiMCAntiProton","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiAntiProton->SetYTitle("#phi");
- fhPhiAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiAntiProton) ;
-
- fhEtaAntiProton = new TH2F
- ("hEtaMCAntiProton","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaAntiProton->SetYTitle("#eta");
- fhEtaAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaAntiProton) ;
-
- 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}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiUnknown->SetYTitle("#phi");
- fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhPhiUnknown) ;
-
- fhEtaUnknown = new TH2F
- ("hEtaMCUnknown","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaUnknown->SetYTitle("#eta");
- fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
- outputContainer->Add(fhEtaUnknown) ;
-
+ 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");
("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) ;
+ 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 ;
}
}
-
//____________________________________________________________________________
void AliAnaPhoton::InitParameters()
{
//Initialize the parameters of the analysis.
AddToHistogramsName("AnaPhoton_");
-
+
fCalorimeter = "EMCAL" ;
fMinDist = 2.;
fMinDist2 = 4.;
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();
}
if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
-
+
//----------------------------------------------------
// Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
//----------------------------------------------------
//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
-
- //Input from second AOD?
- //Int_t input = 0;
- // if (fCalorimeter == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= icalo)
- // input = 1 ;
- // else if(fCalorimeter == "PHOS" && GetReader()->GetPHOSClustersNormalInputEntries() <= icalo)
- // input = 1;
-
- //Get Momentum vector,
- //if (input == 0)
+
+ //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) ;
}
-
- // else if(input == 1)
- // calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
//--------------------------------------
// Cluster selection
//--------------------------------------
- if(GetDebug() > 2)
- printf("AliAnaPhoton::MakeAnalysisFillAOD() 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 pt, skip it
- if(mom.E() < GetMinPt() || mom.E() > GetMaxPt() ) continue ;
- if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",icalo);
-
- //.......................................
- // TOF cut, BE CAREFUL WITH THIS CUT
- Double_t tof = calo->GetTOF()*1e9;
- if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
- if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",icalo);
-
- //.......................................
- if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
- if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",icalo);
-
- //.......................................
- //Check acceptance selection
- if(IsFiducialCutOn()){
- Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
- if(! in ) continue ;
- }
- if(GetDebug() > 2) printf("Fiducial cut passed \n");
-
- //.......................................
- //Skip matched clusters with tracks
- if(fRejectTrackMatch){
- if(IsTrackMatched(calo)) {
- if(GetDebug() > 2) printf("\t Reject matched clusters\n");
- continue ;
- }
- else
- if(GetDebug() > 2) printf(" matching cut passed 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 )
- continue ;
- }
- else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
-
- if(GetDebug() > 0)
- printf("AliAnaPhoton::MakeAnalysisFillAOD() 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());
-
+ if(!ClusterSelected(calo,mom)) continue;
//----------------------------
//Create AOD for analysis
//Set the indeces of the original caloclusters (MC, ID), and calorimeter
Int_t label = calo->GetLabel();
aodph.SetLabel(label);
- //aodph.SetInputFileIndex(input);
aodph.SetCaloLabel(calo->GetID(),-1);
aodph.SetDetector(fCalorimeter);
- //printf("Index %d, Id %d\n",icalo, calo->GetID());
-
- //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
-
+ //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
+
//...............................................
//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());
- //...............................................
- //Set number of cells in this cluster
- //Temporary patch FIXME
- aodph.SetBtag(calo->GetNCells());
- // MEFIX
+ //--------------------------------------------------------------------------------------
+ //Play with the MC stack if available
+ //--------------------------------------------------------------------------------------
+
+ //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());
//-------------------------------------
//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->GetPID(),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()){
//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->GetPID(),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 ;
}
//...............................................
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->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
//--------------------------------------------------------------------------------------
// 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;
//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
+ AliVCluster * calo2 = (AliVCluster*) (pl->At(jcalo)); //Get cluster kinematics
+
//Mixed event, get index of event
Int_t evtIndex2 = 0 ;
calo2->GetMomentum(mom2,vertex) ;
}
- //Check only certain regions
- Bool_t in2 = kTRUE;
- if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
- if(!in2) continue;
+ //--------------------------------------
+ // 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.
aodph.SetTagged(kFALSE);
id2 = calo2->GetID();
indexConverted[icalo]=kTRUE;
- indexConverted[jcalo]=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();
-
- 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());
+ 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());
+ 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
TMath::Abs(dPhi) < fConvDPhiMaxCut &&
TMath::Abs(dPhi) > fConvDPhiMinCut &&
asymmetry < fConvAsymCut ){
- bConverted = kTRUE;
+ bConverted = kTRUE;
}
//printf("Accepted? %d\n",bConverted);
//...........................................
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);
+ 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);
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)){
}
}// Data MC
-
+
break;
}
//Set the indeces of the original caloclusters
aodpair.SetCaloLabel(calo->GetID(),id2);
aodpair.SetDetector(fCalorimeter);
- aodpair.SetPdg(aodph.GetPdg());
+ aodpair.SetIdentifiedParticleType(aodph.GetIdentifiedParticleType());
aodpair.SetTag(aodph.GetTag());
aodpair.SetTagged(kTRUE);
//Add AOD with pair object to aod branch
//printf("\t \t added single cluster %d\n",icalo);
//FIXME, this to MakeAnalysisFillHistograms ...
- fhNCellsPt->Fill(aodph.Pt(),calo->GetNCells());
+ 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);
//Fill histograms
//-------------------------------------------------------------------
- // Access MC information in stack if requested, check that it exists.
- AliStack * stack = 0x0;
- TParticle * primary = 0x0;
- TClonesArray * mcparticles0 = 0x0;
- //TClonesArray * mcparticles1 = 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
- mcparticles0 = GetReader()->GetAODMCParticles(0);
- if(!mcparticles0 && GetDebug() > 0) {
- printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
- }
- // if(GetReader()->GetSecondInputAODTree()){
- // mcparticles1 = GetReader()->GetAODMCParticles(1);
- // if(!mcparticles1 && GetDebug() > 0) {
- // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
- // }
- // }
-
- }
- }// is data and MC
-
+ // 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
Double_t v[3] = {0,0,0}; //vertex ;
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();
- fhNtraNclu->Fill(GetReader()->GetTrackMultiplicity(), naod);
- 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();
-
- if(GetDebug() > 3)
- printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
-
- //If PID used, fill histos with photons in Calorimeter fCalorimeter
- if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
- if(ph->GetDetector() != fCalorimeter) continue;
-
- if(GetDebug() > 2)
- printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
-
+ //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->GetIdentifiedParticleType();
+
+ if(GetDebug() > 3)
+ 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() && pdg != AliCaloPID::kPhoton) continue;
+ if(ph->GetDetector() != fCalorimeter) continue;
+
+ 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 phicluster = ph->Phi();
- Float_t etacluster = ph->Eta();
- Float_t ecluster = ph->E();
-
+ //Fill photon histograms
+ Float_t ptcluster = ph->Pt();
+ Float_t phicluster = ph->Phi();
+ Float_t etacluster = ph->Eta();
+ Float_t ecluster = ph->E();
+
fhEPhoton ->Fill(ecluster);
- fhPtPhoton ->Fill(ptcluster);
- fhPhiPhoton ->Fill(ptcluster,phicluster);
- fhEtaPhoton ->Fill(ptcluster,etacluster);
- if(ecluster > 0.5) fhEtaPhiPhoton ->Fill(etacluster, phicluster);
+ fhPtPhoton ->Fill(ptcluster);
+ fhPhiPhoton ->Fill(ptcluster,phicluster);
+ 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()){
}
//.......................................
- //Play with the MC data if available
- if(IsDataMC()){
-
- Int_t tag =ph->GetTag();
-
- if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
+ //Play with the MC data if available
+ if(IsDataMC()){
+
+ FillAcceptanceHistograms();
+
+ Int_t tag =ph->GetTag();
+
+ if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[mcPhoton])
{
- fhPtMCPhoton ->Fill(ptcluster);
- fhPhiMCPhoton ->Fill(ptcluster,phicluster);
- fhEtaMCPhoton ->Fill(ptcluster,etacluster);
+ fhMCE [mcPhoton] ->Fill(ecluster);
+ fhPtMC [mcPhoton] ->Fill(ptcluster);
+ fhPhiMC[mcPhoton] ->Fill(ecluster,phicluster);
+ fhEtaMC[mcPhoton] ->Fill(ecluster,etacluster);
- if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[mcConversion])
{
- fhPtConversion ->Fill(ptcluster);
- fhPhiConversion ->Fill(ptcluster,phicluster);
- fhEtaConversion ->Fill(ptcluster,etacluster);
- if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
- if(ptcluster > 0.5)fhEtaPhiConversion ->Fill(etacluster,phicluster);
- else fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
+ 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)){
- fhPtPrompt ->Fill(ptcluster);
- fhPhiPrompt ->Fill(ptcluster,phicluster);
- fhEtaPrompt ->Fill(ptcluster,etacluster);
+ 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))
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)&& fhMCE[mcFragmentation])
{
- fhPtFragmentation ->Fill(ptcluster);
- fhPhiFragmentation ->Fill(ptcluster,phicluster);
- fhEtaFragmentation ->Fill(ptcluster,etacluster);
+ 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))
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[mcISR])
{
- fhPtISR ->Fill(ptcluster);
- fhPhiISR ->Fill(ptcluster,phicluster);
- fhEtaISR ->Fill(ptcluster,etacluster);
+ 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))
+ else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
+ !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[mcPi0Decay])
{
- fhPtPi0Decay ->Fill(ptcluster);
- fhPhiPi0Decay ->Fill(ptcluster,phicluster);
- fhEtaPi0Decay ->Fill(ptcluster,etacluster);
+ 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))
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
+ GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) && fhMCE[mcOtherDecay])
{
- fhPtOtherDecay ->Fill(ptcluster);
- fhPhiOtherDecay ->Fill(ptcluster,phicluster);
- fhEtaOtherDecay ->Fill(ptcluster,etacluster);
+ 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(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[mcAntiNeutron])
{
- fhPtAntiNeutron ->Fill(ptcluster);
- fhPhiAntiNeutron ->Fill(ptcluster,phicluster);
- fhEtaAntiNeutron ->Fill(ptcluster,etacluster);
+ 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))
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[mcAntiProton])
{
- fhPtAntiProton ->Fill(ptcluster);
- fhPhiAntiProton ->Fill(ptcluster,phicluster);
- fhEtaAntiProton ->Fill(ptcluster,etacluster);
+ 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{
- fhPtUnknown ->Fill(ptcluster);
- fhPhiUnknown ->Fill(ptcluster,phicluster);
- fhEtaUnknown ->Fill(ptcluster,etacluster);
+ }
+ 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++) {
// }
// printf("\n");
- }
-
- //....................................................................
- // Access MC information in stack if requested, check that it exists.
- Int_t label =ph->GetLabel();
- if(label < 0) {
- 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()) {
+ }
+
+ //....................................................................
+ // 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){
+ }
+
+ 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(!mcparticles0) continue;
- if(label >= mcparticles0->GetEntriesFast()) {
+ }
+ 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, mcparticles0->GetEntriesFast());
+ label, mcparticles->GetEntriesFast());
continue ;
}
//Get the particle
- aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
+ aodprimary = (AliAODMCParticle*) mcparticles->At(label);
- }
-// else {//Second input
-// if(!mcparticles1) continue;
-// if(label >= mcparticles1->GetEntriesFast()) {
-// if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
-// label, mcparticles1->GetEntriesFast());
-// continue ;
-// }
-// //Get the particle
-// aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
-//
-// }//second input
-
- if(!aodprimary){
+ }
+
+ 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
-
+ }
+
+ 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
+
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff