CheckRung/check/gainVsDoca.cxx

Go to the documentation of this file.

//This code can display gains from two files vs doca on one pad. One pad for one entrance angle
 
#include <vector>
#include <string>
 
using namespace std;
 
vector<string>file;
void pushNo()
{
    file.push_back("../../DocaEangle/docaeangle.root");
    file.push_back("docaeangle.root");
}
 
#include "TStyle.h"
#include "TString.h"
 
const int PadNo1 = 0;
const int PadNo2 = 40;
const int TotalNo = 1600;
const int NumSlices = 40;
const int gainmin = 0;
const int gainmax = 3;
TString LegendName = "gain vs doca";
const int Ipadx = 1;
const int Ipady = 1;
const int CanvasX = 800;  //canvas size for ps file
const int CanvasY = 600;
bool NormFlag = false;  //true: nomalize wiregain for each layer
const double histmin = 0.7; //Y axis range of histogram when NormFlag is true
const double histmax = 1.3;
const double delta = 0.15; //Y axis range of histogram when NormFlag is false
 
 
void gainVsDoca()
{
    pushNo();
    const int FileNo = file.size();
    TH2F*** h = new TH2F**[PadNo2];
    TCanvas*** c = new TCanvas**[PadNo2];
    TLegend*** lgd = new TLegend**[PadNo2];
    gStyle->SetOptStat(0);
    gStyle->SetCanvasColor(10);
    for(int i=0;i<PadNo2;i++)
    {   
        h[i] = new TH2F*[FileNo/2];
        c[i] = new TCanvas*[FileNo/2];
        lgd[i] = new TLegend*[FileNo/2];
    }   
 
 
    TFile* f1;
    TFile* f2;
    TTree* tree1;
    TTree* tree2;
    double gain1[TotalNo],gain2[TotalNo],Ip_eangle1[TotalNo],Ip_eangle2[TotalNo],Norm1,Norm2;
 
    stringstream histname;  //histogram name
    stringstream ss_file1;  //file name 
    stringstream ss_file2;
    stringstream ss1;  //picture's drawing content
    stringstream ss2;
    stringstream cut;
    stringstream legendheader;
 
    for(int k=0;k<FileNo;k+=2)
    {
        ss_file1.str("");
        ss_file1<<file[k];
        f1=new TFile(ss_file1.str().c_str());
        tree1=(TTree*)f1->Get("ddgcalib");
        tree1->SetBranchAddress("Out_gain", gain1);
        tree1->SetBranchAddress("Ip_eangle", Ip_eangle1);
        tree1->GetEntry(0);
        tree1->SetMarkerStyle(22);
        tree1->SetMarkerColor(4);
        tree1->SetLineColor(4);
 
        ss_file2.str("");
        ss_file2<<file[k+1];
        f2=new TFile(ss_file2.str().c_str());
        tree2=(TTree*)f2->Get("ddgcalib");
        tree2->SetBranchAddress("Out_gain", gain2);
        tree2->SetBranchAddress("Ip_eangle", Ip_eangle2);
        tree2->GetEntry(0);
        tree2->SetMarkerStyle(20);
        tree2->SetMarkerColor(2);
        tree2->SetLineColor(2);
 
        for(int i=PadNo1;i<PadNo2;i++)
        {
            Norm1=0;
            Norm2=0;
            int k1=0;
            int k2=0;
 
            for(int j=0;j<TotalNo;j++)
            {
                if(Ip_eangle1[j]!=i) continue;
                if(gain1[j]>gainmin && gain1[j]<gainmax) {Norm1 += gain1[j];k1++;}
                if(gain2[j]>gainmin && gain2[j]<gainmax) {Norm2 += gain2[j];k2++;}
            }
            if(k1==0 || k2==0) {Norm1=1; Norm2=1;}
            else {Norm1/= k1; Norm2/= k2;}
            cout<<"Norm1= "<<Norm1<<" Norm2= "<<Norm2<<endl;
 
            legendheader.str("");
            legendheader<<"eangle"<<i<<": "<<LegendName;
            //cout<<"legendheader: "<<legendheader.str().c_str()<<endl;
            lgd[i][k/2] = new TLegend(.3, .75, .89, .89);
            lgd[i][k/2] ->SetFillColor(10);
            lgd[i][k/2] ->SetBorderSize(1);
            lgd[i][k/2] ->SetHeader(legendheader.str().c_str());
            lgd[i][k/2] ->AddEntry(tree1, ss_file1.str().c_str() , "p");
            lgd[i][k/2] ->AddEntry(tree2, ss_file2.str().c_str() , "p");
 
            histname.str("");
            histname<<"eangle"<<i<<"_"<<file[k]<<"-"<<file[k+1];
            //cout<<"histname: "<<histname.str().c_str()<<endl;
            if(NormFlag) h[i][k/2] = new TH2F(histname.str().c_str(),"",100,0,NumSlices,100,histmin,histmax);
            else h[i][k/2] = new TH2F(histname.str().c_str(),"",100,0,NumSlices,100,(Norm1+Norm2)/2-delta,(Norm1+Norm2)/2+delta);
            h[i][k/2] -> GetYaxis()->SetTitle("gain");
            h[i][k/2] -> GetYaxis()->SetTitleOffset(1.);
            h[i][k/2] -> GetXaxis()->SetTitle("doca");
            c[i][k/2] = new TCanvas(histname.str().c_str(),histname.str().c_str(),800,600);
            h[i][k/2] -> Draw();
 
            cut.str("");
            cut<<"Ip_eangle=="<<i;
            cout<<"cut: "<<cut.str().c_str()<<endl;
 
            ss1.str("");
            if(NormFlag) ss1<<"Out_gain/"<<Norm1<<":Id_doca";
            else ss1<<"Out_gain"<<":Id_doca";
            tree1->Draw(ss1.str().c_str(),cut.str().c_str(),"sameLP");
            ss2.str(""); 
            if(NormFlag) ss2<<"Out_gain/"<<Norm2<<":Id_doca";
            else ss2<<"Out_gain"<<":Id_doca";
            tree2->Draw(ss2.str().c_str(),cut.str().c_str(),"sameLP");
            lgd[i][k/2] ->Draw();
        }
    }
 
    TFile* g = new TFile("gainVsDoca.root","RECREATE");
    for(int i=PadNo1;i<PadNo2;i++)
    {
        for(int k=0;k<FileNo;k+=2)
        {
            c[i][k/2] -> Write();
        }
    }
    g->Close();
 
    TCanvas* cc = new TCanvas("cc", "canvas", CanvasX, CanvasY);
    cc->Divide(Ipadx,Ipady);
    cc->Print("gainVsDoca.ps[");
    for(int i=PadNo1;i<PadNo2;i++)
    {
        for(int k=0;k<FileNo;k+=2)
        {
            cc->cd(k/2+1);
            gPad->Clear();
            c[i][k/2]->DrawClonePad();
            c[i][k/2] -> Close();
            lgd[i][k/2] ->Delete();
        }
        cc->Update();
        cc->Print("gainVsDoca.ps");
    }
    cc->Print("gainVsDoca.ps]");
    cc->Close();
 
    delete cc;
    for(int i=PadNo1;i<PadNo2;i++)
    {
        for(int k=0;k<FileNo;k+=2)
        {
            delete h[i][k/2];
            delete c[i][k/2];
        }
    }
}