Fork me on GitHub

Changeset 27 in svn


Ignore:
Timestamp:
Nov 13, 2008, 9:43:35 AM (16 years ago)
Author:
severine ovyn
Message:

Jet resol OK

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/interface/FuncDef.h

    r23 r27  
    7272}
    7373
    74 void GaussValues(TTree * Analyze,string histo,double &rms, double &mean, string min, string max)
     74void GaussValues(TTree * Analyze,string histo,double &rms, double &mean, string min,string max)
    7575{
    7676  string temp = histo;
    … …  
    7979
    8080  string nom = histo.erase(0,histo.find(">>")+2);
    81   TH1F *h = new TH1F(nom.c_str(),"",100,0,2);
    82   Analyze->Draw(temp.c_str(),mintemp.c_str(),maxtemp.c_str());
     81  TH1F *h = new TH1F(nom.c_str(),"",200,0,3);
     82
     83  string all = min + " && " + max;
     84  Analyze->Draw(temp.c_str(),all.c_str());
    8385  h->SetMarkerSize(0.6);
     86//  double MeanFix = ;
    8487  double MeanFix = h->GetMean();
    8588  double RangMin = h->GetMean()-h->GetRMS();
    8689  double RangMax = h->GetMean()+h->GetRMS();
    87   TF1 *Gauss = new TF1("Gauss","gausn",RangMin,RangMax);
     90  TF1 *Gauss = new TF1("Gauss","gaus",RangMin,RangMax);
    8891  Gauss->FixParameter(1,MeanFix);
    8992  h->Fit("Gauss","R");
    … …  
    9598  h->Draw("P");
    9699  h->GetXaxis()->SetTitle("E_{T}^{rec}/E_{T}^{MC} [GeV]");
    97 }
     100Gauss->Delete();
     101 
     102}
     103
     104void GaussValuesETmis(TTree * Analyze,string histo,double &rms, double &mean, string min,string max)
     105{
     106  string temp = histo;
     107  string mintemp = min;
     108  string maxtemp = max;
     109
     110  string nom = histo.erase(0,histo.find(">>")+2);
     111  TH1F *h = new TH1F(nom.c_str(),"",50,-300,300);
     112
     113  string all = min + " && " + max;
     114  Analyze->Draw(temp.c_str(),all.c_str());
     115  h->SetMarkerSize(0.6);
     116  double RangMin = h->GetMean()-h->GetRMS();
     117  double RangMax = h->GetMean()+h->GetRMS();
     118  TF1 *Gauss = new TF1("Gauss","gaus",RangMin,RangMax);
     119  h->Fit("Gauss","R");
     120  h->Fit("Gauss","RI");
     121  h->Fit("Gauss","RI");
     122  Double_t* params = Gauss->GetParameters();
     123  rms=params[2];
     124  mean=params[1];
     125  h->Draw("P");
     126  h->GetXaxis()->SetTitle("E_{T}^{rec}-E_{T}^{MC} [GeV]");
     127}
     128
    98129
    99130TH1F * HiggsMakeNormTH1F(Float_t Scale, Int_t numBin,Float_t minX,Float_t maxX, TTree * Analyze,string histo, Int_t Lcolor, Int_t Fcolor, Int_t Lstyle,Int_t width=2,bool Log=false)

  • trunk/routines/resolutions.C

    r23 r27  
    99#include "interface/FuncDef.h"
    1010
    11 void JetResol(TTree *Analyze)
     11void JetResol()
    1212{
    13    TCanvas *c1 = new TCanvas("c1","JET resol",100,100,600,450);
    14    c1->cd();
     13
     14   setTDRStyle();
     15   gROOT->Reset();
     16
     17   TFile *f1 = new TFile("JETResol.root","read");
     18   TTree *Analyze = (TTree*)f1->Get("Analysis");
    1519 
    1620   const Int_t numBin=7;
    17    double bins[numBin]={0,20,60,100,220,420,860};
     21   double bins[numBin]={0,20,40,60,80,100,220};
    1822   TProfile *ETSoverET = new TProfile("ETSoverET","Jet resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
    19    Analyze->Draw("JetPTResol.SmearePT:JetPTResol.NonSmearePT>>ETSoverET","JetPTResol.NonSmearePT>1");
    2023
    2124   double rms[numBin];
    2225   double mean[numBin];
    2326
    24    TCanvas *c2 = new TCanvas("c2","JET resol",0,0,1000,650);
    25    c2->cd(); int frame=0;
    26    c2->Divide(6,2);
     27   TCanvas *c1 = new TCanvas("c1","JET resol",0,0,1000,650);
     28   c1->cd(); int frame=0;
     29   c1->Divide(6,2);
    2730
    2831   float x[numBin-1];
    2932   float y[numBin-1];
     33   float ex[numBin-1];
     34   float ey[numBin-1];
     35
    3036   float finval=0;//valeur a remplir puis a fitter
    3137
    … …  
    3844        cout<<"binMin "<<binMin<<" binMax "<<binMax<<" bin center "<<binCenter<<endl;
    3945        char tempMin[500];
    40         sprintf(tempMin,"JetPTResol.NonSmearePT>%d",binMin);
     46        if(i==0)binMin=5;
     47        sprintf(tempMin,"JetPTResol.NonSmearePT > %d",binMin);
    4148        string mystringMin(tempMin);
    4249        char tempMax[500];
    43         sprintf(tempMax,"JetPTResol.NonSmearePT<%d",binMax);
     50        sprintf(tempMax,"JetPTResol.NonSmearePT < %d",binMax);
    4451        string mystringMax(tempMax);
    4552        char tempName[500];
    46         sprintf(tempName,"JetPTResol.SmearePT>>hETSoverET%d",i);
     53        sprintf(tempName,"(JetPTResol.SmearePT)>>hETSoverET%d",i);
    4754        string mystringName(tempName);
    4855
    49         c2->cd(++frame);
     56        c1->cd(++frame);
    5057        GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
    5158        x[i]=binCenter;
    5259        finval=rms[i]/mean[i];
    53         y[i]=(finval*2);
     60        y[i]=(finval*100);
     61        ex[i]=0;
     62        ey[i]=0;
    5463        cout<<"finval "<<finval<<" mean val "<<mean[i]<<" rms value "<<rms[i]<<endl;
    55 
    5664      }
    5765
    58    TCanvas *c3 = new TCanvas("c3","JET resol",100,100,600,450);
    59    c3->cd();
     66   TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
     67   c2->cd();
    6068
    61    TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,3000);
     69   TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,170);
    6270
    63    TGraph *gr11 = new TGraph((numBin-1),x,y);
     71   TGraphErrors *gr11 = new TGraphErrors((numBin-1),x,y,ex,ey);
    6472   gr11->Draw("AP");
     73   gr11->SetTitle("");
    6574   gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
    66    gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{MC})/<E_{T}^{rec}/E_{T}^{MC}>");
     75   gr11->GetYaxis()->SetRangeUser(0,50);
     76   gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{MC})_{fit}/<E_{T}^{rec}/E_{T}^{MC}>_{fit}");
    6777
    68    //fitfun->FixParameter(0,0.0541930);
    69    //fitfun->FixParameter(1,0.634908);
    70    gr11->Fit("user","RQ");
    71    gr11->Fit("user","RQ");
    72    gr11->Fit("user","RQ");
     78   Double_t* params = fitfun->GetParameters();
     79
    7380   gr11->Fit("user","R");
     81   gr11->Fit("user","RI");
     82   gr11->Fit("user","RI");
     83   gr11->Fit("user","RI");
     84   char tempResol[500];
     85   sprintf(tempResol,"#frac{#sigma(E_{T}^{rec}/E_{T}^{MC})}{<E_{T}^{rec}/E_{T}^{MC}>} = #frac{%f}{#sqrt{E_{T}^{MC}}} #oplus #frac{%f}{E_{T}^{MC}} #oplus %f",params[1]/100,params[2]/100,params[0]/100);
    7486
     87   TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
     88   leg1->Draw();
    7589
    76 
    77  delete fitfun;
    78 // delete Gauss;
     90   delete fitfun;
    7991}
    8092
    8193
    82 void ETmisResol(TTree *Analyze)
     94void ETmisResol()
    8395{
    8496
    85    TCanvas *c4 = new TCanvas("c4","ETmis resol",100,100,600,450);
    86    c4->cd();
    87  
    88    const Int_t numBin=10;
    89    double bins[numBin]={0,10,20,40,80,160,320,640,1080,2160};
    90    TProfile *ETSoverET = new TProfile("ETSoverET","ETmis resol",(numBin-1),bins,-10,10);
    91    Analyze->Draw("(ETmisResol.EtSmeare-ETmisResol.Et):(ETmisResol.Et)>>ETSoverET","(ETmisResol.Et)>1");
     97   setTDRStyle();
     98   gROOT->Reset();
    9299
    93    double rms[numBin];
    94    double mean[numBin];
     100   TFile *f1 = new TFile("PTMIS.root","read");
     101   TTree *Analyze = (TTree*)f1->Get("Analysis");
     102
     103   const Int_t numBin=6;
     104   double bins[numBin]={10,20,30,50,80,100};
     105   TProfile *ETSoverET = new TProfile("ETSoverET","ETmis resol",(numBin-1),bins,-1000,1000);
     106
     107   double rms[numBin-1];
     108   double mean[numBin-1];
    95109
    96110   TCanvas *c5 = new TCanvas("c5","gerrors2",0,0,1000,650);
    … …  
    98112   c5->Divide(6,2);
    99113
    100    float x[numBin-1];
    101    float y[numBin-1];
     114   double x[numBin-1];
     115   double y[numBin-1];
    102116
    103117   for ( int i=0; i<(numBin-1); i++)  // premiÚre bin : i ==1 et pas i == 0
    … …  
    109123        cout<<"binMin "<<binMin<<" binMax "<<binMax<<" bin center "<<binCenter<<endl;
    110124        char tempMin[500];
    111         sprintf(tempMin,"(ETmisResol.Et)>%d",binMin);
     125        sprintf(tempMin,"ETmisResol.Et>%d",binMin);
    112126        string mystringMin(tempMin);
    113127        char tempMax[500];
    114         sprintf(tempMax,"(ETmisResol.Et)<%d",binMax);
     128        sprintf(tempMax,"ETmisResol.Et<%d",binMax);
    115129        string mystringMax(tempMax);
     130
    116131        char tempName[500];
    117         sprintf(tempName,"(ETmisResol.EtSmeare-ETmisResol.Et)>>hETSoverET%d",i);
     132        sprintf(tempName,"ETmisResol.EtSmeare>>hETSoverET%d",i);
    118133        string mystringName(tempName);
    119134
    120135        c5->cd(++frame);
    121         GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
     136        GaussValuesETmis(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
    122137        x[i]=binCenter;
    123         mean[i]=ETSoverET->GetBinContent(i+1);
    124138        y[i]=rms[i];
    125139        cout<<" mean val "<<mean[i]<<" rms value "<<rms[i]<<endl;
    … …  
    130144   c6->cd();
    131145
    132    TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,3000);
     146   TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",0,200);
    133147
    134    TGraph *gr11 = new TGraph((numBin-1),x,y);
     148   TGraph *gr11 = new TGraph((numBin-1),x,rms);
    135149   gr11->Draw("AP");
    136    gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
    137    gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}>");
    138 
    139    //fitfun->FixParameter(0,0.0541930);
    140    //fitfun->FixParameter(1,0.634908);
     150   gr11->GetXaxis()->SetTitle("MC^{MET} [GeV]");
     151   gr11->GetYaxis()->SetTitle("#sigma(REC^{MET}-MC^{MET})");
    141152   gr11->Fit("user","RQ");
    142153   gr11->Fit("user","RQ");
    … …  
    144155   gr11->Fit("user","R");
    145156
     157   TCanvas *c7 = new TCanvas("c7","ETmis resolution",100,100,600,450);
     158   c7->cd();
     159
     160   TGraph *gr12 = new TGraph((numBin-1),x,mean);
     161   gr12->Draw("AP");
     162   gr12->GetXaxis()->SetTitle("MC^{MET} [GeV]");
     163   gr12->GetYaxis()->SetTitle("<(REC^{MET}-MC^{MET})>");
    146164
    147165
    148166 delete fitfun;
    149 // delete Gauss;
    150167}
    151168
    152169void General()
    153170{
    154   setTDRStyle();
    155   gROOT->Reset();
     171  JetResol();
     172  ETmisResol();
     173}
    156174
    157   TFile *f1 = new TFile("Smearing.root","read");
    158   TTree *Analyze = (TTree*)f1->Get("Analysis");
    159   JetResol(Analyze);
    160175
    161 }
Note: See TracChangeset for help on using the changeset viewer.