- Timestamp:
- Nov 7, 2008, 6:37:38 PM (16 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Resolutions.cpp
r19 r23 172 172 vector<fastjet::PseudoJet> inclusive_jetsS; 173 173 vector<fastjet::PseudoJet> sorted_jetsS; 174 vector<PhysicsTower> towers; 174 175 175 176 fastjet::JetDefinition jet_def; … … 180 181 // set up a CDF midpoint jet definition 181 182 #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); 183 184 jet_def = fastjet::JetDefinition(plugins); 184 185 #else … … 200 201 for(entry = 0; entry < allEntries; ++entry) 201 202 { 203 TLorentzVector PTmisS(0,0,0,0); 202 204 TLorentzVector PTmis(0,0,0,0); 203 205 treeReader->ReadEntry(entry); … … 216 218 inclusive_jetsS.clear(); 217 219 sorted_jetsS.clear(); 220 towers.clear(); 218 221 219 222 // Loop over all particles in event … … 238 241 ) 239 242 ) { 243 if(pid != pMU)PTmis = PTmis + genMomentum;//ptmis 240 244 switch(pid) { 241 245 … … 264 268 // all final particles but muons and neutrinos 265 269 // for calorimetric towers and mission PT 266 if(genMomentum.E()!=0) PTmis = PTmis + genMomentum;//ptmis267 268 if(pid != pMU && genMomentum.Pt() > DET->PT_TRACKS_MIN)270 //if(genMomentum.E()!=0) PTmis = PTmis + genMomentum;//ptmis 271 272 if(pid != pMU) 269 273 { 274 PTmisS = PTmisS + genMomentum; 275 towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()))); 270 276 input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 271 277 } … … 289 295 } // while 290 296 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 311 cout<<"Ptmis "<<PTmis.Et()<<" PTmis smeare "<<PTmisS.Et()<<endl; 312 cout<<"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 291 319 292 320 //***************************** … … 314 342 if(JETSm.Pt()>1) 315 343 { 344 float NonSmeareEt=(JET.E()*sin(JET.Theta())); 345 float SmeareEt=(JETSm.E()*sin(JETSm.Theta())); 346 316 347 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 324 351 } 325 352 -
trunk/interface/FuncDef.h
r17 r23 79 79 80 80 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); 82 82 Analyze->Draw(temp.c_str(),mintemp.c_str(),maxtemp.c_str()); 83 83 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); 85 89 h->Fit("Gauss","R"); 86 90 h->Fit("Gauss","RI"); … … 92 96 h->GetXaxis()->SetTitle("E_{T}^{rec}/E_{T}^{MC} [GeV]"); 93 97 } 94 95 98 96 99 TH1F * 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 9 9 #include "interface/FuncDef.h" 10 10 11 void Plot()11 void JetResol(TTree *Analyze) 12 12 { 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); 25 14 c1->cd(); 26 15 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); 31 19 Analyze->Draw("JetPTResol.SmearePT:JetPTResol.NonSmearePT>>ETSoverET","JetPTResol.NonSmearePT>1"); 32 20 33 21 double rms[numBin]; 34 22 double mean[numBin]; 35 //TF1 *Gauss = new TF1("Gauss","gaus",0,4);36 23 37 TCanvas *c2 = new TCanvas("c2"," gerrors2",0,0,1000,650);24 TCanvas *c2 = new TCanvas("c2","JET resol",0,0,1000,650); 38 25 c2->cd(); int frame=0; 39 26 c2->Divide(6,2); 40 // c2->Divide(6,2);41 27 42 float x[numBin ];43 float y[numBin ];28 float x[numBin-1]; 29 float y[numBin-1]; 44 30 float finval=0;//valeur a remplir puis a fitter 45 31 46 for ( int i=0; i< numBin; i++) // premiÚre bin : i ==1 et pas i == 032 for ( int i=0; i<(numBin-1); i++) // premiÚre bin : i ==1 et pas i == 0 47 33 { 48 34 TAxis *xaxis = ETSoverET->GetXaxis(); … … 64 50 GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax); 65 51 x[i]=binCenter; 66 mean[i]=ETSoverET->GetBinContent(i+1);67 // rms[i]=ETSoverET->GetBinError(i+1);68 52 finval=rms[i]/mean[i]; 69 y[i]= finval;53 y[i]=(finval*2); 70 54 cout<<"finval "<<finval<<" mean val "<<mean[i]<<" rms value "<<rms[i]<<endl; 71 55 72 56 } 73 57 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); 77 59 c3->cd(); 78 60 79 TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1 0,500);61 TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,3000); 80 62 81 TGraph *gr11 = new TGraph((numBin ),x,y);63 TGraph *gr11 = new TGraph((numBin-1),x,y); 82 64 gr11->Draw("AP"); 83 65 gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]"); … … 96 78 // delete Gauss; 97 79 } 80 81 82 void 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 152 void 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.