// --------------------------------------------------------------------- // ----- KrattaEne source file ----- // ----- Created 10/11/13 by J. Lukasik ----- // --------------------------------------------------------------------- #ifndef __CINT__ #include #include "TBox.h" #include "TCanvas.h" #include "TH1F.h" #include "TH1D.h" #include "TH2F.h" #include "TNtuple.h" #include "TProfile.h" #include "TRandom.h" #include "TSystem.h" #include "TTree.h" #include "TVector3.h" #endif #include "TBox.h" #include "TCanvas.h" #include "TH1F.h" #include "TPaveLabel.h" #include "TROOT.h" #include "TMath.h" #include "TRandom3.h" #include "TGraph2D.h" #include "TH3D.h" #include "TSystemDirectory.h" #include "TSystemFile.h" #include "TList.h" #include "KRATTAene.hxx" //////////////////////////////////////////////////////////////////////// ClassImp(KrattaEne) //______________________________________________________________________ KrattaEne::KrattaEne() { /// Initialisation the static arrays: const float local_fct_vert[NMOD]={ 32.83, 27.35, 26.84, 26.73, 36.22, 14.09, 14.53, 33.00, 27.64, 18.55, 26.31, 21.63, 25.63, 16.68, 31.24, 24.80, 21.11, 33.45, 20.61, 22.08, 29.97, 26.72, 73.32, 32.60, 26.70, 18.13, 28.13, 23.35, 28.49, 28.52, 3.48, 30.01, 16.93, 29.64, 17.32}; const float local_thr_sct[KrattaEne::NMOD]={ 17.,16.,15.,16.,14.,27.,38., 16.,16., 0.,15.,16.,26.,36., 16.,15.,18.,13.,15.,31.,26., 16.,14.,18.,16.,15.,35.,26., 15.,14.,24.,21.,16.,26.,42.}; memcpy( fct_vert, local_fct_vert, NMOD*sizeof(float)); memcpy( thr_sct, local_thr_sct, NMOD*sizeof(float)); cpt_fct = 1.01; uball_thck = 7000; run_current =-1; mat_current ="INITIAL"; //mat_current="Au"; TString sAsyeosrootMainDir = gSystem->ExpandPathName("$VMCWORKDIR"); /// TODO: Here is local dependeed variable: FIX THIS: /// TString sKrattaEneDatDir = "calfiles/KRATTA/EnergyCalibr/KrattaEneDat"; TString sKrattaEneDatDir = "calfiles/KRATTA/Energy_calibration/KrattaEneDat"; TString sFullPathToKrattaEneDatDir = sAsyeosrootMainDir + "/" + sKrattaEneDatDir; strcpy ( dat_dir, sFullPathToKrattaEneDatDir.Data() ); strcpy ( TABS,"SRIM"); for(int imod=0;imodClear(); printf("[KrattaEne:] KrattaEne constructed\n"); } //______________________________________________________________________ KrattaEne::~KrattaEne(){} //______________________________________________________________________ void KrattaEne::LoadParameters(Int_t run){ //run_current = run; //printf("KrattaEne::LoadParameters(run=%d)\n",run); float pmin, pact, pmax; int runmin, runmax; char CH; string dummy; for(int imod=0;imodExpandPathName(dat_dir)); ifstream par_file(par_file_name); if(!par_file){ printf("cannot open file %s\n",par_file_name); exit(1); } char CMOD[4]; char cdum10[10]; char caz[2]; int jmod, kmod, kz, ka, iso, kf1, kf2; while(par_file >> CH){ if(CH=='#'){ getline(par_file,dummy); } else{ par_file.putback(CH); //printf("%s\n",dummy.data()); par_file >> CMOD >> jmod; //printf("MOD %02d\n",jmod); for(int ipar=0;ipar> cdum10 >> npar[jmod][ipar]; //printf("%2d %2d %10.6f\n",imod,ipar,par[imod][ipar]); } si0_thck[jmod] = npar[jmod][ 5]; // um of Si active layer for PD0 dl1_thck[jmod] = npar[jmod][ 3]; // um of Si dead layer for PD1 si1_thck[jmod] = npar[jmod][ 6]; // um of Si active layer for PD1 si2_thck[jmod] = npar[jmod][ 7]; // um of Si active layer for PD2 myl0_thck[jmod] = npar[jmod][ 8]; // um of Mylar betw PD1 CSI1 myl1_thck[jmod] = npar[jmod][ 9]; // um of Mylar betw CSI1 CSI2 myl2_thck[jmod] = npar[jmod][10]; // um of Mylar betw CSI2 PD2 csi1_thck[jmod] = npar[jmod][11]; // um CSI1 thickness csi2_thck[jmod] = npar[jmod][12]; // um CSI2 thickness fct_si0[jmod] = npar[jmod][13]; //channels/MeV PD0 fct_si1[jmod] = npar[jmod][14]; //channels/MeV PD1 fct_si2[jmod] = npar[jmod][15]; //channels/MeV PD2 for(int kiso=0;kiso<5;kiso++){ par_file >> CMOD >> kmod >> caz >> kz >> caz >> ka; //printf("Z %d A %d\n",kz,ka); iso = 4; if(kz==1 && ka==1) iso=0; else if(kz==1 && ka==2) iso=1; else if(kz==1 && ka==3) iso=2; else if(kz==2 && ka==3) iso=3; else if(kz==2 && ka==4) iso=4; par_file >> CMOD >> kf1 >> kf2; if(kf1 != 0 || kf2 != 0){ printf("[KrattaEne:] exit(1): check parameter file: %s\n",par_file_name); exit(1); } for(int ipar=0;ipar<10;ipar++){ par_file >> cdum10 >> nparaz[jmod][ipar][iso]; } na11[jmod][iso] = nparaz[jmod][2][iso]; //parameters for CSI1 na21[jmod][iso] = nparaz[jmod][3][iso]; //parameters for CSI1 nof1[jmod][iso] = nparaz[jmod][4][iso]; //parameters for CSI1 na12[jmod][iso] = nparaz[jmod][7][iso]; //parameters for CSI2 na22[jmod][iso] = nparaz[jmod][8][iso]; //parameters for CSI2 nof2[jmod][iso] = nparaz[jmod][9][iso]; //parameters for CSI2 } } } par_file.close(); } //______________________________________________________________________ void KrattaEne::LoadParametersOld(Int_t run){ //run_current = run; //printf("KrattaEne::LoadParametersOld(run=%d)\n",run); float pmin, pact, pmax; int runmin, runmax; for(int imod=0;imodExpandPathName(dat_dir),imod); ifstream par_file(par_file_name); if(!par_file){ for(int ipar=0;ipar> runmin >> runmax){ //printf("runmin,run,runmax %4d %4d %4d\n",runmin,run,runmax); if(run>=runmin && run<=runmax){ for(int ipar=0;ipar> pmin >> par[imod][ipar] >> pmax; //printf("%2d %2d %10.6f\n",imod,ipar,par[imod][ipar]); } par_file.close(); break; } else{ for(int ipar=0;ipar> pmin >> pact >> pmax; par[imod][ipar] = -1; } } } } si0_thck[imod] = par[imod][ 0]; // um of Si active layer for PD0 dl1_thck[imod] = par[imod][ 2]; // um of Si dead layer for PD1 si1_thck[imod] = par[imod][ 3]; // um of Si active layer for PD1 si2_thck[imod] = par[imod][18]; // um of Si active layer for PD2 myl0_thck[imod] = par[imod][ 5]; // um of Mylar betw PD1 CSI1 myl1_thck[imod] = par[imod][11]; // um of Mylar betw CSI1 CSI2 myl2_thck[imod] = par[imod][17]; // um of Mylar betw CSI2 PD2 csi1_thck[imod] = par[imod][ 6]*1000; // um CSI1 thickness csi2_thck[imod] = par[imod][12]*1000; // um CSI2 thickness fct_si0[imod] = par[imod][ 1]; //channels/MeV PD0 fct_si1[imod] = par[imod][ 4]; //channels/MeV PD1 fct_si2[imod] = par[imod][19]; //channels/MeV PD2 a11[imod] = par[imod][ 7]; //Parlog parameters for CSI1 a21[imod] = par[imod][ 8]; //Parlog parameters for CSI1 a31[imod] = par[imod][ 9]; //Parlog parameters for CSI1 a41[imod] = par[imod][10]; //Parlog parameters for CSI1 a12[imod] = par[imod][13]; //Parlog parameters for CSI2 a22[imod] = par[imod][14]; //Parlog parameters for CSI2 a32[imod] = par[imod][15]; //Parlog parameters for CSI2 a42[imod] = par[imod][16]; //Parlog parameters for CSI2 } } //______________________________________________________________________ void KrattaEne::LoadGammaLines(Int_t run){ //run_current = run; char CH; string dummy; for(int imod=0;imodExpandPathName(dat_dir)); ifstream par_file(par_file_name); if(!par_file){ // printf("cannot open file %s\n",par_file_name); //exit(1); } else{ while(par_file >> CH){ if(CH=='#'){ getline(par_file,dummy); } else{ par_file.putback(CH); break; } } while(par_file >> imod >> runmin >> runmax >> gp[0]>> gp[1]>> gp[2]>> gp[3]>> gp[4]>> gp[5]){ if(run>=runmin && run<=runmax){ for(int ipar=0;ipar<3;ipar++){ gamma1[imod][ipar] = gp[ipar]; gamma2[imod][ipar] = gp[ipar+3]; } } } par_file.close(); } // printf("mod run gamma1 | gamma2\n"); // for(imod=0;imodExpandPathName(dat_dir)); ifstream par_file(par_file_name); if(!par_file){ printf("[KrattaEne:] Cannot open file %s\n",par_file_name); //exit(1); } else{ while(par_file >> CH){ if(CH=='#'){ getline(par_file,dummy); } else{ par_file.putback(CH); break; } } while(par_file >> imod){ for(int ipar=0;ipar<5;ipar++){ par_file >> pexp[4-ipar] >> sigma; Pthr_CsI12[imod][0][4-ipar]=fabs(pexp[4-ipar]); Pthr_CsI12_sig[imod][0][4-ipar]=fabs(sigma); } par_file >> CCSI >> pcal[4] >> pcal[3] >> pcal[2] >> pcal[1] >> pcal[0]; for(int ipar=0;ipar<5;ipar++){ //CSI1 if(pexp[ipar]>0){ corr_fct[imod][0][ipar] = pcal[ipar]/pexp[ipar]; } else{ corr_fct[imod][0][ipar] = 1; } } for(int ipar=0;ipar<5;ipar++){ par_file >> pexp[4-ipar] >> sigma; Pthr_CsI12[imod][1][4-ipar]=fabs(pexp[4-ipar]); Pthr_CsI12_sig[imod][1][4-ipar]=fabs(sigma); } par_file >> CCSI >> pcal[4] >> pcal[3] >> pcal[2] >> pcal[1] >> pcal[0]; for(int ipar=0;ipar<5;ipar++){ //CSI1 if(pexp[ipar]>0){ corr_fct[imod][1][ipar] = pcal[ipar]/pexp[ipar]; } else{ corr_fct[imod][1][ipar] = 1; } } } par_file.close(); } float fct4he=0; float nfct4he=0; float fct3he=0; float nfct3he=0; for(imod=0;imodExpandPathName(dat_dir)); ifstream par_file(par_file_name); if(!par_file){ printf("[KrattaEne:] Cannot open file %s\n",par_file_name); //exit(1); } else{ while(par_file >> CH){ if(CH=='#'){ getline(par_file,dummy); } else{ par_file.putback(CH); break; } } while(par_file >> imod >> runmin >> runmax){ for(int ipar=0;ipar<5;ipar++){ par_file >> iz >> ia >> pexp0[ipar] >> sigma0[ipar] >> pexp1[ipar] >> sigma1[ipar]; } if(run >= runmin && run <= runmax){ for(int ipar=0;ipar<5;ipar++){ Pthr_PD01[imod][0][ipar]=pexp0[ipar]; Pthr_PD01_sig[imod][0][ipar]=sigma0[ipar]; Pthr_PD01[imod][1][ipar]=pexp1[ipar]; Pthr_PD01_sig[imod][1][ipar]=sigma1[ipar]; } } } par_file.close(); } // // // // fix missing corrections // for(imod=0;imodGetListOfSpecials(); TObject *obj; TIter nexth(th); TCutG *cg=0; while((obj = nexth())) { if(obj->InheritsFrom(TCutG::Class())){ cg = (TCutG*)obj; //printf("---> %s TCutG found\n",cg->GetName()); if(strstr(cg->GetName(),"CUTGSTAR") || strstr(cg->GetName(),"CUTGPTHR")){ //printf("---> %s TCutG removed\n",cg->GetName()); th->Remove((TCutG*)obj); delete obj; } } } for(int imod=0;imodExpandPathName(dat_dir); TSystemDirectory dir(dirname, dirname); TList *files = dir.GetListOfFiles(); if(files){ TSystemFile *file; TString fname; TIter next(files); while((file=(TSystemFile*)next())){ fname = file->GetName(); if(!file->IsDirectory() && fname.BeginsWith("cal_cutgstar") && fname.EndsWith("root")) { int icsi = atoi(fname(12,1).Data()); int imod = atoi(fname(15+icsi-1,2).Data()); int runmin = atoi(fname(18+icsi-1,4).Data()); int runmax = atoi(fname(23+icsi-1,4).Data()); //printf("%s | mod %d %d %4d %4d\n",fname.Data(),imod,icsi,runmin,runmax); if(run>=runmin && run<=runmax && imod>=0 && imodFindObject(Form("CUTGSTAR1_M%02d",imod)); cutgstar1[imod] = (TCutG*)tcg->Clone();//otherwise memory blowup if(!cutgstar1[imod]){ printf("[KrattaEne::LoadStarCuts:] %s NOT FOUND\n",Form("CUTGSTAR1_M%02d",imod)); } delete tcg; delete tc1; fcutgstar.Close(); } if(icsi==2){ fnames2[imod] = fname; TString full_fname = fname.Prepend("/"); full_fname.Prepend(dirname); TFile fcutgstar(full_fname,"read"); TCanvas *tc1 = (TCanvas*)fcutgstar.Get("c1"); //tc1->GetListOfPrimitives()->ls(); TCutG *tcg = 0; for(int ipar=0;ipar<4;ipar++){ tcg = (TCutG*)tc1->FindObject(Form("CUTGSTAR2Z%d_M%02d",ipar+1,imod)); cutgstar2[imod][ipar] = (TCutG*)tcg->Clone();//otherwise memory blowup } tcg = (TCutG*)tc1->FindObject(Form("CUTGSTAR2_M%02d",imod)); cutgstar2[imod][4] = (TCutG*)tcg->Clone();//otherwise memory blowup for(int ipar=0;ipar<5;ipar++){ if(!cutgstar2[imod][ipar]){ printf("[KrattaEne::LoadStarCuts:] %s NOT FOUND\n",Form("CUTGSTAR2Z%d_M%02d",ipar+1,imod)); } } delete tcg; delete tc1; fcutgstar.Close(); } } } if(!file->IsDirectory() && fname.BeginsWith("cal_cutgpthr") && fname.EndsWith("root")) { int icsi = atoi(fname(12,1).Data()); int imod = atoi(fname(15,2).Data()); int runmin = atoi(fname(18,4).Data()); int runmax = atoi(fname(23,4).Data()); //printf("%s | mod %d %d %4d %4d\n",fname.Data(),imod,icsi,runmin,runmax); if(run>=runmin && run<=runmax && imod>=0 && imodFindObject(Form("CUTGSTAR1_M%02d",imod)); // //fcutgstar.GetObject(Form("CUTGSTAR1_M%02d",imod),cutgstar1[imod]); // if(!cutgstar1[imod]){ // printf("KrattaEne::LoadStarCuts %s NOT FOUND\n",Form("CUTGSTAR1_M%02d",imod)); // } // fcutgstar.Close(); // } if(icsi==2){ fnames2p[imod] = fname; TString full_fname = fname.Prepend("/"); full_fname.Prepend(dirname); TFile fcutgpthr(full_fname,"read"); TCanvas *tc1 = (TCanvas*)fcutgpthr.Get("c1"); TCutG *tcg = (TCutG*)tc1->FindObject(Form("CUTGPTHR2_M%02d",imod)); cutgpthr2[imod] = (TCutG*)tcg->Clone();//otherwise memory blowup if(!cutgpthr2[imod]){ printf("[KrattaEne::LoadStarCuts:] %s NOT FOUND\n",Form("CUTGPTHR2_M%02d",imod)); } delete tcg; delete tc1; fcutgpthr.Close(); } } } } delete files; } for(int imod=0;imodPrint(); } //______________________________________________________________________ void KrattaEne::LoadExtraProtonCuts(Int_t run){ //run_current = run; //printf("KrattaEne::LoadExtraProtonCuts(run=%d)\n",run); TString fnames1[NMOD]; const char *missing=" --------- MISSING --------- "; TSeqCollection *th = gROOT->GetListOfSpecials(); TObject *obj; TIter nexth(th); TCutG *cg=0; while((obj = nexth())) { if(obj->InheritsFrom(TCutG::Class())){ cg = (TCutG*)obj; //printf("---> %s TCutG found\n",cg->GetName()); if(strcmp(cg->GetName(),"CUTGPEXTRA")==0){ //printf("---> %s TCutG removed\n",cg->GetName()); th->Remove((TCutG*)obj); delete obj; } } } for(int imod=0;imodExpandPathName(dat_dir); TSystemDirectory dir(dirname, dirname); TList *files = dir.GetListOfFiles(); if(files){ TSystemFile *file; TString fname; TIter next(files); while((file=(TSystemFile*)next())){ fname = file->GetName(); if(!file->IsDirectory() && fname.BeginsWith("cal_cutgpextra") && fname.EndsWith("root")) { int imod = atoi(fname(16,2).Data()); int runmin = atoi(fname(19,4).Data()); int runmax = atoi(fname(24,4).Data()); //printf("%s | mod %d %4d %4d\n",fname.Data(),imod,runmin,runmax); if(run>=runmin && run<=runmax && imod>=0 && imodFindObject("CUTGPEXTRA"); cutgpextra[imod] = (TCutG*)tcg->Clone();//otherwise memory blowup if(!cutgpextra[imod]){ printf("[KrattaEne::ExtraProtonCuts:] %s NOT FOUND\n",Form("CUTGPEXTRA_M%02d",imod)); } delete tcg; delete tc1; fcutgpextra.Close(); } } } delete files; } for(int imod=0;imodPrint(); } //______________________________________________________________________ void KrattaEne::LoadTime0Intervals(Int_t run){ //run_current = run; //printf("KrattaEne::LoadTime0Intervals(run=%d)\n",run); char CH; string dummy; int runmin, runmax, modc; for(int imod=0;imodExpandPathName(dat_dir),ipd); ifstream par_file(par_file_name); if(!par_file){ for(int imod=0;imod> CH){ if(CH=='#'){ getline(par_file,dummy); } else{ par_file.putback(CH); break; } } while(par_file >> runmin >> runmax){ //printf("runmin,run,runmax %4d %4d %4d\n",runmin,run,runmax); if(run>=runmin && run<=runmax){ for(int imod=0;imod> modc >> time0[imod][ipd][0] >> time0[imod][ipd][1] >> time0[imod][ipd][2] >> time0[imod][ipd][3]; } par_file.close(); break; } else{ for(int imod=0;imod> modc >> time0[imod][ipd][0] >> time0[imod][ipd][1] >> time0[imod][ipd][2] >> time0[imod][ipd][3]; time0[imod][ipd][0] = time0[imod][ipd][1] = time0[imod][ipd][2] = time0[imod][ipd][3] = -1; } } } } } // for(int ipd=0;ipd<3;ipd++){ // for(int imod=0;imod2||imod<0||imod>=NMOD) return false; if(nsigma_minus==0 && nsigma_plus==0){//within default range return (time>=time0[imod][ipd][2] && time<=time0[imod][ipd][3]) ? true : false; } return (time>=time0[imod][ipd][0]-nsigma_minus*time0[imod][ipd][1] && time<=time0[imod][ipd][0]+nsigma_plus*time0[imod][ipd][1]) ? true : false; } //______________________________________________________________________ Float_t KrattaEne::DiffTime(Int_t imod, Int_t ipd, Float_t time){ if(ipd<0||ipd>2||imod<0||imod>=NMOD) return -1000.f; return (time-time0[imod][ipd][0])/time0[imod][ipd][1]; } //______________________________________________________________________ ///void KrattaEne::Calibrate(Int_t run, ASYFadcPeak *peak, Int_t iz, Int_t ia, TRootKRATEvent *kevt, TRootKRATParticle *dat ) void KrattaEne::Calibrate(ASYFadcPeak *peak, TRootKRATEvent *kevt, TRootKRATParticle *dat ) { if (dat != NULL && kevt != NULL) { Calibrate(peak, *kevt, *dat ); } return; } //______________________________________________________________________ ///void KrattaEne::Calibrate(Int_t run, ASYFadcPeak *peak, Int_t iz, Int_t ia, TRootKRATEvent& kevt, TRootKRATParticle& dat ) void KrattaEne::Calibrate( ASYFadcPeak *peak, TRootKRATEvent& kevt, TRootKRATParticle& dat ) { /// cout << "[KrattaEne::Calibrate:] called " << endl; const Double_t u = 931.49386; const Double_t deg2rad = 0.0174533; const Double_t rad2deg = 57.2958; const Double_t pi = 3.1415926536; Int_t run = kevt.Run; Int_t iz = dat.Z; Int_t ia = dat.A; ///SK1213 dat.Run=run; // run No. ///TODO:NEW INTERFACE:Int_t iz = part->Z; ///TODO:NEW INTERFACE:Int_t ia = part->A; //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! correct in final version: if(run!=run_current){ //uncomment //if(run_current==-1){ //comment LoadParametersOld(run); LoadParameters(run); //!!!!! call after LoadParametersOld LoadFineCorrections(run); LoadPD01PthrChannels(run); LoadGammaLines(run); LoadStarCuts(run); LoadExtraProtonCuts(run); LoadTime0Intervals(run); //printf("parameters loaded\n"); //} //comment //if(run!=run_current){ //comment tgt_thickness = tgt->Thickness(run); if(mat_current.CompareTo(tgt->Material(run)) != 0){ if(re_tgt) { delete re_tgt; //printf("tgt deleted\n"); } // re_tgt = new RE_TAB(tgt->Material(run),TABS); //um re_tgt = new RE_TAB(tgt->Material(run),"ATIMA"); //um //printf("|%s| |%s|\n",mat_current.Data(),tgt->Material(run)); //printf("New target material: %s \n",tgt->Material(run)); mat_current = tgt->Material(run); } for(int imod=0;imodtht(imod))); //um } kevt.TgtThck=tgt_thickness; // target thickness sprintf(kevt.TgtMat,tgt->Material(run),"%s"); // target material if(run<=1844){ kevt.BeamZ = 79; kevt.BeamA = 197; sprintf(kevt.Beam,"Au"); // Beam } else if(run<=1974){ kevt.BeamZ = 44; kevt.BeamA = 96; sprintf(kevt.Beam,"Ru"); // Beam } else{ kevt.BeamZ = 40; kevt.BeamA = 96; sprintf(kevt.Beam,"Zr"); // Beam } float e0 = 400.f*kevt.BeamA; float rbeam = re_tgt->RANGE(kevt.BeamZ,kevt.BeamA,e0)-0.5*tgt_thickness; kevt.BeamEnergy = re_tgt->ENERGY(kevt.BeamZ,kevt.BeamA,rbeam); double mbeam = u*Double_t(kevt.BeamA); double mtrgt = u*Double_t(tgt->A(run)); double kbeam = kevt.BeamEnergy; double ebeam = kbeam+mbeam; double pbeam = sqrt(ebeam*ebeam-mbeam*mbeam); kevt.BeamRapidity = float(TMath::ATanH(pbeam/ebeam)); double b_cm = sqrt(kbeam*kbeam+2.*kbeam*mbeam)/(kbeam+mbeam+mtrgt); kevt.CMRapidity = float(TMath::ATanH(b_cm)); printf("[KrattaEne:] Current run %4d tgt %s %4.0f um | Beam %2s E %.3f AMeV y_B %.3f y_CM %.3f\n", run, tgt->Material(run),tgt_thickness, kevt.Beam,kevt.BeamEnergy/kevt.BeamA, kevt.BeamRapidity,kevt.CMRapidity); run_current=run; } //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! correct in final version Int_t MOD = peak->mod; int z = iz; float a = float(ia); dat.Module = MOD; dat.Theta = kgeo->tht(MOD); dat.Phi = kgeo->phi(MOD); kgeo->rndm(MOD); dat.ThetaRndm = kgeo->rndm_tht(); dat.PhiRndm = kgeo->rndm_phi(); /// We trust that the object is already empty /// because of time we will not to repeat any work. // dat.Clear(); /// or: // dat.Energy = -1.f; // dat.EnergyOld = -1.f; // dat.de_csi2 = -1; // dat.de_myl1 = -1; // dat.de_csi1 = -1; // dat.de_myl0 = -1; // dat.de_sil1 = -1; // dat.de_ddl1 = -1; // dat.de_sil0 = -1; // dat.de_ddl0 = -1; // dat.de_win = -1; // dat.de_air = -1; // dat.de_tgt = -1; // dat.WeightCol = -1.f; // dat.WeightInt = -1.f; // dat.WeightGeant = -1.f; // dat.RangeCsI = 0.f; // dat.DeltaECsI = 0.f; // dat.EnergyDE01 = -1.f; // dat.EnergyDE0 = -1.f; // dat.pthr_flag = 0; // dat.Pextra=false; // dat.Gamma=false; // dat.Gamma2[0]=false; // dat.Gamma2[1]=false; // dat.Star=false; // dat.Star2[0]=false; // dat.Star2[1]=false; // dat.Stopped[0]=//dat.Stopped[1]=//dat.Stopped[2]=//dat.Stopped[3]=false; // dat.PunchThrough=false; // dat.InTime=false; // dat.InTime3[0]=dat.InTime3[1]=dat.InTime3[2]=false; dat.ch_csi2 = //CsI2 total light peak->p22+peak->p23; dat.ch_csi1 = //CsI1 total light peak->p12+peak->p13-peak->p10/fct_vert[MOD]; if(dat.ch_csi1<0) dat.ch_csi1 = 0; dat.ch_sil1 = // Si1 amplitude peak->p10*(1.f+1.f/fct_vert[MOD]); if(dat.ch_sil1<0) dat.ch_sil1 = 0; dat.ch_sil0 = peak->p00; // Si0 amplitude dat.ch_slo2 = peak->p23; // Slow for CsI2 dat.ch_fas2 = peak->p22; // Fast for CsI2 dat.ch_sof2 = 0; // Ratio Slow/Fast for CsI2 if(peak->p22>0) dat.ch_sof2 = peak->p23/peak->p22; dat.ch_mod0 = 0; // Mode for Si0 if(peak->p09>0 && peak->p05>0) dat.ch_mod0 = (peak->p09!=peak->p05) ? peak->p05*peak->p09/(peak->p05-peak->p09)*log(peak->p05/peak->p09)*10 : peak->p09*10; dat.ch_sum1 = peak->p12+peak->p13; // CsI1 sum dat.ch_dif1 = peak->p12-peak->p13; // CsI1 difference dat.ch_sum2 = peak->p22+peak->p23; // CsI2 sum dat.ch_dif2 = peak->p22-peak->p23; // CsI2 difference dat.ch_log01 = 0; // log(Si0+Si1) if(peak->p10*(1.f+1.f/fct_vert[MOD])>0) dat.ch_log01 = log(peak->p10*(1.f+1.f/fct_vert[MOD])); dat.ch_slo1 = peak->p13; // Slow for CsI1 dat.ch_fas1 = peak->p12; // Fast for CsI1 dat.ch_sof1 = 0; // Ratio Slow/Fast for CsI1 if(peak->p12-peak->p10/fct_vert[MOD]>0) dat.ch_sof1 = peak->p13/(peak->p12-peak->p10/fct_vert[MOD]); dat.ch_log1 = 0; // log(dat.ch_csi1) if(dat.ch_csi1>0) dat.ch_log1 = log(dat.ch_csi1); dat.ch_log2 = 0; // log(dat.ch_csi1) if(dat.ch_csi2>0) dat.ch_log2 = log(dat.ch_csi2); dat.ch_tm0 = peak->p04; // time0 PD0 dat.ch_tm1 = peak->p14; // time0 PD1 dat.ch_tm2 = peak->p24; // time0 PD2 // check if below gamma-line if(peak->am1>0&&dat.ch_dif1am2>0&&dat.ch_dif2IsInside(dat.ch_sof1,dat.ch_log01)){ dat.Star=true; dat.Star2[0]=true; } int jz = (iz<4) ? iz : 4; if(jz>0&&cutgstar2[MOD][jz-1]) if(cutgstar2[MOD][jz-1]->IsInside(dat.ch_sof2,dat.ch_log1)){ dat.Star=true; dat.Star2[1]=true; } //check if inside "pthr" cut (<==> NOT punching through CsI2) if(cutgpthr2[MOD]) if(!cutgpthr2[MOD]->IsInside(dat.ch_sof2,dat.ch_log2)) dat.PunchThrough=true; //if(!dat.Star2[1]) dat.PunchThrough=true;//Punch Through CSI2 or background //if(!dat.Gamma2[0]&&dat.Star2[0]) printf("dat.Star2[0] %7d %7.1f %7.1f %7.1f %7.1f\n",evt->evt,dat.ch_sil0,dat.ch_sil1,dat.ch_csi1,dat.ch_csi2); //if(!dat.Gamma2[1]&&dat.Star2[1]) printf("dat.Star2[1] %7d %7.1f %7.1f %7.1f %7.1f\n",evt->evt,dat.ch_sil0,dat.ch_sil1,dat.ch_csi1,dat.ch_csi2); if(peak->am2>0&&peak->am1>0&&!dat.PunchThrough) dat.Stopped[3]=true; //stopped in CSI2 if(dat.ch_csi1>=thr_sct[MOD]&&peak->am2<=0) dat.Stopped[2]=true; //stopped in CSI1 if(peak->am1>0&&dat.ch_csi1am2<=0){ dat.Stopped[1]=true; //stopped in PD1 dat.Star=true; } if(peak->am0>0&&peak->am1<=0&&peak->am2<=0){ dat.Stopped[0]=true; //stopped in PD0 dat.Star=true; } if(dat.Stopped[2] //in CSI1 //&& !dat.Gamma2[0] && !dat.Star2[0] && cutgpextra[MOD]->IsInside(dat.ch_csi1,dat.ch_sil1) && iz==1&&ia==1){ dat.Pextra=true; } // float thr2 = (MOD%7 > 4) ? 200.f : 100.f; // //if(peak->am2p12+peak->p13>10*thr2 // && !dat.Star2[0] // && !dat.Gamma2[0] // && ia<=4 // && ia>1 // ){ // dat.Pextra=true; // } float a11CORR = a11[MOD]; if (iz==1&&ia>0&&ia<4) a11CORR *= corr_fct[MOD][0][ia-1]; else if(iz==2&&ia==3) a11CORR *= corr_fct[MOD][0][3]; else a11CORR *= corr_fct[MOD][0][4]; float a12CORR = a12[MOD]; if (iz==1&&ia>0&&ia<4) a12CORR *= corr_fct[MOD][1][ia-1]; else if(iz==2&&ia==3) a12CORR *= corr_fct[MOD][1][3]; else a12CORR *= corr_fct[MOD][1][4]; // if (iz==1&&ia==1){a11CORR *= 0.95;a12CORR *= 0.96;} // if (iz==1&&ia==2){a11CORR *= 0.98;a12CORR *= 0.97;} // if (iz==1&&ia==3){a11CORR *= 0.99;a12CORR *= 0.98;} // if (iz==2&&ia==3){a11CORR *= 0.99;a12CORR *= 0.98;} // calculate residual energy from DeltaE in PD0+PD1 if(dat.ch_csi1>0&&(dat.Stopped[2]||dat.Stopped[3])){ float de_pthr = SI01_PTHR(dat.ch_sil0, dat.ch_sil1,iz, ia, fct_si0[MOD], fct_si1[MOD], si0_thck[MOD],dl1_thck[MOD],si1_thck[MOD],re_si); dat.EneResDE01 = de_pthr; float rng_sil1 = re_si->RANGE(z,a,de_pthr)+si1_thck[MOD]; float ene_sil1 = re_si->ENERGY(z,a,rng_sil1); float rng_ddl1 = re_si->RANGE(z,a,ene_sil1)+dl1_thck[MOD]; float ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); float rng_sil0 = re_si->RANGE(z,a,ene_ddl1)+si0_thck[MOD]; float ene_sil0 = re_si->ENERGY(z,a,rng_sil0); float rng_ddl0 = re_si->RANGE(z,a,ene_sil0)+dl0_thck; float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); float rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; float ene_win = re_cu->ENERGY(z,a,rng_win); float rng_air = re_air->RANGE(z,a,ene_win)+air_thck; float ene_air = re_air->ENERGY(z,a,rng_air); float rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; float ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); dat.EnergyDE01 = ene_tgt; } int iso=4; if(iz==1 && ia==1) iso=0; else if(iz==1 && ia==2) iso=1; else if(iz==1 && ia==3) iso=2; else if(iz==2 && ia==3) iso=3; // calculate residual energy from DeltaE in PD0 if(dat.ch_sil0>0&&(dat.Stopped[1]||dat.Stopped[2]||dat.Stopped[3])){ float de_pthr_si0 = SI0_PTHR(dat.ch_sil0, iz, ia, fct_si0[MOD], si0_thck[MOD],re_si); dat.EneResDE0 = de_pthr_si0; float rng_sil0 = re_si->RANGE(z,a,de_pthr_si0)+si0_thck[MOD]; float ene_sil0 = re_si->ENERGY(z,a,rng_sil0); float rng_ddl0 = re_si->RANGE(z,a,ene_sil0)+dl0_thck; float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); float rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; float ene_win = re_cu->ENERGY(z,a,rng_win); float rng_air = re_air->RANGE(z,a,ene_win)+air_thck; float ene_air = re_air->ENERGY(z,a,rng_air); float rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; float ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); dat.EnergyDE0 = ene_tgt; } if(dat.Pextra){ //dat.de_csi2 = H2E_PTHR(dat.ch_csi1, iz, ia, a11CORR,a21[MOD],a31[MOD],a41[MOD],fct_si1[MOD],csi1_thck[MOD], re_csi); //NEW calibration if(iso<3) dat.de_csi2 = H2E_PTHRI(dat.ch_csi1, iz, ia, na11[MOD][iso],na21[MOD][iso],fct_si1[MOD],csi1_thck[MOD], re_csi,nof1[MOD][iso]); else dat.de_csi2 = H2E_PTHR (dat.ch_csi1, iz, ia, na11[MOD][iso],na21[MOD][iso],fct_si1[MOD],csi1_thck[MOD], re_csi,nof1[MOD][iso]); // float rng_myl1 = re_myl->RANGE(z,a,dat.de_csi2)+myl1_thck[MOD]; // float ene_myl1 = re_myl->ENERGY(z,a,rng_myl1); // dat.de_myl1 = ene_myl1-dat.de_csi2; // float rng_csi1 = re_csi->RANGE(z,a,ene_myl1)+csi1_thck[MOD]; float rng_csi1 = re_csi->RANGE(z,a,dat.de_csi2)+csi1_thck[MOD]; float ene_csi1 = re_csi->ENERGY(z,a,rng_csi1); dat.de_csi1 = ene_csi1-dat.de_csi2; float rng_myl0 = re_myl->RANGE(z,a,ene_csi1)+myl0_thck[MOD]; float ene_myl0 = re_myl->ENERGY(z,a,rng_myl0); dat.de_myl0 = ene_myl0-ene_csi1; dat.de_sil1 = dat.ch_sil1/fct_si1[MOD]; float rng_ddl1 = re_si->RANGE(z,a,dat.de_sil1+ene_myl0)+dl1_thck[MOD]; float ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); dat.de_ddl1 = ene_ddl1-(dat.de_sil1+ene_myl0); dat.de_sil0 = dat.ch_sil0/fct_si0[MOD]; float rng_ddl0 = re_si->RANGE(z,a,dat.de_sil0+ene_ddl1)+dl0_thck; float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); dat.de_ddl0 = ene_ddl0-(dat.de_sil0+ene_ddl1); float rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; float ene_win = re_cu->ENERGY(z,a,rng_win); dat.de_win = ene_win-ene_ddl0; float rng_air = re_air->RANGE(z,a,ene_win)+air_thck; float ene_air = re_air->ENERGY(z,a,rng_air); dat.de_air = ene_air-ene_win; float rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; float ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); dat.de_tgt = ene_tgt-ene_air; dat.EnergyOld = ene_tgt; dat.RangeCsI = rng_csi1/10000.; //[cm] dat.DeltaECsI = re_csi->ENERGY(z,a,rng_csi1); //printf("Pextra: E %d %d %f Ecsi %f Rcsi %f [mm]\n",iz,ia,dat.EnergyOld,dat.DeltaECsI,dat.RangeCsI*10); } else if(dat.Stopped[0]){ // in PD0 if(iso<4 && dat.ch_sil0>cpt_fct*Pthr_PD01[MOD][0][iso]){ dat.pthr_flag = 1; } //printf("PD0 %d %f\n",iso,Pthr_PD01[MOD][1][iso]); dat.de_sil0 = dat.ch_sil0/fct_si0[MOD]; float rng_ddl0 = re_si->RANGE(z,a,dat.de_sil0)+dl0_thck; float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); dat.de_ddl0 = ene_ddl0-dat.de_sil0; float rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; float ene_win = re_cu->ENERGY(z,a,rng_win); dat.de_win = ene_win-ene_ddl0; float rng_air = re_air->RANGE(z,a,ene_win)+air_thck; float ene_air = re_air->ENERGY(z,a,rng_air); dat.de_air = ene_air-ene_win; if(MOD%7==6){//column shadowed by uBall float rng_uball = re_csi->RANGE(z,a,ene_air)+uball_thck; float ene_uball = re_csi->ENERGY(z,a,rng_uball); float de_uball = ene_uball-ene_air; dat.de_air += de_uball; ene_air = ene_uball; } float rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; float ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); dat.de_tgt = ene_tgt-ene_air; dat.EnergyOld = ene_tgt; dat.Energy = ene_tgt; } else if(dat.Stopped[1]){ // in PD1 if(iso<4 && dat.ch_sil1>cpt_fct*Pthr_PD01[MOD][1][iso]){ dat.pthr_flag = 1; } //printf("PD1 %d %f\n",iso,Pthr_PD01[MOD][1][iso]); dat.de_sil1 = dat.ch_sil1/fct_si1[MOD]; float rng_ddl1 = re_si->RANGE(z,a,dat.de_sil1)+dl1_thck[MOD]; float ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); dat.de_ddl1 = ene_ddl1-dat.de_sil1; dat.de_sil0 = dat.ch_sil0/fct_si0[MOD]; float rng_ddl0 = re_si->RANGE(z,a,dat.de_sil0+ene_ddl1)+dl0_thck; float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); dat.de_ddl0 = ene_ddl0-(dat.de_sil0+ene_ddl1); float rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; float ene_win = re_cu->ENERGY(z,a,rng_win); dat.de_win = ene_win-ene_ddl0; float rng_air = re_air->RANGE(z,a,ene_win)+air_thck; float ene_air = re_air->ENERGY(z,a,rng_air); dat.de_air = ene_air-ene_win; if(MOD%7==6){//column shadowed by uBall float rng_uball = re_csi->RANGE(z,a,ene_air)+uball_thck; float ene_uball = re_csi->ENERGY(z,a,rng_uball); float de_uball = ene_uball-ene_air; dat.de_air += de_uball; ene_air = ene_uball; } float rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; float ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); dat.de_tgt = ene_tgt-ene_air; dat.EnergyOld = ene_tgt; dat.Energy = ene_tgt; } else if(dat.Stopped[2]){ // in CSI1 if(iso<4 && dat.ch_csi1>cpt_fct*Pthr_CsI12[MOD][0][iso]){ dat.pthr_flag = 1; } dat.de_csi1 = H2E(dat.ch_csi1, iz, ia, a11CORR,a21[MOD],a31[MOD],a41[MOD],fct_si1[MOD]); float rng_myl0 = re_myl->RANGE(z,a,dat.de_csi1)+myl0_thck[MOD]; float ene_myl0 = re_myl->ENERGY(z,a,rng_myl0); dat.de_myl0 = ene_myl0-dat.de_csi1; // //from amplitudes: // dat.de_sil1 = dat.ch_sil1/fct_si1[MOD]; // float rng_ddl1 = re_si->RANGE(z,a,dat.de_sil1+ene_myl0)+dl1_thck[MOD]; // float ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); // dat.de_ddl1 = ene_ddl1-(dat.de_sil1+ene_myl0); // dat.de_sil0 = dat.ch_sil0/fct_si0[MOD]; // float rng_ddl0 = re_si->RANGE(z,a,dat.de_sil0+ene_ddl1)+dl0_thck; // float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); // dat.de_ddl0 = ene_ddl0-(dat.de_sil0+ene_ddl1); //from tables: float rng_sil1 = re_si->RANGE(z,a,ene_myl0)+si1_thck[MOD]; float ene_sil1 = re_si->ENERGY(z,a,rng_sil1); dat.de_sil1 = ene_sil1-ene_myl0; float rng_ddl1 = re_si->RANGE(z,a,ene_sil1)+dl1_thck[MOD]; float ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); dat.de_ddl1 = ene_ddl1-ene_sil1; float rng_sil0 = re_si->RANGE(z,a,ene_ddl1)+si0_thck[MOD]; float ene_sil0 = re_si->ENERGY(z,a,rng_sil0); dat.de_sil0 = ene_sil0-ene_ddl1; float rng_ddl0 = re_si->RANGE(z,a,ene_sil0)+dl0_thck; float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); dat.de_ddl0 = ene_ddl0-ene_sil0; float rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; float ene_win = re_cu->ENERGY(z,a,rng_win); dat.de_win = ene_win-ene_ddl0; float rng_air = re_air->RANGE(z,a,ene_win)+air_thck; float ene_air = re_air->ENERGY(z,a,rng_air); dat.de_air = ene_air-ene_win; if(MOD%7==6){//column shadowed by uBall float rng_uball = re_csi->RANGE(z,a,ene_air)+uball_thck; float ene_uball = re_csi->ENERGY(z,a,rng_uball); float de_uball = ene_uball-ene_air; dat.de_air += de_uball; ene_air = ene_uball; } float rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; float ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); dat.de_tgt = ene_tgt-ene_air; dat.EnergyOld = ene_tgt; //NEW calibration if(iso<3) dat.de_csi1 = NH2EI(dat.ch_csi1, iz, ia, na11[MOD][iso],na21[MOD][iso],nof1[MOD][iso],fct_si1[MOD]); else dat.de_csi1 = NH2E(dat.ch_csi1, iz, ia, na11[MOD][iso],na21[MOD][iso],nof1[MOD][iso],fct_si1[MOD]); rng_myl0 = re_myl->RANGE(z,a,dat.de_csi1)+myl0_thck[MOD]; ene_myl0 = re_myl->ENERGY(z,a,rng_myl0); dat.de_myl0 = ene_myl0-dat.de_csi1; // //from amplitudes: // dat.de_sil1 = dat.ch_sil1/fct_si1[MOD]; // rng_ddl1 = re_si->RANGE(z,a,dat.de_sil1+ene_myl0)+dl1_thck[MOD]; // ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); // dat.de_ddl1 = ene_ddl1-(dat.de_sil1+ene_myl0); // dat.de_sil0 = dat.ch_sil0/fct_si0[MOD]; // rng_ddl0 = re_si->RANGE(z,a,dat.de_sil0+ene_ddl1)+dl0_thck; // ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); // dat.de_ddl0 = ene_ddl0-(dat.de_sil0+ene_ddl1); //from tables: rng_sil1 = re_si->RANGE(z,a,ene_myl0)+si1_thck[MOD]; ene_sil1 = re_si->ENERGY(z,a,rng_sil1); dat.de_sil1 = ene_sil1-ene_myl0; rng_ddl1 = re_si->RANGE(z,a,ene_sil1)+dl1_thck[MOD]; ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); dat.de_ddl1 = ene_ddl1-ene_sil1; rng_sil0 = re_si->RANGE(z,a,ene_ddl1)+si0_thck[MOD]; ene_sil0 = re_si->ENERGY(z,a,rng_sil0); dat.de_sil0 = ene_sil0-ene_ddl1; rng_ddl0 = re_si->RANGE(z,a,ene_sil0)+dl0_thck; ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); dat.de_ddl0 = ene_ddl0-ene_sil0; rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; ene_win = re_cu->ENERGY(z,a,rng_win); dat.de_win = ene_win-ene_ddl0; rng_air = re_air->RANGE(z,a,ene_win)+air_thck; ene_air = re_air->ENERGY(z,a,rng_air); dat.de_air = ene_air-ene_win; if(MOD%7==6){//column shadowed by uBall float rng_uball = re_csi->RANGE(z,a,ene_air)+uball_thck; float ene_uball = re_csi->ENERGY(z,a,rng_uball); float de_uball = ene_uball-ene_air; dat.de_air += de_uball; ene_air = ene_uball; } rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); dat.de_tgt = ene_tgt-ene_air; dat.Energy = ene_tgt; dat.RangeCsI = re_csi->RANGE(z,a,dat.de_csi1)/10000.; //[cm] dat.DeltaECsI = dat.de_csi1; } else if(dat.Stopped[3]){ // in CSI2 if(iso<4 && dat.ch_csi2>cpt_fct*Pthr_CsI12[MOD][1][iso]){ dat.pthr_flag = 1; } dat.de_csi2 = H2E(dat.ch_csi2, iz, ia, a12CORR,a22[MOD],a32[MOD],a42[MOD],fct_si2[MOD]); float rng_myl1 = re_myl->RANGE(z,a,dat.de_csi2)+myl1_thck[MOD]; float ene_myl1 = re_myl->ENERGY(z,a,rng_myl1); dat.de_myl1 = ene_myl1-dat.de_csi2; float rng_csi1 = re_csi->RANGE(z,a,ene_myl1)+csi1_thck[MOD]; float ene_csi1 = re_csi->ENERGY(z,a,rng_csi1); dat.de_csi1 = ene_csi1-ene_myl1; float rng_myl0 = re_myl->RANGE(z,a,ene_csi1)+myl0_thck[MOD]; float ene_myl0 = re_myl->ENERGY(z,a,rng_myl0); dat.de_myl0 = ene_myl0-ene_csi1; float rng_sil1 = re_si->RANGE(z,a,ene_myl0)+si1_thck[MOD]; float ene_sil1 = re_si->ENERGY(z,a,rng_sil1); dat.de_sil1 = ene_sil1-ene_myl0; float rng_ddl1 = re_si->RANGE(z,a,ene_sil1)+dl1_thck[MOD]; float ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); dat.de_ddl1 = ene_ddl1-ene_sil1; float rng_sil0 = re_si->RANGE(z,a,ene_ddl1)+si0_thck[MOD]; float ene_sil0 = re_si->ENERGY(z,a,rng_sil0); dat.de_sil0 = ene_sil0-ene_ddl1; float rng_ddl0 = re_si->RANGE(z,a,ene_sil0)+dl0_thck; float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); dat.de_ddl0 = ene_ddl0-ene_sil0; // dat.de_sil1 = dat.ch_sil1/fct_si1[MOD]; // float rng_ddl1 = re_si->RANGE(z,a,dat.de_sil1+ene_myl0)+dl1_thck[MOD]; // float ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); // dat.de_ddl1 = ene_ddl1-(dat.de_sil1+ene_myl0); // dat.de_sil0 = dat.ch_sil0/fct_si0[MOD]; // float rng_ddl0 = re_si->RANGE(z,a,dat.de_sil0+ene_ddl1)+dl0_thck; // float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); // dat.de_ddl0 = ene_ddl0-(dat.de_sil0+ene_ddl1); float rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; float ene_win = re_cu->ENERGY(z,a,rng_win); dat.de_win = ene_win-ene_ddl0; float rng_air = re_air->RANGE(z,a,ene_win)+air_thck; float ene_air = re_air->ENERGY(z,a,rng_air); dat.de_air = ene_air-ene_win; if(MOD%7==6){//column shadowed by uBall float rng_uball = re_csi->RANGE(z,a,ene_air)+uball_thck; float ene_uball = re_csi->ENERGY(z,a,rng_uball); float de_uball = ene_uball-ene_air; dat.de_air += de_uball; ene_air = ene_uball; } float rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; float ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); dat.de_tgt = ene_tgt-ene_air; dat.EnergyOld = ene_tgt; // float enemyl1 = re_myl->ENERGY(z,a,myl1_thck[MOD]); // float rngcsi1 = re_csi->RANGE(z,a,enemyl1)+csi1_thck[MOD]; // float de_ptr1= H2E_PTHR(dat.ch_csi1,iz,ia,a11CORR,a21[MOD],a31[MOD],a41[MOD],fct_si1[MOD],rngcsi1, re_csi); //NEW calibration if(iso<3) dat.de_csi2 = NH2EI(dat.ch_csi2, iz, ia, na12[MOD][iso],na22[MOD][iso],nof2[MOD][iso],fct_si2[MOD]); else dat.de_csi2 = NH2E(dat.ch_csi2, iz, ia, na12[MOD][iso],na22[MOD][iso],nof2[MOD][iso],fct_si2[MOD]); rng_myl1 = re_myl->RANGE(z,a,dat.de_csi2)+myl1_thck[MOD]; ene_myl1 = re_myl->ENERGY(z,a,rng_myl1); dat.de_myl1 = ene_myl1-dat.de_csi2; rng_csi1 = re_csi->RANGE(z,a,ene_myl1)+csi1_thck[MOD]; ene_csi1 = re_csi->ENERGY(z,a,rng_csi1); dat.de_csi1 = ene_csi1-ene_myl1; rng_myl0 = re_myl->RANGE(z,a,ene_csi1)+myl0_thck[MOD]; ene_myl0 = re_myl->ENERGY(z,a,rng_myl0); dat.de_myl0 = ene_myl0-ene_csi1; rng_sil1 = re_si->RANGE(z,a,ene_myl0)+si1_thck[MOD]; ene_sil1 = re_si->ENERGY(z,a,rng_sil1); dat.de_sil1 = ene_sil1-ene_myl0; rng_ddl1 = re_si->RANGE(z,a,ene_sil1)+dl1_thck[MOD]; ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); dat.de_ddl1 = ene_ddl1-ene_sil1; rng_sil0 = re_si->RANGE(z,a,ene_ddl1)+si0_thck[MOD]; ene_sil0 = re_si->ENERGY(z,a,rng_sil0); dat.de_sil0 = ene_sil0-ene_ddl1; rng_ddl0 = re_si->RANGE(z,a,ene_sil0)+dl0_thck; ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); dat.de_ddl0 = ene_ddl0-ene_sil0; // dat.de_sil1 = dat.ch_sil1/fct_si1[MOD]; // rng_ddl1 = re_si->RANGE(z,a,dat.de_sil1+ene_myl0)+dl1_thck[MOD]; // ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); // dat.de_ddl1 = ene_ddl1-(dat.de_sil1+ene_myl0); // dat.de_sil0 = dat.ch_sil0/fct_si0[MOD]; // rng_ddl0 = re_si->RANGE(z,a,dat.de_sil0+ene_ddl1)+dl0_thck; // ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); // dat.de_ddl0 = ene_ddl0-(dat.de_sil0+ene_ddl1); rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; ene_win = re_cu->ENERGY(z,a,rng_win); dat.de_win = ene_win-ene_ddl0; rng_air = re_air->RANGE(z,a,ene_win)+air_thck; ene_air = re_air->ENERGY(z,a,rng_air); dat.de_air = ene_air-ene_win; if(MOD%7==6){//column shadowed by uBall float rng_uball = re_csi->RANGE(z,a,ene_air)+uball_thck; float ene_uball = re_csi->ENERGY(z,a,rng_uball); float de_uball = ene_uball-ene_air; dat.de_air += de_uball; ene_air = ene_uball; } rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); dat.de_tgt = ene_tgt-ene_air; dat.Energy = ene_tgt; // enemyl1 = re_myl->ENERGY(z,a,myl1_thck[MOD]); // rngcsi1 = re_csi->RANGE(z,a,enemyl1)+csi1_thck[MOD]; // de_ptr1= H2E_PTHR(dat.ch_csi1,iz,ia,a11CORR,a21[MOD],a31[MOD],a41[MOD],fct_si1[MOD],rngcsi1, re_csi); dat.RangeCsI = rng_csi1/10000.; //[cm] dat.DeltaECsI = re_csi->ENERGY(z,a,rng_csi1); } // else if(dat.PunchThrough){ // de_ptr2= H2E_PTHR(dat.ch_csi2, iz, ia, a12CORR,a22[MOD],a32[MOD],a42[MOD],fct_si2[MOD],csi2_thck[MOD], re_csi); // float rng_csi2 = re_csi->RANGE(z,a,dat.de_ptr2)+csi2_thck[MOD]; // float ene_csi2 = re_csi->ENERGY(z,a,rng_csi2); // dat.de_csi2 = ene_csi2-dat.de_ptr2; //// float rng_myl1 = re_myl->RANGE(z,a,ene_csi2)+myl1_thck[MOD]; //// float ene_myl1 = re_myl->ENERGY(z,a,rng_myl1); //// dat.de_myl1 = ene_myl1-dat.de_csi2; //// float rng_csi1 = re_csi->RANGE(z,a,ene_myl1)+csi1_thck[MOD]; // float ene_csi1 = re_csi->ENERGY(z,a,rng_csi1); // dat.de_csi1 = ene_csi1-ene_myl1; // float rng_myl0 = re_myl->RANGE(z,a,ene_csi1)+myl0_thck[MOD]; // float ene_myl0 = re_myl->ENERGY(z,a,rng_myl0); // dat.de_myl0 = ene_myl0-ene_csi1; // dat.de_sil1 = dat.ch_sil1/fct_si1[MOD]; // float rng_ddl1 = re_si->RANGE(z,a,dat.de_sil1+ene_myl0)+dl1_thck[MOD]; // float ene_ddl1 = re_si->ENERGY(z,a,rng_ddl1); // dat.de_ddl1 = ene_ddl1-(dat.de_sil1+ene_myl0); // dat.de_sil0 = dat.ch_sil0/fct_si0[MOD]; // float rng_ddl0 = re_si->RANGE(z,a,dat.de_sil0+ene_ddl1)+dl0_thck; // float ene_ddl0 = re_si->ENERGY(z,a,rng_ddl0); // dat.de_ddl0 = ene_ddl0-(dat.de_sil0+ene_ddl1); // float rng_win = re_cu->RANGE(z,a,ene_ddl0)+win_thck; // float ene_win = re_cu->ENERGY(z,a,rng_win); // dat.de_win = ene_win-ene_ddl0; // float rng_air = re_air->RANGE(z,a,ene_win)+air_thck; // float ene_air = re_air->ENERGY(z,a,rng_air); // dat.de_air = ene_air-ene_win; // float rng_tgt = re_tgt->RANGE(z,a,ene_air)+tgt_thck[MOD]; // float ene_tgt = re_tgt->ENERGY(z,a,rng_tgt); // dat.de_tgt = ene_tgt-ene_air; // if(de_ptr2>0) hz2a4_ptr2->Fill(ene_tgt); // if(de_ptr2>0) hz2a4_totl->Fill(ene_tgt); // //printf("pthr %8.2f %8.2f\n",de_ptr2/float(ia),ene_tgt/float(ia)); // } // else{ // //printf("I am lost\n"); // } if(dat.Energy<0){ dat.Energy=0; dat.RangeCsI=0; } Double_t theta_i = dat.ThetaRndm; Double_t phi_i = dat.PhiRndm; Double_t m, m2, K, E; Double_t p, p_p, p_t; K = dat.Energy; m = u*Double_t(ia); m2 = m*m; E = K+m; p = (E*E>=m2) ? sqrt(E*E-m2) : 0.; p_p = p*cos(theta_i); p_t = p*sin(theta_i); dat.Pt = p_t; dat.Px = p_t*cos(phi_i); dat.Py = p_t*sin(phi_i); dat.Pz = p_p; dat.RapidityLab = TMath::ATanH(p_p/E); dat.RapidityCM = dat.RapidityLab-kevt.CMRapidity; dat.Xt = p_t/m; dat.Winv = 0.; if(p_t>0){ dat.Winv = 1./(m2*dat.Xt);// *1/d_phi } dat.Gam = E/m; dat.WeightCol = 1; dat.WeightInt = 1; dat.WeightGeant = 1; //Sihver Formula: Phys.Rev.C 47 1225 (1993) if(dat.RangeCsI>0){ float Ap = float(ia); const float At = 0.5*(132.9054519 + 126.90447); // A CsI const float RHO = 4.51; // CsI (* g/cm^3 *) const float Na = 6.0221415*1.e23; //Avogadro (* 1/mol *) float b0 = 1.581 - 0.876*(pow(Ap,-0.33333) + pow(At,-0.33333)); float r0 = 1.166; //(*fm*) giving interaction length of 38.0 cm for protons float SigR = pi*r0*r0*pow(pow(Ap,0.33333) + pow(At,0.33333) - b0*(pow(Ap,-0.33333) + pow(At,-0.33333)),2)*1.e-26;// (*cm^2*) sigma // = sigma_TOT - sigma_elastic - sigma_quasielastic float lambdaTreac = At/(SigR*Na*RHO); float preac = 1.f - exp(-dat.RangeCsI/lambdaTreac);// reaction probability dat.WeightInt = 1.f/(1.f-preac); r0 = 1.522; //(*fm*) giving collision length of 22.3 cm for protons; SigR = pi*r0*r0*pow(pow(Ap,0.33333) + pow(At,0.33333) - b0*(pow(Ap,-0.33333) + pow(At,-0.33333)),2)*1.e-26;// (*cm^2*) sigma // = sigma_TOT float lambdaTcoll = At/(SigR*Na*RHO); float pcoll = 1.f - exp(-dat.RangeCsI/lambdaTcoll);// collision probability dat.WeightCol = 1.f/(1.f-pcoll); // printf("WGHT: Z %2d A %2d E %7.1f R %7.1f (cm) PInt %5.2f PCol %5.2f WInt %5.3f WCol %5.3f\n", // iz,ia,dat.Energy,dat.RangeCsI,preac*100,pcoll*100,dat.WeightInt,dat.WeightCol); } /// Geant efficency correction: dat.WeightGeant = EfficencyCorrection(dat.Z, dat.A, dat.Energy ); //cout << "dat.Z = " << dat.Z; //cout << ", dat.A = " << dat.A; //cout << ", dat.Energy = " << dat.Energy; //cout << ", dat.WeightGeant = " << dat.WeightGeant; //cout << endl; /// Filling the PID flag (See: Swiatowid.cxx) /// EMPTY = 0b00000000, // No signal /// PD0 = 0b00000001, // PD0 fired /// PD1 = 0b00000010, // PD1 fired /// PD2 = 0b00000100, // PD2 fired /// SISI = 0b00001000, // SI0 vs SI1 /// PUNCH = 0b00010000, // CsI2 punched through /// GAMMA = 0b00100000, // Gamma particle /// NEUTRON = 0b01000000, // Neutron /// BCKGRD = 0b10000000, // background /// Commented 2014-02-18 (SK) /// if ( dat.ch_sil0 > 0. ){ dat.PID |= 0b00000001;} /// PD0 /// if ( dat.ch_csi1 > 0. ){ dat.PID |= 0b00000010;} /// PD1 /// if ( dat.ch_csi2 > 0. ){ dat.PID |= 0b00000100;} /// PD2 /// /// if( dat.PunchThrough ){ dat.PID |= 0b00010000;} /// PUNCH /// if( dat.Gamma ){ dat.PID |= 0b00100000;} /// GAMMA /// if( !dat.Star ){ dat.PID |= 0b10000000;} /// BCKGRD } //______________________________________________________________________ float KrattaEne::H2E(float fH, int iz, int ia, float fa1,float fa2,float fa3,float fa4, float fct_si){ /* returns the energy for a given total light H, using Parlog formula */ int iblag; double H; double a1, a2, a3, a4; H = double(fH/fct_si*fa1); a1= double(fa1); a2= double(fa2); a3= double(fa3); a4= double(fa4); // integer mod_etal(10:17)/3*2,4,3,3,2,2/ double E, E0; double a, z; double az2, fbe; bool cond; int niter; double eps=1.e-6; double argx, fx, fpx, dx; E = 0.; iblag = 0; if(H<=0.){ iblag = 5; return float(E); } ia = abs(ia); if(iz <= 0) return float(E); z = double(iz); a = double(ia); //a = (ia > 0) ? double(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } if(a1==-1.){ iblag=15; return float(E); } double a2az2 = a2*az2; double x, a4e; niter=0; E0=a3*a; if(a2 < 1.e-7) x = H; else{ x = H; argx = 1.+ x/a2az2; cond=1; while(cond){ niter++; a4e = (x < E0)? 0. : a4; fx = - H+fbe*(x-a2az2*log(argx)+a2az2*a4e *log((x+a2az2)/(E0+a2az2))); fpx = fbe*(x+a2az2*a4e)/(x+a2az2); dx = -fx/fpx; x = x+dx; argx = 1.+ x/a2az2; cond = fabs(dx/x) >= eps && argx > 0. && niter <= 1000; } if(argx<=0.){ iblag=13; return float(E); } else if(niter > 500){ iblag=14; return float(E); } } E=fbe*x; return float(E); } //______________________________________________________________________ float KrattaEne::H2E_GAMMA(float fH, float fa1,float fa2,float fa3,float fa4, float fct_si){ /* returns the gamma energy for a given total light H assuming gammas are converted into electrons use with caution... not tested... */ int iblag; double H; double a1, a2, a3, a4; H = double(fH/fct_si*fa1); a1= double(fa1); a2= double(fa2); a3= double(fa3); a4= double(fa4); // integer mod_etal(10:17)/3*2,4,3,3,2,2/ double E, E0; double a, z; double az2, fbe; bool cond; int niter; double eps=1.e-6; double argx, fx, fpx, dx; E = 0.; iblag = 0; if(H<=0.){ iblag = 5; return float(E); } z = 1; a = 0.511/931.49386; //Electron; az2 = a*z*z; fbe = 1.; if(a1==-1.){ iblag=15; return float(E); } double a2az2 = a2*az2; double x, a4e; niter=0; E0=a3*a; if(a2 < 1.e-7) x = H; else{ x = H; argx = 1.+ x/a2az2; cond=1; while(cond){ niter++; a4e = (x < E0)? 0. : a4; fx = - H+fbe*(x-a2az2*log(argx)+a2az2*a4e *log((x+a2az2)/(E0+a2az2))); fpx = fbe*(x+a2az2*a4e)/(x+a2az2); dx = -fx/fpx; x = x+dx; argx = 1.+ x/a2az2; cond = fabs(dx/x) >= eps && argx > 0. && niter <= 1000; } if(argx<=0.){ iblag=13; return float(E); } else if(niter > 100){ iblag=14; return float(E); } } E=fbe*x; return float(E); } //______________________________________________________________________ float KrattaEne::E2H(float e_csi, int iz, int ia, float a1,float a2,float a3,float a4, float fct_si){ /* original Parlog returns the total light H for a given energy ene */ // NOTICE: output H is not gain corrected float H; float a, z; float az2, fbe; H = 0.; ia = abs(ia); if(iz <= 0) return H; z = float(iz); a = float(ia); //a = (ia > 0) ? float(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } float a2az2 = a2*az2; H = (a2 != 0.) ? (e_csi-a2az2*log(1.+e_csi/a2az2))+ a4*a2az2*log((e_csi+a2az2)/(a3*a+a2az2)) : e_csi; if(H<0) H=0; H *= fbe*fct_si/a1; return H; } //______________________________________________________________________ float KrattaEne::H2E_PTHR(float fH, int iz, int ia, float fa1,float fa2, float fct_si, float csi2_thcki, RE_TAB *re_csi1,float off){ /* returns the energy for a given total light H */ int iblag; double H; double a1, a2, a3, a4; a2 = double(fa2); a1 = double(fa1); H = double(fH/fct_si*a1); // integer mod_etal(10:17)/3*2,4,3,3,2,2/ double E; double a, z; double az2, fbe; int niter; E = 0.; iblag = 0; if(H<=0.){ iblag = 5; return float(E); } ia = abs(ia); if(iz <= 0) return float(E); z = double(iz); a = double(ia); //a = (ia > 0) ? double(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } if(a1==-1.){ iblag=15; return float(E); } double a2az2 = a2*az2; double x; niter=0; // if(a2 < 1.e-7) x = H; // else{ double eps = 1.e-6; int nsig = 6; int maxfn = 1000; int ier=0; double aa = 0; double bb = 400*a; int ic; double t,fa,fb,c,fc,d,e,tol,rm,s,p,q,r,rone,temp; ier = 0; t = pow(10.0f,-nsig); ic = 2; // float ra = re_csi1->RANGE(iz,float(a),aa); // float ea = re_csi1->ENERGY(iz,float(a),ra+csi2_thcki); // // if(a2>0) fa = ea-aa-a2az2*log((ea+a2az2)/(aa+a2az2))-H; // // else fa = ea-aa-H; // float Hao = aa-a2az2*log(1.f+aa/a2az2); // float Hai = ea-a2az2*log(1.f+ea/a2az2); // fa = Hai-Hao-H; // float rb = re_csi1->RANGE(iz,float(a),bb); // float eb = re_csi1->ENERGY(iz,float(a),rb+csi2_thcki); // // if(a2>0) fb = eb-bb-a2az2*log((eb+a2az2)/(bb+a2az2))-H; // // else fb = eb-bb-H; // float Hbo = bb-a2az2*log(1.f+bb/a2az2); // float Hbi = eb-a2az2*log(1.f+eb/a2az2); // fb = Hbi-Hbo-H; float ra = re_csi1->RANGE(iz,float(a),aa); float ea = re_csi1->ENERGY(iz,float(a),ra+csi2_thcki); float Hao = 0; if(aa>off) Hao = aa-off-a2az2*log(1.f+(aa-off)/a2az2); float Hai = 0; if(ea>off) Hai = ea-off-a2az2*log(1.f+(ea-off)/a2az2); fa = Hai-Hao-H; float rb = re_csi1->RANGE(iz,float(a),bb); float eb = re_csi1->ENERGY(iz,float(a),rb+csi2_thcki); float Hbo = 0; if(bb>off) Hbo = bb-off-a2az2*log(1.f+(bb-off)/a2az2); float Hbi = 0; if(eb>off) Hbi = eb-off-a2az2*log(1.f+(eb-off)/a2az2); fb = Hbi-Hbo-H; // if(fH>930){ // //printf("O aa bb %f %f %f %f %f %f %f %f\n",ea,eb,fa,fb,ra,ea,rb,eb); // printf("O %f %f %f %f %f %f %f %f %f\n",fH,ea,eb,fa,fb,Hao,Hai,Hbo,Hbi); // } //printf("aa bb %f %f %f %f\n",ea,eb,fa,fb); // test for same sign if(fa*fb>0.f){//terminal error - f(a) and f(b) have the same sign ier = 130; maxfn = ic; x = aa; goto end; } lab5: c = aa; fc = fa; d = bb-c; e = d; lab10: if(fabs(fc)>=fabs(fb)) goto lab15; aa = bb; bb = c; c = aa; fa = fb; fb = fc; fc = fa; lab15: tol = (fabs(bb)>=0.1) ? t*fabs(bb) : 0.1*t; rm = (c-bb)*0.5f; // test for convergence criteria if(fabs(fb)<=eps || fabs(c-bb)<=tol){//convergence of b aa = c; maxfn = ic; //printf("a,b %f %f\n",aa,bb); x = bb; goto end; } // check evaluation counter if(ic>=maxfn){// maxfn evaluations ier = 129; aa = c; maxfn = ic; x = bb; goto end; } // is bisection forced if(fabs(e)0.f) q = -q; if(p<0.f) p = -p; s = e; e = d; // if abs(p/q).ge.75*abs(c-b) then // force bisection if(p+p>=3.0f*rm*q) goto lab30; // if abs(p/q).ge..5*abs(s) then force // bisection. s = the value of p/q // on the step before the last one if(p+p>=fabs(s*q)) goto lab30; d = p/q; goto lab35; // bisection lab30: e = rm; d = e; // increment b lab35: aa = bb; fa = fb; temp = d; if(fabs(temp)<=0.5f*tol){ temp = 0.f; if(rm!=0.f) temp = 0.5f*tol*rm/fabs(rm); } bb = bb+temp; rb = re_csi1->RANGE(iz,float(a),bb); eb = re_csi1->ENERGY(iz,float(a),rb+csi2_thcki); // if(a2>0) fb = eb-bb-a2az2*log((eb+a2az2)/(bb+a2az2))-H; // else fb = eb-bb-H; Hbo = 0; if(bb>off) Hbo = bb-off-a2az2*log(1.f+(bb-off)/a2az2); Hbi = 0; if(eb>off) Hbi = eb-off-a2az2*log(1.f+(eb-off)/a2az2); fb = Hbi-Hbo-H; ic = ic+1; if(fb*fc<=0.f) goto lab10; goto lab5; // } end: E=fbe*x; return float(E); } //______________________________________________________________________ float KrattaEne::H2E_PTHRI(float fH, int iz, int ia, float fa1, float fa2, float fct_si, float csi2_thcki, RE_TAB *re_csi1,float off){ /* returns the energy for a given total light H */ int iblag; double H; double a1, a2, a3, a4; a2 = double(fa2); a1 = double(fa1); H = double(fH/fct_si*a1); // integer mod_etal(10:17)/3*2,4,3,3,2,2/ double E; double a, z; double az2, fbe; int niter; E = 0.; iblag = 0; if(H<=0.){ iblag = 5; return float(E); } ia = abs(ia); if(iz <= 0) return float(E); z = double(iz); a = double(ia); //a = (ia > 0) ? double(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } if(a1==-1.){ iblag=15; return float(E); } double a2az2 = a2*az2; double x; niter=0; // if(a2 < 1.e-7) x = H; // else{ double eps = 1.e-6; int nsig = 6; int maxfn = 1000; int ier=0; double aa = 0; double bb = 400*a; int ic; double t,fa,fb,c,fc,d,e,tol,rm,s,p,q,r,rone,temp; ier = 0; t = pow(10.0f,-nsig); ic = 2; float ra = re_csi1->RANGE(iz,float(a),aa); float ea = re_csi1->ENERGY(iz,float(a),ra+csi2_thcki); float Hao = 0; if(aa>off) Hao = aa-off+a2az2*log(1.f+(aa-off)/a2az2); float Hai = 0; if(ea>off) Hai = ea-off+a2az2*log(1.f+(ea-off)/a2az2); fa = Hai-Hao-H; float rb = re_csi1->RANGE(iz,float(a),bb); float eb = re_csi1->ENERGY(iz,float(a),rb+csi2_thcki); float Hbo = 0; if(bb>off) Hbo = bb-off+a2az2*log(1.f+(bb-off)/a2az2); float Hbi = 0; if(eb>off) Hbi = eb-off+a2az2*log(1.f+(eb-off)/a2az2); fb = Hbi-Hbo-H; // if(fH>925){ // printf("I %f %f %f %f %f %f %f %f %f\n",fH,ea,eb,fa,fb,Hao,Hai,Hbo,Hbi); // } // test for same sign if(fa*fb>0.f){//terminal error - f(a) and f(b) have the same sign ier = 130; maxfn = ic; x = aa; goto end; } lab5: c = aa; fc = fa; d = bb-c; e = d; lab10: if(fabs(fc)>=fabs(fb)) goto lab15; aa = bb; bb = c; c = aa; fa = fb; fb = fc; fc = fa; lab15: tol = (fabs(bb)>=0.1) ? t*fabs(bb) : 0.1*t; rm = (c-bb)*0.5f; // test for convergence criteria if(fabs(fb)<=eps || fabs(c-bb)<=tol){//convergence of b aa = c; maxfn = ic; //printf("a,b %f %f\n",aa,bb); x = bb; goto end; } // check evaluation counter if(ic>=maxfn){// maxfn evaluations ier = 129; aa = c; maxfn = ic; x = bb; goto end; } // is bisection forced if(fabs(e)0.f) q = -q; if(p<0.f) p = -p; s = e; e = d; // if abs(p/q).ge.75*abs(c-b) then // force bisection if(p+p>=3.0f*rm*q) goto lab30; // if abs(p/q).ge..5*abs(s) then force // bisection. s = the value of p/q // on the step before the last one if(p+p>=fabs(s*q)) goto lab30; d = p/q; goto lab35; // bisection lab30: e = rm; d = e; // increment b lab35: aa = bb; fa = fb; temp = d; if(fabs(temp)<=0.5f*tol){ temp = 0.f; if(rm!=0.f) temp = 0.5f*tol*rm/fabs(rm); } bb = bb+temp; rb = re_csi1->RANGE(iz,float(a),bb); eb = re_csi1->ENERGY(iz,float(a),rb+csi2_thcki); Hbo = 0; if(bb>off) Hbo = bb-off+a2az2*log(1.f+(bb-off)/a2az2); Hbi = 0; if(eb>off) Hbi = eb-off+a2az2*log(1.f+(eb-off)/a2az2); fb = Hbi-Hbo-H; ic = ic+1; if(fb*fc<=0.f) goto lab10; goto lab5; // } end: E=fbe*x; return float(E); } //______________________________________________________________________ float KrattaEne::SI01_PTHR(float ch_sil0, float ch_sil1,int iz, int ia, float fct_sil0, float fct_sil1, float si0_thcki,float dl1_thcki,float si1_thcki, RE_TAB *re_si1){ /* returns the residual energy for a given loss in PD0+PD1 */ int iblag; double E; double a, z; double az2, fbe; int niter; double dE0 = ch_sil0/fct_sil0; double dE1 = ch_sil1/fct_sil1; E = 0.; iblag = 0; ia = abs(ia); if(iz <= 0) return float(E); z = double(iz); a = double(ia); //a = (ia > 0) ? double(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } double x; niter=0; double eps = 1.e-6; int nsig = 6; int maxfn = 1000; int ier=0; double aa = 0; double bb = 600*a; int ic; double t,fa,fb,c,fc,d,e,tol,rm,s,p,q,r,rone,temp; ier = 0; t = pow(10.0f,-nsig); ic = 2; // float ra = re_si1->RANGE(iz,float(a),aa); // float ea1 = re_si1->ENERGY(iz,float(a),ra+si1_thcki); // float ed1 = re_si1->ENERGY(iz,float(a),ra+si1_thcki+dl1_thcki); // float ea0 = re_si1->ENERGY(iz,float(a),ra+si1_thcki+dl1_thcki+si0_thcki); // fa = ea0-ed1+ea1-aa-dE1-dE0; // float rb = re_si1->RANGE(iz,float(a),bb); // float eb1 = re_si1->ENERGY(iz,float(a),rb+si1_thcki); // ed1 = re_si1->ENERGY(iz,float(a),rb+si1_thcki+dl1_thcki); // float eb0 = re_si1->ENERGY(iz,float(a),rb+si1_thcki+dl1_thcki+si0_thcki); // fb = eb0-ed1+eb1-bb-dE1-dE0; float ra = re_si1->RANGE(iz,float(a),aa+dE1); float ed1 = re_si1->ENERGY(iz,float(a),ra+dl1_thcki); float ea0 = re_si1->ENERGY(iz,float(a),ra+dl1_thcki+si0_thcki); fa = ea0-ed1-dE0; float rb = re_si1->RANGE(iz,float(a),bb+dE1); ed1 = re_si1->ENERGY(iz,float(a),rb+dl1_thcki); float eb0 = re_si1->ENERGY(iz,float(a),rb+dl1_thcki+si0_thcki); fb = eb0-ed1-dE0; //printf("aa bb %f %f %f %f\n",ea,eb,fa,fb); // test for same sign if(fa*fb>0.f){//terminal error - f(a) and f(b) have the same sign ier = 130; maxfn = ic; x = aa; goto end; } lab5: c = aa; fc = fa; d = bb-c; e = d; lab10: if(fabs(fc)>=fabs(fb)) goto lab15; aa = bb; bb = c; c = aa; fa = fb; fb = fc; fc = fa; lab15: tol = (fabs(bb)>=0.1) ? t*fabs(bb) : 0.1*t; rm = (c-bb)*0.5f; // test for convergence criteria if(fabs(fb)<=eps || fabs(c-bb)<=tol){//convergence of b aa = c; maxfn = ic; //printf("a,b %f %f\n",aa,bb); x = bb; goto end; } // check evaluation counter if(ic>=maxfn){// maxfn evaluations ier = 129; aa = c; maxfn = ic; x = bb; goto end; } // is bisection forced if(fabs(e)0.f) q = -q; if(p<0.f) p = -p; s = e; e = d; // if abs(p/q).ge.75*abs(c-b) then // force bisection if(p+p>=3.0f*rm*q) goto lab30; // if abs(p/q).ge..5*abs(s) then force // bisection. s = the value of p/q // on the step before the last one if(p+p>=fabs(s*q)) goto lab30; d = p/q; goto lab35; // bisection lab30: e = rm; d = e; // increment b lab35: aa = bb; fa = fb; temp = d; if(fabs(temp)<=0.5f*tol){ temp = 0.f; if(rm!=0.f) temp = 0.5f*tol*rm/fabs(rm); } bb = bb+temp; // rb = re_si1->RANGE(iz,float(a),bb); // eb1 = re_si1->ENERGY(iz,float(a),rb+si1_thcki); // ed1 = re_si1->ENERGY(iz,float(a),rb+si1_thcki+dl1_thcki); // eb0 = re_si1->ENERGY(iz,float(a),rb+si1_thcki+dl1_thcki+si0_thcki); // fb = eb0-ed1+eb1-bb-dE1-dE0; rb = re_si1->RANGE(iz,float(a),bb+dE1); ed1 = re_si1->ENERGY(iz,float(a),rb+dl1_thcki); eb0 = re_si1->ENERGY(iz,float(a),rb+dl1_thcki+si0_thcki); fb = eb0-ed1-dE0; ic = ic+1; if(fb*fc<=0.f) goto lab10; goto lab5; // } end: E=fbe*x; return float(E); } //______________________________________________________________________ float KrattaEne::SI0_PTHR(float ch_sil0, int iz, int ia, float fct_sil0, float si0_thcki, RE_TAB *re_si1){ /* returns the residual energy for a given loss in PD0 */ int iblag; double E; double a, z; double az2, fbe; int niter; double dE0 = ch_sil0/fct_sil0; E = 0.; iblag = 0; ia = abs(ia); if(iz <= 0) return float(E); z = double(iz); a = double(ia); //a = (ia > 0) ? double(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } double x; niter=0; double eps = 1.e-6; int nsig = 6; int maxfn = 1000; int ier=0; double aa = 0; double bb = 800*a; int ic; double t,fa,fb,c,fc,d,e,tol,rm,s,p,q,r,rone,temp; ier = 0; t = pow(10.0f,-nsig); ic = 2; float ra = re_si1->RANGE(iz,float(a),aa); float ea0 = re_si1->ENERGY(iz,float(a),ra+si0_thcki); fa = ea0-aa-dE0; float rb = re_si1->RANGE(iz,float(a),bb); float eb0 = re_si1->ENERGY(iz,float(a),rb+si0_thcki); fb = eb0-bb-dE0; //printf("aa bb %f %f %f %f\n",ea,eb,fa,fb); // test for same sign if(fa*fb>0.f){//terminal error - f(a) and f(b) have the same sign ier = 130; maxfn = ic; x = aa; goto end; } lab5: c = aa; fc = fa; d = bb-c; e = d; lab10: if(fabs(fc)>=fabs(fb)) goto lab15; aa = bb; bb = c; c = aa; fa = fb; fb = fc; fc = fa; lab15: tol = (fabs(bb)>=0.1) ? t*fabs(bb) : 0.1*t; rm = (c-bb)*0.5f; // test for convergence criteria if(fabs(fb)<=eps || fabs(c-bb)<=tol){//convergence of b aa = c; maxfn = ic; //printf("a,b %f %f\n",aa,bb); x = bb; goto end; } // check evaluation counter if(ic>=maxfn){// maxfn evaluations ier = 129; aa = c; maxfn = ic; x = bb; goto end; } // is bisection forced if(fabs(e)0.f) q = -q; if(p<0.f) p = -p; s = e; e = d; // if abs(p/q).ge.75*abs(c-b) then // force bisection if(p+p>=3.0f*rm*q) goto lab30; // if abs(p/q).ge..5*abs(s) then force // bisection. s = the value of p/q // on the step before the last one if(p+p>=fabs(s*q)) goto lab30; d = p/q; goto lab35; // bisection lab30: e = rm; d = e; // increment b lab35: aa = bb; fa = fb; temp = d; if(fabs(temp)<=0.5f*tol){ temp = 0.f; if(rm!=0.f) temp = 0.5f*tol*rm/fabs(rm); } bb = bb+temp; rb = re_si1->RANGE(iz,float(a),bb); eb0 = re_si1->ENERGY(iz,float(a),rb+si0_thcki); fb = eb0-bb-dE0; ic = ic+1; if(fb*fc<=0.f) goto lab10; goto lab5; // } end: E=fbe*x; return float(E); } //______________________________________________________________________ float KrattaEne::NH2E(float xx, int iz, int ia, float a1f,float a2f,float offf, float fct_sii){ /* returns the energy for a given total light H, using Parlog formula */ double fH = double(xx);// double fa2 = double(a2f); double fct_si= double(fct_sii); double off = double(offf); double fa1 = double(a1f); int iblag; double H; double a1, a2; H = double(fH/fct_si*fa1); a1= double(fa1); a2= double(fa2); // integer mod_etal(10:17)/3*2,4,3,3,2,2/ double E; double a, z; double az2, fbe; bool cond; int niter; double eps=1.e-6; double argx, fx, fpx, dx; E = 0.; iblag = 0; if(H<=0.){ iblag = 5; return float(E); } ia = abs(ia); if(iz <= 0) return float(E); z = double(iz); a = double(ia); //a = (ia > 0) ? double(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } if(a1==-1.){ iblag=15; return float(E); } double a2az2 = a2*az2; double x; niter=0; if(a2 < 1.e-7) x = H; else{ x = H; argx = 1.+ x/a2az2; cond=1; while(cond){ niter++; fx = - H+fbe*(x-a2az2*log(argx)); fpx = fbe*x/(x+a2az2); dx = -fx/fpx; x = x+dx; argx = 1.+ x/a2az2; cond = fabs(dx/x) >= eps && argx > 0. && niter <= 1000; } if(argx<=0.){ iblag=13; return float(E); } else if(niter > 100){ iblag=14; return float(E); } } E=fbe*x+off; return float(E); } //______________________________________________________________________ float KrattaEne::NH2EI(float xx, int iz, int ia, float a1f,float a2f,float offf, float fct_sii){ /* returns the energy for a given total light H, using "inverse" Parlog with +Log() formula */ double fH = double(xx);// double fa2 = double(a2f); double fct_si= double(fct_sii); double off = double(offf); double fa1 = double(a1f); int iblag; double H; double a1, a2; H = double(fH/fct_si*fa1); a1= double(fa1); a2= double(fa2); // integer mod_etal(10:17)/3*2,4,3,3,2,2/ double E; double a, z; double az2, fbe; bool cond; int niter; double eps=1.e-6; double argx, fx, fpx, dx; E = 0.; iblag = 0; // if(H<=0.){ // iblag = 5; // return double(E); // } ia = abs(ia); if(iz <= 0) return float(E); z = double(iz); a = double(ia); //a = (ia > 0) ? double(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } if(a1==-1.){ iblag=15; return float(E); } double a2az2 = a2*az2; double x; niter=0; if(a2 < 1.e-7) x = H; else{ x = H; argx = 1.+ x/a2az2; cond=1; while(cond){ niter++; fx = - H+fbe*(x+a2az2*log(argx)); fpx = fbe*(x+2.*a2az2)/(x+a2az2); dx = -fx/fpx; x = x+dx; argx = 1.+ x/a2az2; cond = fabs(dx/x) >= eps && argx > 0. && niter <= 1000; } if(argx<=0.){ iblag=13; printf("iblag %d %f %f\n",iblag,argx,fbe*x+off); return float(E); } else if(niter > 500){ iblag=14; printf("iblag %d %f %f\n",iblag,argx,fbe*x+off); return float(E); } } E=fbe*x+off; return float(E); } //______________________________________________________________________ float KrattaEne::NE2H(float xx, int iz, int ia, float a1f,float a2f,float offf, float fct_sii){ /* original Parlog returns the total light H for a given energy ene */ // NOTICE: output H is not gain corrected double e_csi = double(xx);// float a2 = a2f; float fct_si= fct_sii; float off = offf; float a1 = a1f; double H; double a, z; double az2, fbe; H = 0.; ia = abs(ia); if(iz <= 0) return float(H); z = double(iz); a = double(ia); //a = (ia > 0) ? float(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } double a2az2 = a2*az2; H = (a2 != 0.) ? (e_csi-off-a2az2*log(1.+(e_csi-off)/a2az2)) : e_csi-off; if(H<0) H=0; H *= fbe*fct_si/a1; return float(H); } //______________________________________________________________________ float KrattaEne::NE2HI(float xx, int iz, int ia, float a1f,float a2f,float offf, float fct_sii){ /* "inverse" Parlog with + Log() returns the total light H for a given energy ene */ // NOTICE: output H is not gain corrected double e_csi = double(xx);// float a2 = a2f; float fct_si= fct_sii; float off = offf; float a1 = a1f; double H; double a, z; double az2, fbe; H = 0.; ia = abs(ia); if(iz <= 0) return float(H); z = double(iz); a = double(ia); //a = (ia > 0) ? float(ia) : 2.08*z + 0.0029*z*z; if(iz==1 && ia==-1) a= 0.511/931.49386; //Gamma if(iz==1 && ia==-2) a= 0.511/931.49386; //Electron if(iz==1 && ia==-3) a=139.56995/931.49386; //Pion if(iz==1 && ia==-4) a=105.65841/931.49386; //Muon if(iz==1 && ia==-5) a=493.67700/931.49386; //Kaon if(ia == 8 && iz == 4){ // Be8 = 2* He4 az2 = a*z*z/8.; fbe = 2.; } else{ // tout le reste az2 = a*z*z; fbe = 1.; } double a2az2 = a2*az2; H = (a2 != 0.) ? (e_csi-off+a2az2*log(1.+(e_csi-off)/a2az2)) : e_csi-off; if(H<0) H=0; H *= fbe*fct_si/a1; return float(H); } //______________________________________________________________________ float KrattaEne::EfficencyCorrection( int iz, int ia, float E_kin_in_MeV) { float efficencyCorrection = -1.; float p0 = 0.; float p1 = 0.; float p2 = 0.; float E = E_kin_in_MeV/1000.; /// E [GeV] see http://kratta.ifj.edu.pl/ float reactionProb = 0.; /// Reaction probability /// Energy limitation for this Geant correction in [GeV] /// proton E < 0.250 /// deuteron E < 0.320 /// tritium E < 0.380 /// 3He E < 0.880 /// 4He E < 0.980 /// Z=1 if(iz==1 && ia==1 && E < 0.250 ){ p0=-0.00022196; p1=0.0944907; p2=7.3126; }; //=proton if(iz==1 && ia==2 && E < 0.320 ){ p0=-0.0145808; p1=0.620002; p2=2.36021; }; //=deuteron if(iz==1 && ia==3 && E < 0.380 ){ p0=-0.02462; p1=0.751908; p2=1.35404; }; //=tritium /// Z=2 if(iz==2 && ia==3 && E < 0.880 ){ p0=-0.0194343; p1=0.287043; p2=0.281509; }; //=3He if(iz==2 && ia==4 && E < 0.980 ){ p0=-0.0253448; p1=0.312913; p2=0.198148; }; //=4He //if(iz==2 && ia==6 && E < ){ p0=0 ;p1=0 ;p2=0; }; ; //=6He reactionProb = p0+p1*E+p2*E*E; if( 0. <= reactionProb && reactionProb < 1.){ efficencyCorrection = 1. / (1. - reactionProb ); } //if( ia > 0 && iz > 0 ){ // cout << "iz = " << iz; // cout << ", ia = " << ia; // cout << ", E = " << E; // cout << ", reactionProb = " << reactionProb; // cout << ", efficencyCorrection = " << efficencyCorrection; // cout << endl; //} return efficencyCorrection; } //______________________________________________________________________ //void KrattaEne::EfficencyCorrection( TRootKRATParticle* dat ) //{ // float rp = EfficencyCorrection(dat->A, dat->Z, dat->Energy); // dat->??? *= rp; // return; //} //////////////////////// GLOBAL FUNCTIONS ////////////////////////////// //______________________________________________________________________ void KRATTAene(Int_t mod=7){ printf("module %d\n",mod); KrattaEne *kene = new KrattaEne(); // kene->LoadParametersOld(1555); // kene->LoadGammaLines(1601); // kene->LoadStarCuts(1601); // kene->LoadStarCuts(3601); kene->LoadParameters(1555); kene->LoadFineCorrections(1555); kene->LoadPD01PthrChannels(1555); // // TH1F *h4[3]; // TH1F *h3[3]; // for(int i=0;i<3;i++){ // h4[i] = new TH1F(Form("h4%d",i),Form("h4%d",i),100,0.8,1); // h3[i] = new TH1F(Form("h3%d",i),Form("h3%d",i),100,0.8,1); // } // float fct4he=0; // float nfct4he=0; // float fct3he=0; // float nfct3he=0; // int imod; // int NMOD=35; // for(imod=0;imodcorr_fct[imod][1][4]!=1){ // printf("%2d 4He %10.6f %10.6f %10.6f\n",imod, // kene->corr_fct[imod][1][4]/kene->corr_fct[imod][1][0], // kene->corr_fct[imod][1][4]/kene->corr_fct[imod][1][1], // kene->corr_fct[imod][1][4]/kene->corr_fct[imod][1][2]); // fct4he += kene->corr_fct[imod][1][4]/kene->corr_fct[imod][1][2]; // nfct4he += 1; // h4[0]->Fill(kene->corr_fct[imod][1][4]/kene->corr_fct[imod][1][0]); // h4[1]->Fill(kene->corr_fct[imod][1][4]/kene->corr_fct[imod][1][1]); // h4[2]->Fill(kene->corr_fct[imod][1][4]/kene->corr_fct[imod][1][2]); // } // } // for(imod=0;imodcorr_fct[imod][1][3]!=1){ // printf("%2d 3He %10.6f %10.6f %10.6f\n",imod, // kene->corr_fct[imod][1][3]/kene->corr_fct[imod][1][0], // kene->corr_fct[imod][1][3]/kene->corr_fct[imod][1][1], // kene->corr_fct[imod][1][3]/kene->corr_fct[imod][1][2]); // fct3he += kene->corr_fct[imod][1][3]/kene->corr_fct[imod][1][2]; // nfct3he += 1; // h3[0]->Fill(kene->corr_fct[imod][1][3]/kene->corr_fct[imod][1][0]); // h3[1]->Fill(kene->corr_fct[imod][1][3]/kene->corr_fct[imod][1][1]); // h3[2]->Fill(kene->corr_fct[imod][1][3]/kene->corr_fct[imod][1][2]); // } // } // printf("%f %f\n",fct4he/nfct4he,fct3he/nfct3he); // TCanvas *tc = new TCanvas("tc","tc",0,0,1200,1200); // tc->Divide(3,2,0.01,0.01,10); // tc->cd(1); // h4[0]->Draw(); // tc->cd(2); // h4[1]->Draw(); // tc->cd(3); // h4[2]->Draw(); // tc->cd(4); // h3[0]->Draw(); // tc->cd(5); // h3[1]->Draw(); // tc->cd(6); // h3[2]->Draw(); } ////////////////////////////////////////////////////////////////////////