add plot for origin and distribution of MC particles, other fixes
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 27 Nov 2013 14:51:28 +0000 (14:51 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 27 Nov 2013 14:51:28 +0000 (14:51 +0000)
PWGGA/CaloTrackCorrelations/macros/QA/DrawAnaCaloTrackQA.C

index a630777..15ab424 100644 (file)
@@ -53,6 +53,8 @@ void DrawAnaCaloTrackQA(TString listName = "Pi0IM_GammaTrackCorr_EMCAL_default",
   //Plot basic correlation QA
   CorrelationQA();
   
+  // MC basic QA plots, cluster origins (only if it run on MC)
+  MCQA();
 }
 
 //___________
@@ -126,7 +128,7 @@ void CaloQA()
   
   TLegend l2(0.3,0.7,0.83,0.85);
   l2.SetTextSize(0.04);
-  l2.AddEntry(hCorr,"No Exotics + non linearity","P");
+  l2.AddEntry(hCorr,"No Exotics + non lin. ++","P");
   l2.AddEntry(hTM,  "+ Track matching","P");
   l2.AddEntry(hShSh,"+ #lambda^{2}_{0} < 0.4","P");
   l2.SetBorderSize(0);
@@ -134,7 +136,88 @@ void CaloQA()
   l2.Draw();
   
   
-  ccalo->cd(3);
+  // Plot track-matching residuals
+  // first test did not have this histogram, add protection
+  TH2F* hTrackMatchResEtaPhi = (TH2F*) GetHisto("QA_hTrackMatchedDEtaDPhi");
+  if(hTrackMatchResEtaPhi)
+  {
+    ccalo->cd(3);
+    gPad->SetLogz();
+    
+    hTrackMatchResEtaPhi->SetAxisRange(-0.025,0.025,"X");
+    hTrackMatchResEtaPhi->SetAxisRange(-0.025,0.025,"Y");
+    hTrackMatchResEtaPhi->SetTitleOffset(1.5,"Y");
+    hTrackMatchResEtaPhi->SetTitle("Track-cluster residual #Delta #phi vs #Delta #eta, E > 0.5 GeV");
+    hTrackMatchResEtaPhi->Draw("colz");
+    
+    ccalo->cd(4);
+    gPad->SetLogy();
+    
+    TH2F* h2TrackMatchResEtaNeg = (TH2F*) GetHisto("QA_hTrackMatchedDEta");
+    TH2F* h2TrackMatchResEtaPos = (TH2F*) GetHisto("QA_hTrackMatchedDEtaPos");
+    TH2F* h2TrackMatchResPhiNeg = (TH2F*) GetHisto("QA_hTrackMatchedDPhi");
+    TH2F* h2TrackMatchResPhiPos = (TH2F*) GetHisto("QA_hTrackMatchedDPhiPos");
+    
+    h2TrackMatchResEtaNeg->Add(h2TrackMatchResEtaPos,-1);
+    h2TrackMatchResPhiNeg->Add(h2TrackMatchResPhiPos,-1);
+    
+    Float_t binMin = hCorr->FindBin(0.5);
+    TH1F* hTrackMatchResEtaNeg = (TH1F*) h2TrackMatchResEtaNeg->ProjectionY("TMProjEtaNeg",binMin, 1000);
+    TH1F* hTrackMatchResEtaPos = (TH1F*) h2TrackMatchResEtaPos->ProjectionY("TMProjEtaPos",binMin, 1000);
+    TH1F* hTrackMatchResPhiNeg = (TH1F*) h2TrackMatchResPhiNeg->ProjectionY("TMProjPhiNeg",binMin, 1000);
+    TH1F* hTrackMatchResPhiPos = (TH1F*) h2TrackMatchResPhiPos->ProjectionY("TMProjPhiPos",binMin, 1000);
+    
+    hTrackMatchResEtaNeg->SetTitle("Track-cluster residuals, E > 1 GeV");
+    hTrackMatchResEtaNeg->SetAxisRange(-0.05,0.05,"X");
+    hTrackMatchResEtaNeg->Sumw2();
+    hTrackMatchResEtaNeg->SetMarkerStyle(25);
+    hTrackMatchResEtaNeg->SetMarkerColor(2);
+    hTrackMatchResEtaNeg->Draw("");
+    
+    hTrackMatchResEtaPos->Sumw2();
+    hTrackMatchResEtaPos->SetMarkerStyle(25);
+    hTrackMatchResEtaPos->SetMarkerColor(4);
+    hTrackMatchResEtaPos->Draw("same");
+    
+    hTrackMatchResPhiNeg->Sumw2();
+    hTrackMatchResPhiNeg->SetMarkerStyle(24);
+    hTrackMatchResPhiNeg->SetMarkerColor(2);
+    hTrackMatchResPhiNeg->Draw("same");
+    
+    hTrackMatchResPhiPos->Sumw2();
+    hTrackMatchResPhiPos->SetMarkerStyle(24);
+    hTrackMatchResPhiPos->SetMarkerColor(4);
+    hTrackMatchResPhiPos->Draw("same");
+    
+    TLine l0(0,hTrackMatchResEtaNeg->GetMinimum(),0,hTrackMatchResEtaNeg->GetMaximum()*1.2);
+    l0.Draw("same");
+    
+    TLegend l3(0.55,0.7,0.83,0.85);
+    l3.SetTextSize(0.04);
+    l3.AddEntry(hTrackMatchResEtaNeg,"#Delta #eta, Negative","P");
+    l3.AddEntry(hTrackMatchResEtaPos,"#Delta #eta, Positive","P");
+    l3.AddEntry(hTrackMatchResPhiNeg,"#Delta #phi, Negative","P");
+    l3.AddEntry(hTrackMatchResPhiPos,"#Delta #phi, Positive","P");
+    l3.SetBorderSize(0);
+    l3.SetFillColor(0);
+    l3.Draw();
+  }
+  ccalo->Print(Form("%s_CaloHisto.eps",histoTag.Data()));
+  
+  
+  TCanvas * ccalo2 = new TCanvas(Form("CaloHisto2_%s",histoTag.Data()),"",500,500);
+  ccalo2->Divide(2,2);
+  
+  ccalo2->cd(3);
+//  gPad->SetLogz();
+//  TH2F* hCellAmpId   = (TH2F*) GetHisto("QA_hAmpId");
+//  hCellAmpId->SetTitle("Cell Id vs energy");
+//  hCellAmpId->SetYTitle("Cell Id");
+//  //hCellAmpId->SetAxisRange(300.,900.,"Y");
+//  hCellAmpId->SetAxisRange(0.,30.,"X");
+//  hCellAmpId->SetTitleOffset(1.5,"Y");
+//  hCellAmpId->Draw("colz");
+
   gPad->SetLogz();
   
   TH2F* hClusterTime   = (TH2F*) GetHisto("QA_hClusterTimeEnergy");
@@ -143,43 +226,16 @@ void CaloQA()
   hClusterTime->SetAxisRange(300.,900.,"Y");
   hClusterTime->SetAxisRange(0.,30.,"X");
   hClusterTime->SetTitleOffset(1.5,"Y");
-  
   hClusterTime->Draw("colz");
   
-  ccalo->cd(4);
-  gPad->SetLogz();
-  
-  TH2F* hClusterM02   = (TH2F*) GetHisto("AnaPhoton_hLam0E");
-  hClusterM02->SetTitle("Cluster energy vs #lambda^{2}_{0}");
-  //hClusterM02->SetAxisRange(300.,900.,"Y");
-  hClusterM02->SetAxisRange(0.,30.,"X");
-  hClusterM02->SetTitleOffset(1.5,"Y");
-  hClusterM02->Draw("colz");
-  
-  ccalo->Print(Form("%s_CaloHisto.eps",histoTag.Data()));
-  
-  
-  TCanvas * ccalo2 = new TCanvas(Form("CaloHisto2_%s",histoTag.Data()),"",500,500);
-  ccalo2->Divide(2,2);
-  
   ccalo2->cd(1);
-  gPad->SetLogz();
-  TH2F* hCellAmpId   = (TH2F*) GetHisto("QA_hAmpId");
-  hCellAmpId->SetTitle("Cell Id vs energy");
-  hCellAmpId->SetYTitle("Cell Id");
-  //hCellAmpId->SetAxisRange(300.,900.,"Y");
-  hCellAmpId->SetAxisRange(0.,30.,"X");
-  hCellAmpId->SetTitleOffset(1.5,"Y");
-  hCellAmpId->Draw("colz");
-  
-  ccalo2->cd(2);
   
   TH2F* hCellActivity  = (TH2F*) GetHisto("QA_hGridCells");
   hCellActivity->SetTitle("Hits per cell (E > 0.2 GeV)");
   hCellActivity->SetTitleOffset(1.5,"Y");
   hCellActivity->Draw("colz");
   
-  ccalo2->cd(3);
+  ccalo2->cd(2);
   
   TH2F* hCellActivity  = (TH2F*) GetHisto("QA_hGridCells");
   TH2F* hCellActivityE = (TH2F*) GetHisto("QA_hGridCellsE");
@@ -232,9 +288,9 @@ void TrackQA()
   
   hPhi     ->SetMinimum(1);
   
-  hPhi     ->Draw();
-  hPhiSPD  ->Draw("same");
-  hPhiNoSPD->Draw("same");
+  hPhi     ->Draw("H");
+  hPhiSPD  ->Draw("Hsame");
+  hPhiNoSPD->Draw("Hsame");
   
   TLegend l(0.12,0.8,0.4,0.9);
   l.SetTextSize(0.04);
@@ -293,8 +349,8 @@ void Pi0QA()
   }
   else
   {
+    h2DMass = (TH2F*) hMassE[0]->Clone("hMassProj");
     h2DMass->SetTitle("Invariant mass vs cluster pair E");
-    h2DMass = hMassE[0];
   }
   
   h2DMass->SetTitleOffset(1.6,"Y");
@@ -304,18 +360,29 @@ void Pi0QA()
   
   // Pi0 Invariant mass projection, in PbPb 6 centrality bins from 0 to 50%, all in pp
   cpi0->cd(2);
-  TH1F* hMass[10];
-  TH1F* hMix [10];
+  TH1F* hMass   [10];
+  TH1F* hMix    [10];
   TH1F* hMassEta[10];
   TH1F* hMassPi0[10];
   
+  //Init to 0
+  for(Int_t icen=0; icen<10; icen++ )
+  {
+    hMass   [icen] = 0;
+    hMix    [icen] = 0;
+    hMassEta[icen] = 0;
+    hMassPi0[icen] = 0;
+  }
+  
   TH1F * hX = (TH1F*) hMassE[0]->ProjectionX("hEPairCen0",0,10000);
   Int_t binmin = hX->FindBin(2);  // Project histo from 2 GeV pairs
   Int_t binmax = hX->FindBin(10); // Project histo up to 10 GeV pairs
+  Float_t maxPi0 = 0;
+  Float_t maxEta = 0;
   for(Int_t icen = 0; icen < 6; icen++)
   {
     if(!hMassE[icen]) continue;
-    
+
     hMass[icen] = (TH1F*) hMassE   [icen]->ProjectionY(Form("hMassCen%d",icen),binmin,binmax);
     hMix [icen] = (TH1F*) hMixMassE[icen]->ProjectionY(Form("hMixCen%d" ,icen),binmin,binmax);
     hMass[icen]->Sumw2();
@@ -345,11 +412,14 @@ void Pi0QA()
     hMassEta[icen]->SetAxisRange(0.4,0.9);
     hMassEta[icen]->SetMarkerSize(0.5);
     hMassEta[icen]->Scale(1./scale);
+    
+    if(maxEta < hMassEta[icen]->GetMaximum()) maxEta = hMassEta[icen]->GetMaximum();
+    if(maxPi0 < hMassPi0[icen]->GetMaximum()) maxPi0 = hMassPi0[icen]->GetMaximum();
   }
-  
+
   //gPad->SetLogy();
   //gPad->SetGridy();
-  hMassPi0[0]->SetMinimum(0.9);
+  hMassPi0[0]->SetMinimum(0.8);
   hMassPi0[0]->SetTitleOffset(1.6,"Y");
   hMassPi0[0]->SetYTitle("Real / Mixed");
   hMassPi0[0]->SetTitle("#pi^{0} peak, 2 < E_{pair}< 10 GeV");
@@ -357,7 +427,7 @@ void Pi0QA()
   
   if(hMass[1]) // PbPb
   {
-    hMassPi0[0]->SetMaximum(hMassPi0[5]->GetMaximum()*1.2);
+    hMassPi0[0]->SetMaximum(maxPi0*1.2);
     hMassPi0[5]->Draw("Hsame");
     hMassPi0[4]->Draw("Hsame");
     hMassPi0[3]->Draw("Hsame");
@@ -381,19 +451,24 @@ void Pi0QA()
     l.SetFillColor(0);
     l.Draw();
   }
-  
+
   TLine l1(0.04,1,0.24,1);
   l1.Draw("same");
   
   // Pi0 invariant mass per EMCal super module
   cpi0->cd(3);
   
-  TH1F* hSM[10];
+  TH1F* hSM   [10];
+  TH1F* hMixSM[10];
   binmin = hX->FindBin(4);  // Project histo from 3 GeV pairs
   binmax = hX->FindBin(20); // Project histo up to 20 GeV pairs
+  Float_t maxSM = 0;
+
   for(Int_t ism = 0; ism < 10; ism++)
   {
-    TH2F* hTmpSM = (TH2F*) GetHisto(Form("QA_hIM_Mod%d",ism));
+    TH2F* hTmpSM = (TH2F*) GetHisto(Form("AnaPi0_hReMod_%d",ism));
+    if(!hTmpSM) hTmpSM = (TH2F*) GetHisto(Form("QA_hIM_Mod%d",ism));
+    
     hSM[ism] = (TH1F*) hTmpSM->ProjectionY(Form("hMassSM%d",ism),binmin,binmax);
     hSM[ism]->Sumw2();
     hSM[ism]->SetMarkerStyle(26);
@@ -401,14 +476,35 @@ void Pi0QA()
     //hSM[ism]->Scale(1./hSM[ism]->Integral(0,10000));
     hSM[ism]->SetMarkerColor(color[ism]);
     hSM[ism]->SetLineColor(color[ism]);
-    hSM[ism]->SetAxisRange(0.04,0.24);
     hSM[ism]->SetMarkerSize(0.5);
+
+    TH2F* hTmpMixSM = (TH2F*) GetHisto(Form("AnaPi0_hMiMod_%d",ism));
+    if(hTmpMixSM)
+    {
+      hMixSM[ism] = (TH1F*) hTmpMixSM->ProjectionY(Form("hMassMixSM%d",ism),binmin,binmax);
+      hMixSM[ism]->Sumw2();
+      hMixSM[ism]->Rebin(2);
+      hSM[ism]->Divide(hMixSM[ism]);
+      hSM[ism]->Fit("pol0","Q","",0.25,0.35);
+      Float_t scale = 1;
+      if(hSM[ism]->GetFunction("pol0")) scale = hSM[ism]->GetFunction("pol0")->GetParameter(0);
+      //printf("Scale factor %f for cen %d\n",scale,icen);
+      hSM[ism]->Scale(1./scale);
+      hSM[ism]->SetYTitle("Real / Mixed");
+
+    }
+    
+    if(maxSM < hSM[ism]->GetMaximum()) maxSM = hSM[ism]->GetMaximum();
+
   }
   
   hSM[0]->SetTitle("#pi^{0} peak in Modules, 4 < E_{pair}< 10 GeV");
   hSM[0]->SetTitleOffset(1.6,"Y");
+  hSM[0]->SetAxisRange(0.04,0.24);
+  hSM[0]->SetMaximum(maxSM*1.2);
+  hSM[0]->SetMinimum(0.8);
+
   hSM[0]->Draw("H");
-  
   TLegend lsm(0.12,0.5,0.35,0.85);
   lsm.SetTextSize(0.04);
   lsm.AddEntry(hSM[0],Form("Mod %d",0),"P");
@@ -423,12 +519,14 @@ void Pi0QA()
   lsm.SetFillColor(0);
   lsm.Draw();
   
+  l1.Draw("same");
+  
   // Pi0 Invariant mass projection, in PbPb 6 centrality bins from 0 to 50%, all in pp
   cpi0->cd(4);
   
   //gPad->SetLogy();
   //gPad->SetGridy();
-  hMassEta[0]->SetMinimum(0.9);
+  hMassEta[0]->SetMinimum(0.8);
   hMassEta[0]->SetTitleOffset(1.6,"Y");
   hMassEta[0]->SetYTitle("Real / Mixed");
   hMassEta[0]->SetTitle("#eta peak, 2 < E_{pair}< 10 GeV");
@@ -436,7 +534,7 @@ void Pi0QA()
   
   if(hMass[1]) // PbPb
   {
-    hMassEta[0]->SetMaximum(hMassEta[5]->GetMaximum()*1.2);
+    hMassEta[0]->SetMaximum(maxEta*1.2);
     hMassEta[5]->Draw("Hsame");
     hMassEta[4]->Draw("Hsame");
     hMassEta[3]->Draw("Hsame");
@@ -457,7 +555,7 @@ void Pi0QA()
     l2.SetFillColor(0);
     l2.Draw();
   }
-  
+
   cpi0->Print(Form("%s_Pi0Histo.eps",histoTag.Data()));
   
 }
@@ -479,7 +577,7 @@ void CorrelationQA()
   
   //Azimuthal correlation
   cCorrelation->cd(1);
-  //gPad->SetLogy();
+  gPad->SetLogy();
   TH1F* hDeltaPhi[4];
   
   TLegend l(0.35,0.6,0.83,0.85);
@@ -506,7 +604,7 @@ void CorrelationQA()
     hDeltaPhi[ibin]->Scale(1./scale);
     //printf("ibin %d, scale %f\n",ibin,scale);
 
-    hDeltaPhi[ibin]->SetAxisRange(-1.8,4.1);
+    hDeltaPhi[ibin]->SetAxisRange(-1.6,4.7);
 
     hDeltaPhi[ibin]->SetMarkerStyle(24);
     hDeltaPhi[ibin]->SetMarkerColor(color[ibin]);
@@ -519,8 +617,8 @@ void CorrelationQA()
   }
   
   
-  //hDeltaPhi[2]->SetMaximum(hDeltaPhi[2]->GetMaximum()*1.1);
-  //hDeltaPhi[2]->SetMinimum(0.0001);
+  hDeltaPhi[2]->SetMaximum(hDeltaPhi[2]->GetMaximum()*10);
+  hDeltaPhi[2]->SetMinimum(0.8);
   
    hDeltaPhi[2]->Draw("H");
    hDeltaPhi[1]->Draw("Hsame");
@@ -576,6 +674,212 @@ void CorrelationQA()
 
 }
 
+//___________
+void MCQA()
+{
+  // Basic calorimeter QA histograms
+  
+  TH2F* h2ClusterPho = (TH2F*) GetHisto("QA_hRecoMCE_Photon_Match0");    // not track-matched
+  TH2F* h2ClusterPi0 = (TH2F*) GetHisto("QA_hRecoMCE_Pi0_Match0");       // not track-matched
+  TH2F* h2ClusterEta = (TH2F*) GetHisto("QA_hRecoMCE_Eta_Match0");       // not track-matched
+  TH2F* h2ClusterEle = (TH2F*) GetHisto("QA_hRecoMCE_Electron_Match1");  // Track-matched
+  
+  if(!h2ClusterPho) return;
+  
+  
+  TH1F* hPrimPho = (TH1F*) GetHisto("QA_hGenMCAccE_Photon");
+  TH1F* hPrimPi0 = (TH1F*) GetHisto("QA_hGenMCAccE_Pi0");
+  TH1F* hPrimEta = (TH1F*) GetHisto("QA_hGenMCAccE_Eta");
+  
+  TCanvas * cmc = new TCanvas(Form("MCHisto_%s",histoTag.Data()),"",1000,1000);
+  cmc->Divide(2,2);
+  
+  cmc->cd(1);
+  gPad->SetLogy();
+  
+  TH1F* hClusterPho = (TH1F*) h2ClusterPho->ProjectionX("ClusterPho",0,1000);
+  TH1F* hClusterPi0 = (TH1F*) h2ClusterPi0->ProjectionX("ClusterPi0",0,1000);
+  TH1F* hClusterEta = (TH1F*) h2ClusterEta->ProjectionX("ClusterEta",0,1000);
+  
+  hClusterPho->SetTitle("Cluster origin spectra, primary spectra in Calo acceptance");
+  hClusterPho->Sumw2();
+  hClusterPho->SetMarkerColor(1);
+  hClusterPho->SetMarkerStyle(20);
+  hClusterPho->SetAxisRange(0.,50.,"X");
+  hClusterPho->SetXTitle("E_{rec,gen} (GeV)");
+  hClusterPho->Draw("");
+
+  hClusterPi0->Sumw2();
+  hClusterPi0->SetMarkerColor(4);
+  hClusterPi0->SetMarkerStyle(21);
+  hClusterPi0->Draw("same");
+
+  hClusterEta->Sumw2();
+  hClusterEta->SetMarkerColor(2);
+  hClusterEta->SetMarkerStyle(22);
+  hClusterEta->Draw("same");
+
+  hPrimPho->Sumw2();
+  hPrimPho->SetMarkerColor(1);
+  hPrimPho->SetMarkerStyle(24);
+  hPrimPho->Draw("same");
+  
+  hPrimPi0->Sumw2();
+  hPrimPi0->SetMarkerColor(4);
+  hPrimPi0->SetMarkerStyle(25);
+  hPrimPi0->Draw("same");
+  
+  hPrimEta->Sumw2();
+  hPrimEta->SetMarkerColor(2);
+  hPrimEta->SetMarkerStyle(26);
+  hPrimEta->Draw("same");
+  
+  TLegend l(0.45,0.6,0.83,0.89);
+  l.SetTextSize(0.04);
+  l.AddEntry(hClusterPho,"#gamma cluster","P");
+  l.AddEntry(hClusterPi0,"#pi^{0} (merged) cluster","P");
+  l.AddEntry(hClusterEta,"#eta (merged) cluster","P");
+  l.AddEntry(hPrimPho,"#gamma generated","P");
+  l.AddEntry(hPrimPi0,"#pi^{0} generated","P");
+  l.AddEntry(hPrimEta,"#eta generated","P");
+  l.SetBorderSize(0);
+  l.SetFillColor(0);
+  l.Draw();
+  
+  
+  cmc->cd(2);
+  gPad->SetLogy();
+  TH1F* hRatPho = (TH1F*) hClusterPho->Clone("hGenRecoPho");
+  TH1F* hRatPi0 = (TH1F*) hClusterPi0->Clone("hGenRecoPi0");
+  TH1F* hRatEta = (TH1F*) hClusterEta->Clone("hGenRecoEta");
+  
+  hRatPho->Divide(hPrimPho);
+  hRatPi0->Divide(hPrimPi0);
+  hRatEta->Divide(hPrimEta);
+  
+  hRatPho->SetTitle("Reconstructed cluster / Generated particle in Calo acc.");
+  hRatPho->SetYTitle("Ratio");
+  hRatPho->SetXTitle("E(GeV)");
+  hRatPho->SetMinimum(1e-3);
+  hRatPho->SetMaximum(10);
+  hRatPho->Draw("");
+  hRatPi0->Draw("same");
+  hRatEta->Draw("same");
+
+  TLegend l2(0.15,0.7,0.3,0.85);
+  l2.SetTextSize(0.04);
+  l2.AddEntry(hRatPho,"#gamma","P");
+  l2.AddEntry(hRatPi0,"#pi^{0}","P");
+  l2.AddEntry(hRatEta,"#eta","P");
+  l2.SetBorderSize(0);
+  l2.SetFillColor(0);
+  l2.Draw();
+
+  cmc->cd(3);
+  //gPad->SetLogy();
+
+  TH2F* h2PrimPhoPhi = (TH2F*) GetHisto("AnaPhoton_hPhiPrim_MCPhoton");
+  TH2F* h2PrimPi0Phi = (TH2F*) GetHisto("AnaPi0_hPrimPi0Phi");
+  TH2F* h2PrimEtaPhi = (TH2F*) GetHisto("AnaPi0_hPrimEtaPhi");
+  
+  Int_t binMin = hPrimPho->FindBin(3);
+  
+  TH1F* hPrimPhoPhi = (TH1F*) h2PrimPhoPhi->ProjectionY("PrimPhoPhi",binMin,1000);
+  TH1F* hPrimPi0Phi = (TH1F*) h2PrimPi0Phi->ProjectionY("PrimPi0Phi",binMin,1000);
+  TH1F* hPrimEtaPhi = (TH1F*) h2PrimEtaPhi->ProjectionY("PrimEtaPhi",binMin,1000);
+
+  hPrimPhoPhi->Scale(1./hPrimPhoPhi->Integral(0,1000));
+  hPrimPi0Phi->Scale(1./hPrimPi0Phi->Integral(0,1000));
+  hPrimEtaPhi->Scale(1./hPrimEtaPhi->Integral(0,1000));
+
+  Float_t maxPhi = hPrimPhoPhi->GetMaximum();
+  if(maxPhi < hPrimPi0Phi->GetMaximum()) maxPhi =  hPrimPi0Phi->GetMaximum();
+  if(maxPhi < hPrimEtaPhi->GetMaximum()) maxPhi =  hPrimEtaPhi->GetMaximum();
+
+  Float_t minPhi = hPrimPhoPhi->GetMinimum();
+  if(minPhi > hPrimPi0Phi->GetMinimum()) minPhi =  hPrimPi0Phi->GetMinimum();
+  if(minPhi > hPrimEtaPhi->GetMinimum()) minPhi =  hPrimEtaPhi->GetMinimum();
+
+  hPrimPi0Phi->SetMaximum(maxPhi*1.1);
+  hPrimPi0Phi->SetMinimum(minPhi);
+  TGaxis::SetMaxDigits(3);
+
+  hPrimPi0Phi->SetYTitle("1/total entries dN/d#phi");
+  hPrimPi0Phi->SetTitle("Generated particles #phi for E > 3 GeV");
+  hPrimPi0Phi->SetTitleOffset(1.6,"Y");
+  hPrimPi0Phi->Sumw2();
+  hPrimPi0Phi->SetMarkerColor(4);
+  hPrimPi0Phi->SetMarkerStyle(21);
+  hPrimPi0Phi->Draw("");
+  
+  hPrimPhoPhi->Sumw2();
+  hPrimPhoPhi->SetMarkerColor(1);
+  hPrimPhoPhi->SetMarkerStyle(20);
+  Float_t scale = TMath::RadToDeg();
+  ScaleXaxis(hPrimPhoPhi, TMath::RadToDeg());
+  hPrimPhoPhi->Draw("same");
+
+
+  hPrimEtaPhi->Sumw2();
+  hPrimEtaPhi->SetMarkerColor(2);
+  hPrimEtaPhi->SetMarkerStyle(22);
+  hPrimEtaPhi->Draw("same");
+
+  cmc->cd(4);
+  //gPad->SetLogy();
+  
+  TH2F* h2PrimPhoEta = (TH2F*) GetHisto("AnaPhoton_hEtaPrim_MCPhoton");
+  TH2F* h2PrimPi0Eta = (TH2F*) GetHisto("AnaPi0_hPrimPi0Rapidity");
+  TH2F* h2PrimEtaEta = (TH2F*) GetHisto("AnaPi0_hPrimEtaRapidity");
+  
+  Int_t binMin = hPrimPho->FindBin(3);
+  
+  TH1F* hPrimPhoEta = (TH1F*) h2PrimPhoEta->ProjectionY("PrimPhoEta",binMin,1000);
+  TH1F* hPrimPi0Eta = (TH1F*) h2PrimPi0Eta->ProjectionY("PrimPi0Eta",binMin,1000);
+  TH1F* hPrimEtaEta = (TH1F*) h2PrimEtaEta->ProjectionY("PrimEtaEta",binMin,1000);
+  
+  hPrimPhoEta->Scale(1./hPrimPhoEta->Integral(0,1000));
+  hPrimPi0Eta->Scale(1./hPrimPi0Eta->Integral(0,1000));
+  hPrimEtaEta->Scale(1./hPrimEtaEta->Integral(0,1000));
+  
+  Float_t maxEta = hPrimPhoEta->GetMaximum();
+  if(maxEta < hPrimPi0Eta->GetMaximum()) maxEta =  hPrimPi0Eta->GetMaximum();
+  if(maxEta < hPrimEtaEta->GetMaximum()) maxEta =  hPrimEtaEta->GetMaximum();
+  
+  Float_t minEta = hPrimPhoEta->GetMinimum();
+  if(minEta > hPrimPi0Eta->GetMinimum()) minEta =  hPrimPi0Eta->GetMinimum();
+  if(minEta > hPrimEtaEta->GetMinimum()) minEta =  hPrimEtaEta->GetMinimum();
+  
+  hPrimPi0Eta->SetMaximum(maxEta*1.1);
+  hPrimPi0Eta->SetMinimum(minEta);
+  TGaxis::SetMaxDigits(3);
+  
+  hPrimPi0Eta->SetYTitle("1/total entries dN/d#eta");
+  hPrimPi0Eta->SetTitle("Generated particles #eta for E > 3 GeV");
+  hPrimPi0Eta->SetTitleOffset(1.6,"Y");
+  hPrimPi0Eta->Sumw2();
+  hPrimPi0Eta->SetMarkerColor(4);
+  hPrimPi0Eta->SetMarkerStyle(21);
+  hPrimPi0Eta->Draw("");
+  
+  hPrimPhoEta->Sumw2();
+  hPrimPhoEta->SetMarkerColor(1);
+  hPrimPhoEta->SetMarkerStyle(20);
+  Float_t scale = TMath::RadToDeg();
+  hPrimPhoEta->Draw("same");
+  
+  hPrimEtaEta->Sumw2();
+  hPrimEtaEta->SetMarkerColor(2);
+  hPrimEtaEta->SetMarkerStyle(22);
+  hPrimEtaEta->Draw("same");
+  
+  cmc->Print(Form("%s_MCHisto.eps",histoTag.Data()));
+
+
+
+  
+ }
+
 
 //____________________________________________________________________
 void GetFileAndList(TString fileName, TString listName, Bool_t export)
@@ -604,3 +908,37 @@ TObject * GetHisto(TString histoName)
   else     return file->Get       (histoName);
 }
 
+//___________________________________________________
+void ScaleAxis(TAxis *a, Double_t scale)
+{
+  if (!a) return; // just a precaution
+  if (a->GetXbins()->GetSize())
+  {
+    // an axis with variable bins
+    // note: bins must remain in increasing order, hence the "Scale"
+    // function must be strictly (monotonically) increasing
+    TArrayD X(*(a->GetXbins()));
+    for(Int_t i = 0; i < X.GetSize(); i++) X[i] = scale*X[i];
+    a->Set((X.GetSize() - 1), X.GetArray()); // new Xbins
+  }
+  else
+  {
+    // an axis with fix bins
+    // note: we modify Xmin and Xmax only, hence the "Scale" function
+    // must be linear (and Xmax must remain greater than Xmin)
+    a->Set(a->GetNbins(),
+           scale*a->GetXmin(), // new Xmin
+           scale*a->GetXmax()); // new Xmax
+  }
+  return;
+}
+
+//___________________________________________________
+void ScaleXaxis(TH1 *h, Double_t scale)
+{
+  if (!h) return; // just a precaution
+  ScaleAxis(h->GetXaxis(), scale);
+  return;
+}
+
+