ShwfitAna Class Reference

#include <ShwfitAna.h>

Inheritance diagram for ShwfitAna:
NueAnaBase

List of all members.

Public Member Functions

 ShwfitAna (Shwfit &sf)
virtual ~ShwfitAna ()
void Analyze (int evtn, RecRecordImp< RecCandHeader > *srobj)
void doSlopes (NtpSREvent *event, RecRecordImp< RecCandHeader > *srobj)
void FitLShower (Float_t pulseheight)
void FitLShower_Dan (Float_t pulseheight)
void TransVar (TH1F *h, PlaneView::EPlaneView pv)
void Reset (int snarl, int event)
void SetCutParams (int planes, float striphcut, float planephcut, float contplanephcut)
Bool_t PassCuts (int PhNStrips, int PhNPlanes)
void FitTShower (Float_t pulseheight)

Public Attributes

int sfDPlaneCut
float sfContPhPlaneCut
float sfPhStripCut
float sfPhPlaneCut

Private Member Functions

Double_t GetMaximumX (TF1 *efit, Double_t xmin=0, Double_t xmax=0)
float BuildUVVar (float u, float v)

Private Attributes

ShwfitfShwfit
float asym_peak
float asym_vert
float molrad_peak
float molrad_vert
float mean
float rms
float skew
float kurt

Detailed Description

Definition at line 11 of file ShwfitAna.h.


Constructor & Destructor Documentation

ShwfitAna::ShwfitAna ( Shwfit sf  ) 

Definition at line 77 of file ShwfitAna.cxx.

00077                               :
00078    fShwfit(sf)
00079 {}

ShwfitAna::~ShwfitAna (  )  [virtual]

Definition at line 81 of file ShwfitAna.cxx.

00082 {}


Member Function Documentation

void ShwfitAna::Analyze ( int  evtn,
RecRecordImp< RecCandHeader > *  srobj 
) [virtual]

Implements NueAnaBase.

Definition at line 86 of file ShwfitAna.cxx.

References asym_peak, asym_vert, BuildUVVar(), Shwfit::contPlaneCount, Shwfit::contPlaneCount015, Shwfit::contPlaneCount030, Shwfit::contPlaneCount050, Shwfit::contPlaneCount075, Shwfit::contPlaneCount100, Shwfit::contPlaneCount200, doSlopes(), Shwfit::energyPlane0, Shwfit::energyPlane1, Shwfit::energyPlane2, NueAnaBase::evtstp0mip, NueAnaBase::evtstp1mip, NtpVtxFinder::FindVertex(), FitLShower(), FitLShower_Dan(), FitTShower(), fShwfit, SntpHelpers::GetEvent(), RecRecordImp< T >::GetHeader(), RecPhysicsHeader::GetSnarl(), VHS::GetStrip(), SntpHelpers::GetStripIndex(), SntpHelpers::GetTrack(), SntpHelpers::GetTrackIndex(), Shwfit::hiPhPlaneCount, Shwfit::hiPhPlaneCountM2, Shwfit::hiPhPlaneCountM4, Shwfit::hiPhPlaneCountP2, Shwfit::hiPhPlaneCountP4, Shwfit::hiPhStripCount, Shwfit::hiPhStripCountM2, Shwfit::hiPhStripCountM4, Shwfit::hiPhStripCountP2, Shwfit::hiPhStripCountP4, ReleaseType::IsCedar(), Msg::kDebug, Msg::kError, ANtpDefaultValue::kFloat, PlaneView::kU, kurt, PlaneView::kV, Shwfit::lenepl, mean, molrad_peak, molrad_vert, MSG, NtpSRPlane::n, Shwfit::par_a, Shwfit::par_b, Shwfit::par_e0, NtpSRPulseHeight::pe, NtpSRStrip::ph0, NtpSRStrip::ph1, Shwfit::ph_hist, NtpSRStrip::plane, NtpSRTrack::plane, NtpSRStrip::planeview, NueAnaBase::release, Reset(), rms, sfContPhPlaneCut, sfDPlaneCut, sfPhPlaneCut, sfPhStripCut, skew, CalDetConstants::STRIPWIDTH, Shwfit::tenestu, Shwfit::tenestu_9s_2pe, Shwfit::tenestu_9s_2pe_dw, Shwfit::tenestv, Shwfit::tenestv_9s_2pe, Shwfit::tenestv_9s_2pe_dw, NtpSRStrip::tpos, TransVar(), Shwfit::u_asym_peak, Shwfit::u_asym_peak_9s_2pe, Shwfit::u_asym_peak_9s_2pe_dw, Shwfit::u_asym_vert, Shwfit::u_asym_vert_9s_2pe, Shwfit::u_asym_vert_9s_2pe_dw, Shwfit::u_kurt, Shwfit::u_kurt_9s_2pe, Shwfit::u_kurt_9s_2pe_dw, Shwfit::u_mean, Shwfit::u_mean_9s_2pe, Shwfit::u_mean_9s_2pe_dw, Shwfit::u_molrad_peak, Shwfit::u_molrad_peak_9s_2pe, Shwfit::u_molrad_peak_9s_2pe_dw, Shwfit::u_molrad_vert, Shwfit::u_molrad_vert_9s_2pe, Shwfit::u_molrad_vert_9s_2pe_dw, Shwfit::u_rms, Shwfit::u_rms_9s_2pe, Shwfit::u_rms_9s_2pe_dw, Shwfit::u_skew, Shwfit::u_skew_9s_2pe, Shwfit::u_skew_9s_2pe_dw, Shwfit::uv_asym_peak, Shwfit::uv_asym_peak_9s_2pe, Shwfit::uv_asym_peak_9s_2pe_dw, Shwfit::uv_asym_vert, Shwfit::uv_asym_vert_9s_2pe, Shwfit::uv_asym_vert_9s_2pe_dw, Shwfit::uv_kurt, Shwfit::uv_kurt_9s_2pe, Shwfit::uv_kurt_9s_2pe_dw, Shwfit::uv_mean, Shwfit::uv_mean_9s_2pe, Shwfit::uv_mean_9s_2pe_dw, Shwfit::uv_molrad_peak, Shwfit::uv_molrad_peak_9s_2pe, Shwfit::uv_molrad_peak_9s_2pe_dw, Shwfit::uv_molrad_vert, Shwfit::uv_molrad_vert_9s_2pe, Shwfit::uv_molrad_vert_9s_2pe_dw, Shwfit::uv_ratio, Shwfit::uv_ratio_9s_2pe, Shwfit::uv_ratio_9s_2pe_dw, Shwfit::uv_rms, Shwfit::uv_rms_9s_2pe, Shwfit::uv_rms_9s_2pe_dw, Shwfit::uv_skew, Shwfit::uv_skew_9s_2pe, Shwfit::uv_skew_9s_2pe_dw, Shwfit::v_asym_peak, Shwfit::v_asym_peak_9s_2pe, Shwfit::v_asym_peak_9s_2pe_dw, Shwfit::v_asym_vert, Shwfit::v_asym_vert_9s_2pe, Shwfit::v_asym_vert_9s_2pe_dw, Shwfit::v_kurt, Shwfit::v_kurt_9s_2pe, Shwfit::v_kurt_9s_2pe_dw, Shwfit::v_mean, Shwfit::v_mean_9s_2pe, Shwfit::v_mean_9s_2pe_dw, Shwfit::v_molrad_peak, Shwfit::v_molrad_peak_9s_2pe, Shwfit::v_molrad_peak_9s_2pe_dw, Shwfit::v_molrad_vert, Shwfit::v_molrad_vert_9s_2pe, Shwfit::v_molrad_vert_9s_2pe_dw, Shwfit::v_rms, Shwfit::v_rms_9s_2pe, Shwfit::v_rms_9s_2pe_dw, Shwfit::v_skew, Shwfit::v_skew_9s_2pe, Shwfit::v_skew_9s_2pe_dw, Shwfit::vtxEnergy, NtpVtxFinder::VtxPlane(), NtpVtxFinder::VtxU(), and NtpVtxFinder::VtxV().

Referenced by NueRecordAna::Analyze(), NueDisplayModule::Analyze(), and NueDisplayModule::UpdateDisplay().

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    //loop over strips to fill histograms
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   // new contPlaneCount - Minerba
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;//Munits::meters
00169       MSG("ShwfitAna",Msg::kDebug)<< "plane " << deltaplanes << " stripPh " << stripPh <<endl;
00170       //fill longitudnal energy deposition histogram
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       //fill transverse energy deposition histograms
00184       //u view
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       //v view
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       //unknown view
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   //  contPlaneCount 
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; // ~0.75PE
00250   bool PlaneCut030=true; // ~1.50PE
00251   bool PlaneCut050=true; // ~2.00PE
00252   bool PlaneCut075=true; // ~3.00PE
00253   bool PlaneCut100=true; // ~4.00PE
00254   bool PlaneCut200=true; // ~8.00PE
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    //fit EM shower
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    //fill transverse variables (the following variables filled inside
00328    //   function and then stored
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    //fill uv variables
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 }

float ShwfitAna::BuildUVVar ( float  u,
float  v 
) [private]

Definition at line 584 of file ShwfitAna.cxx.

References ANtpDefaultValue::IsDefault(), and ANtpDefaultValue::kFloat.

Referenced by Analyze().

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 }                                                         

void ShwfitAna::doSlopes ( NtpSREvent event,
RecRecordImp< RecCandHeader > *  srobj 
)

Definition at line 810 of file ShwfitAna.cxx.

References Shwfit::complexity, NueAnaBase::evtstp0mip, NueAnaBase::evtstp1mip, fShwfit, VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), VHS::GetStrip(), SntpHelpers::GetStripIndex(), RecHeader::GetVldContext(), Msg::kError, Detector::kFar, Detector::kNear, PlaneView::kU, PlaneView::kV, Shwfit::LongE, MSG, NtpSRStrip::plane, NtpSRStrip::planeview, Shwfit::slopefix, NtpSRStrip::tpos, NtpSRVertex::u, Shwfit::UBeamLike, Shwfit::UFitquality, Shwfit::ULongE, Shwfit::UOffset, Shwfit::USlope, Shwfit::USlopeMom, Shwfit::UVSlope, Shwfit::UWDiff, NtpSRVertex::v, Shwfit::VBeamLike, Shwfit::VFitquality, Shwfit::VLongE, Shwfit::VOffset, Shwfit::VSlope, Shwfit::VSlopeMom, NtpSREvent::vtx, Shwfit::VWDiff, Shwfit::wcomplexity, NtpSRStrip::z, and NtpSRVertex::z.

Referenced by Analyze().

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    //the slope of the beam 3deg inc in Y for far and 3deg dec in Y for near
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          // else
00913          //  ULongE+=stripPh;
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)//only want hits in front of vertex
00921          VLongE+=stripPh*TMath::Cos(TMath::ATan((strip->tpos-event->vtx.v)/(strip->z-event->vtx.z)));
00922          // else
00923          //  VLongE+=stripPh;
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]);  //den is the energy in each plane
00935        }else{
00936          num[i]=0;
00937          wcenter[i]=0;
00938        }
00939 
00940      }
00941   
00942    if(USlopeMom>0){
00943    UBeamLikeOffset=UBeamLikeOffset/USlopeMom; //USlopeMom currently contains total energy in UPlane
00944    }
00945 
00946    for(int i=vmin;i<=vmax;i+=2)
00947      {
00948        if(den[i]!=0){
00949          wcenter[i]=num[i]/den[i];
00950          VBeamLikeOffset+=den[i]*(wcenter[i]-beamslope*zpos[i]);  //den is the energy in each plane
00951        }else{
00952          num[i]=0;
00953          wcenter[i]=0;
00954        }
00955      }
00956    
00957    if(VSlopeMom>0){
00958    VBeamLikeOffset=VBeamLikeOffset/VSlopeMom; //VSlopeMom currently contains total energy in UPlane
00959    }
00960 
00961 
00962 
00963 Int_t stripmin=umin; if(vmin<umin)stripmin=vmin;
00964 Int_t stripmax=umax; if(vmax>umax)stripmax=vmax;
00965 
00966 Float_t complex=0;
00967 Float_t wcomplex=0;
00968  Float_t ncomplex=0;
00969 
00970 for (int i=stripmin;i<=stripmax;i++)
00971 {
00972         complex+=stripcount[i]*stripcount[i+1];
00973         wcomplex+=stripcount[i]*stripcount[i+1]*stripe[i]*stripe[i+1];
00974         ncomplex+=stripe[i]*stripe[i+1];
00975 }   
00976 
00977  if(ncomplex>0) wcomplex=wcomplex/ncomplex;
00978  else wcomplex=0;
00979 
00980 
00981 
00982    Float_t VSw=0;
00983    Float_t VSwxy=0;
00984    Float_t VSwx=0;
00985    Float_t VSwy=0;
00986    Float_t VSwxx=0;
00987    
00988    Float_t USw=0;
00989    Float_t USwxy=0;
00990    Float_t USwx=0;
00991    Float_t USwy=0;
00992    Float_t USwxx=0;
00993 
00994 
00995 
00996 
00997    Float_t pos;
00998 
00999 
01000 
01001 
01002    for(int i=umin;i<=umax;i+=2){
01003 
01004      if (i>999)continue;
01005 
01006      UBeamLike+=den[i]*den[i]*(wcenter[i]-UBeamLikeOffset-beamslope*zpos[i])*(wcenter[i]-UBeamLikeOffset-beamslope*zpos[i]);
01007 
01008              
01009      if(den[i]>0){
01010        
01011        pos=num[i]/(den[i]);
01012        weight[i]=den[i]/tote;
01013        
01014        USw=USw+weight[i];
01015        USwxy=USwxy+weight[i]*zpos[i]*pos;
01016        USwx=USwx+weight[i]*zpos[i];
01017        USwy=USwy+weight[i]*pos;
01018        USwxx=USwxx+weight[i]*zpos[i]*zpos[i];
01019        
01020 
01021      }
01022 
01023          
01024    }
01025 
01026    if(USlopeMom>0){
01027      UBeamLike=TMath::Sqrt(UBeamLike)/USlopeMom;    //USlopeMom still contains total energy in U plane
01028    } 
01029 
01030  
01031    for(int i=vmin;i<=vmax;i+=2){
01032      
01033      
01034      if (i>999)continue;
01035      
01036       VBeamLike+=den[i]*den[i]*(wcenter[i]-VBeamLikeOffset-beamslope*zpos[i])*(wcenter[i]-VBeamLikeOffset-beamslope*zpos[i]);
01037 
01038      
01039      if(den[i]>0){
01040        
01041        pos=num[i]/(den[i]);
01042        weight[i]=den[i]/tote;
01043        
01044        VSw=VSw+weight[i];
01045        VSwxy=VSwxy+weight[i]*zpos[i]*pos;
01046        VSwx=VSwx+weight[i]*zpos[i];
01047        VSwy=VSwy+weight[i]*pos;
01048        VSwxx=VSwxx+weight[i]*zpos[i]*zpos[i];
01049        
01050 
01051        
01052      }
01053 
01054        
01055        
01056    }
01057 
01058    if(VSlopeMom>0){
01059     VBeamLike=TMath::Sqrt(VBeamLike)/VSlopeMom;    //VSlopeMom still contains total energy in V plane
01060    }
01061 
01062    Float_t Udelta=(USw*USwxx-USwx*USwx);
01063 
01064    if (Udelta>0){
01065      Uslope=(USw*USwxy-USwx*USwy)/Udelta;
01066      Uuncer=TMath::Sqrt(USw/Udelta);
01067      
01068      
01069      
01070      UOffset=(USwxx*USwy-USwx*USwxy)/Udelta;
01071     
01072      
01073      Float_t Uvtx=event->vtx.u;
01074 
01075     
01076      
01077       Ugoodfit=0;
01078      Float_t SumWeight=0;
01079      
01080      for(int i=umin;i<=umax;i+=2){
01081        if (i>999)continue;
01082        if(den[i]>0){
01083          Ugoodfit=Ugoodfit+(UOffset+Uslope*zpos[i])*(UOffset+Uslope*zpos[i])/(weight[i]*weight[i]);
01084          SumWeight=SumWeight+1/weight[i];
01085 
01086           UWd=UWd+(UOffset+Uslope*zpos[i]-Uvtx)*(UOffset+Uslope*zpos[i]-Uvtx)/(weight[i]*weight[i]);
01087          
01088        }
01089      }
01090      
01091      Ugoodfit=TMath::Sqrt(Ugoodfit)/SumWeight;
01092      
01093      UWd=TMath::Sqrt(UWd)/SumWeight;
01094      
01095      
01096      
01097 
01098    }
01099    
01100    
01101    
01102    
01103    Float_t Vdelta=(VSw*VSwxx-VSwx*VSwx);
01104 
01105    if (Vdelta>0){
01106      Vslope=(VSw*VSwxy-VSwx*VSwy)/Vdelta;
01107      Vuncer=TMath::Sqrt(VSw/Vdelta);
01108      
01109      
01110      
01111      VOffset=(VSwxx*VSwy-VSwx*VSwxy)/Vdelta;
01112 
01113      
01114      Float_t Vvtx=event->vtx.v;
01115      
01116      Vgoodfit=0;
01117      Float_t SumWeight=0;
01118      for(int i=vmin;i<=vmax;i+=2){
01119        if (i>999)continue;
01120        if(den[i]>0){
01121          Vgoodfit=Vgoodfit+(VOffset+Vslope*zpos[i])*(VOffset+Vslope*zpos[i])/(weight[i]*weight[i]);
01122          SumWeight=SumWeight+1/weight[i];
01123 
01124          VWd=VWd+(VOffset+Vslope*zpos[i]-Vvtx)*(VOffset+Vslope*zpos[i]-Vvtx)/(weight[i]*weight[i]);
01125 
01126        }
01127      }
01128      
01129      Vgoodfit=TMath::Sqrt(Vgoodfit)/SumWeight;
01130      
01131           VWd=TMath::Sqrt(VWd)/SumWeight;
01132    }
01133    
01134 
01135    if(USlopeMom>0){
01136      USlopeMom=1/USlopeMom;
01137      USlopeMom=USlopeMom*Uslope;
01138    }
01139 
01140    if(VSlopeMom>0){
01141      VSlopeMom=1/VSlopeMom;
01142      VSlopeMom=VSlopeMom*Vslope;
01143    }
01144    
01145    
01146    fShwfit.UBeamLike=UBeamLike;
01147    fShwfit.VBeamLike=VBeamLike;
01148 
01149    fShwfit.slopefix=beamslope;
01150    fShwfit.USlope=Uslope;
01151    fShwfit.VSlope=Vslope;
01152    fShwfit.UOffset=UOffset;
01153    fShwfit.VOffset=VOffset;
01154    fShwfit.UFitquality=Ugoodfit;
01155    fShwfit.VFitquality=Vgoodfit;
01156    fShwfit.UWDiff=UWd;
01157    fShwfit.VWDiff=VWd;
01158    fShwfit.USlopeMom=USlopeMom;
01159    fShwfit.VSlopeMom=VSlopeMom;
01160    fShwfit.UVSlope=Uslope*Uslope+Vslope*Vslope;
01161 
01162    fShwfit.ULongE=ULongE;
01163    fShwfit.VLongE=VLongE;
01164    fShwfit.LongE = ULongE + VLongE;
01165    
01166 
01167    fShwfit.complexity=complex;
01168 
01169    fShwfit.wcomplexity=wcomplex;
01170 
01171    return;
01172  }

void ShwfitAna::FitLShower ( Float_t  pulseheight  ) 

Definition at line 445 of file ShwfitAna.cxx.

References Shwfit::caldet_comp, Shwfit::chisq, Shwfit::chisq_ndf, Shwfit::conv, Shwfit::e0_pe_ratio, Shwfit::efit, fShwfit, GetMaximumX(), Msg::kDebug, ANtpDefaultValue::kFloat, Shwfit::lenepl, Shwfit::max_pe_plane, MSG, Shwfit::par_a, Shwfit::par_b, Shwfit::par_e0, Shwfit::shwmax, Shwfit::shwmaxplane, and Shwfit::shwmaxplane_diff.

Referenced by Analyze().

00446 {
00447    fShwfit.efit->SetParameters(3.,0.5,pulseheight);
00448    fShwfit.lenepl->Fit(fShwfit.efit,"RLLQ0+");
00449    //   fShwfit.lenepl->Print("all");
00450    //   fShwfit.efit->Print("all");
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       //      fShwfit.shwmax=1.*fShwfit.lenepl->GetFunction(fShwfit.efit->GetName())->GetMaximumX();
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 }

void ShwfitAna::FitLShower_Dan ( Float_t  pulseheight  ) 

Definition at line 490 of file ShwfitAna.cxx.

References Shwfit::Beta_Maxwell, Shwfit::Beta_Maxwell3, Shwfit::chisq_Maxwell, Shwfit::chisq_Maxwell3, MuELoss::e, Shwfit::E_ratio_2, Shwfit::E_ratio_half, Shwfit::efit_maxwell, Shwfit::efit_maxwell3, Shwfit::Energy_Maxwell, Shwfit::Energy_Maxwell3, fShwfit, Shwfit::lenepl, Shwfit::n_ratio_2, Shwfit::n_ratio_half, Shwfit::ndf_Maxwell, Shwfit::ndf_Maxwell3, Shwfit::ph_hist, Shwfit::pos_E_split, and Shwfit::pos_n_split.

Referenced by Analyze().

00491 {
00492    //DDC Mawell fit Code
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      //cout<<"n_split_point = "<<n_split_point<<" n_split_val = "<<n_split_val<<" 1/2 Integral = "<<n_integral<<endl;
00561      //cout<<"E_split_point = "<<E_split_point<<" E_split_val = "<<E_split_val<<" 1/2 Integral = "<<E_integral<<endl;
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 }

void ShwfitAna::FitTShower ( Float_t  pulseheight  ) 

Definition at line 597 of file ShwfitAna.cxx.

References fShwfit, Shwfit::tenestu, Shwfit::tenestv, Shwfit::trans_u_chisq, Shwfit::trans_u_mean, Shwfit::trans_u_ndf, Shwfit::trans_u_sigma, Shwfit::trans_v_chisq, Shwfit::trans_v_mean, Shwfit::trans_v_ndf, Shwfit::trans_v_sigma, Shwfit::ufit, and Shwfit::vfit.

Referenced by Analyze().

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 }

Double_t ShwfitAna::GetMaximumX ( TF1 *  efit,
Double_t  xmin = 0,
Double_t  xmax = 0 
) [private]

Definition at line 789 of file ShwfitAna.cxx.

References MuELoss::e.

Referenced by FitLShower().

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 }

Bool_t ShwfitAna::PassCuts ( int  PhNStrips,
int  PhNPlanes 
)

Definition at line 778 of file ShwfitAna.cxx.

References fShwfit, Shwfit::hiPhPlaneCount, and Shwfit::hiPhStripCount.

Referenced by NueDisplayModule::UpdateDisplay().

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 }

void ShwfitAna::Reset ( int  snarl,
int  event 
)

Definition at line 683 of file ShwfitAna.cxx.

References Shwfit::efit, Shwfit::efit_maxwell, Shwfit::efit_maxwell3, fShwfit, hadfunc(), Shwfit::hfit, Shwfit::lenepl, maxwell(), maxwell3(), Shwfit::ph_hist, Shwfit::Reset(), shwfunc(), Shwfit::tenestu, Shwfit::tenestu_9s_2pe, Shwfit::tenestu_9s_2pe_dw, Shwfit::tenestv, Shwfit::tenestv_9s_2pe, Shwfit::tenestv_9s_2pe_dw, Shwfit::ufit, and Shwfit::vfit.

Referenced by Analyze().

00684 {
00685 
00686 //putting histogram creators here so we can name them according to
00687 //event and snarl number, and we won't get errors about replacing 
00688 //existing histograms when we have multiple events in a snarl
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    //DDC Mawell fit Code
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 }

void ShwfitAna::SetCutParams ( int  planes,
float  striphcut,
float  planephcut,
float  contplanephcut 
) [inline]

Definition at line 26 of file ShwfitAna.h.

References sfContPhPlaneCut, sfDPlaneCut, sfPhPlaneCut, and sfPhStripCut.

Referenced by NueDisplayModule::Ana(), and NueModule::Analyze().

00026 {sfDPlaneCut=planes; sfPhStripCut=striphcut; sfPhPlaneCut=planephcut; sfContPhPlaneCut=contplanephcut; }

void ShwfitAna::TransVar ( TH1F *  h,
PlaneView::EPlaneView  pv 
)

Definition at line 614 of file ShwfitAna.cxx.

References asym_peak, asym_vert, Msg::kDebug, ANtpDefaultValue::kFloat, kurt, mean, molrad_peak, molrad_vert, MSG, rms, and skew.

Referenced by Analyze().

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 }


Member Data Documentation

float ShwfitAna::asym_peak [private]

Definition at line 44 of file ShwfitAna.h.

Referenced by Analyze(), and TransVar().

float ShwfitAna::asym_vert [private]

Definition at line 45 of file ShwfitAna.h.

Referenced by Analyze(), and TransVar().

Definition at line 40 of file ShwfitAna.h.

Referenced by Analyze(), doSlopes(), FitLShower(), FitLShower_Dan(), FitTShower(), PassCuts(), and Reset().

float ShwfitAna::kurt [private]

Definition at line 51 of file ShwfitAna.h.

Referenced by Analyze(), and TransVar().

float ShwfitAna::mean [private]

Definition at line 48 of file ShwfitAna.h.

Referenced by Analyze(), and TransVar().

float ShwfitAna::molrad_peak [private]

Definition at line 46 of file ShwfitAna.h.

Referenced by Analyze(), and TransVar().

float ShwfitAna::molrad_vert [private]

Definition at line 47 of file ShwfitAna.h.

Referenced by Analyze(), and TransVar().

float ShwfitAna::rms [private]

Definition at line 49 of file ShwfitAna.h.

Referenced by Analyze(), and TransVar().

Definition at line 35 of file ShwfitAna.h.

Referenced by Analyze(), and SetCutParams().

Definition at line 28 of file ShwfitAna.h.

Referenced by Analyze(), and SetCutParams().

Definition at line 38 of file ShwfitAna.h.

Referenced by Analyze(), and SetCutParams().

Definition at line 37 of file ShwfitAna.h.

Referenced by Analyze(), and SetCutParams().

float ShwfitAna::skew [private]

Definition at line 50 of file ShwfitAna.h.

Referenced by Analyze(), and TransVar().


The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1