bump_analys.C

Go to the documentation of this file.

{
        // Plot opening angle between 2 bumps and invariant mass of 2 bumps
        // It is assumed that high enrgy pi0 was generated which lead to 1 cluster with 2 bumps

        Double_t threshold=0.02;

  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  rootlogon();
  basiclibs();

        TFile* f = new TFile("cluster_emc.root"); //file you want to analyse
        TTree *t=(TTree *) f->Get("pndsim") ;
        TClonesArray* bump_array=new TClonesArray("PndEmcBump");
        t->SetBranchAddress("EmcBump",&bump_array);

        TFile* fsim = new TFile("sim_emc.root"); //file you want to analyse
        TTree *tsim=(TTree *) fsim->Get("pndsim") ;
        TClonesArray* mctrack_array=new TClonesArray("PndMCTrack");
        tsim->SetBranchAddress("MCTrack",&mctrack_array);

        TFile* fpar = new TFile("simparams.root");
        fpar->Get("FairBaseParSet");

        PndEmcMapper *emcMap=PndEmcMapper::Instance(2);

        TVector3 mc_momentum1, mc_momentum2;
        TVector3 mc_vertex1, mc_vertex2;

        TH1F *h_MC_angle= new TH1F("h_MC_angle","MC truth: opening angle of #pi^{0}",200,0.,15.);
        TH1F *h_angle= new TH1F("h_angle","reconstructed opening angle of #pi^{0}",200,0.,15.);
        TH1F *h_egamma= new TH1F("h_egamma","Energy of #gamma 's (MC truth)",100,0.0,10.0);
        TH1F *h_ebump= new TH1F("h_ebump","Bump energy",100,0.0,10.0);
        TH1F *h_MC= new TH1F("h_MC","Invariant mass",100,0.0,0.3);
        TH1F *h_MC_mass= new TH1F("h_MC_mass","Invariant mass MC",100,0.0,0.3);

        TH2F *h_energy_vs_angle_MC=new TH2F("h_energy_vs_angle_MC", "Energy vs. opening angle", 100,0.0,10.0,100,0.,15.0);


        TVector3 v1,v2; // Bump position
        Double_t bump_energy1=0;
        Double_t bump_energy2=0;

        // Bump energy
        cout<<"Number of events : "<<t->GetEntriesFast()<<endl;
        for (Int_t j=0; j< t->GetEntriesFast(); j++)
        {
                t->GetEntry(j);
                for (Int_t i=0; i<bump_array->GetEntriesFast(); i++)
                {
                        PndEmcBump *bump1=(PndEmcBump*)bump_array->At(i);
                        bump_energy1=bump1->energy();
                        if (bump_energy1<threshold) continue;
                        h_ebump->Fill(bump_energy1);
                        v1=bump1->where();

                        for (Int_t k=i+1; k<bump_array->GetEntriesFast(); k++)
                        {
                                PndEmcBump *bump2=(PndEmcBump*)bump_array->At(k);
                                bump_energy2=bump2->energy();
                                if (bump_energy2<threshold) continue;
                                v2=bump2->where();
                                double alpha=v1.Angle(v2);
                                double invMass=sqrt(2*bump_energy1*bump_energy2*(1-cos(alpha)));
                                h_angle->Fill(alpha*TMath::RadToDeg());
                                h_MC->Fill(invMass);
                        }

                }
        }

        TVector3 pi0_momentum;
        for (Int_t j=0; j< t->GetEntriesFast(); j++)
        {
                tsim->GetEntry(j);

                for (Int_t i=0; i<mctrack_array->GetEntriesFast(); i++)
                {
                        PndMCTrack *mctrack1=(PndMCTrack *) mctrack_array->At(i);
                        if (mctrack1->GetPdgCode()==22) // photons
                        {
                                mc_momentum1=mctrack1->GetMomentum();
                                mc_vertex1=mctrack1->GetStartVertex();
                                if (mc_vertex1.Mag()>2.) continue;
                                h_egamma->Fill(mc_momentum1.Mag());

                                for (Int_t k=i+1; k<mctrack_array->GetEntriesFast(); k++)
                                {
                                        PndMCTrack *mctrack2=(PndMCTrack *) mctrack_array->At(k);
                                        if (mctrack2->GetPdgCode()==22) // photons
                                        {
                                                mc_momentum2=mctrack2->GetMomentum();
                                                mc_vertex2=mctrack2->GetStartVertex();
                                                if (mc_vertex2.Mag()>2.) continue;
                                                double alphaMC=mc_momentum1.Angle(mc_momentum2);
                                                double invMassMC=sqrt(2*mc_momentum1.Mag()*mc_momentum2.Mag()*(1-cos(alphaMC)));
                                                h_MC_angle->Fill(alphaMC*TMath::RadToDeg());
                                                h_MC_mass->Fill(invMassMC);
                                                pi0_momentum=mc_momentum1+mc_momentum2;
                                                h_energy_vs_angle_MC->Fill(pi0_momentum.Mag(),alphaMC*TMath::RadToDeg());

                                        }
                                }
                        }
                }
        }

        TCanvas* c1 = new TCanvas("c1", "c1", 100, 100, 800, 800);
        c1->Divide(1,2);
        c1->cd(1);
        h_MC_angle->SetTitle("MC truth: opening angle of #pi^{0}");
        h_MC_angle->GetXaxis()->SetTitle("Opening angle, degree");
        h_MC_angle->Draw();

        c1->cd(2);
        h_angle->SetTitle("reconstructed opening angle of #pi^{0}");
        h_angle->GetXaxis()->SetTitle("Opening angle, degree");
        h_angle->Draw();

        TCanvas* c2 = new TCanvas("c2", "c2", 100, 100, 800, 800);
        c2->Divide(1,2);
        c2->cd(1);
        h_egamma->SetTitle("Energy of #gamma 's (MC truth)");
        h_egamma->GetXaxis()->SetTitle("Energy, GeV");
        h_egamma->Draw();

        c2->cd(2);
        h_ebump->SetTitle("Energy of bumps");
        h_ebump->GetXaxis()->SetTitle("Energy, GeV");
        h_ebump->Draw();

        TCanvas* c3 = new TCanvas("c3", "c3", 100, 100, 800, 800);
        c3->Divide(1,2);
        c3->cd(1);
        h_MC->SetTitle("Invariant mass");
        h_MC->GetXaxis()->SetTitle("M, GeV");
        h_MC->Draw();

        c3->cd(2);
        h_MC_mass->SetTitle("Invariant mass MC");
        h_MC_mass->GetXaxis()->SetTitle("M, GeV");
        h_MC_mass->Draw();

        TCanvas* c4 = new TCanvas("c4", "c4", 100, 100, 800, 800);
        h_energy_vs_angle_MC->SetTitle("Energy vs. opening angle (MC Truth)");
        h_energy_vs_angle_MC->GetXaxis()->SetTitle("Energy, GeV");
        h_energy_vs_angle_MC->GetYaxis()->SetTitle("Angle, degree");
        h_energy_vs_angle_MC->Draw();
}