00001 #include "TH1F.h"
00002 #include "TF1.h"
00003 #include "TCanvas.h"
00004 #include "TClonesArray.h"
00005 #include "TMinuit.h"
00006 #include "StandardNtuple/NtpStRecord.h"
00007 #include "CandNtupleSR/NtpSRRecord.h"
00008 #include "CandNtupleSR/NtpSREvent.h"
00009 #include "CandNtupleSR/NtpSRTrack.h"
00010 #include "CandNtupleSR/NtpSRShower.h"
00011 #include "CandNtupleSR/NtpSRStrip.h"
00012 #include "MessageService/MsgService.h"
00013 #include "NueAna/NueAnaTools/SntpHelpers.h"
00014 #include "AnalysisNtuples/Module/ANtpInfoObjectFillerNue.h"
00015 #include "AnalysisNtuples/ANtpNueInfo.h"
00016 #include "AnalysisNtuples/ANtpDefaultValue.h"
00017 #include "VertexFinder/NtpVtxFinder/NtpVtxFinder.h"
00018
00019 CVSID("$Id: ANtpInfoObjectFillerNue.cxx,v 1.3 2008/11/10 15:03:30 boehm Exp $");
00020
00021 static Double_t shwfunc(Double_t *x, Double_t *par)
00022 {
00023
00024
00025 Float_t R = 1.46676;
00026 Float_t xx=R*x[0];
00027
00028
00029 Double_t lnf = TMath::Log(R*par[2]*par[1])+(par[0]-1)*TMath::Log(xx*par[1])-
00030 par[1]*xx-TMath::LnGamma(par[0]);
00031
00032 Double_t f = exp(lnf);
00033
00034
00035
00036
00037 return f;
00038 }
00039
00040
00041
00042 ANtpInfoObjectFillerNue::ANtpInfoObjectFillerNue():
00043 lenepl(0),
00044 tenestu(0),
00045 tenestv(0),
00046 tenestu_9s_2pe_dw(0),
00047 tenestv_9s_2pe_dw(0),
00048 efit(0)
00049 {
00050 }
00051
00052 ANtpInfoObjectFillerNue::~ANtpInfoObjectFillerNue()
00053 {
00054 if(lenepl!=0){ delete lenepl; lenepl=0; }
00055 if(tenestu!=0){ delete tenestu; tenestu=0; }
00056 if(tenestv!=0){ delete tenestv; tenestv=0; }
00057 if(efit!=0){ delete efit; efit=0; }
00058
00059 if(tenestv_9s_2pe_dw!=0){ delete tenestv_9s_2pe_dw; tenestv_9s_2pe_dw=0;}
00060 if(tenestu_9s_2pe_dw!=0){ delete tenestu_9s_2pe_dw; tenestu_9s_2pe_dw=0;}
00061 }
00062
00063
00064
00065 void ANtpInfoObjectFillerNue::Analyze(int evtn, NtpStRecord * st, ANtpNueInfo* NueVars)
00066 {
00067 if(st==0){ return; }
00068 Reset(st->GetHeader().GetSnarl(), evtn, NueVars);
00069
00070 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<"In ANtpInfoObjectFillerNue::Analyze"<<endl;
00071 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<"On Snarl "<<st->GetHeader().GetSnarl()
00072 <<" event "<<evtn<<endl;
00073
00074 NtpSREvent *event = SntpHelpers::GetEvent(evtn,st);
00075
00076 if(!event){
00077 MSG("ANtpInfoObjectFillerNue",Msg::kError)<<"Couldn't get event "<<evtn
00078 <<" from Snarl "<<st->GetHeader().GetSnarl()<<endl;
00079 return;
00080 }
00081
00082 Int_t vtxPlane = event->vtx.plane;
00083 Float_t vtxU = event->vtx.u;
00084 Float_t vtxV = event->vtx.v;
00085
00086 NtpVtxFinder vtxf(evtn, st);
00087 if(vtxf.FindVertex() > 0){
00088 vtxPlane = vtxf.VtxPlane();
00089 vtxU = vtxf.VtxU();
00090 vtxV = vtxf.VtxV();
00091 }
00092
00093 float* evtstp0mip = new float[5500];
00094 float* evtstp1mip = new float[5500];
00095
00096 for(int i = 0; i < 5500; i++){
00097 evtstp0mip[i] = evtstp1mip[i] = -1;
00098 }
00099
00100 FillEventEnergy(evtstp0mip, evtstp1mip, evtn, st, 5500);
00101
00102 for(int i=0;i<event->nstrip;i++){
00103 Int_t index = SntpHelpers::GetStripIndex(i,event);
00104 NtpSRStrip *strip = SntpHelpers::GetStrip(index,st);
00105 if(!strip) continue;
00106
00107 if(!evtstp0mip){
00108 MSG("ANtpInfoObjectFillerNue",Msg::kError)<<"No mip strip information"<<endl;
00109 continue;
00110 }
00111
00112 Float_t stripPh = evtstp0mip[index] + evtstp1mip[index];
00113 double strippe = strip->ph0.pe+strip->ph1.pe;
00114
00115 Float_t deltaplanes = strip->plane-event->vtx.plane;
00116 if(deltaplanes>=0.)
00117 lenepl->Fill(deltaplanes-0.5,stripPh);
00118
00119 const Float_t STRIPWIDTH=0.041;
00120
00121
00122 if(strip->planeview==PlaneView::kU){
00123 double dist = (strip->tpos-event->vtx.u)/STRIPWIDTH;
00124
00125 tenestu->Fill(dist,stripPh);
00126 if(dist < 9.0 && strippe > 2.0){
00127 tenestu_9s_2pe_dw->Fill(dist,stripPh*stripPh);
00128 }
00129 }
00130
00131
00132 else if(strip->planeview==PlaneView::kV){
00133 double dist = (strip->tpos-event->vtx.v)/STRIPWIDTH;
00134
00135 tenestv->Fill(dist,stripPh);
00136 if(dist < 9.0 && strippe > 2.0){
00137 tenestv_9s_2pe_dw->Fill(dist,stripPh*stripPh);
00138 }
00139 }
00140
00141 else{
00142 MSG("ANtpInfoObjectFillerNue",Msg::kError)<<"Don't know what to do with a PlaneView "
00143 <<strip->planeview<<" skipping"<<endl;
00144 continue;
00145 }
00146 }
00147
00148
00149 if(event->plane.n > 4){
00150 FitLShower(event->ph.mip, NueVars);
00151 }
00152 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<"Fit shower "<<NueVars->par_a<<" "
00153 <<NueVars->par_b<<" "<<NueVars->par_e0<<endl;
00154
00155
00156
00157
00158 TransVar(tenestu, PlaneView::kU);
00159 NueVars->u_asym_peak=asym_peak;
00160 NueVars->u_asym_vert=asym_vert;
00161 NueVars->u_molrad_peak=molrad_peak;
00162 NueVars->u_molrad_vert=molrad_vert;
00163 NueVars->u_mean=mean;
00164 NueVars->u_rms=rms;
00165 NueVars->u_skew=skew;
00166 NueVars->u_kurt=kurt;
00167
00168 TransVar(tenestv, PlaneView::kV);
00169 NueVars->v_asym_peak=asym_peak;
00170 NueVars->v_asym_vert=asym_vert;
00171 NueVars->v_molrad_peak=molrad_peak;
00172 NueVars->v_molrad_vert=molrad_vert;
00173 NueVars->v_mean=mean;
00174 NueVars->v_rms=rms;
00175 NueVars->v_skew=skew;
00176 NueVars->v_kurt=kurt;
00177
00178
00179
00180 NueVars->uv_asym_peak = BuildUVVar(NueVars->u_asym_peak, NueVars->v_asym_peak);
00181 NueVars->uv_asym_vert = BuildUVVar(NueVars->u_asym_vert, NueVars->v_asym_vert);
00182 NueVars->uv_molrad_peak = BuildUVVar(NueVars->u_molrad_peak, NueVars->v_molrad_peak);
00183 NueVars->uv_molrad_vert = BuildUVVar(NueVars->u_molrad_vert, NueVars->v_molrad_vert);
00184 NueVars->uv_mean = BuildUVVar(NueVars->u_mean, NueVars->v_mean);
00185 NueVars->uv_rms = BuildUVVar(NueVars->u_rms, NueVars->v_rms);
00186 NueVars->uv_skew = BuildUVVar(NueVars->u_skew, NueVars->v_skew);
00187 NueVars->uv_kurt = BuildUVVar(NueVars->u_kurt, NueVars->v_kurt);
00188
00189 if(tenestu->Integral()>0&&tenestv->Integral()>0){
00190 NueVars->uv_ratio=tenestu->Integral()/tenestv->Integral();
00191 if(NueVars->uv_ratio>100.){
00192 NueVars->uv_ratio= ANtpDefVal::kFloat;
00193 }
00194 }
00195
00196 TransVar(tenestu_9s_2pe_dw, PlaneView::kU);
00197 NueVars->u_asym_peak_9s_2pe_dw=asym_peak;
00198 NueVars->u_asym_vert_9s_2pe_dw=asym_vert;
00199 NueVars->u_molrad_peak_9s_2pe_dw=molrad_peak;
00200 NueVars->u_molrad_vert_9s_2pe_dw=molrad_vert;
00201 NueVars->u_mean_9s_2pe_dw=mean;
00202 NueVars->u_rms_9s_2pe_dw=rms;
00203 NueVars->u_skew_9s_2pe_dw=skew;
00204 NueVars->u_kurt_9s_2pe_dw=kurt;
00205
00206
00207 TransVar(tenestv_9s_2pe_dw, PlaneView::kV);
00208 NueVars->v_asym_peak_9s_2pe_dw=asym_peak;
00209 NueVars->v_asym_vert_9s_2pe_dw=asym_vert;
00210 NueVars->v_molrad_peak_9s_2pe_dw=molrad_peak;
00211 NueVars->v_molrad_vert_9s_2pe_dw=molrad_vert;
00212 NueVars->v_mean_9s_2pe_dw=mean;
00213 NueVars->v_rms_9s_2pe_dw=rms;
00214 NueVars->v_skew_9s_2pe_dw=skew;
00215 NueVars->v_kurt_9s_2pe_dw=kurt;
00216
00217
00218 NueVars->uv_asym_peak_9s_2pe_dw = BuildUVVar(NueVars->u_asym_peak_9s_2pe_dw, NueVars->v_asym_peak_9s_2pe_dw);
00219 NueVars->uv_asym_vert_9s_2pe_dw = BuildUVVar(NueVars->u_asym_vert_9s_2pe_dw, NueVars->v_asym_vert_9s_2pe_dw);
00220 NueVars->uv_molrad_peak_9s_2pe_dw = BuildUVVar(NueVars->u_molrad_peak_9s_2pe_dw, NueVars->v_molrad_peak_9s_2pe_dw);
00221 NueVars->uv_molrad_vert_9s_2pe_dw = BuildUVVar(NueVars->u_molrad_vert_9s_2pe_dw, NueVars->v_molrad_vert_9s_2pe_dw);
00222 NueVars->uv_mean_9s_2pe_dw = BuildUVVar(NueVars->u_mean_9s_2pe_dw, NueVars->v_mean_9s_2pe_dw);
00223 NueVars->uv_rms_9s_2pe_dw = BuildUVVar(NueVars->u_rms_9s_2pe_dw, NueVars->v_rms_9s_2pe_dw);
00224 NueVars->uv_skew_9s_2pe_dw = BuildUVVar(NueVars->u_skew_9s_2pe_dw, NueVars->v_skew_9s_2pe_dw);
00225 NueVars->uv_kurt_9s_2pe_dw = BuildUVVar(NueVars->u_kurt_9s_2pe_dw, NueVars->v_kurt_9s_2pe_dw);
00226
00227
00228 if(tenestu_9s_2pe_dw->Integral()>0&&tenestv_9s_2pe_dw->Integral()>0){
00229 NueVars->uv_ratio_9s_2pe_dw=tenestu_9s_2pe_dw->Integral()/tenestv_9s_2pe_dw->Integral();
00230 if(NueVars->uv_ratio_9s_2pe_dw>100.){
00231 NueVars->uv_ratio_9s_2pe_dw= ANtpDefVal::kFloat;
00232 }
00233 }
00234
00235
00236 if (event->plane.beg>event->plane.end) return;
00237
00238 Float_t plane[14];
00239 for (int i = 0; i<14; i++){
00240 plane[i] = 0;
00241 }
00242
00243 vector<Float_t> stpph;
00244 vector<Float_t> stpph2;
00245 vector<Float_t>::iterator p;
00246 vector<Int_t> stppl;
00247 vector<Int_t> stpstp;
00248 vector<Float_t> stpt;
00249 vector<Float_t> stpz;
00250 vector<Int_t> planev;
00251 Float_t total_ph = 0;
00252 Int_t vtxpl = event->vtx.plane;
00253
00254
00255 for (Int_t istp = 0; istp<event->nstrip; istp++){
00256 Int_t index = event->stp[istp];
00257 NtpSRStrip *strip = SntpHelpers::GetStrip(index,st);
00258 if(!strip) continue;
00259 float charge = strip->ph0.sigcor+strip->ph1.sigcor;
00260
00261
00262
00263 stpph.push_back(charge);
00264 stpph2.push_back(charge);
00265 stppl.push_back(strip->plane);
00266 stpstp.push_back(strip->strip);
00267 stpt.push_back(strip->tpos);
00268 stpz.push_back(strip->z);
00269 planev.push_back(int(strip->planeview));
00270 total_ph += charge;
00271 }
00272
00273 for(unsigned int istp=0; istp<stpph.size(); istp++){
00274 if (stppl[istp]>=vtxpl-2&&stppl[istp]<vtxpl+12)
00275 plane[stppl[istp]-vtxpl+2] += stpph[istp];
00276 }
00277
00278 sort(stpph2.begin(), stpph2.end());
00279
00280 Float_t sum_1_plane = 0, sum_2_planes = 0,sum_3_planes = 0;
00281 Float_t sum_4_planes = 0, sum_5_planes = 0,sum_6_planes = 0;
00282 Float_t sum_2_counts = 0, sum_4_counts = 0, sum_6_counts = 0;
00283 Float_t sum_8_counts = 0, sum_10_counts = 0, sum_12_counts = 0;
00284
00285
00286 if (total_ph>0.){
00287 for (Int_t i = 0; i<14; i++){
00288 if (i<14){
00289 sum_1_plane = plane[i];
00290 if (sum_1_plane/total_ph > NueVars->fract_1_plane) {
00291 NueVars->fract_1_plane = sum_1_plane/total_ph;
00292 }
00293 }
00294 if (i<13){
00295 sum_2_planes = plane[i]+plane[i+1];
00296 if (sum_2_planes/total_ph > NueVars->fract_2_planes)
00297 NueVars->fract_2_planes = sum_2_planes/total_ph;
00298 }
00299 if (i<12){
00300 sum_3_planes = plane[i]+plane[i+1]+plane[i+2];
00301 if (sum_3_planes/total_ph > NueVars->fract_3_planes)
00302 NueVars->fract_3_planes = sum_3_planes/total_ph;
00303 }
00304 if (i<11){
00305 sum_4_planes = plane[i]+plane[i+1]+plane[i+2]+plane[i+3];
00306 if (sum_4_planes/total_ph > NueVars->fract_4_planes)
00307 NueVars->fract_4_planes = sum_4_planes/total_ph;
00308 }
00309 if (i<10) {
00310 sum_5_planes = plane[i]+plane[i+1]+plane[i+2]+plane[i+3]+plane[i+4];
00311 if (sum_5_planes/total_ph > NueVars->fract_5_planes)
00312 NueVars->fract_5_planes = sum_5_planes/total_ph;
00313 }
00314 if (i<9) {
00315 sum_6_planes = plane[i]+plane[i+1]+plane[i+2]+plane[i+3]+plane[i+4]+plane[i+5];
00316 if (sum_6_planes/total_ph > NueVars->fract_6_planes)
00317 NueVars->fract_6_planes = sum_6_planes/total_ph;
00318 }
00319 }
00320 }
00321
00322 p = stpph2.end();
00323 while (p!=stpph2.begin()&&p>stpph2.end()-12){
00324 p--;
00325 if (p>=stpph2.end()-2) sum_2_counts += *p;
00326 if (p>=stpph2.end()-4) sum_4_counts += *p;
00327 if (p>=stpph2.end()-6) sum_6_counts += *p;
00328 if (p>=stpph2.end()-8) sum_8_counts += *p;
00329 if (p>=stpph2.end()-10) sum_10_counts += *p;
00330 if (p>=stpph2.end()-12) sum_12_counts += *p;
00331 }
00332
00333 if (total_ph>0) {
00334 NueVars->fract_2_counters = sum_2_counts/total_ph;
00335 NueVars->fract_4_counters = sum_4_counts/total_ph;
00336 NueVars->fract_6_counters = sum_6_counts/total_ph;
00337 NueVars->fract_8_counters = sum_8_counts/total_ph;
00338 NueVars->fract_10_counters = sum_10_counts/total_ph;
00339 NueVars->fract_12_counters = sum_12_counts/total_ph;
00340 }
00341
00342 delete [] evtstp0mip;
00343 delete [] evtstp1mip;
00344
00345
00346 }
00347
00348
00349 void ANtpInfoObjectFillerNue::FitLShower(Float_t pulseheight, ANtpNueInfo *NueVars)
00350 {
00351 efit->SetParameters(3.,0.5,pulseheight);
00352 lenepl->Fit(efit,"RLQ0+");
00353
00354
00355
00356 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<" STATUS: "<<gMinuit->fCstatu
00357 <<" "<<efit->GetParameter(0)
00358 <<" "<<efit->GetNDF()<<" "<<pulseheight<<endl;
00359
00360 string fitstatus = (string)(gMinuit->fCstatu);
00361 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<" STATUS: "<<fitstatus
00362 <<", "<<efit->GetParameter(0)
00363 <<", "<<efit->GetNDF()<<" "<<pulseheight<<endl;
00364
00365 if(fitstatus=="CONVERGED "&&efit->GetParameter(0)<29.9&&
00366 efit->GetNDF()>0&&pulseheight>0){
00367 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<" filling vars"<<endl;
00368
00369 NueVars->par_a= efit->GetParameter(0);
00370 NueVars->par_b= efit->GetParameter(1);
00371 NueVars->par_e0 = efit->GetParameter(2);
00372 NueVars->chisq= efit->GetChisquare();
00373 NueVars->shwmax=1.*GetMaximumX(efit);
00374
00375
00376
00377 NueVars->shwmaxplane=(int)(lenepl->
00378 GetBinContent(lenepl->FindBin(NueVars->shwmax)));
00379 NueVars->conv=1;
00380 if(efit->GetNDF()>0){
00381 NueVars->chisq_ndf=NueVars->chisq/efit->GetNDF();
00382 }
00383 if(pulseheight!=0){
00384 NueVars->e0_pe_ratio=NueVars->par_e0/pulseheight;
00385 }
00386 NueVars->caldet_comp=ANtpDefVal::kFloat;
00387 NueVars->max_pe_plane=lenepl->GetMaximumBin();
00388 NueVars->shwmaxplane_diff=NueVars->shwmaxplane-NueVars->max_pe_plane;
00389 }
00390
00391 return;
00392 }
00393
00394 float ANtpInfoObjectFillerNue::BuildUVVar(float u, float v)
00395 {
00396 float uv = ANtpDefVal::kFloat;
00397 if(!ANtpDefVal::IsDefault(u)&& !ANtpDefVal::IsDefault(v)&&
00398 u<1.e15 && v<1.e15){
00399 uv=sqrt(u*u+v*v)/sqrt(2.0);
00400 }
00401
00402 return uv;
00403 }
00404
00405
00406 void ANtpInfoObjectFillerNue::TransVar(TH1F *histo, PlaneView::EPlaneView pv)
00407 {
00408 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<"In ANtpInfoObjectFillerNue::TransVar"<<endl;
00409 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<"plane view is "<<pv<<endl;
00410 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<"Hist name is "<<histo->GetName()<<endl;
00411 MSG("ANtpInfoObjectFillerNue",Msg::kDebug)<<"Hist has "<<histo->GetEntries()<<" entries"<<endl;
00412
00413 Int_t THISTBINS = (Int_t)(histo->GetNbinsX());
00414
00415 Int_t binmax=histo->GetMaximumBin();
00416 Int_t binvert=Int_t(((THISTBINS-1)/2)+1);
00417 Float_t peak_bin=histo->Integral(binmax,binmax);
00418 Float_t vert_bin=histo->Integral(binvert,binvert);
00419 Float_t tot=histo->Integral(1,THISTBINS);
00420
00421 molrad_peak = molrad_vert = 0;
00422 mean = rms = skew = kurt = 0;
00423
00424 asym_peak= ANtpDefVal::kFloat;
00425 asym_vert= ANtpDefVal::kFloat;
00426
00427 if(tot-peak_bin){
00428 asym_peak=(TMath::Abs(histo->Integral(binmax+1,THISTBINS)-histo->Integral(1,binmax-1)))
00429 /(tot-peak_bin);
00430 }
00431
00432 if(tot-vert_bin){
00433 asym_vert=(TMath::Abs(histo->Integral(binvert+1,THISTBINS)-histo->Integral(1,binvert-1)))
00434 /(tot-vert_bin);
00435 }
00436
00437 Float_t ratio;
00438
00439 if(tot){
00440 for(Int_t i=0; i<=binvert-1;i++){
00441 ratio=histo->Integral(binmax-i>0?binmax-i:1,
00442 binmax+i<THISTBINS?binmax+i:THISTBINS)/tot;
00443 if(ratio>0.90){molrad_peak=i+1; break;}
00444 }
00445 for(Int_t i=0; i<=binvert-1;i++){
00446 ratio=histo->Integral(binvert-i,binvert+i)/tot;
00447 if(ratio>0.90){molrad_vert=i+1; break;}
00448 }
00449
00450 }
00451
00452 mean=histo->GetMean();
00453 rms=histo->GetRMS();
00454 Int_t n_count=0;
00455
00456 for(Int_t i=1;i<=THISTBINS;i++){
00457
00458 if(histo->GetBinContent(i)){
00459 skew=skew+(TMath::Power((histo->GetBinCenter(i)-mean),3)
00460 *histo->GetBinContent(i));
00461 kurt=kurt+(TMath::Power((histo->GetBinCenter(i)-mean),4)
00462 *histo->GetBinContent(i));
00463 n_count++;
00464 }
00465
00466 }
00467 if(rms>0 && n_count>1){
00468 skew=skew/((Float_t)(n_count-1)*TMath::Power((rms),3));
00469 kurt=(kurt/((Float_t)(n_count-1)*TMath::Power((rms),4)))-3;
00470 }
00471 else{skew= ANtpDefVal::kFloat; kurt= ANtpDefVal::kFloat;}
00472
00473
00474 }
00475
00476 void ANtpInfoObjectFillerNue::Reset(int snarl, int event, ANtpNueInfo *NueVars)
00477 {
00478
00479
00480
00481
00482 if(lenepl!=0){ delete lenepl; lenepl=0; }
00483 if(tenestu!=0){ delete tenestu; tenestu=0; }
00484 if(tenestv!=0){ delete tenestv; tenestv=0; }
00485 if(efit!=0){ delete efit; efit=0; }
00486
00487 if(tenestv_9s_2pe_dw!=0){ delete tenestv_9s_2pe_dw; tenestv_9s_2pe_dw=0;}
00488 if(tenestu_9s_2pe_dw!=0){ delete tenestu_9s_2pe_dw; tenestu_9s_2pe_dw=0;}
00489
00490
00491 const Int_t LHISTBINS=30;
00492 const Int_t THISTBINS=41;
00493 NueVars->Reset();
00494
00495 char ln[100];
00496 char tun[100];
00497 char tvn[100];
00498 sprintf(ln,"lenepl_%d_%d",snarl,event);
00499 sprintf(tun,"tenestu_%d_%d",snarl,event);
00500 sprintf(tvn,"tenestv_%d_%d",snarl,event);
00501 lenepl = new TH1F(ln,"longitudinal energy by plane",
00502 LHISTBINS,0.0,LHISTBINS);
00503
00504 tenestu = new TH1F(tun,"trasverse energy by strip (U view)",
00505 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00506 tenestv = new TH1F(tvn,"trasverse energy by strip (V view)",
00507 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00508
00509 char tun_9s_2pe_dw[100];
00510 char tvn_9s_2pe_dw[100];
00511
00512 sprintf(tun_9s_2pe_dw,"tenestu_9s_2pe_dw_%d_%d",snarl,event);
00513 sprintf(tvn_9s_2pe_dw,"tenestv_9s_2pe_dw_%d_%d",snarl,event);
00514
00515 tenestu_9s_2pe_dw = new TH1F(tun_9s_2pe_dw,"trasverse energy by strip (U view)",
00516 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00517 tenestv_9s_2pe_dw = new TH1F(tvn_9s_2pe_dw,"trasverse energy by strip (V view)",
00518 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00519
00520
00521
00522 int hmin=0;
00523 int hmax=30;
00524 int npare=3;
00525 char efn[100];
00526 sprintf(efn,"efit_%d_%d",snarl,event);
00527 efit = new TF1(efn,shwfunc,hmin+0.001,hmax,npare);
00528 efit->SetParNames("a","b","e0");
00529 efit->SetParLimits(0,hmin+0.001,hmax);
00530 efit->SetParLimits(1,0.001,20000);
00531 efit->SetParLimits(2,0+0.001,1000000);
00532
00533 }
00534
00535
00536 Double_t ANtpInfoObjectFillerNue::GetMaximumX(TF1* efit, Double_t xmin, Double_t xmax)
00537 {
00538 Double_t fXmin = 0.001;
00539 Double_t fXmax = 30;
00540 Double_t fNpx = 30;
00541
00542 Double_t x,y;
00543 if (xmin >= xmax) {xmin = fXmin; xmax = fXmax;}
00544 Double_t dx = (xmax-xmin)/fNpx;
00545 Double_t xxmax = xmin;
00546 Double_t yymax = efit->Eval(xmin+dx);
00547 for (Int_t i=0;i<fNpx;i++) {
00548 x = xmin + (i+0.5)*dx;
00549 y = efit->Eval(x);
00550 if (y > yymax) {xxmax = x; yymax = y;}
00551 }
00552 if (dx < 1.e-9*(fXmax-fXmin)) return TMath::Min(xmax,xxmax);
00553 else return GetMaximumX(efit, TMath::Max(xmin,xxmax-dx), TMath::Min(xmax,xxmax+dx));
00554 }
00555
00556
00557
00558 void ANtpInfoObjectFillerNue::FillEventEnergy(float* ph0, float* ph1, int evtn, NtpStRecord* st, const int SIZE)
00559 {
00560 NtpSREvent *event = 0;
00561 event = SntpHelpers::GetEvent(evtn,st);
00562 if(event == 0){
00563 cout<<"No event"<<endl;
00564 return;
00565 }
00566
00567 for(int j=0;j<event->ntrack;j++){
00568 int tindex = SntpHelpers::GetTrackIndex(j,event);
00569 NtpSRTrack *track = SntpHelpers::GetTrack(tindex,st);
00570 for(int k = 0; k < track->nstrip; k++){
00571 int index = SntpHelpers::GetStripIndex(k, track);
00572 if(index > SIZE){
00573 MSG("SntpHelper",Msg::kFatal)<<" evt mip array insufficiently large: "
00574 <<"index "<<index<<"/"<<SIZE<<" requested"<<std::endl;
00575 }
00576
00577 float mips1 = track->stpph1mip[k];
00578 float mips0 = track->stpph0mip[k];
00579
00580 if(ph0[index] < 0) ph0[index] = mips0;
00581 if(ph1[index] < 0) ph1[index] = mips1;
00582 }
00583 }
00584 for(int j=0;j<event->nshower;j++){
00585 int sindex = SntpHelpers::GetShowerIndex(j,event);
00586 NtpSRShower *shw = SntpHelpers::GetShower(sindex,st);
00587 for(int k = 0; k < shw->nstrip; k++){
00588 int index = SntpHelpers::GetStripIndex(k, shw);
00589 if(index > SIZE){
00590 MSG("SntpHelper",Msg::kFatal)<<" evt mip array insufficiently large: "
00591 <<"index "<<index<<"/"<<SIZE<<" requested"<<std::endl;
00592 }
00593
00594 float mips1 = shw->stpph1mip[k];
00595 float mips0 = shw->stpph0mip[k];
00596
00597 if(ph0[index] < 0) ph0[index] = mips0;
00598 if(ph1[index] < 0) ph1[index] = mips1;
00599 }
00600 }
00601 return;
00602 }
00603
