emc_correction_QA.C

Go to the documentation of this file.

// QA macro for EMC cluster energy correction
// Depending on "test" input parametr checks
// 1 - energy ratio of photon vs theta and E
// 2 - energy distribution for single 1 GeV gammas
// 3 - pi0 invariant mass
// Two methods can be tested depending on input parametr (method=1 or 2)
int emc_correction_QA(TString InputFile="emc_complete_QA.root", Int_t test=1, Int_t version=1)
{
        // test=1
        TH1F *h_energy= new TH1F("h_energy","Cluster energy",100,0.8,1.2);
        TH1F *h_energy_corr1= new TH1F("h_energy_corr1","Cluster energy (corrected), (hist)",100,0.8,1.2);
        TH1F *h_energy_corr2= new TH1F("h_energy_corr2","Cluster energy (corrected), (parametrization)",100,0.8,1.2);
        TH1F *h_energy_corr_old= new TH1F("h_energy_corr_old","Cluster energy (corrected), (old)",100,0.8,1.2);
        TH1F *h_energy_shashlyk= new TH1F("h_energy_shashlyk","Cluster energy (shashlyk)",100,0.8,1.2);
        TH1F *h_energy_shashlyk_corr1= new TH1F("h_energy_shashlyk_corr1","Cluster energy (shashlyk) (corrected), (hist)",100,0.8,1.2);
        TH1F *h_energy_shashlyk_corr2= new TH1F("h_energy_shashlyk_corr2","Cluster energy (shashlyk) (corrected), (parametrization)",100,0.8,1.2);
        TH1F *h_energy_shashlyk_corr_old= new TH1F("h_energy_shashlyk_corr_old","Cluster energy (shashlyk) (corrected)",100,0.8,1.2);

        // test=2
        TH2F *hEnergyRatioTheta=new TH2F("hEnergyRatioTheta","Energy_reco/Energy_MC",100, 0., 180., 100, 0.8, 1.2);
        hEnergyRatioTheta->GetXaxis()->SetTitle("Cluster Reconstructed #theta Angle (#circ)");
        TH2F *hEnergyRatioEnergy=new TH2F("hEnergyRatioEnergy","Energy_reco/Energy_MC",100, 0., 10., 100, 0.8, 1.2);
        hEnergyRatioEnergy->GetXaxis()->SetTitle("Cluster Reconstructed Energy (GeV)");
        TH2F *hEnergyRatioCorrEnergy1=new TH2F("hEnergyRatioCorrEnergy1","Energy_reco (corrected)/Energy_MC, (hist)",100, 0., 10., 100, 0.8, 1.2);
        hEnergyRatioCorrEnergy1->GetXaxis()->SetTitle("Cluster Reconstructed Energy (GeV)");
        TH2F *hEnergyRatioCorrTheta1=new TH2F("hEnergyRatioCorrTheta1","Energy_reco (corrected)/Energy_MC, (hist)",100, 0., 180., 100, 0.8, 1.2);
        hEnergyRatioCorrTheta1->GetXaxis()->SetTitle("Cluster Reconstructed #theta Angle (#circ)");
        TH2F *hEnergyRatioCorrEnergy2=new TH2F("hEnergyRatioCorrEnergy2","Energy_reco (corrected)/Energy_MC, (parametrization)",100, 0., 10., 100, 0.8, 1.2);
        hEnergyRatioCorrEnergy2->GetXaxis()->SetTitle("Cluster Reconstructed Energy (GeV)");
        TH2F *hEnergyRatioCorrTheta2=new TH2F("hEnergyRatioCorrTheta2","Energy_reco (corrected)/Energy_MC, (parametrization)",100, 0., 180., 100, 0.8, 1.2);
        hEnergyRatioCorrTheta2->GetXaxis()->SetTitle("Cluster Reconstructed #theta Angle (#circ)");
        TH2F *hEnergyRatioCorrEnergyOld=new TH2F("hEnergyRatioCorrEnergyOld","Energy_reco (corrected)/Energy_MC, (old)",100, 0., 10., 100, 0.8, 1.2);
        hEnergyRatioCorrEnergyOld->GetXaxis()->SetTitle("Cluster Reconstructed Energy (GeV)");
        TH2F *hEnergyRatioCorrThetaOld=new TH2F("hEnergyRatioCorrThetaOld","Energy_reco (corrected)/Energy_MC, (old)",100, 0., 180., 100, 0.8, 1.2);
        hEnergyRatioCorrThetaOld->GetXaxis()->SetTitle("Cluster Reconstructed #theta Angle (#circ)");

        // test=3
        TH1F *h_mpi0= new TH1F("h_mpi0","pi0 invariant mass",100,0.05,0.2);
        TH1F *h_mpi0_corr1= new TH1F("h_mpi0_corr1","pi0 invariant mass (corrected) (hist)",100,0.05,0.2);
        TH1F *h_mpi0_corr2= new TH1F("h_mpi0_corr2","pi0 invariant mass (corrected) (parametrization)",100,0.05,0.2);
        TH1F *h_mpi0_corr_old= new TH1F("h_mpi0_corr_old","pi0 invariant mass (corrected) (old)",100,0.05,0.2);

        TFile *f=new TFile(InputFile);
        TTree *t=(TTree *) f->Get("pndsim") ;

        TClonesArray* mctrack_array=new TClonesArray("PndMCTrack");
        t->SetBranchAddress("MCTrack",&mctrack_array);

        TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
        t->SetBranchAddress("EmcCluster",&cluster_array);

        Double_t cluster_energy, cluster_theta, cluster_phi; //position of the cluster
        TVector3 mc_momentum;
        Double_t thetaMC, phiMC, energyMC;
        Double_t max_energy;

        // Here the cluster calibrator is set
        // first parameter is method (1-histogram, 2-parametrization)
        // second parameter is version (defined by geometry)
        PndEmcAbsClusterCalibrator * calibrator1= PndEmcClusterCalibrator::MakeEmcClusterCalibrator(1, version);
        PndEmcAbsClusterCalibrator * calibrator2= PndEmcClusterCalibrator::MakeEmcClusterCalibrator(2, version);

        if ((test==1)||(test==2))
        {
                for (Int_t j=0; j<t->GetEntriesFast(); j++)
                {
                        t->GetEntry(j);
                        PndMCTrack *mctrack=(PndMCTrack *) mctrack_array->At(0);
                        mc_momentum=mctrack->GetMomentum();
                        energyMC=mc_momentum.Mag();
                        thetaMC=mc_momentum.Theta()*180./TMath::Pi();
                        max_energy=0;
                        for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
                        {
                                PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(i);
                                cluster_energy=cluster->energy();

                                if (cluster_energy>max_energy)
                                {
                                        max_energy=cluster_energy;
                                        TVector3 cluster_pos=cluster->where();
                                        cluster_theta=cluster_pos.Theta()*180./TMath::Pi();
                                        Double_t energy_corr1 = calibrator1->Energy(cluster);
                                        Double_t energy_corr2 = calibrator2->Energy(cluster);
                                        Double_t energy_corr_old = cluster->GetEnergyCorrected();
                                }

                                Double_t e_ratio=max_energy/energyMC;
                                Double_t e_ratio_corr1=energy_corr1/energyMC;
                                Double_t e_ratio_corr2=energy_corr2/energyMC;
                                Double_t e_ratio_corr_old=energy_corr_old/energyMC;

                                if (cluster_theta<5.)
                                {
                                        h_energy_shashlyk->Fill(max_energy);
                                        h_energy_shashlyk_corr1->Fill(energy_corr1);
                                        h_energy_shashlyk_corr2->Fill(energy_corr2);
                                        h_energy_shashlyk_corr_old->Fill(energy_corr_old);
                                } else {
                                        h_energy->Fill(max_energy);
                                        h_energy_corr1->Fill(energy_corr1);
                                        h_energy_corr2->Fill(energy_corr2);
                                        h_energy_corr_old->Fill(energy_corr_old);
                                }

                                hEnergyRatioEnergy->Fill( energyMC, e_ratio);
                                hEnergyRatioTheta->Fill( thetaMC, e_ratio);
                                hEnergyRatioCorrEnergy1->Fill(energyMC, e_ratio_corr1);
                                hEnergyRatioCorrTheta1->Fill(thetaMC, e_ratio_corr1);
                                hEnergyRatioCorrEnergy2->Fill(energyMC, e_ratio_corr2);
                                hEnergyRatioCorrTheta2->Fill(thetaMC, e_ratio_corr2);
                                hEnergyRatioCorrEnergyOld->Fill(energyMC, e_ratio_corr_old);
                                hEnergyRatioCorrThetaOld->Fill(thetaMC, e_ratio_corr_old);
                        }
                }
        }

        if (test==3)
        {
                double threshold=0.02;
                TVector3 v1,v2;
                TVector3 v1_corr1,v2_corr1;
                TVector3 v1_corr2,v2_corr2;
                Double_t energy1, energy2, energy1_corr1, energy2_corr1;
                Double_t energy1_corr2, energy2_corr2, energy1_corr_old, energy2_corr_old;
                // Cluster energy
                for (Int_t j=0; j< t->GetEntriesFast(); j++)
                {
                        t->GetEntry(j);
                        for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
                        {
                                PndEmcCluster *cluster1=(PndEmcCluster*)cluster_array->At(i);
                                energy1=cluster1->energy();
                                if (energy1<threshold) continue;
                                v1=cluster1->where();
                                energy1_corr1 = calibrator1->Energy(cluster1);
                                v1_corr1 = calibrator1->Where(cluster1);
                                energy1_corr2 = calibrator2->Energy(cluster1);
                                v1_corr2 = calibrator2->Where(cluster1);
                                energy1_corr_old = cluster1->GetEnergyCorrected();

                                for (Int_t k=i; k<cluster_array->GetEntriesFast(); k++)
                                {
                                        PndEmcCluster *cluster2=(PndEmcCluster*)cluster_array->At(k);
                                        if (cluster2->energy()<threshold) continue;
                                        energy2=cluster2->energy();
                                        energy2_corr1 = calibrator1->Energy(cluster2);
                                        energy2_corr2 = calibrator2->Energy(cluster2);
                                        energy2_corr_old = cluster2->GetEnergyCorrected();
                                        v2=cluster2->where();
                                        v2_corr1 = calibrator1->Where(cluster2);
                                        v2_corr2 = calibrator2->Where(cluster2);
                                        double alpha=v1.Angle(v2);
                                        double invMass=sqrt(2*energy1*energy2*(1-cos(alpha)));
                                        h_mpi0->Fill(invMass);
                                        double alpha_corr1=v1_corr1.Angle(v2_corr1);
                                        double invMass_corr1=sqrt(2*energy1_corr1*energy2_corr1*(1-cos(alpha_corr1)));
                                        h_mpi0_corr1->Fill(invMass_corr1);
                                        double alpha_corr2=v1_corr2.Angle(v2_corr2);
                                        double invMass_corr2=sqrt(2*energy1_corr2*energy2_corr2*(1-cos(alpha_corr2)));
                                        h_mpi0_corr2->Fill(invMass_corr2);
                                        double invMass_corr_old=sqrt(2*energy1_corr_old*energy2_corr_old*(1-cos(alpha)));
                                        h_mpi0_corr_old->Fill(invMass_corr_old);
                                }
                        }
                }
        }


if (test==1)
{
        TCanvas* c1 = new TCanvas("c1", "Energy reconstruction vs Energy", 100, 100, 800, 800);
        c1->Divide(2,2);
        c1->cd(1);
        hEnergyRatioEnergy->Draw("contz");
        TLine *l=new TLine(0.,1.,10.,1.0);
        l->SetLineColor(2);
        l->SetLineWidth(2);
        l->Draw();

        c1->cd(2);
        hEnergyRatioCorrEnergy1->Draw("contz");
        TLine *l=new TLine(0.,1.,10.,1.0);
        l->SetLineColor(2);
        l->SetLineWidth(2);
        l->Draw();

        c1->cd(3);
        hEnergyRatioCorrEnergy2->Draw("contz");
        TLine *l=new TLine(0.,1.,10.,1.0);
        l->SetLineColor(2);
        l->SetLineWidth(2);
        l->Draw();

//      c1->cd(4);
//      hEnergyRatioCorrEnergyOld->Draw("contz");
//      TLine *l=new TLine(0.,1.,10.,1.0);
//      l->SetLineColor(2);
//      l->SetLineWidth(2);
//      l->Draw();

        TCanvas* c2 = new TCanvas("c2", "Energy reconstruction vs Theta", 100, 100, 800, 800);
        c2->Divide(2,2);
        c2->cd(1);
        hEnergyRatioTheta->Draw("contz");
        TLine *l=new TLine(0.,1.,180.,1.0);
        l->SetLineColor(2);
        l->SetLineWidth(2);
        l->Draw();
        c2->cd(2);
        hEnergyRatioCorrTheta1->Draw("contz");
        TLine *l=new TLine(0.,1.,180.,1.0);
        l->SetLineColor(2);
        l->SetLineWidth(2);
        l->Draw();
        c2->cd(3);
        hEnergyRatioCorrTheta2->Draw("contz");
        TLine *l=new TLine(0.,1.,180.,1.0);
        l->SetLineColor(2);
        l->SetLineWidth(2);
        l->Draw();
//      c2->cd(4);
//      hEnergyRatioCorrThetaOld->Draw("contz");
//      TLine *l=new TLine(0.,1.,180.,1.0);
//      l->SetLineColor(2);
//      l->SetLineWidth(2);
//      l->Draw();
}

if (test==2)
{
        TCanvas* c3 = new TCanvas("c3", "Energy (TS)", 100, 100, 800, 800);
        c3->Divide(2,2);
        c3->cd(1);
        h_energy->Draw();
        c3->cd(2);
        h_energy_corr1->Draw();
        c3->cd(3);
        h_energy_corr2->Draw();
//      c3->cd(4);
//      h_energy_corr_old->Draw();

        TCanvas* c4 = new TCanvas("c4", "Energy (shashlyk)", 100, 100, 800, 800);
        c4->Divide(2,2);
        c4->cd(1);
        h_energy_shashlyk->Draw();
        c4->cd(2);
        h_energy_shashlyk_corr1->Draw();
        c4->cd(3);
        h_energy_shashlyk_corr2->Draw();
//      c4->cd(4);
//      h_energy_shashlyk_corr_old->Draw();
}
if (test==3)
{
        gStyle->SetOptFit(kTRUE);

        double pars[5]={800,0.13,0.01,100,-500};
        TF1 *f1 = new TF1("f1","gaus(0)+pol1(3)",0.08,0.18);
        TF1 *f2 = new TF1("f2","gaus(0)+pol1(3)",0.08,0.18);
        TF1 *f3 = new TF1("f3","gaus(0)+pol1(3)",0.08,0.18);
        TF1 *f4 = new TF1("f4","gaus(0)+pol1(3)",0.08,0.18);
        f1->SetParameters(pars);
        f2->SetParameters(pars);
        f3->SetParameters(pars);
        f4->SetParameters(pars);

        TCanvas* c5 = new TCanvas("c5", "pi0 mass", 100, 100, 800, 800);
        c5->Divide(2,2);
        c5->cd(1);
        h_mpi0->Draw();
        h_mpi0->Fit(f1,"R");
        c5->cd(2);
        h_mpi0_corr1->Draw();
        h_mpi0_corr1->Fit(f2,"R");
        c5->cd(3);
        h_mpi0_corr2->Draw();
        h_mpi0_corr2->Fit(f3,"R");
//      c5->cd(4);
//      h_mpi0_corr_old->Draw();
//      h_mpi0_corr_old->Fit(f4,"R");
}
  return 0;
}