UsefullCode/Comparison/MakeChargeDistribution.C

0001
0002
0003 #include "TROOT.h"
0004 #include "TFile.h"
0005 #include "TDirectory.h"
0006 #include "TChain.h"
0007 #include "TObject.h"
0008 #include "TCanvas.h"
0009 #include "TMath.h"
0010 #include "TLegend.h"
0011 #include "TGraph.h"
0012 #include "TH1.h"
0013 #include "TH2.h"
0014 #include "TH3.h"
0015 #include "TTree.h"
0016 #include "TF1.h"
0017 #include "TPaveText.h"
0018 #include "PlotFunction.h"
0019
0020
0021 #include<vector>
0022 #include<tdrstyle.C>
0023
0024 void PlotMacro_Core(string input, string input2, string moduleName, string output);
0025
0026 int DataType = 2;
0027
0028 void MakeChargeDistribution(){
0029    gROOT->Reset();
0030    setTDRStyle();
0031    gStyle->SetPadTopMargin   (0.05);
0032    gStyle->SetPadBottomMargin(0.10);
0033    gStyle->SetPadRightMargin (0.18);
0034    gStyle->SetPadLeftMargin  (0.13);
0035    gStyle->SetTitleSize(0.04, "XYZ");
0036    gStyle->SetTitleXOffset(1.1);
0037    gStyle->SetTitleYOffset(1.35);
0038    gStyle->SetPalette(1);
0039    gStyle->SetCanvasColor(0);
0040    gStyle->SetBarOffset(0);
0041
0042
0043    system("mkdir -p Pictures/Charge");
0044    PlotMacro_Core("file:../../Data_Run_247252_to_247990_PCL/Gains_Tree.root", "file:../../Data_Run_247992_to_247992_PCL/Gains_Tree.root", "SiStripCalib"          , "Checks");
0045 }
0046
0047
0048 TF1* getPeakOfLandau(TH1* InputHisto, char* name, double LowRange=50, double HighRange=5400)
0049 {
0050    // perform fit with standard landau
0051    TF1* MyLandau = new TF1(name,"landau",LowRange, HighRange);
0052    MyLandau->SetParameter(1,300);
0053    InputHisto->Fit(MyLandau,"0QR WW");
0054    return MyLandau;
0055 }
0056
0057
0058 void PlotMacro_Core(string input, string input2, string moduleName, string output)
0059 {
0060    FILE* pFile;
0061    TCanvas* c1;
0062    TObject** Histos = new TObject*[10];
0063    std::vector<string> legend;
0064
0065    unsigned int  tree1_Index;
0066    unsigned int  tree1_DetId;
0067    unsigned char tree1_APVId;
0068    unsigned char tree1_SubDet;
0069    float         tree1_x;
0070    float         tree1_y;
0071    float         tree1_z;
0072    float         tree1_Eta;
0073    float         tree1_R;
0074    float         tree1_Phi;
0075    float         tree1_Thickness;
0076    float         tree1_FitMPV;
0077    float         tree1_FitMPVErr;
0078    float         tree1_FitWidth;
0079    float         tree1_FitWidthErr;
0080    float         tree1_FitChi2NDF;
0081    double        tree1_Gain;
0082    double        tree1_PrevGain;
0083    double        tree1_NEntries;
0084    bool          tree1_isMasked;
0085
0086
0087    unsigned int  tree2_Index;
0088    unsigned int  tree2_DetId;
0089    unsigned char tree2_APVId;
0090    unsigned char tree2_SubDet;
0091    float         tree2_x;
0092    float         tree2_y;
0093    float         tree2_z;
0094    float         tree2_Eta;
0095    float         tree2_R;
0096    float         tree2_Phi;
0097    float         tree2_Thickness;
0098    float         tree2_FitMPV;
0099    float         tree2_FitMPVErr;
0100    float         tree2_FitWidth;
0101    float         tree2_FitWidthErr;
0102    float         tree2_FitChi2NDF;
0103    double        tree2_Gain;
0104    double        tree2_PrevGain;
0105    double        tree2_NEntries;
0106    bool          tree2_isMasked;
0107
0108
0109    TFile* f1     = new TFile(input.c_str());
0110    TTree *t1     = (TTree*)GetObjectFromPath(f1,moduleName+"/APVGain");
0111    TH2D* ChargeDistrib1  = (TH2D*)GetObjectFromPath(f1,moduleName+"/Charge_Vs_Index");
0112    t1->SetBranchAddress("Index"             ,&tree1_Index      );
0113    t1->SetBranchAddress("DetId"             ,&tree1_DetId      );
0114    t1->SetBranchAddress("APVId"             ,&tree1_APVId      );
0115    t1->SetBranchAddress("SubDet"            ,&tree1_SubDet     );
0116    t1->SetBranchAddress("x"                 ,&tree1_x          );
0117    t1->SetBranchAddress("y"                 ,&tree1_y          );
0118    t1->SetBranchAddress("z"                 ,&tree1_z          );
0119    t1->SetBranchAddress("Eta"               ,&tree1_Eta        );
0120    t1->SetBranchAddress("R"                 ,&tree1_R          );
0121    t1->SetBranchAddress("Phi"               ,&tree1_Phi        );
0122    t1->SetBranchAddress("Thickness"         ,&tree1_Thickness  );
0123    t1->SetBranchAddress("FitMPV"            ,&tree1_FitMPV     );
0124    t1->SetBranchAddress("FitMPVErr"         ,&tree1_FitMPVErr  );
0125    t1->SetBranchAddress("FitWidth"          ,&tree1_FitWidth   );
0126    t1->SetBranchAddress("FitWidthErr"       ,&tree1_FitWidthErr);
0127    t1->SetBranchAddress("FitChi2NDF"        ,&tree1_FitChi2NDF );
0128    t1->SetBranchAddress("Gain"              ,&tree1_Gain       );
0129    t1->SetBranchAddress("PrevGain"          ,&tree1_PrevGain   );
0130    t1->SetBranchAddress("NEntries"          ,&tree1_NEntries   );
0131    t1->SetBranchAddress("isMasked"          ,&tree1_isMasked   );
0132
0133
0134    TFile* f2     = new TFile(input2.c_str());
0135    TTree *t2     = (TTree*)GetObjectFromPath(f2,moduleName+"/APVGain");
0136    TH2D* ChargeDistrib2  = (TH2D*)GetObjectFromPath(f2,moduleName+"/Charge_Vs_Index");
0137    t2->SetBranchAddress("Index"             ,&tree2_Index      );
0138    t2->SetBranchAddress("DetId"             ,&tree2_DetId      );
0139    t2->SetBranchAddress("APVId"             ,&tree2_APVId      );
0140    t2->SetBranchAddress("SubDet"            ,&tree2_SubDet     );
0141    t2->SetBranchAddress("x"                 ,&tree2_x          );
0142    t2->SetBranchAddress("y"                 ,&tree2_y          );
0143    t2->SetBranchAddress("z"                 ,&tree2_z          );
0144    t2->SetBranchAddress("Eta"               ,&tree2_Eta        );
0145    t2->SetBranchAddress("R"                 ,&tree2_R          );
0146    t2->SetBranchAddress("Phi"               ,&tree2_Phi        );
0147    t2->SetBranchAddress("Thickness"         ,&tree2_Thickness  );
0148    t2->SetBranchAddress("FitMPV"            ,&tree2_FitMPV     );
0149    t2->SetBranchAddress("FitMPVErr"         ,&tree2_FitMPVErr  );
0150    t2->SetBranchAddress("FitWidth"          ,&tree2_FitWidth   );
0151    t2->SetBranchAddress("FitWidthErr"       ,&tree2_FitWidthErr);
0152    t2->SetBranchAddress("FitChi2NDF"        ,&tree2_FitChi2NDF );
0153    t2->SetBranchAddress("Gain"              ,&tree2_Gain       );
0154    t2->SetBranchAddress("PrevGain"          ,&tree2_PrevGain   );
0155    t2->SetBranchAddress("NEntries"          ,&tree2_NEntries   );
0156    t2->SetBranchAddress("isMasked"          ,&tree2_isMasked   );
0157
0158    pFile = fopen("GainDiff.txt","w");
0159
0160    unsigned int PreviousId = 0;
0161    unsigned int NAPV       = 0;
0162    double       MPV1       = 0;
0163    double       MPV2       = 0;
0164
0165    double Min = 10000;
0166    double Max =-10000;
0167
0168    printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
0169    printf("Looping on the Tree          :");
0170    int TreeStep = t1->GetEntries()/50;if(TreeStep==0)TreeStep=1;
0171    for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
0172       if(ientry%TreeStep==0){printf(".");fflush(stdout);}
0173       t1->GetEntry(ientry);
0174       t2->GetEntry(ientry);
0175
0176
0177        if(tree2_FitMPV>0 and tree2_FitMPV<150){
0178           TH1D* Proj1         = ChargeDistrib1 ->ProjectionY("proj1" ,tree1_Index, tree1_Index, "e");
0179           TH1D* Proj2         = ChargeDistrib2 ->ProjectionY("proj2" ,tree2_Index, tree2_Index, "e");
0180
0181           TCanvas* c1 = new TCanvas("c1", "c1", 600,600);
0182 //          c1->SetLogy(true);
0183
0184
0185           TH1D* frame = new TH1D("frame", "frame", 1,0, 1000);
0186           frame->GetXaxis()->SetNdivisions(505);
0187           frame->SetTitle("");
0188           frame->SetStats(kFALSE);
0189           frame->GetXaxis()->SetTitle("charge/path-length (ADC/mm)");
0190           frame->GetYaxis()->SetTitle("#clusters");
0191           frame->SetMaximum(std::max(tree1_NEntries, tree2_NEntries)*0.04);
0192           frame->SetMinimum(0.0);
0193           frame->GetYaxis()->SetTitleOffset(1.50);
0194           frame->Draw();
0195
0196           //Proj1->Scale(1.0/Proj1->Integral());
0197           Proj1->SetLineColor(2);
0198           Proj1->SetLineWidth(2);
0199           Proj1->Draw("H same");
0200
0201           TF1* Fit1 = new TF1("MyLandau1","landau", 0, 2000);
0202           Fit1->SetParameters(tree1_NEntries/tree1_FitWidth, tree1_FitMPV, tree1_FitWidth);
0203           Fit1->SetLineColor(2);
0204           Fit1->SetLineWidth(2);
0205           Fit1->Draw("L same");
0206
0207           TF1* Fit1b = getPeakOfLandau(Proj1, "MyLandau1bis");
0208           Fit1b->SetLineColor(8);
0209           Fit1b->SetLineWidth(2);
0210           Fit1b->Draw("L same");
0211
0212           printf("period1 MPV vs new MPV --> %6.2f  vs %6.2f\n", tree1_FitMPV, Fit1b->GetParameter(1));
0213
0214
0215           //Proj2->Scale(1.0/Proj2->Integral());
0216           Proj2->SetLineColor(4);
0217           Proj2->SetLineWidth(2);
0218           Proj2->Draw("H same");
0219
0220           TF1* Fit2 = new TF1("MyLandau2","landau", 0, 2000);
0221           Fit2->SetParameters(tree2_NEntries/tree2_FitWidth, tree2_FitMPV, tree2_FitWidth);
0222           Fit2->SetLineColor(4);
0223           Fit2->SetLineWidth(2);
0224           Fit2->Draw("L same");
0225
0226
0227           TF1* Fit2b = getPeakOfLandau(Proj2, "MyLandau2bis");
0228           Fit2b->SetLineColor(8);
0229           Fit2b->SetLineWidth(2);
0230           Fit2b->Draw("L same");
0231
0232           printf("period2 MPV vs new MPV --> %6.2f  vs %6.2f\n", tree2_FitMPV, Fit2b->GetParameter(1));
0233
0234           char buffer[256];
0235           sprintf(buffer, "Pictures/Charge/SubDet%i_Id%i_Apv%i_N%i.png",tree1_SubDet, tree1_DetId, tree1_APVId, (int)tree1_NEntries);
0236           c1->SaveAs(buffer);
0237        }
0238
0239
0240
0241
0242       if(tree1_APVId==0 && PreviousId>0){
0243          double Mean = (MPV2/NAPV) / (MPV1/NAPV);
0244          if(Mean<Min)Min=Mean;
0245          if(Mean>Max)Max=Mean;
0246          if(NAPV>0) fprintf(pFile,"%i  %f\n",PreviousId,Mean);
0247          NAPV=0; MPV1=0; MPV2=0;
0248       }
0249       PreviousId = tree1_DetId;
0250       if(tree1_FitMPV<=0 || tree2_FitMPV<=0)continue;
0251       NAPV++;
0252       MPV1+=tree1_FitMPV;
0253       MPV2+=tree2_FitMPV;
0254    }printf("\n");
0255    fclose(pFile);
0256
0257    printf("Min=%f - Max=%f\n",Min,Max);
0258 }