Fork me on GitHub

Changeset 23 in svn


Ignore:
Timestamp:
Nov 7, 2008, 6:37:38 PM (16 years ago)
Author:
severine ovyn
Message:

Working on the resolution

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/Resolutions.cpp

    r19 r23  
    172172  vector<fastjet::PseudoJet> inclusive_jetsS;
    173173  vector<fastjet::PseudoJet> sorted_jetsS;
     174  vector<PhysicsTower> towers;
    174175
    175176  fastjet::JetDefinition jet_def;
    … …  
    180181  // set up a CDF midpoint jet definition
    181182  #ifdef ENABLE_PLUGIN_CDFCONES
    182   plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
     183  plugins = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
    183184  jet_def = fastjet::JetDefinition(plugins);
    184185  #else
    … …  
    200201  for(entry = 0; entry < allEntries; ++entry)
    201202    {
     203      TLorentzVector PTmisS(0,0,0,0);
    202204      TLorentzVector PTmis(0,0,0,0);
    203205      treeReader->ReadEntry(entry);
    … …  
    216218      inclusive_jetsS.clear();
    217219      sorted_jetsS.clear();
     220      towers.clear();
    218221
    219222      // Loop over all particles in event
    … …  
    238241              )
    239242            ) {
     243                if(pid != pMU)PTmis = PTmis + genMomentum;//ptmis
    240244                switch(pid) {
    241245                   
    … …  
    264268                // all final particles but muons and neutrinos 
    265269                // for calorimetric towers and mission PT
    266                 if(genMomentum.E()!=0) PTmis = PTmis + genMomentum;//ptmis
    267 
    268                 if(pid != pMU && genMomentum.Pt() > DET->PT_TRACKS_MIN)
     270                //if(genMomentum.E()!=0) PTmis = PTmis + genMomentum;//ptmis
     271
     272                if(pid != pMU)
    269273                  {
     274                     PTmisS = PTmisS + genMomentum;
     275                     towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
    270276                     input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
    271277                  }
    … …  
    289295        } // while
    290296
     297  TLorentzVector toWerS(0,0,0,0);
     298  //calcul de ETmis au depart des tours calorimetriques..
     299/*  for(unsigned int i=0; i < towers.size(); i++) {
     300      if(towers[i].fourVector.pt() < 0.5) continue;
     301          toWerS.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
     302          PTmisS = PTmisS + toWerS;
     303          }
     304*/
     305//PTmis = (0,0,0,0)-PTmis;
     306//PTmisS = (0,0,0,0)-PTmisS;
     307
     308//float ET=PTmis.Et()<<endl;
     309//float ETS=PTmisS.Et()<<endl;
     310
     311cout<<"Ptmis "<<PTmis.Et()<<" PTmis smeare "<<PTmisS.Et()<<endl;
     312cout<<"Ptmis "<<(-PTmis).Pt()<<" PTmis smeare "<<(-PTmisS).Pt()<<endl;
     313
     314  elementEtmis= (ETMIS*) branchetmis->NewEntry();
     315  elementEtmis->Et = (PTmis).Et();
     316  elementEtmis->EtSmeare = (PTmisS).Et();
     317
     318
    291319
    292320      //*****************************
    … …  
    314342        if(JETSm.Pt()>1)
    315343          {
     344            float NonSmeareEt=(JET.E()*sin(JET.Theta()));
     345            float SmeareEt=(JETSm.E()*sin(JETSm.Theta()));
     346
    316347            elementJet= (RESOLJET*) branchjet->NewEntry();
    317             elementJet->NonSmearePT = JET.Pt();
    318             elementJet->SmearePT = (JETSm.Pt()/JET.Pt());
    319             /*cout<<"valeur obtenue "<<JETSm.Pt()/JET.Pt()<<endl;
    320             cout<<" pt non smeare "<<JET.Pt()<<" phi "<<JET.Phi()<<" eta "<<JET.Eta()<<endl;
    321             cout<<"pt smeare "<<JETSm.Pt()<<" phi "<<JETSm.Phi()<<" eta "<<JETSm.Eta()<<endl;
    322             cout<<"************"<<endl;
    323 */
     348            elementJet->NonSmearePT = JET.Et();
     349            elementJet->SmearePT = JETSm.Et()/JET.Et();
     350         
    324351          }
    325352

  • trunk/interface/FuncDef.h

    r17 r23  
    7979
    8080  string nom = histo.erase(0,histo.find(">>")+2);
    81   TH1F *h = new TH1F(nom.c_str(),"",50,0.5,1.5);
     81  TH1F *h = new TH1F(nom.c_str(),"",100,0,2);
    8282  Analyze->Draw(temp.c_str(),mintemp.c_str(),maxtemp.c_str());
    8383  h->SetMarkerSize(0.6);
    84   TF1 *Gauss = new TF1("Gauss","gausn",0.8,1.3);
     84  double MeanFix = h->GetMean();
     85  double RangMin = h->GetMean()-h->GetRMS();
     86  double RangMax = h->GetMean()+h->GetRMS();
     87  TF1 *Gauss = new TF1("Gauss","gausn",RangMin,RangMax);
     88  Gauss->FixParameter(1,MeanFix);
    8589  h->Fit("Gauss","R");
    8690  h->Fit("Gauss","RI");
    … …  
    9296  h->GetXaxis()->SetTitle("E_{T}^{rec}/E_{T}^{MC} [GeV]");
    9397}
    94 
    9598
    9699TH1F * 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

    r9 r23  
    99#include "interface/FuncDef.h"
    1010
    11 void Plot()
     11void JetResol(TTree *Analyze)
    1212{
    13 
    14   setTDRStyle();
    15   gROOT->Reset();
    16 
    17   TFile *f1 = new TFile("Smearing.root","read");
    18   TTree *Analyze = (TTree*)f1->Get("Analysis");
    19 
    20 
    21 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    22 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    23 
    24    TCanvas *c1 = new TCanvas("c1","gerrors2",100,100,600,450);
     13   TCanvas *c1 = new TCanvas("c1","JET resol",100,100,600,450);
    2514   c1->cd();
    2615 
    27    const Int_t numBin=5;
    28 //   double bins[numBin]={5,10,25,55,85,115,205,400};
    29    double bins[5]={10,30,50,80,120};
    30    TProfile *ETSoverET = new TProfile("ETSoverET","Resolution des electrons",(numBin-1),bins,-10,10);
     16   const Int_t numBin=7;
     17   double bins[numBin]={0,20,60,100,220,420,860};
     18   TProfile *ETSoverET = new TProfile("ETSoverET","Jet resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
    3119   Analyze->Draw("JetPTResol.SmearePT:JetPTResol.NonSmearePT>>ETSoverET","JetPTResol.NonSmearePT>1");
    3220
    3321   double rms[numBin];
    3422   double mean[numBin];
    35    //TF1 *Gauss = new TF1("Gauss","gaus",0,4);
    3623
    37    TCanvas *c2 = new TCanvas("c2","gerrors2",0,0,1000,650);
     24   TCanvas *c2 = new TCanvas("c2","JET resol",0,0,1000,650);
    3825   c2->cd(); int frame=0;
    3926   c2->Divide(6,2);
    40 //   c2->Divide(6,2);
    4127
    42    float x[numBin];
    43    float y[numBin];
     28   float x[numBin-1];
     29   float y[numBin-1];
    4430   float finval=0;//valeur a remplir puis a fitter
    4531
    46    for ( int i=0; i<numBin; i++)  // premiÚre bin : i ==1 et pas i == 0
     32   for ( int i=0; i<(numBin-1); i++)  // premiÚre bin : i ==1 et pas i == 0
    4733      {
    4834        TAxis *xaxis = ETSoverET->GetXaxis();
    … …  
    6450        GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
    6551        x[i]=binCenter;
    66         mean[i]=ETSoverET->GetBinContent(i+1);
    67    //     rms[i]=ETSoverET->GetBinError(i+1);
    6852        finval=rms[i]/mean[i];
    69         y[i]=finval;
     53        y[i]=(finval*2);
    7054        cout<<"finval "<<finval<<" mean val "<<mean[i]<<" rms value "<<rms[i]<<endl;
    7155
    7256      }
    7357
    74 y[0]=0.2;
    75 
    76    TCanvas *c3 = new TCanvas("c3","gerrors2",100,100,600,450);
     58   TCanvas *c3 = new TCanvas("c3","JET resol",100,100,600,450);
    7759   c3->cd();
    7860
    79    TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",10,500);
     61   TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,3000);
    8062
    81    TGraph *gr11 = new TGraph((numBin),x,y);
     63   TGraph *gr11 = new TGraph((numBin-1),x,y);
    8264   gr11->Draw("AP");
    8365   gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
    … …  
    9678// delete Gauss;
    9779}
     80
     81
     82void ETmisResol(TTree *Analyze)
     83{
     84
     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");
     92
     93   double rms[numBin];
     94   double mean[numBin];
     95
     96   TCanvas *c5 = new TCanvas("c5","gerrors2",0,0,1000,650);
     97   c5->cd(); int frame=0;
     98   c5->Divide(6,2);
     99
     100   float x[numBin-1];
     101   float y[numBin-1];
     102
     103   for ( int i=0; i<(numBin-1); i++)  // premiÚre bin : i ==1 et pas i == 0
     104      {
     105        TAxis *xaxis = ETSoverET->GetXaxis();
     106        float binCenter = xaxis->GetBinCenter(i+1);
     107        int binMin = (int)xaxis->GetBinLowEdge(i+1);
     108        int binMax = (int)xaxis->GetBinUpEdge(i+1);
     109        cout<<"binMin "<<binMin<<" binMax "<<binMax<<" bin center "<<binCenter<<endl;
     110        char tempMin[500];
     111        sprintf(tempMin,"(ETmisResol.Et)>%d",binMin);
     112        string mystringMin(tempMin);
     113        char tempMax[500];
     114        sprintf(tempMax,"(ETmisResol.Et)<%d",binMax);
     115        string mystringMax(tempMax);
     116        char tempName[500];
     117        sprintf(tempName,"(ETmisResol.EtSmeare-ETmisResol.Et)>>hETSoverET%d",i);
     118        string mystringName(tempName);
     119
     120        c5->cd(++frame);
     121        GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
     122        x[i]=binCenter;
     123        mean[i]=ETSoverET->GetBinContent(i+1);
     124        y[i]=rms[i];
     125        cout<<" mean val "<<mean[i]<<" rms value "<<rms[i]<<endl;
     126
     127      }
     128
     129   TCanvas *c6 = new TCanvas("c6","ETmis resolution",100,100,600,450);
     130   c6->cd();
     131
     132   TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,3000);
     133
     134   TGraph *gr11 = new TGraph((numBin-1),x,y);
     135   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);
     141   gr11->Fit("user","RQ");
     142   gr11->Fit("user","RQ");
     143   gr11->Fit("user","RQ");
     144   gr11->Fit("user","R");
     145
     146
     147
     148 delete fitfun;
     149// delete Gauss;
     150}
     151
     152void General()
     153{
     154  setTDRStyle();
     155  gROOT->Reset();
     156
     157  TFile *f1 = new TFile("Smearing.root","read");
     158  TTree *Analyze = (TTree*)f1->Get("Analysis");
     159  JetResol(Analyze);
     160
     161}
Note: See TracChangeset for help on using the changeset viewer.