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/NtpSRStrip.h"
00011 #include "MessageService/MsgService.h"
00012 #include "NueAna/NueAnaTools/SntpHelpers.h"
00013 #include "NueAna/ShwfitAna.h"
00014 #include "AnalysisNtuples/ANtpDefaultValue.h"
00015 #include "VertexFinder/NtpVtxFinder/NtpVtxFinder.h"
00016
00017 CVSID("$Id: ShwfitAna.cxx,v 1.38 2008/11/19 18:22:51 rhatcher Exp $");
00018
00019 static Double_t shwfunc(Double_t *x, Double_t *par)
00020 {
00021
00022
00023 Float_t R = 1.46676;
00024 Float_t xx=R*x[0];
00025
00026
00027 Double_t lnf = TMath::Log(R*par[2]*par[1])+(par[0]-1)*TMath::Log(xx*par[1])-
00028 par[1]*xx-TMath::LnGamma(par[0]);
00029
00030 Double_t f = exp(lnf);
00031
00032
00033
00034
00035 return f;
00036 }
00037
00038 static Double_t hadfunc(Double_t *x, Double_t *par)
00039 {
00040
00041
00042 Float_t xx=x[0];
00043 Double_t f1 = par[2]*TMath::Power(par[1],par[0])*
00044 TMath::Power(xx*1.46676,par[0]-1)*
00045 TMath::Exp(-par[1]*(xx*1.46676))/TMath::Gamma(par[0]);
00046 Double_t f2 = par[2]*TMath::Power(par[4],par[3])*
00047 TMath::Power(xx*0.1694,par[3]-1)*
00048 TMath::Exp(-par[4]*(xx*0.1694))/TMath::Gamma(par[3]);
00049 Double_t f=(par[5]*f1+(1-par[5])*0.1154*f2)*1.46676;
00050 return f;
00051 }
00052
00053
00054
00055 static Double_t maxwell(Double_t *x, Double_t *par)
00056 {
00057
00058
00059 Float_t xx = x[0];
00060
00061 Double_t f = par[1]*(4/TMath::Sqrt(TMath::Pi()))*TMath::Power(par[0],(3/2))*TMath::Power(xx,2)*TMath::Exp(-par[0]*TMath::Power(xx,2));
00062
00063 return f;
00064 }
00065 static Double_t maxwell3(Double_t *x, Double_t *par)
00066 {
00067
00068
00069 Float_t xx = x[0];
00070
00071 Double_t f = par[1]*2*TMath::Power(par[0],2)*TMath::Power(xx,3)*TMath::Exp(-par[0]*TMath::Power(xx,2));
00072
00073 return f;
00074 }
00075
00076
00077 ShwfitAna::ShwfitAna(Shwfit &sf):
00078 fShwfit(sf)
00079 {}
00080
00081 ShwfitAna::~ShwfitAna()
00082 {}
00083
00084
00085
00086 void ShwfitAna::Analyze(int evtn, RecRecordImp<RecCandHeader> *srobj)
00087 {
00088 if(srobj==0){
00089 return;
00090 }
00091 if(((dynamic_cast<NtpStRecord *>(srobj))==0)&&
00092 ((dynamic_cast<NtpSRRecord *>(srobj))==0)){
00093 return;
00094 }
00095
00096 MSG("ShwfitAna",Msg::kDebug)<<"In ShwfitAna::Analyze"<<endl;
00097 MSG("ShwfitAna",Msg::kDebug)<<"On Snarl "<<srobj->GetHeader().GetSnarl()
00098 <<" event "<<evtn<<endl;
00099
00100 Reset(srobj->GetHeader().GetSnarl(),evtn);
00101
00102 NtpSREvent *event = SntpHelpers::GetEvent(evtn,srobj);
00103
00104 if(!event){
00105 MSG("ShwfitAna",Msg::kError)<<"Couldn't get event "<<evtn
00106 <<" from Snarl "<<srobj->GetHeader().GetSnarl()<<endl;
00107 return;
00108 }
00109
00110 Int_t vtxPlane = event->vtx.plane;
00111 Float_t vtxU = event->vtx.u;
00112 Float_t vtxV = event->vtx.v;
00113
00114 if(ReleaseType::IsCedar(release)){
00115 NtpStRecord* st = dynamic_cast<NtpStRecord *>(srobj);
00116 NtpVtxFinder vtxf(evtn, st);
00117 if(vtxf.FindVertex() > 0){
00118 vtxPlane = vtxf.VtxPlane();
00119 vtxU = vtxf.VtxU();
00120 vtxV = vtxf.VtxV();
00121 }
00122 }
00123
00124
00125 fShwfit.hiPhStripCountM4=0;
00126 fShwfit.hiPhPlaneCountM4=0;
00127 fShwfit.hiPhStripCountM2=0;
00128 fShwfit.hiPhPlaneCountM2=0;
00129 fShwfit.hiPhStripCount=0;
00130 fShwfit.hiPhPlaneCount=0;
00131 fShwfit.hiPhStripCountP2=0;
00132 fShwfit.hiPhPlaneCountP2=0;
00133 fShwfit.hiPhStripCountP4=0;
00134 fShwfit.hiPhPlaneCountP4=0;
00135
00136
00137 fShwfit.contPlaneCount=0;
00138
00139 float VertexEnergy = 0.0;
00140 Float_t vtxEnergy0 = 0.0;
00141 Float_t vtxEnergy1 = 0.0;
00142 Float_t vtxEnergy2 = 0.0;
00143
00144
00145 doSlopes(event,srobj);
00146
00147
00148 for(int i=0;i<event->nstrip;i++){
00149 Int_t index = SntpHelpers::GetStripIndex(i,event);
00150 NtpSRStrip *strip = SntpHelpers::GetStrip(index,srobj);
00151 if(!strip){
00152 continue;
00153 }
00154 if(!evtstp0mip){
00155 MSG("ShwfitAna",Msg::kError)<<"No mip strip information"<<endl;
00156 continue;
00157 }
00158
00159 Float_t deltaplanes = strip->plane-event->vtx.plane;
00160 Float_t stripPh = evtstp0mip[index] + evtstp1mip[index];
00161 double strippe = strip->ph0.pe+strip->ph1.pe;
00162
00163 if(deltaplanes == 0 || deltaplanes == 1) VertexEnergy += stripPh;
00164 if(deltaplanes == 0) vtxEnergy0 += stripPh;
00165 if(deltaplanes == 1) vtxEnergy1 += stripPh;
00166 if(deltaplanes == 2) vtxEnergy2 += stripPh;
00167
00168 const Float_t STRIPWIDTH=0.041;
00169 MSG("ShwfitAna",Msg::kDebug)<< "plane " << deltaplanes << " stripPh " << stripPh <<endl;
00170
00171 if(deltaplanes>=0.){
00172 fShwfit.lenepl->Fill(deltaplanes-0.5,stripPh);
00173 fShwfit.ph_hist->Fill(stripPh);
00174
00175 if(sfDPlaneCut>0&&deltaplanes<sfDPlaneCut-4&&stripPh>sfPhStripCut) fShwfit.hiPhStripCountM4++;
00176 if(sfDPlaneCut>0&&deltaplanes<sfDPlaneCut-2&&stripPh>sfPhStripCut) fShwfit.hiPhStripCountM2++;
00177 if(sfDPlaneCut>0&&deltaplanes<sfDPlaneCut&&stripPh>sfPhStripCut) fShwfit.hiPhStripCount++;
00178 if(sfDPlaneCut>0&&deltaplanes<sfDPlaneCut+2&&stripPh>sfPhStripCut) fShwfit.hiPhStripCountP2++;
00179 if(sfDPlaneCut>0&&deltaplanes<sfDPlaneCut+4&&stripPh>sfPhStripCut) fShwfit.hiPhStripCountP4++;
00180
00181 }
00182
00183
00184
00185 if(strip->planeview==PlaneView::kU){
00186 double dist = (strip->tpos-event->vtx.u)/STRIPWIDTH;
00187
00188 fShwfit.tenestu->Fill(dist,stripPh);
00189 if(dist < 9.0 && strippe > 2.0){
00190 fShwfit.tenestu_9s_2pe->Fill(dist,stripPh);
00191 fShwfit.tenestu_9s_2pe_dw->Fill(dist,stripPh*stripPh);
00192 }
00193 }
00194
00195
00196 else if(strip->planeview==PlaneView::kV){
00197 double dist = (strip->tpos-event->vtx.v)/STRIPWIDTH;
00198
00199 fShwfit.tenestv->Fill(dist,stripPh);
00200 if(dist < 9.0 && strippe > 2.0){
00201 fShwfit.tenestv_9s_2pe->Fill(dist,stripPh);
00202 fShwfit.tenestv_9s_2pe_dw->Fill(dist,stripPh*stripPh);
00203 }
00204 }
00205
00206 else{
00207 MSG("ShwfitAna",Msg::kError)<<"Don't know what to do with a PlaneView "
00208 <<strip->planeview<<" skipping"<<endl;
00209 continue;
00210 }
00211 }
00212
00213
00214 fShwfit.vtxEnergy = VertexEnergy;
00215 fShwfit.energyPlane0 = vtxEnergy0;
00216 fShwfit.energyPlane1 = vtxEnergy1;
00217 fShwfit.energyPlane2 = vtxEnergy2;
00218
00219
00220 int planebins=fShwfit.lenepl->GetNbinsX();
00221
00222 if(sfDPlaneCut<planebins) planebins=sfDPlaneCut;
00223 for(int i=1;i<=planebins+4;i++){
00224 if(fShwfit.lenepl->GetBinContent(i)>sfPhPlaneCut&&i<=planebins-4) fShwfit.hiPhPlaneCountM4++;
00225 if(fShwfit.lenepl->GetBinContent(i)>sfPhPlaneCut&&i<=planebins-2) fShwfit.hiPhPlaneCountM2++;
00226 if(fShwfit.lenepl->GetBinContent(i)>sfPhPlaneCut&&i<=planebins) fShwfit.hiPhPlaneCount++;
00227 if(fShwfit.lenepl->GetBinContent(i)>sfPhPlaneCut&&i<=planebins+2) fShwfit.hiPhPlaneCountP2++;
00228 if(fShwfit.lenepl->GetBinContent(i)>sfPhPlaneCut&&i<=planebins+4) fShwfit.hiPhPlaneCountP4++;
00229 }
00230
00231
00232
00233 if(sfContPhPlaneCut>0){
00234 bool PlaneCut=true;
00235 fShwfit.contPlaneCount=0;
00236 if(fShwfit.lenepl->GetBinContent(1)>sfContPhPlaneCut){
00237 fShwfit.contPlaneCount++;
00238 }
00239 for(int i=2;i<=30; i++){
00240 if(fShwfit.lenepl->GetBinContent(i)>sfContPhPlaneCut&&PlaneCut){
00241 fShwfit.contPlaneCount++;
00242 }else {
00243 PlaneCut=false;
00244 }
00245 }
00246 }
00247
00248
00249 bool PlaneCut015=true;
00250 bool PlaneCut030=true;
00251 bool PlaneCut050=true;
00252 bool PlaneCut075=true;
00253 bool PlaneCut100=true;
00254 bool PlaneCut200=true;
00255
00256 fShwfit.contPlaneCount015=0;
00257 fShwfit.contPlaneCount030=0;
00258 fShwfit.contPlaneCount050=0;
00259 fShwfit.contPlaneCount075=0;
00260 fShwfit.contPlaneCount100=0;
00261 fShwfit.contPlaneCount200=0;
00262
00263 if(fShwfit.lenepl->GetBinContent(1)>0.15) fShwfit.contPlaneCount015++;
00264 if(fShwfit.lenepl->GetBinContent(1)>0.30) fShwfit.contPlaneCount030++;
00265 if(fShwfit.lenepl->GetBinContent(1)>0.50) fShwfit.contPlaneCount050++;
00266 if(fShwfit.lenepl->GetBinContent(1)>0.75) fShwfit.contPlaneCount075++;
00267 if(fShwfit.lenepl->GetBinContent(1)>1.00) fShwfit.contPlaneCount100++;
00268 if(fShwfit.lenepl->GetBinContent(1)>2.00) fShwfit.contPlaneCount200++;
00269
00270 for(int i=2;i<=30; i++){
00271 if(fShwfit.lenepl->GetBinContent(i)>0.15&&PlaneCut015){
00272 fShwfit.contPlaneCount015++;
00273 }else {
00274 PlaneCut015=false;
00275 }
00276 if(fShwfit.lenepl->GetBinContent(i)>0.30&&PlaneCut030){
00277 fShwfit.contPlaneCount030++;
00278 }else {
00279 PlaneCut030=false;
00280 }
00281 if(fShwfit.lenepl->GetBinContent(i)>0.50&&PlaneCut050){
00282 fShwfit.contPlaneCount050++;
00283 }else {
00284 PlaneCut050=false;
00285 }
00286 if(fShwfit.lenepl->GetBinContent(i)>0.75&&PlaneCut075){
00287 fShwfit.contPlaneCount075++;
00288 }else {
00289 PlaneCut075=false;
00290 }
00291 if(fShwfit.lenepl->GetBinContent(i)>1.00&&PlaneCut100){
00292 fShwfit.contPlaneCount100++;
00293 }else {
00294 PlaneCut100=false;
00295 }
00296 if(fShwfit.lenepl->GetBinContent(i)>2.00&&PlaneCut200){
00297 fShwfit.contPlaneCount200++;
00298 }else {
00299 PlaneCut200=false;
00300 }
00301 }
00302
00303
00304 MSG("ShwfitAna",Msg::kDebug)<<"StripCount "<< fShwfit.hiPhStripCount
00305 << " PlaneCount " << fShwfit.hiPhPlaneCount
00306 << " contPlaneCount "<<fShwfit.contPlaneCount<< endl;
00308
00309
00310 MSG("ShwfitAna",Msg::kDebug)<<"In ShwfitAna, trying to fit"<<endl;
00311
00312 Int_t trk_plane_num=0;
00313 Int_t ntrack=event->ntrack;
00314 if(ntrack>0){
00315 Int_t trkNum=SntpHelpers::GetTrackIndex(0,event);
00316 NtpSRTrack *track = SntpHelpers::GetTrack(trkNum,srobj);
00317 trk_plane_num=track->plane.n;
00318 }
00319
00320 if(event->plane.n > 4){
00321 FitLShower(event->ph.mip);
00322 FitLShower_Dan(event->ph.mip);
00323 }
00324 MSG("ShwfitAna",Msg::kDebug)<<"Fit shower "<<fShwfit.par_a<<" "
00325 <<fShwfit.par_b<<" "<<fShwfit.par_e0<<endl;
00326
00327
00328
00329 FitTShower(event->ph.mip);
00330
00331 TransVar(fShwfit.tenestu, PlaneView::kU);
00332 fShwfit.u_asym_peak=asym_peak;
00333 fShwfit.u_asym_vert=asym_vert;
00334 fShwfit.u_molrad_peak=molrad_peak;
00335 fShwfit.u_molrad_vert=molrad_vert;
00336 fShwfit.u_mean=mean;
00337 fShwfit.u_rms=rms;
00338 fShwfit.u_skew=skew;
00339 fShwfit.u_kurt=kurt;
00340
00341 TransVar(fShwfit.tenestv, PlaneView::kV);
00342 fShwfit.v_asym_peak=asym_peak;
00343 fShwfit.v_asym_vert=asym_vert;
00344 fShwfit.v_molrad_peak=molrad_peak;
00345 fShwfit.v_molrad_vert=molrad_vert;
00346 fShwfit.v_mean=mean;
00347 fShwfit.v_rms=rms;
00348 fShwfit.v_skew=skew;
00349 fShwfit.v_kurt=kurt;
00350
00351
00352
00353 fShwfit.uv_asym_peak = BuildUVVar(fShwfit.u_asym_peak, fShwfit.v_asym_peak);
00354 fShwfit.uv_asym_vert = BuildUVVar(fShwfit.u_asym_vert, fShwfit.v_asym_vert);
00355 fShwfit.uv_molrad_peak = BuildUVVar(fShwfit.u_molrad_peak, fShwfit.v_molrad_peak);
00356 fShwfit.uv_molrad_vert = BuildUVVar(fShwfit.u_molrad_vert, fShwfit.v_molrad_vert);
00357 fShwfit.uv_mean = BuildUVVar(fShwfit.u_mean, fShwfit.v_mean);
00358 fShwfit.uv_rms = BuildUVVar(fShwfit.u_rms, fShwfit.v_rms);
00359 fShwfit.uv_skew = BuildUVVar(fShwfit.u_skew, fShwfit.v_skew);
00360 fShwfit.uv_kurt = BuildUVVar(fShwfit.u_kurt, fShwfit.v_kurt);
00361
00362 if(fShwfit.tenestu->Integral()>0&&fShwfit.tenestv->Integral()>0){
00363 fShwfit.uv_ratio=fShwfit.tenestu->Integral()/fShwfit.tenestv->Integral();
00364 if(fShwfit.uv_ratio>100.){
00365 fShwfit.uv_ratio= ANtpDefVal::kFloat;
00366 }
00367 }
00368
00369 TransVar(fShwfit.tenestu_9s_2pe, PlaneView::kU);
00370 fShwfit.u_asym_peak_9s_2pe=asym_peak;
00371 fShwfit.u_asym_vert_9s_2pe=asym_vert;
00372 fShwfit.u_molrad_peak_9s_2pe=molrad_peak;
00373 fShwfit.u_molrad_vert_9s_2pe=molrad_vert;
00374 fShwfit.u_mean_9s_2pe=mean;
00375 fShwfit.u_rms_9s_2pe=rms;
00376 fShwfit.u_skew_9s_2pe=skew;
00377 fShwfit.u_kurt_9s_2pe=kurt;
00378
00379 TransVar(fShwfit.tenestv_9s_2pe, PlaneView::kV);
00380 fShwfit.v_asym_peak_9s_2pe=asym_peak;
00381 fShwfit.v_asym_vert_9s_2pe=asym_vert;
00382 fShwfit.v_molrad_peak_9s_2pe=molrad_peak;
00383 fShwfit.v_molrad_vert_9s_2pe=molrad_vert;
00384 fShwfit.v_mean_9s_2pe=mean;
00385 fShwfit.v_rms_9s_2pe=rms;
00386 fShwfit.v_skew_9s_2pe=skew;
00387 fShwfit.v_kurt_9s_2pe=kurt;
00388
00389 fShwfit.uv_asym_peak_9s_2pe = BuildUVVar(fShwfit.u_asym_peak_9s_2pe, fShwfit.v_asym_peak_9s_2pe);
00390 fShwfit.uv_asym_vert_9s_2pe = BuildUVVar(fShwfit.u_asym_vert_9s_2pe, fShwfit.v_asym_vert_9s_2pe);
00391 fShwfit.uv_molrad_peak_9s_2pe = BuildUVVar(fShwfit.u_molrad_peak_9s_2pe, fShwfit.v_molrad_peak_9s_2pe);
00392 fShwfit.uv_molrad_vert_9s_2pe = BuildUVVar(fShwfit.u_molrad_vert_9s_2pe, fShwfit.v_molrad_vert_9s_2pe);
00393 fShwfit.uv_mean_9s_2pe = BuildUVVar(fShwfit.u_mean_9s_2pe, fShwfit.v_mean_9s_2pe);
00394 fShwfit.uv_rms_9s_2pe = BuildUVVar(fShwfit.u_rms_9s_2pe, fShwfit.v_rms_9s_2pe);
00395 fShwfit.uv_skew_9s_2pe = BuildUVVar(fShwfit.u_skew_9s_2pe, fShwfit.v_skew_9s_2pe);
00396 fShwfit.uv_kurt_9s_2pe = BuildUVVar(fShwfit.u_kurt_9s_2pe, fShwfit.v_kurt_9s_2pe);
00397 if(fShwfit.tenestu_9s_2pe->Integral()>0&&fShwfit.tenestv_9s_2pe->Integral()>0){
00398 fShwfit.uv_ratio_9s_2pe=fShwfit.tenestu_9s_2pe->Integral()/fShwfit.tenestv_9s_2pe->Integral();
00399 if(fShwfit.uv_ratio_9s_2pe>100.){
00400 fShwfit.uv_ratio_9s_2pe= ANtpDefVal::kFloat;
00401 }
00402 }
00403
00404 TransVar(fShwfit.tenestu_9s_2pe_dw, PlaneView::kU);
00405 fShwfit.u_asym_peak_9s_2pe_dw=asym_peak;
00406 fShwfit.u_asym_vert_9s_2pe_dw=asym_vert;
00407 fShwfit.u_molrad_peak_9s_2pe_dw=molrad_peak;
00408 fShwfit.u_molrad_vert_9s_2pe_dw=molrad_vert;
00409 fShwfit.u_mean_9s_2pe_dw=mean;
00410 fShwfit.u_rms_9s_2pe_dw=rms;
00411 fShwfit.u_skew_9s_2pe_dw=skew;
00412 fShwfit.u_kurt_9s_2pe_dw=kurt;
00413
00414
00415 TransVar(fShwfit.tenestv_9s_2pe_dw, PlaneView::kV);
00416 fShwfit.v_asym_peak_9s_2pe_dw=asym_peak;
00417 fShwfit.v_asym_vert_9s_2pe_dw=asym_vert;
00418 fShwfit.v_molrad_peak_9s_2pe_dw=molrad_peak;
00419 fShwfit.v_molrad_vert_9s_2pe_dw=molrad_vert;
00420 fShwfit.v_mean_9s_2pe_dw=mean;
00421 fShwfit.v_rms_9s_2pe_dw=rms;
00422 fShwfit.v_skew_9s_2pe_dw=skew;
00423 fShwfit.v_kurt_9s_2pe_dw=kurt;
00424
00425
00426 fShwfit.uv_asym_peak_9s_2pe_dw = BuildUVVar(fShwfit.u_asym_peak_9s_2pe_dw, fShwfit.v_asym_peak_9s_2pe_dw);
00427 fShwfit.uv_asym_vert_9s_2pe_dw = BuildUVVar(fShwfit.u_asym_vert_9s_2pe_dw, fShwfit.v_asym_vert_9s_2pe_dw);
00428 fShwfit.uv_molrad_peak_9s_2pe_dw = BuildUVVar(fShwfit.u_molrad_peak_9s_2pe_dw, fShwfit.v_molrad_peak_9s_2pe_dw);
00429 fShwfit.uv_molrad_vert_9s_2pe_dw = BuildUVVar(fShwfit.u_molrad_vert_9s_2pe_dw, fShwfit.v_molrad_vert_9s_2pe_dw);
00430 fShwfit.uv_mean_9s_2pe_dw = BuildUVVar(fShwfit.u_mean_9s_2pe_dw, fShwfit.v_mean_9s_2pe_dw);
00431 fShwfit.uv_rms_9s_2pe_dw = BuildUVVar(fShwfit.u_rms_9s_2pe_dw, fShwfit.v_rms_9s_2pe_dw);
00432 fShwfit.uv_skew_9s_2pe_dw = BuildUVVar(fShwfit.u_skew_9s_2pe_dw, fShwfit.v_skew_9s_2pe_dw);
00433 fShwfit.uv_kurt_9s_2pe_dw = BuildUVVar(fShwfit.u_kurt_9s_2pe_dw, fShwfit.v_kurt_9s_2pe_dw);
00434
00435
00436 if(fShwfit.tenestu_9s_2pe_dw->Integral()>0&&fShwfit.tenestv_9s_2pe_dw->Integral()>0){
00437 fShwfit.uv_ratio_9s_2pe_dw=fShwfit.tenestu_9s_2pe_dw->Integral()/fShwfit.tenestv_9s_2pe_dw->Integral();
00438 if(fShwfit.uv_ratio_9s_2pe_dw>100.){
00439 fShwfit.uv_ratio_9s_2pe_dw= ANtpDefVal::kFloat;
00440 }
00441 }
00442 }
00443
00444
00445 void ShwfitAna::FitLShower(Float_t pulseheight)
00446 {
00447 fShwfit.efit->SetParameters(3.,0.5,pulseheight);
00448 fShwfit.lenepl->Fit(fShwfit.efit,"RLLQ0+");
00449
00450
00451
00452 MSG("ShwfitAna",Msg::kDebug)<<" STATUS: "<<gMinuit->fCstatu
00453 <<" "<<fShwfit.efit->GetParameter(0)
00454 <<" "<<fShwfit.efit->GetNDF()<<" "<<pulseheight<<endl;
00455
00456 string fitstatus = (string)(gMinuit->fCstatu);
00457 MSG("ShwfitAna",Msg::kDebug)<<" STATUS: "<<fitstatus
00458 <<", "<<fShwfit.efit->GetParameter(0)
00459 <<", "<<fShwfit.efit->GetNDF()<<" "<<pulseheight<<endl;
00460
00461 if(fitstatus=="CONVERGED "&&fShwfit.efit->GetParameter(0)<29.9&&
00462 fShwfit.efit->GetNDF()>0&&pulseheight>0){
00463 MSG("ShwfitAna",Msg::kDebug)<<" filling vars"<<endl;
00464
00465 fShwfit.par_a=fShwfit.efit->GetParameter(0);
00466 fShwfit.par_b=fShwfit.efit->GetParameter(1);
00467 fShwfit.par_e0=fShwfit.efit->GetParameter(2);
00468 fShwfit.chisq=fShwfit.efit->GetChisquare();
00469 fShwfit.shwmax=1.*GetMaximumX(fShwfit.efit);
00470
00471
00472
00473 fShwfit.shwmaxplane=(int)(fShwfit.lenepl->
00474 GetBinContent(fShwfit.lenepl->FindBin(fShwfit.shwmax)));
00475 fShwfit.conv=1;
00476 if(fShwfit.efit->GetNDF()>0){
00477 fShwfit.chisq_ndf=fShwfit.chisq/fShwfit.efit->GetNDF();
00478 }
00479 if(pulseheight!=0){
00480 fShwfit.e0_pe_ratio=fShwfit.par_e0/pulseheight;
00481 }
00482 fShwfit.caldet_comp=ANtpDefVal::kFloat;
00483 fShwfit.max_pe_plane=fShwfit.lenepl->GetMaximumBin();
00484 fShwfit.shwmaxplane_diff=fShwfit.shwmaxplane-fShwfit.max_pe_plane;
00485 }
00486
00487 return;
00488 }
00489
00490 void ShwfitAna::FitLShower_Dan(Float_t pulseheight)
00491 {
00492
00493 fShwfit.efit_maxwell->SetParameters(1.0,pulseheight);
00494 fShwfit.lenepl->Fit(fShwfit.efit_maxwell,"RLQ0+");
00495 fShwfit.Beta_Maxwell = fShwfit.efit_maxwell->GetParameter(0);
00496 int np = 1000;
00497 double *x = new double[np];
00498 double *w = new double[np];
00499 fShwfit.efit_maxwell->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
00500 fShwfit.Energy_Maxwell = fShwfit.efit_maxwell->IntegralFast(np,x,w,0,1000);
00501 fShwfit.chisq_Maxwell=fShwfit.efit_maxwell->GetChisquare();
00502 fShwfit.ndf_Maxwell=fShwfit.efit_maxwell->GetNDF();
00503
00504 fShwfit.efit_maxwell3->SetParameters(1.0,pulseheight);
00505 fShwfit.lenepl->Fit(fShwfit.efit_maxwell3,"RLQ0+");
00506 fShwfit.Beta_Maxwell3 = fShwfit.efit_maxwell3->GetParameter(0);
00507 int np3 = 1000;
00508 double *x3 = new double[np3];
00509 double *w3 = new double[np3];
00510 fShwfit.efit_maxwell3->CalcGaussLegendreSamplingPoints(np3,x3,w3,1e-15);
00511 fShwfit.Energy_Maxwell3 = fShwfit.efit_maxwell3->IntegralFast(np3,x3,w3,0,1000);
00512 fShwfit.chisq_Maxwell3=fShwfit.efit_maxwell3->GetChisquare();
00513 fShwfit.ndf_Maxwell3=fShwfit.efit_maxwell3->GetNDF();
00514
00515
00516 float E_half_num = 0;
00517 float E_half_den = 0;
00518 float n_half_num = 0;
00519 float n_half_den = 0;
00520 float E_2_num = 0;
00521 float E_2_den = 0;
00522 float n_2_num = 0;
00523 float n_2_den = 0;
00524 float E_split_val = 0;
00525 float E_split_point = 0;
00526 float n_split_val = 0;
00527 float n_split_point = 0;
00528 float n_integral = 0;
00529 float E_integral = 0;
00530 float ph_hist_length = 0;
00531 float ph_hist_counter = 0;
00532 for (int q = 1; q<=fShwfit.ph_hist->GetNbinsX(); q++){
00533 n_integral += fShwfit.ph_hist->GetBinContent(q);
00534 E_integral += fShwfit.ph_hist->GetBinContent(q)*(fShwfit.ph_hist->GetBinLowEdge(q)+(1.0/2.0)*fShwfit.ph_hist->GetBinWidth(q));
00535 ph_hist_counter++;
00536 if(fShwfit.ph_hist->GetBinContent(q)>0){
00537 ph_hist_length += ph_hist_counter;
00538 ph_hist_counter = 0;
00539 }
00540 }
00541 for (int q = 1; q<=fShwfit.ph_hist->GetNbinsX(); q++){
00542 if (q<=2){
00543 n_2_num += fShwfit.ph_hist->GetBinContent(q);
00544 E_2_num += fShwfit.ph_hist->GetBinContent(q)*(fShwfit.ph_hist->GetBinLowEdge(q)+(1.0/2.0)*fShwfit.ph_hist->GetBinWidth(q));
00545 }else{
00546 n_2_den += fShwfit.ph_hist->GetBinContent(q);
00547 E_2_den += fShwfit.ph_hist->GetBinContent(q)*(fShwfit.ph_hist->GetBinLowEdge(q)+(1.0/2.0)*fShwfit.ph_hist->GetBinWidth(q));
00548 }
00549 if (q<=ph_hist_length/2){
00550 n_half_num += fShwfit.ph_hist->GetBinContent(q);
00551 E_half_num += fShwfit.ph_hist->GetBinContent(q)*(fShwfit.ph_hist->GetBinLowEdge(q)+(1.0/2.0)*fShwfit.ph_hist->GetBinWidth(q));
00552 }else{
00553 n_half_den += fShwfit.ph_hist->GetBinContent(q);
00554 E_half_den += fShwfit.ph_hist->GetBinContent(q)*(fShwfit.ph_hist->GetBinLowEdge(q)+(1.0/2.0)*fShwfit.ph_hist->GetBinWidth(q));
00555 }
00556 if(n_split_val<(1.0/2.0)*n_integral) n_split_val += fShwfit.ph_hist->GetBinContent(q);
00557 else if (n_split_point==0) n_split_point = (fShwfit.ph_hist->GetBinLowEdge(q)+fShwfit.ph_hist->GetBinWidth(q));
00558 if(E_split_val<(1.0/2.0)*E_integral) E_split_val += fShwfit.ph_hist->GetBinContent(q)*(fShwfit.ph_hist->GetBinLowEdge(q)+(1.0/2.0)*fShwfit.ph_hist->GetBinWidth(q));
00559 else if (E_split_point==0) E_split_point = (fShwfit.ph_hist->GetBinLowEdge(q)+fShwfit.ph_hist->GetBinWidth(q));
00560
00561
00562 }
00563 if(E_half_den!=0) fShwfit.E_ratio_half = E_half_num/E_half_den;
00564 else fShwfit.E_ratio_half = 0;
00565 if(n_half_den!=0) fShwfit.n_ratio_half = n_half_num/n_half_den;
00566 else fShwfit.n_ratio_half = 0;
00567 if(E_2_den!=0) fShwfit.E_ratio_2 = E_2_num/E_2_den;
00568 else fShwfit.E_ratio_2 = 0;
00569 if(n_2_den!=0) fShwfit.n_ratio_2 = n_2_num/n_2_den;
00570 else fShwfit.n_ratio_2 = 0;
00571 fShwfit.pos_E_split = E_split_point;
00572 fShwfit.pos_n_split = n_split_point;
00573
00574
00575 delete [] x;
00576 delete [] x3;
00577 delete [] w;
00578 delete [] w3;
00579
00580 return;
00581 }
00582
00583
00584 float ShwfitAna::BuildUVVar(float u, float v)
00585 {
00586 float uv = ANtpDefVal::kFloat;
00587 if(!ANtpDefVal::IsDefault(u)&& !ANtpDefVal::IsDefault(v)&&
00588 u<1.e15 && v<1.e15){
00589 uv=sqrt(u*u+v*v)/sqrt(2.0);
00590 }
00591
00592 return uv;
00593 }
00594
00595
00596
00597 void ShwfitAna::FitTShower(Float_t )
00598 {
00599 fShwfit.tenestu->Fit(fShwfit.ufit,"RLQ0+");
00600 fShwfit.trans_u_mean = fShwfit.ufit->GetParameter(1);
00601 fShwfit.trans_u_sigma = fShwfit.ufit->GetParameter(2);
00602 fShwfit.trans_u_chisq = fShwfit.ufit->GetChisquare();
00603 fShwfit.trans_u_ndf = fShwfit.ufit->GetNDF();
00604
00605 fShwfit.tenestv->Fit(fShwfit.vfit,"RLQ0+");
00606 fShwfit.trans_v_mean = fShwfit.vfit->GetParameter(1);
00607 fShwfit.trans_v_sigma = fShwfit.vfit->GetParameter(2);
00608 fShwfit.trans_v_chisq = fShwfit.vfit->GetChisquare();
00609 fShwfit.trans_v_ndf = fShwfit.vfit->GetNDF();
00610
00611 }
00612
00613
00614 void ShwfitAna::TransVar(TH1F *histo, PlaneView::EPlaneView pv)
00615 {
00616 MSG("ShwfitAna",Msg::kDebug)<<"In ShwfitAna::TransVar"<<endl;
00617 MSG("ShwfitAna",Msg::kDebug)<<"plane view is "<<pv<<endl;
00618 MSG("ShwfitAna",Msg::kDebug)<<"Hist name is "<<histo->GetName()<<endl;
00619 MSG("ShwfitAna",Msg::kDebug)<<"Hist has "<<histo->GetEntries()<<" entries"<<endl;
00620
00621 Int_t THISTBINS = (Int_t)(histo->GetNbinsX());
00622
00623 Int_t binmax=histo->GetMaximumBin();
00624 Int_t binvert=Int_t(((THISTBINS-1)/2)+1);
00625 Float_t peak_bin=histo->Integral(binmax,binmax);
00626 Float_t vert_bin=histo->Integral(binvert,binvert);
00627 Float_t tot=histo->Integral(1,THISTBINS);
00628
00629 molrad_peak = molrad_vert = 0;
00630 mean = rms = skew = kurt = 0;
00631
00632 asym_peak= ANtpDefVal::kFloat;
00633 asym_vert= ANtpDefVal::kFloat;
00634
00635 if(tot-peak_bin){
00636 asym_peak=(TMath::Abs(histo->Integral(binmax+1,THISTBINS)-histo->Integral(1,binmax-1)))
00637 /(tot-peak_bin);
00638 }
00639
00640 if(tot-vert_bin){
00641 asym_vert=(TMath::Abs(histo->Integral(binvert+1,THISTBINS)-histo->Integral(1,binvert-1)))
00642 /(tot-vert_bin);
00643 }
00644
00645 Float_t ratio;
00646
00647 if(tot){
00648 for(Int_t i=0; i<=binvert-1;i++){
00649 ratio=histo->Integral(binmax-i>0?binmax-i:1,
00650 binmax+i<THISTBINS?binmax+i:THISTBINS)/tot;
00651 if(ratio>0.90){molrad_peak=i+1; break;}
00652 }
00653 for(Int_t i=0; i<=binvert-1;i++){
00654 ratio=histo->Integral(binvert-i,binvert+i)/tot;
00655 if(ratio>0.90){molrad_vert=i+1; break;}
00656 }
00657
00658 }
00659
00660 mean=histo->GetMean();
00661 rms=histo->GetRMS();
00662 Int_t n_count=0;
00663
00664 for(Int_t i=1;i<=THISTBINS;i++){
00665
00666 if(histo->GetBinContent(i)){
00667 skew=skew+(TMath::Power((histo->GetBinCenter(i)-mean),3)
00668 *histo->GetBinContent(i));
00669 kurt=kurt+(TMath::Power((histo->GetBinCenter(i)-mean),4)
00670 *histo->GetBinContent(i));
00671 n_count++;
00672 }
00673
00674 }
00675 if(rms>0 && n_count>1){
00676 skew=skew/((Float_t)(n_count-1)*TMath::Power((rms),3));
00677 kurt=(kurt/((Float_t)(n_count-1)*TMath::Power((rms),4)))-3;
00678 }
00679 else{skew= ANtpDefVal::kFloat; kurt= ANtpDefVal::kFloat;}
00680
00681 }
00682
00683 void ShwfitAna::Reset(int snarl, int event)
00684 {
00685
00686
00687
00688
00689 const Int_t LHISTBINS=30;
00690 const Int_t THISTBINS=41;
00691 fShwfit.Reset();
00692
00693 char ln[100];
00694 char tun[100];
00695 char tvn[100];
00696 sprintf(ln,"lenepl_%d_%d",snarl,event);
00697 sprintf(tun,"tenestu_%d_%d",snarl,event);
00698 sprintf(tvn,"tenestv_%d_%d",snarl,event);
00699 fShwfit.lenepl = new TH1F(ln,"longitudinal energy by plane",
00700 LHISTBINS,0.0,LHISTBINS);
00701
00702 fShwfit.tenestu = new TH1F(tun,"trasverse energy by strip (U view)",
00703 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00704 fShwfit.tenestv = new TH1F(tvn,"trasverse energy by strip (V view)",
00705 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00706
00707 char tun_9s_2pe[100];
00708 char tvn_9s_2pe[100];
00709 char tun_9s_2pe_dw[100];
00710 char tvn_9s_2pe_dw[100];
00711
00712 sprintf(tun_9s_2pe,"tenestu_9s_2pe_%d_%d",snarl,event);
00713 sprintf(tvn_9s_2pe,"tenestv_9s_2pe_%d_%d",snarl,event);
00714 sprintf(tun_9s_2pe_dw,"tenestu_9s_2pe_dw_%d_%d",snarl,event);
00715 sprintf(tvn_9s_2pe_dw,"tenestv_9s_2pe_dw_%d_%d",snarl,event);
00716
00717 fShwfit.tenestu_9s_2pe = new TH1F(tun_9s_2pe,"trasverse energy by strip (U view)",
00718 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00719 fShwfit.tenestv_9s_2pe = new TH1F(tvn_9s_2pe,"trasverse energy by strip (V view)",
00720 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00721 fShwfit.tenestu_9s_2pe_dw = new TH1F(tun_9s_2pe_dw,"trasverse energy by strip (U view)",
00722 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00723 fShwfit.tenestv_9s_2pe_dw = new TH1F(tvn_9s_2pe_dw,"trasverse energy by strip (V view)",
00724 THISTBINS,-THISTBINS/2.,THISTBINS/2.);
00725
00726
00727
00728 int hmin=0;
00729 int hmax=30;
00730 int npare=3;
00731 int nparh=6;
00732 char efn[100];
00733 char hfn[100];
00734 sprintf(efn,"efit_%d_%d",snarl,event);
00735 sprintf(hfn,"hfit_%d_%d",snarl,event);
00736 fShwfit.efit = new TF1(efn,shwfunc,hmin+0.001,hmax,npare);
00737 fShwfit.hfit = new TF1(hfn,hadfunc,hmin+0.001,hmax,nparh);
00738 fShwfit.efit->SetParNames("a","b","e0");
00739 fShwfit.efit->SetParLimits(0,hmin+0.001,hmax);
00740 fShwfit.efit->SetParLimits(1,0.001,20000);
00741 fShwfit.efit->SetParLimits(2,0+0.001,1000000);
00742
00743 fShwfit.hfit->SetParNames("a","b","e0");
00744 fShwfit.hfit->SetParLimits(0,hmin+0.001,hmax);
00745 fShwfit.hfit->SetParLimits(1,0.001,20000);
00746 fShwfit.hfit->SetParLimits(2,0.001,1000000);
00747
00748
00749
00750 char phh[100];
00751
00752 sprintf(phh,"ph_hist_%d_%d",snarl,event);
00753 fShwfit.ph_hist = new TH1F(phh,"energy per srip", 100,0.5,10.5);
00754
00755 int npar_maxwell=2;
00756 char mfn[100];
00757 sprintf(mfn,"efit_maxwell_%d_%d",snarl,event);
00758 fShwfit.efit_maxwell = new TF1(mfn,maxwell,hmin+0.001,hmax,npar_maxwell);
00759 fShwfit.efit_maxwell->SetParNames("Beta","Energy");
00760 fShwfit.efit_maxwell->SetParLimits(0,hmin+0.001,hmax);
00761
00762 int npar_maxwell3=2;
00763 char mfn3[100];
00764 sprintf(mfn3,"efit_maxwell3_%d_%d",snarl,event);
00765 fShwfit.efit_maxwell3 = new TF1(mfn3,maxwell3,hmin+0.001,hmax,npar_maxwell3);
00766 fShwfit.efit_maxwell3->SetParNames("Beta","Energy");
00767 fShwfit.efit_maxwell3->SetParLimits(0,hmin+0.001,hmax);
00768
00769 char ufn[100];
00770 sprintf(ufn,"ufit_%d_%d",snarl,event);
00771 fShwfit.ufit = new TF1(ufn,"gaus",-20,20);
00772
00773 char vfn[100];
00774 sprintf(vfn,"vfit_%d_%d",snarl,event);
00775 fShwfit.vfit = new TF1(vfn,"gaus",-20,20);
00776
00777 }
00778 Bool_t ShwfitAna::PassCuts(int PhNStrips, int PhNPlanes )
00779 {
00780 bool passstrips=false;
00781 bool passplanes=false;
00782 if(PhNStrips>0&&fShwfit.hiPhStripCount>PhNStrips||PhNStrips<=0) passstrips=true;
00783 if(PhNPlanes>0&&fShwfit.hiPhPlaneCount>PhNPlanes||PhNPlanes<=0) passplanes=true;
00784 cout << "StripCount " << fShwfit.hiPhStripCount << " PlaneCount " << fShwfit.hiPhPlaneCount << endl;
00785 if(passstrips&&passplanes) return true;
00786 else return false;
00787 }
00788
00789 Double_t ShwfitAna::GetMaximumX(TF1* efit, Double_t xmin, Double_t xmax)
00790 {
00791 Double_t fXmin = 0.001;
00792 Double_t fXmax = 30;
00793 Double_t fNpx = 30;
00794
00795 Double_t x,y;
00796 if (xmin >= xmax) {xmin = fXmin; xmax = fXmax;}
00797 Double_t dx = (xmax-xmin)/fNpx;
00798 Double_t xxmax = xmin;
00799 Double_t yymax = efit->Eval(xmin+dx);
00800 for (Int_t i=0;i<fNpx;i++) {
00801 x = xmin + (i+0.5)*dx;
00802 y = efit->Eval(x);
00803 if (y > yymax) {xxmax = x; yymax = y;}
00804 }
00805 if (dx < 1.e-9*(fXmax-fXmin)) return TMath::Min(xmax,xxmax);
00806 else return GetMaximumX(efit, TMath::Max(xmin,xxmax-dx), TMath::Min(xmax,xxmax+dx));
00807 }
00808
00809
00810 void ShwfitAna::doSlopes( NtpSREvent *event , RecRecordImp<RecCandHeader> *srobj)
00811 {
00812
00813
00814
00815 Float_t Uslope=0;
00816 Float_t Uuncer=0;
00817 Float_t UOffset=0;
00818 Float_t Ugoodfit=0;
00819
00820 Float_t Vslope=0;
00821 Float_t Vuncer=0;
00822 Float_t VOffset=0;
00823 Float_t Vgoodfit=0;
00824
00825 Float_t USlopeMom=0;
00826 Float_t VSlopeMom=0;
00827
00828 Float_t UBeamLike=0;
00829 Float_t VBeamLike=0;
00830
00831 Float_t UBeamLikeOffset=0;
00832 Float_t VBeamLikeOffset=0;
00833
00834 Float_t VWd=0;
00835
00836 Float_t UWd=0;
00837
00838
00839 Float_t num[1000];
00840 Float_t den[1000];
00841 Float_t zpos[1000];
00842 Float_t wcenter[1000];
00843
00844 Int_t stripcount[1000];
00845 Float_t stripe[1000];
00846
00847 Int_t umin=10000;
00848 Int_t umax=0;
00849 Int_t vmin=10000;
00850 Int_t vmax=0;
00851 Int_t curplane;
00852 Float_t curstrip;
00853
00854 Float_t ULongE=0;
00855 Float_t VLongE=0;
00856
00857 Float_t tote=0.0;
00858
00859 Float_t weight[1000];
00860 for(int i=0;i<1000;i++){
00861 num[i]=0;
00862 den[i]=0;
00863 wcenter[i]=0;
00864 weight[i]=0;
00865 zpos[i]=0.0;
00866 stripcount[i]=0;
00867 stripe[i]=0;
00868 }
00869
00870
00871
00872
00873
00874 Float_t beamslope=0;
00875 if(srobj->GetHeader().GetVldContext().GetDetector()== Detector::kFar)
00876 beamslope=0.037058;
00877 else if(srobj->GetHeader().GetVldContext().GetDetector()== Detector::kNear)
00878 beamslope=-0.037058;
00879
00880
00881 for(int i=0;i<event->nstrip;i++){
00882 Int_t index = SntpHelpers::GetStripIndex(i,event);
00883 NtpSRStrip *strip = SntpHelpers::GetStrip(index,srobj);
00884 if(!strip) continue;
00885 if(!evtstp0mip){
00886 MSG("ShwfitAna",Msg::kError)<<"No mip strip information"<<endl;
00887 continue;
00888 }
00889
00890 Float_t stripPh = evtstp0mip[index] + evtstp1mip[index];
00891
00892 curplane=strip->plane;
00893
00894 stripcount[curplane]++;
00895 stripe[curplane]+=stripPh;
00896
00897 curstrip=strip->tpos;
00898 zpos[curplane]=strip->z;
00899
00900 num[curplane]+=curstrip*stripPh;
00901 den[curplane]+=stripPh;
00902
00903 tote+=stripPh;
00904
00905 if(strip->planeview==PlaneView::kU)
00906 {
00907 if(umin>curplane)umin=curplane;
00908 if(umax<curplane)umax=curplane;
00909 USlopeMom+=stripPh;
00910 if((strip->z-event->vtx.z)>0.0)
00911 ULongE+=stripPh*TMath::Cos(TMath::ATan((strip->tpos-event->vtx.u)/(strip->z-event->vtx.z)));
00912
00913
00914 }
00915 else if(strip->planeview==PlaneView::kV)
00916 {
00917 if(vmin>curplane)vmin=curplane;
00918 if(vmax<curplane)vmax=curplane;
00919 VSlopeMom+=stripPh;
00920 if((strip->z-event->vtx.z)>0.0)
00921 VLongE+=stripPh*TMath::Cos(TMath::ATan((strip->tpos-event->vtx.v)/(strip->z-event->vtx.z)));
00922
00923
00924 }
00925
00926 }
00927
00928
00929
00930 for(int i=umin;i<=umax;i+=2)
00931 {
00932 if(den[i]!=0){
00933 wcenter[i]=num[i]/den[i];
00934 UBeamLikeOffset+=den[i]*(wcenter[i]-beamslope*zpos[i]);
00935 }else{
00936 num[i]=0;
00937 wcenter[i]=0;
00938 }
00939
00940 }
00941
00942 UBeamLikeOffset=UBeamLikeOffset/USlopeMom;
00943
00944
00945 for(int i=vmin;i<=vmax;i+=2)
00946 {
00947 if(den[i]!=0){
00948 wcenter[i]=num[i]/den[i];
00949 VBeamLikeOffset+=den[i]*(wcenter[i]-beamslope*zpos[i]);
00950 }else{
00951 num[i]=0;
00952 wcenter[i]=0;
00953 }
00954 }
00955
00956 VBeamLikeOffset=VBeamLikeOffset/VSlopeMom;
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966 Int_t stripmin=umin; if(vmin<umin)stripmin=vmin;
00967 Int_t stripmax=umax; if(vmax>umax)stripmax=vmax;
00968
00969 Float_t complex=0;
00970 Float_t wcomplex=0;
00971 Float_t ncomplex=0;
00972
00973 for (int i=stripmin;i<=stripmax;i++)
00974 {
00975 complex+=stripcount[i]*stripcount[i+1];
00976 wcomplex+=stripcount[i]*stripcount[i+1]*stripe[i]*stripe[i+1];
00977 ncomplex+=stripe[i]*stripe[i+1];
00978 }
00979
00980 if(ncomplex>0) wcomplex=wcomplex/ncomplex;
00981 else wcomplex=0;
00982
00983
00984
00985 Float_t VSw=0;
00986 Float_t VSwxy=0;
00987 Float_t VSwx=0;
00988 Float_t VSwy=0;
00989 Float_t VSwxx=0;
00990
00991 Float_t USw=0;
00992 Float_t USwxy=0;
00993 Float_t USwx=0;
00994 Float_t USwy=0;
00995 Float_t USwxx=0;
00996
00997
00998
00999
01000 Float_t pos;
01001
01002
01003
01004
01005 for(int i=umin;i<=umax;i+=2){
01006
01007 if (i>999)continue;
01008
01009 UBeamLike+=den[i]*den[i]*(wcenter[i]-UBeamLikeOffset-beamslope*zpos[i])*(wcenter[i]-UBeamLikeOffset-beamslope*zpos[i]);
01010
01011
01012 if(den[i]>0){
01013
01014 pos=num[i]/(den[i]);
01015 weight[i]=den[i]/tote;
01016
01017 USw=USw+weight[i];
01018 USwxy=USwxy+weight[i]*zpos[i]*pos;
01019 USwx=USwx+weight[i]*zpos[i];
01020 USwy=USwy+weight[i]*pos;
01021 USwxx=USwxx+weight[i]*zpos[i]*zpos[i];
01022
01023
01024 }
01025
01026
01027 }
01028
01029 UBeamLike=TMath::Sqrt(UBeamLike)/USlopeMom;
01030
01031
01032 for(int i=vmin;i<=vmax;i+=2){
01033
01034
01035 if (i>999)continue;
01036
01037 VBeamLike+=den[i]*den[i]*(wcenter[i]-VBeamLikeOffset-beamslope*zpos[i])*(wcenter[i]-VBeamLikeOffset-beamslope*zpos[i]);
01038
01039
01040 if(den[i]>0){
01041
01042 pos=num[i]/(den[i]);
01043 weight[i]=den[i]/tote;
01044
01045 VSw=VSw+weight[i];
01046 VSwxy=VSwxy+weight[i]*zpos[i]*pos;
01047 VSwx=VSwx+weight[i]*zpos[i];
01048 VSwy=VSwy+weight[i]*pos;
01049 VSwxx=VSwxx+weight[i]*zpos[i]*zpos[i];
01050
01051
01052
01053 }
01054
01055
01056
01057 }
01058
01059 VBeamLike=TMath::Sqrt(VBeamLike)/VSlopeMom;
01060
01061 Float_t Udelta=(USw*USwxx-USwx*USwx);
01062
01063 if (Udelta>0){
01064 Uslope=(USw*USwxy-USwx*USwy)/Udelta;
01065 Uuncer=TMath::Sqrt(USw/Udelta);
01066
01067
01068
01069 UOffset=(USwxx*USwy-USwx*USwxy)/Udelta;
01070
01071
01072 Float_t Uvtx=event->vtx.u;
01073
01074
01075
01076 Ugoodfit=0;
01077 Float_t SumWeight=0;
01078
01079 for(int i=umin;i<=umax;i+=2){
01080 if (i>999)continue;
01081 if(den[i]>0){
01082 Ugoodfit=Ugoodfit+(UOffset+Uslope*zpos[i])*(UOffset+Uslope*zpos[i])/(weight[i]*weight[i]);
01083 SumWeight=SumWeight+1/weight[i];
01084
01085 UWd=UWd+(UOffset+Uslope*zpos[i]-Uvtx)*(UOffset+Uslope*zpos[i]-Uvtx)/(weight[i]*weight[i]);
01086
01087 }
01088 }
01089
01090 Ugoodfit=TMath::Sqrt(Ugoodfit)/SumWeight;
01091
01092 UWd=TMath::Sqrt(UWd)/SumWeight;
01093
01094
01095
01096
01097 }
01098
01099
01100
01101
01102 Float_t Vdelta=(VSw*VSwxx-VSwx*VSwx);
01103
01104 if (Vdelta>0){
01105 Vslope=(VSw*VSwxy-VSwx*VSwy)/Vdelta;
01106 Vuncer=TMath::Sqrt(VSw/Vdelta);
01107
01108
01109
01110 VOffset=(VSwxx*VSwy-VSwx*VSwxy)/Vdelta;
01111
01112
01113 Float_t Vvtx=event->vtx.v;
01114
01115 Vgoodfit=0;
01116 Float_t SumWeight=0;
01117 for(int i=vmin;i<=vmax;i+=2){
01118 if (i>999)continue;
01119 if(den[i]>0){
01120 Vgoodfit=Vgoodfit+(VOffset+Vslope*zpos[i])*(VOffset+Vslope*zpos[i])/(weight[i]*weight[i]);
01121 SumWeight=SumWeight+1/weight[i];
01122
01123 VWd=VWd+(VOffset+Vslope*zpos[i]-Vvtx)*(VOffset+Vslope*zpos[i]-Vvtx)/(weight[i]*weight[i]);
01124
01125 }
01126 }
01127
01128 Vgoodfit=TMath::Sqrt(Vgoodfit)/SumWeight;
01129
01130 VWd=TMath::Sqrt(VWd)/SumWeight;
01131 }
01132
01133
01134 if(USlopeMom>0){
01135 USlopeMom=1/USlopeMom;
01136 USlopeMom=USlopeMom*Uslope;
01137 }
01138
01139 if(VSlopeMom>0){
01140 VSlopeMom=1/VSlopeMom;
01141 VSlopeMom=VSlopeMom*Vslope;
01142 }
01143
01144
01145 fShwfit.UBeamLike=UBeamLike;
01146 fShwfit.VBeamLike=VBeamLike;
01147
01148 fShwfit.slopefix=beamslope;
01149 fShwfit.USlope=Uslope;
01150 fShwfit.VSlope=Vslope;
01151 fShwfit.UOffset=UOffset;
01152 fShwfit.VOffset=VOffset;
01153 fShwfit.UFitquality=Ugoodfit;
01154 fShwfit.VFitquality=Vgoodfit;
01155 fShwfit.UWDiff=UWd;
01156 fShwfit.VWDiff=VWd;
01157 fShwfit.USlopeMom=USlopeMom;
01158 fShwfit.VSlopeMom=VSlopeMom;
01159 fShwfit.UVSlope=Uslope*Uslope+Vslope*Vslope;
01160
01161 fShwfit.ULongE=ULongE;
01162 fShwfit.VLongE=VLongE;
01163 fShwfit.LongE = ULongE + VLongE;
01164
01165
01166 fShwfit.complexity=complex;
01167
01168 fShwfit.wcomplexity=wcomplex;
01169
01170 return;
01171 }
01172
01173
01174
01175
