macro/detectors/emc/dedicated/reco_analys.C

Go to the documentation of this file.

{
  // Macro loads a file after reconstruction and plots difference between initial direction of particle and angular position of cluster
        TFile* f = new TFile("cluster_emc.root"); //file you want to analyse
        TTree *t=(TTree *) f->Get("pndsim") ;
        TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
        t->SetBranchAddress("EmcCluster",&cluster_array);

        t->AddFriend("pndsim", "digi_emc.root");
        TClonesArray* digi_array=new TClonesArray("PndEmcDigi");
        t->SetBranchAddress("EmcDigi",&digi_array);

        TFile* fsim = new TFile("sim_emc.root"); //file you want to analyse
        TTree *tsim=(TTree *) fsim->Get("pndsim") ;

        PndEmcMapper::Init(1);

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

        TVector3 photon_momentum;

        double cluster_energy;
        double cluster_theta, cluster_phi; //position of the cluster
        double theta, phi; // angular position of the initial particle
        double theta_diff, phi_diff;
        int ndigi, npoint;
        double max_energy=0;

        TH1F *h1= new TH1F("h1","Theta difference",200,-5.,5.);
        TH1F *h2= new TH1F("h2","Phi difference",200,-5.,5.);
        TH1F *h3= new TH1F("h3","Cluster energy",100,0.85,1.05);
        TH2F *h2theta= new TH2F("h2theta","Theta difference",200,0.,180.,200,-5.,5.);
        TH2F *h2phi= new TH2F("h2phi","Phi difference",200,0.,180.,200,-5.,5.);
        TH1F *hE1= new TH1F("hE1","E1",200,0.,1.05);
        TH1F *hE1E9= new TH1F("hE1E9","E1 / E9",200,0.,1.05);
        TH1F *hE9E25= new TH1F("hE9E25","E9 / E25",200,0.,1.05);

        // Cluster angular position
        // Entrance point is determined by minimal time

        // 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 *cluster=(PndEmcCluster*)cluster_array->At(i);
                        cluster_energy=cluster->energy();
                        if ((cluster->NumberOfDigis()>1)&&(cluster_energy>0.02))
                                h3->Fill(cluster_energy);
                        PndEmcClusterEnergySums esum(*cluster, digi_array);
                        hE1->Fill(esum.E1());
                        hE1E9->Fill(esum.E1E9());
                        hE9E25->Fill(esum.E9E25());

                }
        }

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

                PndMCTrack *mctrack=(PndMCTrack *) mctrack_array->At(0);
                photon_momentum=mctrack->GetMomentum();
                theta=photon_momentum.Theta();
                phi=photon_momentum.Phi();


                // Loop over clusters
                // If we have 1 initial particle and several cluster
                // we can separate cluster from the first interaction by maximum energy

                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();
                                cluster_phi=cluster_pos.Phi();
                        }

                }

                if (max_energy>0.6)
                {
                        theta_diff=(cluster_theta-theta)*180./TMath::Pi();
                        h1->Fill(theta_diff);
                        h2theta->Fill(theta*TMath::RadToDeg(),theta_diff);

                        phi_diff=(cluster_phi-phi)*180./TMath::Pi();
                        h2->Fill(phi_diff);
                        h2phi->Fill(phi*TMath::RadToDeg(),phi_diff);
                }

        }

        TCanvas* c1 = new TCanvas("c1", "Cluster Energy", 100, 100, 800, 800);
        h3->SetTitle("Cluster energy of 1 GeV photon");
        h3->GetXaxis()->SetTitle("Energy, GeV");
        h3->Draw();


        TCanvas* c2 = new TCanvas("c2", "#theta_{reco} - #theta_{truth}", 100, 100, 800, 800);
        h1->SetTitle("Difference between cluster and initial photon theta angle");
        h1->GetXaxis()->SetTitle("#theta_{reco} - #theta_{truth}, degree");
        h1->Draw();

//      TF1 *f1 = new TF1("f1","gaus(0)",-2.,2.);
//      Double_t par1[3]={60,0,0.5};
//      f1->SetParameters(par1);
//      f1->SetLineColor(2);
        //
//      h1->Fit("f1","RB");
//      double mu1=f1->GetParameter(1);
//      double sigma1=f1->GetParameter(2);

        TCanvas* c3 = new TCanvas("c3", "#phi_{reco} - #phi_{truth}", 100, 100, 800, 800);
        h2->SetTitle("Difference between cluster and initial photon phi angle");
        h2->GetXaxis()->SetTitle("#phi_{reco} - #phi_{truth}, degree");
        h2->Draw();

        TCanvas* c4 = new TCanvas("c4", "#theta_{reco} - #theta_{truth} vs #theta_{truth}", 100, 100, 800, 800);
        h2theta->SetTitle("Difference between cluster and initial photon theta angle");
        h2theta->GetXaxis()->SetTitle("#theta_{truth}, degree");
        h2theta->GetYaxis()->SetTitle("#theta_{reco} - #theta_{truth}, degree");
        h2theta->Draw();

        TCanvas* c5 = new TCanvas("c5", "#phi_{reco} - #phi_{truth} vs #phi_{truth}", 100, 100, 800, 800);
        h2phi->SetTitle("Difference between cluster and initial photon phi angle");
        h2phi->GetXaxis()->SetTitle("#phi_{truth}, degree");
        h2phi->GetYaxis()->SetTitle("#phi_{reco} - #phi_{truth}, degree");
        h2phi->Draw();

        TCanvas* c6 = new TCanvas("c6", "Cluster Properties", 100, 100, 800, 800);
        c6->Divide(2,2);
        c6->cd(1); hE1->Draw();
        c6->cd(2); hE1E9->Draw();
        c6->cd(3); hE9E25->Draw();


}