#include <sys/stat.h>
#include <stdlib.h>
#include <errno.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <stdexcept>
#include "TROOT.h"
#include "TH2.h"
#include "TFile.h"
#include "TRandom3.h"
#include "TTree.h"
#include "TStopwatch.h"
#include "TRadCor.h"
#include "TEPCrossPart.h"
#include "TMuCrossPart.h"
#include "TPiCrossPart.h"
#include "TKnCrossPart.h"
#include "TKcCrossPart.h"
#include "TGGCrossPart.h"
#include "TRadGlobal.h"
#include "TKinemCut.h"
Go to the source code of this file.
Classes | |
class | EVENT |
Functions | |
int | GetOption (unsigned int, char **) |
void | hist_sample (Int_t i=0) |
void | hist_sample_angle (Int_t i=0) |
void | hist_sample_comp () |
int | main (int argc, char **argv) |
Variables | |
Double_t | mass [] |
char * | proc_name [] |
TRadCor * | gRad |
TRadGlobal * | gGlobal |
TKinemCut * | gCut |
TConstants * | gConst |
TStopwatch * | timer |
Double_t | EBeam |
Int_t | NSim |
Int_t | NRad |
Int_t | IsTuple |
Int_t | IsSmear |
Int_t | IsBrem |
Int_t | IsHardPhoton |
Int_t | IsNoVacPol |
Int_t | IsVerbose |
Int_t | IsFSR |
Int_t | IsAngleHist |
Int_t | IsDTSmear |
string | fname |
string | RandomFileName |
double | smear_par [9] |
int | type |
const int | MaxList = 22 |
bool | InList [MaxList] |
TH1D * | gH [MaxList] |
TH1D * | gH2 [MaxList] |
int GetOption | ( | unsigned | int, | |
char ** | ||||
) |
Definition at line 210 of file hist.C.
References al, at(), TKinemCut::AverageTheta(), cos(), TKinemCut::CosPsi(), TKinemCut::DeltaPhi(), TKinemCut::DeltaTheta(), dt, EBeam, TKinemCut::EMin(), fname, gConst, gCut, TRadGlobal::Get_RelativeError(), TRadGlobal::Get_TotalError(), gGlobal, genRecEmupikp::i, TKinemCut::Init(), TRadGlobal::Init(), InList, IsAngleHist, IsBrem, IsDTSmear, IsFSR, IsHardPhoton, IsNoVacPol, IsSmear, IsTuple, IsVerbose, ganga-rec::j, M_PI, MaxList, NRad, NSim, TKinemCut::PAverage(), TKinemCut::PCross(), TKinemCut::Print(), TRadGlobal::Print(), TConstants::Print(), RandomFileName, TRadGlobal::Set_dE_Abs(), TRadGlobal::Set_E(), TRadGlobal::Set_RelativeError(), TRadGlobal::Set_Theta0_Rel(), TRadGlobal::Set_ThetaInt(), TRadGlobal::Set_ThetaMin(), TRadGlobal::Set_TotalError(), TRadGlobal::Set_Type(), TConstants::SetAlphaScale(), smear_par, TKinemCut::ThetaRangeM(), TKinemCut::ThetaRangeP(), and type.
Referenced by main().
00210 { 00211 EBeam = 185.; 00212 NRad = 1000; 00213 NSim = 0; 00214 IsTuple = 0; 00215 IsSmear = 0; 00216 IsBrem = 0; 00217 IsHardPhoton = 1; 00218 IsNoVacPol = 0; 00219 IsVerbose = 0; 00220 IsFSR = 1; 00221 IsAngleHist = 0; 00222 IsDTSmear = 0; 00223 type = 0; 00224 fname = ""; 00225 RandomFileName = ".random_stat.root"; 00226 int proc = 0; 00227 double pc = 0.; 00228 double de = -1.; 00229 double nt0 = -1.; 00230 double dt = 0.25; 00231 double dp = 0.15; 00232 double at = 1.1; 00233 double td = -1.; 00234 double am = 90.; 00235 double cm = 90.; 00236 double em = 50.; 00237 double ti = 0.473; 00238 double al = 1.; 00239 double thm_min = -1.; 00240 double thp_min = -1.; 00241 double thm_max = -1.; 00242 double thp_max = -1.; 00243 double te = 0.1; 00244 double re = 0.1; 00245 for(int j = 0; j<18; j++) InList[j] = true; 00246 for(int j = 18; j<MaxList; j++) InList[j] = false; 00247 00248 for (unsigned int i = 1; i < argc; i++){ 00249 if(strlen(argv[i])>1) 00250 if (!strcmp(argv[i], "-e")){ 00251 EBeam = atof(argv[++i]); 00252 00253 }else if (!strcmp(argv[i], "-ns")){ 00254 NSim = atoi(argv[++i]); 00255 00256 }else if (!strcmp(argv[i], "-type")){ 00257 type = atoi(argv[++i]); 00258 00259 }else if (!strcmp(argv[i], "-v")){ 00260 IsVerbose = 1; 00261 00262 }else if (!strcmp(argv[i], "-nr")){ 00263 NRad = atoi(argv[++i]); 00264 00265 }else if (!strcmp(argv[i], "-de")){ 00266 de = atof(argv[++i]); 00267 00268 }else if (!strcmp(argv[i], "-nt0")){ 00269 nt0 = atof(argv[++i]); 00270 00271 }else if (!strcmp(argv[i], "-dt")){ 00272 dt = atof(argv[++i]); 00273 00274 }else if (!strcmp(argv[i], "-tc")){ 00275 at = atof(argv[++i]); 00276 00277 }else if (!strcmp(argv[i], "-dp")){ 00278 dp = atof(argv[++i]); 00279 00280 }else if (!strcmp(argv[i], "-td")){ 00281 td = atof(argv[++i]); 00282 00283 }else if (!strcmp(argv[i], "-ti")){ 00284 ti = atof(argv[++i]); 00285 00286 }else if (!strcmp(argv[i], "-al")){ 00287 al = atof(argv[++i]); 00288 00289 }else if (!strcmp(argv[i], "-am")){ 00290 am = atof(argv[++i]); 00291 00292 }else if (!strcmp(argv[i], "-cm")){ 00293 cm = atof(argv[++i]); 00294 00295 }else if (!strcmp(argv[i], "-em")){ 00296 em = atof(argv[++i]); 00297 00298 }else if (!strcmp(argv[i], "-pc")){ 00299 pc = atof(argv[++i]); 00300 00301 }else if (!strcmp(argv[i], "-thm")){ 00302 thm_min = atof(argv[++i]); 00303 thm_max = atof(argv[++i]); 00304 00305 }else if (!strcmp(argv[i], "-thp")){ 00306 thp_min = atof(argv[++i]); 00307 thp_max = atof(argv[++i]); 00308 00309 }else if (!strcmp(argv[i], "-p")){ 00310 proc = atoi(argv[++i]); 00311 00312 }else if (!strcmp(argv[i], "-er")){ 00313 te = atof(argv[++i]); 00314 00315 }else if (!strcmp(argv[i], "-rer")){ 00316 re = atof(argv[++i]); 00317 00318 }else if (!strcmp(argv[i], "-tup")){ 00319 IsTuple = 1; 00320 00321 }else if (!strcmp(argv[i], "-brem")){ 00322 IsBrem = 1; 00323 00324 }else if (!strcmp(argv[i], "-fname")){ 00325 fname = argv[++i]; 00326 00327 }else if (!strcmp(argv[i], "-rname")){ 00328 RandomFileName = argv[++i]; 00329 00330 }else if (!strcmp(argv[i], "-hsm")){ 00331 IsSmear = 1; 00332 00333 }else if (!strcmp(argv[i], "-ah")){ 00334 IsAngleHist = 1; 00335 00336 }else if (!strcmp(argv[i], "-hp")){ 00337 IsHardPhoton = 1; 00338 00339 }else if (!strcmp(argv[i], "-nhp")){ 00340 IsHardPhoton = 0; 00341 00342 }else if (!strcmp(argv[i], "-nfsr")){ 00343 IsFSR = 0; 00344 00345 }else if (!strcmp(argv[i], "-nvp")){ 00346 IsNoVacPol = 1; 00347 00348 }else if (!strcmp(argv[i], "-exclude")){ 00349 for(int j = 0; j<MaxList; j++) InList[j] = true; 00350 while((i+1 < argc) && bool(strncmp(argv[i+1],"-",1))){ 00351 int j = atoi(argv[++i]); 00352 if( j<MaxList ) InList[j] = false; 00353 } 00354 00355 }else if (!strcmp(argv[i], "-include")){ 00356 for(int j = 0; j<MaxList; j++) InList[j] = false; 00357 while((i+1 < argc) && bool(strncmp(argv[i+1],"-",1))){ 00358 int j = atoi(argv[++i]); 00359 if( j<MaxList ) InList[j] = true; 00360 } 00361 00362 }else if (!strcmp(argv[i], "-asm")){ 00363 IsDTSmear = 1; 00364 int k = 0; 00365 while((i+1 < argc) && bool(strncmp(argv[i+1],"-",1))){ 00366 smear_par[k] = atof(argv[++i]); 00367 k++; 00368 } 00369 if ( k != 6 && k != 9 ) throw logic_error("Number of parameters for smearing must be 6 or 9!"); 00370 if ( k<9 ) { 00371 smear_par[6] = 1.; 00372 smear_par[7] = 0.; 00373 smear_par[8] = 0.; 00374 } 00375 }else { 00376 cerr<<"Unknown option: "<<argv[i]<<endl; 00377 } 00378 } 00379 00380 if( ( EBeam < 100 || EBeam > 25000 ) && !IsNoVacPol ){ 00381 cerr<<"Invalid value of beam energy:"<<EBeam<<endl; 00382 exit(1); 00383 } 00384 00385 gGlobal->Set_TotalError(te); 00386 00387 gGlobal->Set_RelativeError(re); 00388 00389 gGlobal->Set_Type(proc); 00390 00391 gGlobal->Set_E(EBeam); 00392 00393 gGlobal->Set_ThetaInt(ti); 00394 00395 if(0>td) 00396 gGlobal->Set_ThetaMin(at-dt/2); 00397 else 00398 gGlobal->Set_ThetaMin(td); 00399 00400 if(0>de){ 00401 if(proc==0) 00402 gGlobal->Set_dE_Abs(0.015*EBeam); 00403 else 00404 gGlobal->Set_dE_Abs(0.0003*EBeam); 00405 } 00406 else 00407 gGlobal->Set_dE_Abs(de); 00408 00409 if(0>nt0) 00410 gGlobal->Set_Theta0_Rel(1.5); 00411 else 00412 gGlobal->Set_Theta0_Rel(nt0); 00413 00414 if(0<thm_min){ 00415 gCut->ThetaRangeM(thm_min, thm_max); 00416 gGlobal->Set_ThetaMin((thm_min<M_PI-thm_max)?thm_min:M_PI-thm_max); 00417 cout<<"Miminal angle = "<<((thm_min<M_PI-thm_max)?thm_min:M_PI-thm_max)<<endl; 00418 // gGlobal->Set_ThetaMin(1.1); 00419 } 00420 00421 if(0<thp_min){ 00422 gCut->ThetaRangeP(thp_min, thp_max); 00423 } 00424 gCut->CosPsi(cos(pc)); 00425 00426 gCut->DeltaTheta(dt); 00427 00428 gCut->AverageTheta(at); 00429 00430 gCut->DeltaPhi(dp); 00431 00432 gCut->PAverage(am); 00433 00434 gCut->PCross(cm); 00435 00436 gCut->EMin(em); 00437 00438 gConst->SetAlphaScale(al); 00439 00440 if(IsVerbose) gConst->Print(); 00441 00442 gGlobal->Init(); 00443 if(IsVerbose) gGlobal->Print(); 00444 00445 gCut->Init(); 00446 if(IsVerbose) gCut->Print(); 00447 00448 if(NSim>0) 00449 cout<<"NSim = "<<NSim<<endl; 00450 00451 if(fname.size()>0) 00452 cout<<"fname = "<<fname<<endl; 00453 00454 cout<<"Cross-section statistical precision will be better than " 00455 <<gGlobal->Get_TotalError()<<" nb and " 00456 <<gGlobal->Get_RelativeError()*100<<"%"<<endl; 00457 00458 if(IsBrem) 00459 cout<<"Bremsstrahlung on beam pipe is switched on!"<<endl; 00460 00461 if(IsSmear) 00462 cout<<"Angles will be smeared!"<<endl; 00463 00464 if(IsTuple && NSim>0) 00465 cout<<"Ntuple will be sampled!"<<endl; 00466 00467 if(!IsHardPhoton) 00468 cout<<"Hard photon on big angle is not included!"<<endl; 00469 00470 if(IsNoVacPol) 00471 cout<<"Vacuum polarization is not included!"<<endl; 00472 00473 if(!IsFSR) 00474 cout<<"Final state radiation is not included!"<<endl; 00475 00476 cout<<flush; 00477 00478 return proc; 00479 }
void hist_sample | ( | Int_t | i = 0 |
) |
Definition at line 481 of file hist.C.
References alpha, EBeam, f1, f2, fname, TRadCor::GetCrossPart(), TVCrossPart::GetEvent(), gRad, genRecEmupikp::i, IsBrem, IsSmear, IsTuple, ganga-rec::j, M_PI, TRadCor::MakeEvent(), mass, nbins, NSim, EvtCyclic3::prev(), proc_name, deljobs::string, timer, EVENT::Tphi, EVENT::Tphit, EVENT::Tptot, EVENT::Tptott, EVENT::Tth, EVENT::Ttht, EVENT::Type, type, and EVENT::Vz.
Referenced by main().
00481 { 00482 ostringstream ost; 00483 00484 // Read angular resolution data 00485 double theta_res[3],phi_res[3],d; 00486 if(IsSmear){ 00487 ost<<"angular"<<Int_t(EBeam+0.5)<<".dat"<<ends; 00488 cout<<"Try to open file with angular resolution."<<endl<<flush; 00489 cout<<ost.str()<<endl<<flush; 00490 ifstream IN(ost.str().c_str()); 00491 if(IN){ 00492 for(int i=0;i<3;i++)IN>>theta_res[i]>>d; 00493 for(int i=0;i<3;i++)IN>>phi_res[i]>>d; 00494 IN.close(); 00495 }else{ 00496 cout<<"Unable to open file. Assume that detector is ideal."<<endl; 00497 IsSmear = 0; 00498 } 00499 } 00500 // Open file. If it does not exist it will be created. 00501 TFile *f = NULL; 00502 string file_name; 00503 if(fname == ""){ 00504 ost.str(""); ost<<""; 00505 if(!IsTuple){ 00506 ost<<proc_name[iproc]<<Int_t(EBeam+0.5)<<"hist_new_a2.root"<<ends; 00507 }else{ 00508 ost<<proc_name[iproc]<<Int_t(EBeam+0.5)<<"tuple.root"<<ends; 00509 } 00510 file_name = ost.str(); 00511 } else { 00512 file_name = fname; 00513 } 00514 cout<<"Try to open output file."<<endl<<flush; 00515 cout<<file_name<<endl<<flush; 00516 f = new TFile(file_name.c_str(),"RECREATE"); 00517 cout<<"file opened."<<endl<<flush; 00518 00519 Double_t width = 0.5; 00520 Double_t pmin = 90.; 00521 Double_t pmax = TMath::Sqrt(EBeam*EBeam-mass[iproc]*mass[iproc]); 00522 pmax *= 1. + 1.e-5; // to make sure that all events will be inside histogram ranges 00523 Int_t nbins = int((pmax-pmin)/width); 00524 // to use this method bin width should be much less than momentum resolution 00525 00526 const Int_t ncth = 4, ndth = 5; 00527 TH2D *h[ncth][ndth]; 00528 TTree *tree=NULL; 00529 static EVENT event; 00530 if(!IsTuple){ 00531 // Create new histogram 00532 for(int i=0; i<ncth; i++){ 00533 for(int j=0; j<ndth; j++){ 00534 ost.str(""); ost<<proc_name[iproc]<<"_"<<i<<"_"<<j<<ends; 00535 string name = ost.str(); 00536 00537 ost.str(""); ost<<setw(3)<<setprecision(2)<<1.0+i*0.1; 00538 string hth = ost.str(); 00539 00540 ost.str(""); ost<<setw(4)<<setprecision(2)<<0.25-j*0.05; 00541 string dth = ost.str(); 00542 00543 ost.str(""); ost<<"Energy="<<Int_t(EBeam+0.5)<<" MeV, Theta_Cut=" 00544 <<hth<<", Delta_Theta="<<dth<<ends; 00545 00546 h[i][j] = new TH2D(name.c_str(), ost.str().c_str(),nbins,pmin,pmax,nbins,pmin,pmax); 00547 h[i][j]->SetDirectory(f); 00548 } 00549 } 00550 }else{ 00551 // Create new tree 00552 ost.str("");ost<<proc_name[iproc]<<"tuple"<<ends; 00553 tree = new TTree(ost.str().c_str(),"Simulation"); 00554 tree->Branch("event",&event,"Tptot[3]/F:Tth[3]/F:Tphi[3]/F:Tptott[3]/F:Ttht[3]/F:Tphit[3]/F:Vz:Type/I"); 00555 tree->SetDirectory(f); 00556 } 00557 00558 const Double_t alpha = 0.365073; // from 260 MeV point 00559 Double_t ar;// = TMath::Power(Ebeam/(Ebeam-pmin),alpha); 00560 const Double_t br = 105.76; // from 260 MeV point 00561 Double_t r; 00562 00563 Double_t p1,t1,f1,p2,t2,f2; 00564 Double_t bp1,bp2; 00565 Double_t smt1,smt2,smf1,smf2,th_a; 00566 Double_t prev, now, f_shift; 00567 UInt_t type; 00568 prev = timer->RealTime(); 00569 for(int i=0; i<NSim; i++){ 00570 00571 gRad->MakeEvent(p1,t1,f1,p2,t2,f2,type); 00572 p1 *= EBeam; 00573 p2 *= EBeam; 00574 00575 //simulation of bremsstrahlung photons parameterization from global simulation 00576 if(iproc == 0 && IsBrem ){ 00577 r = gRandom->Rndm(); 00578 ar = TMath::Power(p1/(p1-pmin),alpha); 00579 bp1 = p1*(1.-TMath::Power(1./(ar*r-br*(r-1.)),1/alpha)); 00580 00581 r = gRandom->Rndm(); 00582 ar = TMath::Power(p2/(p2-pmin),alpha); 00583 bp2 = p2*(1.-TMath::Power(1./(ar*r-br*(r-1.)),1/alpha)); 00584 00585 }else{ 00586 bp1 = p1; 00587 00588 bp2 = p2; 00589 00590 } 00591 00592 if(IsSmear){ 00593 //smearing of particle angles 00594 f_shift = 2*M_PI*gRandom->Rndm(); 00595 00596 f1 += f_shift; 00597 f2 += f_shift; 00598 00599 if(gRandom->Rndm()<theta_res[0]) 00600 smt1 = t1 + gRandom->Gaus(0.,theta_res[1]); 00601 else 00602 smt1 = t1 + gRandom->Gaus(0.,theta_res[2]); 00603 00604 if(gRandom->Rndm()<theta_res[0]) 00605 smt2 = t2 + gRandom->Gaus(0.,theta_res[1]); 00606 else 00607 smt2 = t2 + gRandom->Gaus(0.,theta_res[2]); 00608 00609 if(gRandom->Rndm()<phi_res[0]) 00610 smf1 = f1 + gRandom->Gaus(0.,phi_res[1]); 00611 else 00612 smf1 = f1 + gRandom->Gaus(0.,phi_res[2]); 00613 00614 if(gRandom->Rndm()<phi_res[0]) 00615 smf2 = f2 + gRandom->Gaus(0.,phi_res[1]); 00616 else 00617 smf2 = f2 + gRandom->Gaus(0.,phi_res[2]); 00618 00619 }else{ 00620 00621 smt1 = t1; 00622 smt2 = t2; 00623 smf1 = f1; 00624 smf2 = f2; 00625 00626 } 00627 00628 if(smf1<0) smf1 += 2*M_PI; 00629 if(smf1>2*M_PI) smf1 -= 2*M_PI; 00630 if(smf2<0) smf2 += 2*M_PI; 00631 if(smf2>2*M_PI) smf2 -= 2*M_PI; 00632 if(f1<0) f1 += 2*M_PI; 00633 if(f1>2*M_PI) f1 -= 2*M_PI; 00634 if(f2<0) f2 += 2*M_PI; 00635 if(f2>2*M_PI) f2 -= 2*M_PI; 00636 00637 if(IsTuple){ 00638 event.Tptot[0] = bp1; 00639 event.Tptot[1] = bp2; 00640 event.Tth[0] = smt1; 00641 event.Tth[1] = smt2; 00642 event.Tphi[0] = smf1; 00643 event.Tphi[1] = smf2; 00644 event.Tptott[0]= p1; 00645 event.Tptott[1]= p2; 00646 event.Ttht[0] = t1; 00647 event.Ttht[1] = t2; 00648 event.Tphit[0] = f1; 00649 event.Tphit[1] = f2; 00650 event.Type = type; 00651 00652 double ttt[4*4]; 00653 int npart; 00654 gRad->GetCrossPart()->GetEvent(ttt,npart); 00655 00656 int j=2; 00657 if(npart>3) 00658 if(ttt[2*4+3]<ttt[3*4+3]) j = 3; 00659 00660 event.Tptot[2] = 1000*ttt[j*4+3]; 00661 event.Tth[2] = acos(ttt[j*4+2]); 00662 event.Tphi[2] = atan2(ttt[j*4+1],ttt[j*4+0]); 00663 00664 } 00665 00666 th_a=0.5*(smt1+M_PI-smt2); 00667 Double_t delta_th = TMath::Abs(smt1+smt2-M_PI); 00668 if(th_a > 1.0 && th_a < M_PI-1.0 && 00669 delta_th<0.25 && 00670 TMath::Abs(TMath::Abs(smf1-smf2)-M_PI)<0.15){ 00671 if(!IsTuple){ 00672 for(int i=0; i<ncth; i++) 00673 for(int j=0; j<ndth; j++){ 00674 if(th_a > 1.+i*0.1 && th_a < M_PI - (1.+i*0.1) && 00675 delta_th<0.25-j*0.05 ) h[i][j]->Fill(bp1,bp2); 00676 } 00677 } 00678 if(!(i%(NSim/100))){ 00679 timer->Stop(); 00680 now = timer->RealTime(); 00681 if(now-prev>1){ 00682 cout<<"REALTIME = "<<now<<", " 00683 <<"CPUTIME = "<<timer->CpuTime()<<", " 00684 <<"EVENTNUM = "<<i<<", " 00685 <<"p1 = "<<p1<<", " 00686 <<"p2 = "<<p2 00687 <<endl 00688 <<flush; 00689 prev = now; 00690 } 00691 timer->Start(kFALSE); 00692 } 00693 if(IsTuple) event.Vz = 1.; 00694 }else{ 00695 if(IsTuple) event.Vz = 0.; 00696 } 00697 if(IsTuple) tree->Fill(); 00698 } 00699 cout<<"Histogram has been sampled!"<<timer->RealTime()<<endl<<flush; 00700 00701 if(!IsTuple){ 00702 for(int i=0; i<ncth; i++) 00703 for(int j=0; j<ndth; j++) 00704 h[i][j]->Write(); 00705 }else{ 00706 tree->Write(); 00707 } 00708 // Delete previous stored histograms 00709 f->Purge(); 00710 cout<<"Histogram has been written!"<<timer->RealTime()<<endl<<flush; 00711 if(!IsTuple) 00712 for(int i=0; i<ncth; i++) 00713 for(int j=0; j<ndth; j++) 00714 delete h[i][j]; 00715 if(IsTuple) delete tree; 00716 delete f; 00717 cout<<"File has been closed!"<<timer->RealTime()<<endl<<flush; 00718 }
void hist_sample_angle | ( | Int_t | i = 0 |
) |
Definition at line 720 of file hist.C.
References EBeam, f1, f2, fname, gRad, genRecEmupikp::i, ganga-rec::j, M_PI, TRadCor::MakeEvent(), nbins, NSim, EvtCyclic3::prev(), proc_name, deljobs::string, timer, and type.
Referenced by main().
00720 { 00721 ostringstream ost; 00722 00723 // Open file. If it does not exist it will be created. 00724 TFile *f = NULL; 00725 string file_name; 00726 if(fname == ""){ 00727 ost.str(""); ost<<""; 00728 ost<<proc_name[iproc]<<Int_t(EBeam+0.5)<<".root"<<ends; 00729 file_name = ost.str(); 00730 } else { 00731 file_name = fname; 00732 } 00733 cout<<"Try to open output file."<<endl<<flush; 00734 cout<<file_name<<endl<<flush; 00735 f = new TFile(file_name.c_str(),"RECREATE"); 00736 00737 const double MAX_DPHI = 0.3, MAX_DTH = 0.5; 00738 const int ndphi = 6, ndth = 10; 00739 TH1D *hdth[ndphi]; 00740 TH1D *hdphi[ndth]; 00741 00742 // Create new histogram 00743 for(int i=0; i<ndphi; i++){ 00744 ost.str(""); ost<<proc_name[iproc]<<"_dth_"<<i; 00745 string name = ost.str(); 00746 00747 ost.str(""); ost<<setw(3)<<setprecision(2)<<MAX_DPHI - i*0.05; 00748 string dphi = ost.str(); 00749 00750 ost.str(""); ost<<"#Delta#theta, Energy = "<<Int_t(EBeam+0.5)<<" MeV, " 00751 <<"#Delta#phi = "<<dphi; 00752 int nbins = 600; 00753 double dth_min = 0; 00754 double dth_max = MAX_DTH; 00755 hdth[i] = new TH1D(name.c_str(), ost.str().c_str(), nbins, dth_min, dth_max); 00756 hdth[i]->SetDirectory(f); 00757 } 00758 00759 for(int i=0; i<ndth; i++){ 00760 ost.str(""); ost<<proc_name[iproc]<<"_dphi_"<<i; 00761 string name = ost.str(); 00762 00763 ost.str(""); ost<<setw(3)<<setprecision(2)<<MAX_DTH - i*0.05; 00764 string dth = ost.str(); 00765 00766 ost.str(""); ost<<"#Delta#phi, Energy = "<<Int_t(EBeam+0.5)<<" MeV, " 00767 <<"#Delta#theta = "<<dth; 00768 int nbins = 600; 00769 double dphi_min = 0; 00770 double dphi_max = MAX_DPHI; 00771 hdphi[i] = new TH1D(name.c_str(), ost.str().c_str(), nbins, dphi_min, dphi_max); 00772 hdphi[i]->SetDirectory(f); 00773 } 00774 00775 Double_t p1,t1,f1,p2,t2,f2; 00776 UInt_t type; 00777 double prev = timer->RealTime(), now; 00778 for(int i=0; i<NSim; i++){ 00779 00780 gRad->MakeEvent(p1,t1,f1,p2,t2,f2,type); 00781 double a_delta_th = fabs(t1 + t2 - M_PI); 00782 double a_delta_phi = fabs(fabs(f1 - f2) - M_PI); 00783 00784 for(int j=0; j<ndphi; j++) 00785 if ( a_delta_phi < MAX_DPHI - j*0.05 ) hdth[j]->Fill(a_delta_th); 00786 00787 for(int j=0; j<ndth; j++) 00788 if ( a_delta_th < MAX_DTH - j*0.05 ) hdphi[j]->Fill(a_delta_phi); 00789 00790 if(!(i%(NSim/100))){ 00791 timer->Stop(); 00792 now = timer->RealTime(); 00793 if(now-prev>1){ 00794 cout<<"REALTIME = "<<now<<", " 00795 <<"CPUTIME = "<<timer->CpuTime()<<", " 00796 <<"EVENTNUM = "<<i<<", " 00797 <<"p1 = "<<p1<<", " 00798 <<"p2 = "<<p2 00799 <<endl 00800 <<flush; 00801 prev = now; 00802 } 00803 timer->Start(kFALSE); 00804 } 00805 } 00806 cout<<"Histogram has been sampled!"<<timer->RealTime()<<endl<<flush; 00807 00808 for(int i=0; i<ndphi; i++) 00809 hdth[i]->Write(); 00810 00811 for(int i=0; i<ndth; i++) 00812 hdphi[i]->Write(); 00813 00814 // Delete previous stored histograms 00815 f->Purge(); 00816 cout<<"Histogram has been written!"<<timer->RealTime()<<endl<<flush; 00817 00818 for(int i=0; i<ndphi; i++) 00819 delete hdth[i]; 00820 00821 for(int i=0; i<ndth; i++) 00822 delete hdphi[i]; 00823 00824 delete f; 00825 cout<<"File has been closed!"<<timer->RealTime()<<endl<<flush; 00826 }
void hist_sample_comp | ( | ) |
Definition at line 828 of file hist.C.
References EBeam, ebeam, f1, f2, fname, gConst, TRadGlobal::Get_E(), TRadGlobal::Get_ThetaMin(), gGlobal, gRad, genRecEmupikp::i, M_PI, TRadCor::MakeEvent(), me, TConstants::Me(), NSim, EvtCyclic3::prev(), timer, and type.
00828 { 00829 // Open file. If it does not exist it will be created. 00830 TFile *f = NULL; 00831 if(fname == "") return; 00832 cout<<"Try to open output file."<<endl<<flush; 00833 cout<<fname<<endl<<flush; 00834 f = new TFile(fname.c_str(),"RECREATE"); 00835 00836 double elemass = gConst->Me(); 00837 double ebeam = gGlobal->Get_E(); 00838 double emin = 0.8*ebeam; 00839 double thmin = gGlobal->Get_ThetaMin()*180/M_PI; 00840 double thmax = 180 - thmin; 00841 TH1D *an = new TH1D("an","Electron angle", 200, thmin - 1, thmax + 1); 00842 TH1D *en = new TH1D("en","Electron energy", 200, emin - 10, ebeam + 10); 00843 TH1D *ac = new TH1D("ac","Accolinearity", 200, -1, 11); 00844 TH1D *me = new TH1D("me","Missing energy", 200, -10, 0.4*ebeam); 00845 TH1D *ma = new TH1D("ma","Angle of missing energy",200,-1,181); 00846 00847 Double_t prev, now; 00848 prev = timer->RealTime(); 00849 for(int i=0; i<NSim; i++){ 00850 Double_t p1,t1,f1,p2,t2,f2; 00851 UInt_t type; 00852 gRad->MakeEvent(p1,t1,f1,p2,t2,f2,type); 00853 p1 *= EBeam; 00854 p2 *= EBeam; 00855 00856 an->Fill(t1*180/M_PI); 00857 double ele_en = sqrt(p1*p1 + elemass*elemass); 00858 double pos_en = sqrt(p2*p2 + elemass*elemass); 00859 if(ele_en>ebeam) ele_en = ebeam; 00860 if(pos_en>ebeam) pos_en = ebeam; 00861 en->Fill(ele_en-1e-12); 00862 ac->Fill(fabs(t1 + t2 - M_PI)*180/M_PI); 00863 double miss_en = 2*ebeam - ele_en - pos_en; 00864 if(miss_en>0.02*ebeam){ 00865 me->Fill(miss_en); 00866 double cth0,sth0,cphi0,sphi0,cth1,sth1,cphi1,sphi1; 00867 sincos(t1, &sth0, &cth0); 00868 sincos(t2, &sth1, &cth1); 00869 sincos(f1,&sphi0,&cphi0); 00870 sincos(f2,&sphi1,&cphi1); 00871 double px = p1*sth0*cphi0 + p2*sth1*cphi1; 00872 double py = p1*sth0*sphi0 + p2*sth1*sphi1; 00873 double pz = p1*cth0 + p2*cth1; 00874 ma->Fill(acos(pz/sqrt(px*px+py*py+pz*pz))*180/M_PI); 00875 } 00876 00877 if(!(i%(NSim/100))){ 00878 timer->Stop(); 00879 now = timer->RealTime(); 00880 if(now-prev>1){ 00881 cout<<"REALTIME = "<<now<<", " 00882 <<"CPUTIME = "<<timer->CpuTime()<<", " 00883 <<"EVENTNUM = "<<i<<", " 00884 <<"p1 = "<<p1<<", " 00885 <<"p2 = "<<p2 00886 <<endl 00887 <<flush; 00888 prev = now; 00889 } 00890 timer->Start(kFALSE); 00891 } 00892 } 00893 cout<<"Histogram has been sampled!"<<timer->RealTime()<<endl<<flush; 00894 f->Write(); 00895 delete f; 00896 cout<<"File has been closed!"<<timer->RealTime()<<endl<<flush; 00897 }
int main | ( | int | argc, | |
char ** | argv | |||
) |
Definition at line 87 of file hist.C.
References TKinemCut::AverageThetaMin(), TVCrossPart::BornCrossSection(), fname, gConst, gCut, TRadGlobal::Get_Norm(), TVCrossPart::GetEvent(), TRadCor::GetMax(), TRadCor::GetN(), TVCrossPart::GetNPart(), GetOption(), TRadCor::GetWN(), gGlobal, gH, gH2, gRad, hist_sample(), hist_sample_angle(), genRecEmupikp::i, TRadCor::Init(), InList, IsAngleHist, IsDTSmear, IsFSR, IsHardPhoton, IsNoVacPol, M_PI, TRadCor::MakeCrossSection(), MatrixElements, NRad, NSim, RandomFileName, TVCrossPart::SetHardPhoton(), TRadCor::SetNEvents(), TVCrossPart::SetNoFSR(), TRadCor::SetPartList(), TRadCor::SetSmear(), TVCrossPart::SetZeroVP(), smear_par, and timer.
00087 { 00088 int proc; 00089 timer = new TStopwatch(); 00090 timer->Start(kFALSE); 00091 00092 gConst = new TConstants(); 00093 gGlobal = new TRadGlobal(); 00094 gCut = new TKinemCut(); 00095 00096 try{ 00097 proc = GetOption(argc, argv); 00098 }catch(std::logic_error lge){ 00099 cout<<lge.what()<<endl; 00100 return 1; 00101 } 00102 00103 if(gRandom) delete gRandom; 00104 gRandom = new TRandom3(); 00105 struct stat stbuf; 00106 if( stat(RandomFileName.c_str(),&stbuf) == 0 ){ 00107 cout<<"Reading random seed from file \""<<RandomFileName<<"\"."<<endl; 00108 TFile *frndm = new TFile(RandomFileName.c_str(),"read"); 00109 frndm->ReadTObject(gRandom,gRandom->GetName()); 00110 delete frndm; 00111 // gRandom->ReadRandom(RandomFileName.c_str()); 00112 } else { 00113 cout<<strerror(errno)<<endl; 00114 } 00115 00116 TVCrossPart *MatrixElements; 00117 if(proc==0) 00118 MatrixElements = new TEPCrossPart(); 00119 else if(proc==1 || proc==5) 00120 MatrixElements = new TMuCrossPart(); 00121 else if(proc==2) 00122 MatrixElements = new TPiCrossPart(); 00123 else if(proc==3) 00124 MatrixElements = new TKnCrossPart(); 00125 else if(proc==4) 00126 MatrixElements = new TKcCrossPart(); 00127 else if(proc==6) 00128 MatrixElements = new TGGCrossPart(); 00129 else 00130 return 1; 00131 00132 if(IsNoVacPol) 00133 MatrixElements->SetZeroVP(); 00134 00135 if(!IsFSR) 00136 MatrixElements->SetNoFSR(); 00137 00138 gRad = new TRadCor(MatrixElements); 00139 gRad->SetNEvents(NRad); 00140 gRad->SetPartList(InList); 00141 gRad->Init(); 00142 00143 MatrixElements->SetHardPhoton(IsHardPhoton); 00144 if(IsDTSmear){ 00145 gRad->SetSmear(true); 00146 MatrixElements->GetEvent()->SetdPar(smear_par); 00147 } 00148 00149 #ifdef HIST 00150 TFile *f = NULL; 00151 if(fname == "") return -1; 00152 cout<<"Try to open output file."<<endl<<flush; 00153 cout<<fname<<endl<<flush; 00154 f = new TFile(fname.c_str(),"RECREATE"); 00155 00156 int nelem = MatrixElements->GetNPart(); 00157 for(int i=0;i<nelem+1;i++) 00158 if(InList[i]||i==nelem){ 00159 ostringstream hname; 00160 hname.str(""); hname<<"dth_"<<i; 00161 gH[i] = new TH1D(hname.str().c_str(),"#Delta#theta",34*5,0,34*0.03); 00162 gH[i]->Sumw2(); 00163 hname.str(""); hname<<"dth2_"<<i; 00164 gH2[i] = new TH1D(hname.str().c_str(),"#Delta#theta",34*5,0,34*0.03); 00165 gH2[i]->Sumw2(); 00166 } 00167 #endif 00168 gRad->MakeCrossSection(); 00169 #ifdef HIST 00170 for(int i=0;i<nelem;i++) 00171 if(InList[i]){ 00172 gH[nelem]->Add(gH[i],1./gRad->GetWN(i)); 00173 gH2[nelem]->Add(gH2[i],gRad->GetMax(i)/gRad->GetN(i)); 00174 } 00175 00176 cout<<"Histogram has been sampled!"<<timer->RealTime()<<endl<<flush; 00177 f->Write(); 00178 delete f; 00179 cout<<"File has been closed!"<<timer->RealTime()<<endl<<flush; 00180 #endif 00181 if(proc==2)((TPiCrossPart*)MatrixElements)->GetFormFactor()->Print(); 00182 cout<<"Born Cross Section is "<<2*M_PI*gGlobal->Get_Norm()*MatrixElements->BornCrossSection(gCut->AverageThetaMin())<<" nb"<<endl; 00183 if ( fname != "" && NSim > 0 ) { 00184 // hist_sample_comp(); 00185 if(IsAngleHist) 00186 hist_sample_angle(proc); 00187 else 00188 hist_sample(proc); 00189 } 00190 00191 delete MatrixElements; 00192 delete gRad; 00193 delete gCut; 00194 delete gGlobal; 00195 if( stat(RandomFileName.c_str(),&stbuf) == 0 ){ 00196 remove(RandomFileName.c_str()); 00197 } 00198 cout<<"Writing random seed to file \""<<RandomFileName<<"\"."<<endl; 00199 // gRandom->WriteRandom(RandomFileName.c_str()); 00200 TFile *frndm = new TFile(RandomFileName.c_str(),"recreate"); 00201 frndm->WriteTObject(gRandom,gRandom->GetName()); 00202 delete frndm; 00203 00204 delete gRandom; 00205 delete timer; 00206 00207 return 0; 00208 }
Definition at line 76 of file hist.C.
Referenced by MdcCalibFunSvc::checkConst(), DimRpcReader::currentFile(), MixerAlg::execute(), TofCalibManager::fillRoot(), TofCalibManager::fillTxt(), getFnamesAtBkk(), getHistList(), GetOption(), hist_sample(), hist_sample_angle(), hist_sample_comp(), TMuKFun::Init(), raw_ofstream::init_fstream(), MixerAlg::initialize(), MdcCosGeom::initPream(), EvtJetSet::jetSetInit(), main(), xmlBase::XmlParser::parse(), EvtPythia::pythiaInit(), raw_ofstream::real_fname(), ReaderRpc< Reader >::rpcHandler(), SqliteInterface::select_db(), and RawFileTools::wildcard_correct().
Definition at line 85 of file hist.C.
Referenced by main(), TRadCor::MakeIntegral(), and TRadCor::MakeMaximum().
Definition at line 67 of file hist.C.
Referenced by Mcgpj::execute(), Mcgpj::finalize(), hist_sample(), hist_sample_angle(), hist_sample_comp(), init_prime_gen_(), Mcgpj::initialize(), main(), and makeevent_().
Definition at line 84 of file hist.C.
Referenced by GetOption(), init_prime_gen_(), Mcgpj::initialize(), and main().
Int_t IsAngleHist |
Int_t IsBrem |
Int_t IsDTSmear |
Int_t IsFSR |
Int_t IsHardPhoton |
Int_t IsNoVacPol |
Int_t IsSmear |
Int_t IsTuple |
Int_t IsVerbose |
Double_t mass[] |
const int MaxList = 22 |
Definition at line 83 of file hist.C.
Referenced by GetOption(), init_prime_gen_(), and Mcgpj::initialize().
Int_t NRad |
Int_t NSim |
Definition at line 74 of file hist.C.
Referenced by GetOption(), hist_sample(), hist_sample_angle(), hist_sample_comp(), and main().
char* proc_name[] |
Initial value:
{ "ee", "mm", "pp", "kc", "kn", "tt", "gg" }
Definition at line 57 of file hist.C.
Referenced by hist_sample(), and hist_sample_angle().
Definition at line 76 of file hist.C.
Referenced by genfinish_(), GetOption(), init_prime_gen_(), and main().
double smear_par[9] |
TStopwatch* timer |
Definition at line 71 of file hist.C.
Referenced by VertexFit::BuildVirtualParticle(), VertexFit::Fit(), KinematicFit::Fit(), KinematicFit::fit(), KalmanKinematicFit::Fit(), KalmanKinematicFit::fit(), VertexFit::fitVertex(), hist_sample(), hist_sample_angle(), hist_sample_comp(), main(), TRadCor::MakeCrossSection(), VertexFit::swimVertex(), and VertexFit::vertexCovMatrix().
int type |
Definition at line 78 of file hist.C.
Referenced by MdcUtilitySvc::cellTrackPassed(), JobOptionsMgr::clientOptsTemplate(), RecMucTrackCnv::DataObjectToTObject(), HltInfCnv::DataObjectToTObject(), MucTrackCnv::DataObjectToTObject(), EvtTauola::decay(), EvtPythia::decay(), EvtPyGaGa::decay(), EvtPycont::decay(), EvtMultibody::decay(), EvtLundCharm::decay(), EvtLunda::decay(), EvtJscont::decay(), EvtJetSet::decay(), BesGeoTrack::Draw(), EvtMHelAmp::EvtMHelAmp(), CalibEventSelect::execute(), EFtoTDS::execute(), EFTest::execute(), HltEventMaker::execute(), TConformalFinder0::findClusters2(), G__RootEventData_rootcint_457_0_5(), G__RootEventData_rootcint_480_0_5(), G__RootEventData_rootcint_520_0_3(), KinematicFit::gda(), KalmanKinematicFit::gda(), GetOption(), hist_sample(), hist_sample_angle(), hist_sample_comp(), MucPID::init(), MdcID::is_axial(), KalFitAlg::kalman_fitting_anal(), KalFitAlg::kalman_fitting_calib(), KalFitAlg::kalman_fitting_csmalign(), KalFitAlg::kalman_fitting_MdcxReco_Csmc_Sew(), MdcxFittedHel::OriginIncluded(), EvtMTree::parsenode(), DedxCorrecSvc::PathL(), EvtPto3PAmpFactory::processAmp(), SqliteInterface::query(), MysqlInterface::query(), McCnvSvc::storageType(), RecMucTrackCnv::TObjectToDataObject(), MucTrackCnv::TObjectToDataObject(), VertexFit::UpdateConstraints(), KinematicFit::updateConstraints(), KalmanKinematicFit::updateConstraints(), and utf8_toUtf16().