ShwFit Class Reference

#include <ShwFit.h>

List of all members.

Public Member Functions

 ShwFit ()
 ~ShwFit ()
void Reset ()
void Insert (double ph, double z)
void Fit3d (int psf=0)
void Insert3d (double ph, double u, double v, double z)
void SetDetector (int mydet)
Double_t GetMaximumX (TF1 *efit, Double_t xmin=0, Double_t xmax=0)
void BuildPlaneMap ()
int GetPlaneFromZ (double z)

Public Attributes

double par_a
double par_b
double par_e0
double par_a_err
double par_b_err
double par_e0_err
double chisq
double ndf
double prob
double shwmax
double shwmaxplane
double conv
double pp_chisq
int pp_ndf
int pp_igood
double pp_p
double cmp_chisq
int cmp_ndf
double peakdiff
double pred_e_a
double pred_g_a
double pred_b
double pred_e0
double pred_e_chisq
double pred_e_ndf
double pred_g_chisq
double pred_g_ndf
double zshift
double pre_over
double pre_under
double post_over
double post_under
std::vector< double > v_ph
std::vector< double > v_u
std::vector< double > v_v
std::vector< double > v_z
int detector
int debug

Static Public Attributes

static TH1F * lenepl = 0
static TF1 * efit = 0
static TH1F * lenepl3d = 0
static std::map< double, int > * planemap = 0

Detailed Description

Definition at line 11 of file ShwFit.h.


Constructor & Destructor Documentation

ShwFit::ShwFit (  ) 

Definition at line 91 of file ShwFit.cxx.

References BuildPlaneMap(), debug, MsgService::Instance(), Msg::kDebug, Msg::kError, MSG, planemap, and Reset().

00092 {
00093         debug=0;
00094         if(MsgService::Instance()->IsActive("ShwFit",Msg::kDebug))debug=1;
00095 
00096 
00097         Reset();
00098    //const Int_t LHISTBINS=100;
00099   // if(lenepl)delete lenepl;
00100    //    lenepl = new TH1F("lenepl","longitudinal energy by plane",LHISTBINS,0.0,LHISTBINS);
00101         //lenepl->SetDirectory(0);
00102         if(!planemap)BuildPlaneMap();
00103 
00104 
00105 
00106         
00107         
00108         if(debug)MSG("ShwFit",Msg::kError)<<"ParticleFinder ShwFit is in debug mode... there will be a purposeful memory leak\n";
00109         
00110         
00111 }

ShwFit::~ShwFit (  ) 

Definition at line 184 of file ShwFit.cxx.

References debug, efit, lenepl, and lenepl3d.

00185 {
00186 
00187         if(!debug)
00188         {
00189                 if(lenepl){delete lenepl;lenepl=0;}
00190                 if(lenepl3d){delete lenepl3d;lenepl3d=0;}
00191                 if(efit){delete efit;efit=0;}
00192         }
00193 
00194 }


Member Function Documentation

void ShwFit::BuildPlaneMap (  ) 

Definition at line 113 of file ShwFit.cxx.

References UgliGeomHandle::GetSteelPlnHandleVector(), infile, SimFlag::kData, Msg::kDebug, Detector::kFar, MSG, outfile, and planemap.

Referenced by ShwFit().

00114 {
00115         //try to open plane location file ... for faster running!
00116         
00117         ifstream infile("planemap.txt");
00118         if(infile.is_open())
00119         {
00120                 MSG("ShwFit",Msg::kDebug)<<"loading shwfit planemap from file!....  ";
00121         
00122                 planemap=new std::map<double,int>;
00123                 string line;
00124                 int cnt=0;
00125                 while(!infile.eof())
00126                 {
00127                         getline(infile,line);
00128                         int planei=0;
00129                         double planez=0;
00130                         sscanf(line.c_str(),"%lf %d",&planez, &planei);
00131                         planemap->insert(make_pair(planez,planei));     
00132                         cnt++;
00133                 }
00134                 
00135                 MSG("ShwFit",Msg::kDebug)<<"loaded "<<cnt<<" planes....";
00136         
00137                 infile.close();
00138         }else{
00139 
00140 
00141                 VldContext vldc(Detector::kFar , SimFlag::kData , VldTimeStamp());
00142                 UgliGeomHandle geo(vldc);       
00143 
00144                 planemap=new std::map<double,int>;
00145 
00146                 //iterate over the planes, and insert the upper z and plane number into the map
00147                 std::vector<UgliSteelPlnHandle>  plns = geo.GetSteelPlnHandleVector();
00148  
00149                 ofstream outfile("planemap.txt");
00150                 
00151                 
00152                 for(unsigned int i=0;i<plns.size();i++)
00153                 {
00154                         planemap->insert(make_pair(plns[i].GetZ0(),i));
00155         
00156                         if(outfile.is_open())
00157                         {
00158                                 outfile << (double)plns[i].GetZ0() << " "<<i<<"\n";
00159                         }
00160                 }
00161                 
00162                 outfile.close();
00163         
00164         }
00165         
00166 
00167 
00168 }

void ShwFit::Fit3d ( int  psf = 0  ) 

see if we need a z shift!

fast version //N do not draw

detail version

detail version

Definition at line 365 of file ShwFit.cxx.

References MuELoss::a, chisq, cmp_chisq, cmp_ndf, conv, debug, MuELoss::e, efit, GetMaximumX(), GetPlaneFromZ(), Msg::kDebug, Msg::kError, lenepl3d, Munits::m, MSG, ndf, par_a, par_a_err, par_b, par_b_err, par_e0, par_e0_err, peakdiff, post_over, post_under, pp_chisq, pp_igood, pp_ndf, pp_p, pre_over, pre_under, pred_b, pred_e0, pred_e_a, pred_e_chisq, pred_e_ndf, pred_g_a, pred_g_chisq, pred_g_ndf, prob, MuELoss::R, shwfunc3d(), shwmax, shwmaxplane, v_ph, v_u, v_v, and v_z.

Referenced by Finder::FinalizeParticles3D(), and PrimaryShowerFinder::FindPrimaryShower().

00366 {
00367 
00368         double Rsteel = 1.445; // for steel
00369         double Rscint=0.0241; //for scint only
00370         double R = Rscint+Rsteel;
00371         
00372         double offset=Rsteel / 2; //offset from front of scint plane to count as 0 radiation lengths...
00373 
00374 /*
00375         //build a vector measuring distance in radiation lengths from first hit
00376         VldContext vldc(Detector::kFar , SimFlag::kData , VldTimeStamp());
00377         UgliGeomHandle geo(vldc);       
00378 */
00379         std::vector<int> planes;
00380         
00381         for(unsigned int i=0;i<v_z.size();i++)
00382         {
00383 /*              PlexPlaneId p;
00384                 p = geo.GetPlaneIdFromZ(v_z[i]);        
00385                 planes.push_back(p.GetPlane());
00386 */
00387                 planes.push_back(GetPlaneFromZ(v_z[i]));        
00388         }
00389         
00390 
00391         std::vector<double> rdlen;      
00392         double rddist=offset;   
00393         rdlen.push_back(rddist);
00394                 
00395         int sm1=0;
00396         int sm2=0;
00397         
00398         if(v_z[0]<15.5)sm1=1;
00399         if(v_z[0]>=15.5)sm2=1;
00400         
00401         for(unsigned int i=1;i<v_z.size();i++)
00402         {
00403                 double dz = v_z[i]-v_z[i-1];
00404                 double du = v_u[i]-v_u[i-1];
00405                 double dv = v_v[i]-v_v[i-1];
00406         
00407                 double tdist = sqrt(dz*dz +du*du +dv*dv);
00408 
00409                 if(planes[i]==planes[i-1]) //in the same scint plane
00410                 {
00411                         rddist+=tdist * Rscint;
00412                 }else{
00413                         int dp = planes[i]-planes[i-1];
00414                         double dt = sqrt(du*du+dv*dv);
00415                         rddist += sqrt( R*R*dp*dp +  dt*dt);
00416                 }
00417                 rdlen.push_back(rddist);
00418                 
00419                 if(v_z[i]<15.5)sm1=1;
00420                 if(v_z[i]>=15.5)sm2=1;  
00421         }
00422         
00423         
00424         if(sm1 && sm2)
00425         {
00426                 MSG("ShwFit",Msg::kError)<<"Attempting to fit between super modules! exiting\n";
00427                 return; 
00428         }
00429         
00430         
00431         //for(unsigned int i=0;i<v_z.size();i++)
00432                 //std::cout << "z " << v_z[i] << " plane " <<planes[i] <<" rdlen "<<rdlen[i] <<endl;    
00433         
00434 
00435         //decide what size histo to make, and make it   
00436         int nbins = v_z.size()*2;       
00437         double rlmax = rdlen[rdlen.size()-1]*2;
00438 
00439         //don't delete the histo if we want to be viewing it later
00440         //this will leak... so just don't run this module in debug mode
00441         //except when debugging!
00442         if(!debug) if (lenepl3d){delete lenepl3d;lenepl3d=0;}
00443         
00444         if(debug && psf)
00445                 lenepl3d = new TH1F("lenepl3dPSF","longitudinal energy by plane",nbins,0.0,rlmax);
00446         else
00447                 lenepl3d = new TH1F("lenepl3d","longitudinal energy by plane",nbins,0.0,rlmax);
00448         //hide the histo unless we are debugging the primary shower
00449         if(!debug)      lenepl3d->SetDirectory(0);
00450 
00451 
00452         //try a smear fill...            
00453         for(unsigned int i=0;i<rdlen.size();i++) 
00454         {
00455                 double ph = v_ph[i];
00456                 ph/=R;
00457                 //cout<<"adding to hist "<<rdlen[i]<<" e "<<ph<<"\n";   
00458                 if(i==0)
00459                 {
00460                         lenepl3d->Fill(rdlen[i],ph);//*0.8);
00461                         //lenepl3d->Fill(rdlen[i]+lenepl3d->GetBinWidth(1),ph*0.2);
00462                         //lenepl3d->Fill(rdlen[i]+lenepl3d->GetBinWidth(1)*2,ph*0.1);
00463                 }else{
00464                         //lenepl3d->Fill(rdlen[i]+lenepl3d->GetBinWidth(1),ph*0.2);
00465                         lenepl3d->Fill(rdlen[i],ph);//*0.60);
00466                         //lenepl3d->Fill(rdlen[i]-lenepl3d->GetBinWidth(1),ph*0.2);
00467                         
00468                         //lenepl3d->Fill(rdlen[i]+lenepl3d->GetBinWidth(1),ph*0.05);
00469                         //lenepl3d->Fill(rdlen[i]-lenepl3d->GetBinWidth(1),ph*0.05);
00470                 }
00471         }
00472         
00474         //calculate expected values here....
00476 
00477 
00478 
00479         
00480         
00481         double chisq_elec=-1;
00482         double chisq_gamma=-1;
00483 
00484         //what should a be for this energy?
00485         //a = 1 +b(lnE/Ec -0.5)  (or+0.5 for gamma)  b=0.5
00486         //Ec is about 36.64 MeV assuming steel has z 26 and scint is about z=7
00487 
00488         //calculate total e
00489         double etot=0;
00490         for(int i=0;i<(int)v_ph.size();i++)
00491                 etot+=v_ph[i];
00492         
00493         //make it MeV with roughly 25 mip/gev
00494 
00495         etot=etot/23*1000;// *2; 
00496 
00497         if(etot<1e-5)return;
00498 
00499         pred_b=0.55;
00500         pred_e0=etot;
00501 
00502 
00503         double a=1+0.5*(log(etot/36.64)-0.5);
00504         TF1* predElec = 0;
00505         if(psf)predElec = new TF1("predElecPSF",shwfunc3d,0.1,rlmax,3);
00506         else predElec = new TF1("predElec",shwfunc3d,0.1,rlmax,3);
00507         predElec->SetParNames("a","b","e0");
00508         predElec->SetParameter(0,a);
00509         predElec->SetParameter(1,0.5);  
00510         predElec->SetParameter(2,etot);
00511         predElec->SetLineColor(2);
00512 
00513         pred_e_a=a;
00514         
00515         MSG("ShwFit",Msg::kDebug)<<"elec a b e0 "<<a<<" "<<0.5<<" "<<etot<<"\n";
00516 
00517 
00518         double maxX = GetMaximumX(predElec,0,rlmax);
00519         
00520         int maxbin = lenepl3d->GetMaximumBin();
00521         double maxH = lenepl3d->GetBinCenter(maxbin);
00522         
00523         MSG("ShwFit",Msg::kDebug)<<"pred max "<<maxX<<" hist max "<<maxH<<"\n";
00524         
00525 /*      
00527         //do we need to reset the offset of the histogram?
00528         if(fabs(maxH-maxX)>0.05)
00529         {
00530                 zshift = maxX-maxH;
00531 
00532                 //lenepl3d->Reset();
00533                 if(lenepl3d)delete lenepl3d;lenepl3d=0;
00534                 lenepl3d = new TH1F("lenepl3d","longitudinal energy by plane",nbins,0.0,rlmax);
00535                 lenepl3d->SetDirectory(0);
00536                 
00537                 //try a smear fill...            
00538                 for(unsigned int i=0;i<rdlen.size();i++) 
00539                 {
00540                         double ph = v_ph[i];
00541                         ph/=R;
00542                 
00543                         if(i==0)
00544                         {
00545                                 lenepl3d->Fill(rdlen[i]+zshift,ph);// *0.8);
00546                                 //lenepl3d->Fill(rdlen[i]+zshift+lenepl3d->GetBinWidth(1),ph*0.2);
00547                                 //lenepl3d->Fill(rdlen[i]+zshift+lenepl3d->GetBinWidth(1)*2,ph*0.1);
00548                         }else{
00549                         //      lenepl3d->Fill(rdlen[i]+zshift+lenepl3d->GetBinWidth(1),ph*0.2);
00550                                 lenepl3d->Fill(rdlen[i]+zshift,ph);// *0.60);
00551                         //      lenepl3d->Fill(rdlen[i]+zshift-lenepl3d->GetBinWidth(1),ph*0.2);
00552                         
00553                                 //lenepl3d->Fill(rdlen[i]+zshift+lenepl3d->GetBinWidth(1),ph*0.05);
00554                                 //lenepl3d->Fill(rdlen[i]+zshift-lenepl3d->GetBinWidth(1),ph*0.05);
00555                         }
00556                 }       
00557                 
00558                 //store zshift in terms of actual distance, not radiation length
00559                 zshift/=R; //to get number of planes
00560                 zshift*=0.0354;//average plane width
00561         
00562                 int maxbin = lenepl3d->GetMaximumBin();
00563                 double maxH = lenepl3d->GetBinCenter(maxbin);
00564                 printf("new hist max %f\n",maxH);
00565         }
00567 */
00568 
00569 
00570         a=1+0.5*(log(etot/36.64)+0.5); //gamma
00571         TF1* predGamma = 0;
00572         if(psf) predGamma = new TF1("predGammaPSF",shwfunc3d,0.1,rlmax,3);
00573         else predGamma = new TF1("predGamma",shwfunc3d,0.1,rlmax,3);
00574         predGamma->SetParNames("a","b","e0");
00575         predGamma->SetParameter(0,a);
00576         predGamma->SetParameter(1,0.5);         
00577         predGamma->SetParameter(2,etot);
00578         predGamma->SetLineColor(3);
00579 
00580         pred_g_a=a;
00581         MSG("ShwFit",Msg::kDebug)<<"gamma a b e0 "<<a<<" "<<0.5<<" "<<etot<<"\n";
00582   
00583         //calculate chi sq's ....
00584 
00585   
00586         chisq_elec=0;
00587         chisq_gamma=0;
00588         int dc_e=0;
00589         int dc_g=0;
00590         for(int i=1;i<lenepl3d->GetNbinsX()+1;i++)
00591         {
00592                 double val = lenepl3d->GetBinContent(i);
00593                 double x = lenepl3d->GetBinCenter(i);
00594                 MSG("ShwFit",Msg::kDebug)<<"v "<<val<<" x "<<x<<" eval "<<predElec->Eval(x)<<" "<<predGamma->Eval(x)<<"\n";
00595   
00596                 //      double predval = predElec->Integral( lenepl3d->GetBinLowEdge(i),lenepl3d->GetBinLowEdge(i+1),(double *)0,1e-3);
00597   
00598                 double bincenter = lenepl3d->GetBinCenter(i);
00599                 double predval = predElec->Eval(bincenter);
00600         
00601         //      if(predval>1e-4 && val>0)
00602                 if(predval)chisq_elec+=(val-predval)*(val-predval)/(predval);
00603                 //else dc_e++;
00604         
00605                 //predval = predGamma->Integral( lenepl3d->GetBinLowEdge(i),lenepl3d->GetBinLowEdge(i+1),(double *)0,1e-3);
00606                 predval = predGamma->Eval(bincenter);
00607         
00608                 //if(predval>1e-4 && val>0)
00609                 if(predval)
00610                         chisq_gamma+=(val-predval)*(val-predval)/(predval);
00611                 //else dc_g++;
00612   
00613         }
00614 
00615         pred_e_chisq=chisq_elec;
00616         pred_e_ndf=lenepl3d->GetNbinsX()-dc_e;
00617         pred_g_chisq=chisq_gamma;
00618         pred_g_ndf=lenepl3d->GetNbinsX()-dc_g;  
00619   
00620         chisq_elec/=lenepl3d->GetNbinsX()-dc_e;
00621         chisq_gamma/=lenepl3d->GetNbinsX()-dc_g;
00622   
00623         post_over=0;
00624         post_under=0;
00625         pre_over=0;
00626         pre_under=0;
00627         int m = lenepl3d->GetMaximumBin();
00628  
00629 
00630 
00631 
00633 
00634         //do the fit
00635 
00636         int hmin=0;
00637         int npare=3;
00638    
00639         double pulseheight = lenepl3d->Integral();
00640         if(!debug)      if(efit){delete efit;efit=0;}
00641 
00642         efit = new TF1("efit",shwfunc3d,hmin+0.001,rlmax,npare);
00643         efit->SetParNames("a","b","e0","zshift");
00644         
00645         //efit->SetParLimits(0,hmin+0.001,rlmax);
00646         //efit->SetParLimits(1,0.001,20000);    
00647         //efit->SetParLimits(2,0+0.001,1000000);
00648   
00649         efit->SetParameters(3,0.5,300); 
00650         efit->SetParLimits(0,1.5,5);
00651         efit->SetParLimits(1,0.01,1.5);         
00652         efit->SetParLimits(2,0,1000000);
00653   
00654        
00655 
00656         //efit->SetParameters(3,-5,5);
00657    
00658         //efit->SetParameters(3.,0.5,pulseheight *60*25,0); //gues e as sum of ph in projection....
00659    
00660         TCanvas * c=0; 
00661         
00662         if(debug){
00663                 MSG("ShwFit",Msg::kDebug)<<(psf!=0)<<" "<<(c!=0)<<"\n";
00664                 if(psf){
00665                         c = new TCanvas("psfcan","psfcan");
00666                         c->cd();
00667                         lenepl3d->Draw();
00668                 }
00669                 MSG("ShwFit",Msg::kDebug)<<(psf!=0)<<" "<<(c!=0)<<"\n";
00670         }
00671 
00672         if(!debug){             
00673                 lenepl3d->Fit(efit,"NWLLQ"); 
00674         }else{
00675             if(psf)
00676                         lenepl3d->Fit(efit,"WLLV+"); 
00677                 else
00678                         lenepl3d->Fit(efit,"NWLLV"); 
00679         }       
00680 
00681 
00682         TH1F * predhist= 0;
00683 
00684         if(psf){
00685                 predhist = new TH1F("predhistPSF","longitudinal energy by plane",nbins,0.0,rlmax);
00686         //      predhist->SetDirectory(0);
00687                 predhist->SetLineColor(2);
00688         }else{
00689                 predhist = new TH1F("predhist","longitudinal energy by plane",nbins,0.0,rlmax);
00690                 predhist->SetDirectory(0);
00691                 predhist->SetLineColor(2);
00692         }
00693         
00694         
00695         for(int i=1;i<lenepl3d->GetNbinsX()+1;i++)
00696         {
00697 
00698                 double x = lenepl3d->GetBinCenter(i);
00699                 //      double pred = predElec->Integral(lenepl3d->GetBinLowEdge(i),lenepl3d->GetBinLowEdge(i+1),(double *)0,1e-3);
00700                 double pred = predElec->Eval(lenepl3d->GetBinCenter(i));
00701    
00702                 double val = lenepl3d->GetBinContent(i);
00703         
00704                 if(predhist)predhist->Fill(x,pred);
00705         
00706                 if(i>m)
00707                 {
00708                         if(pred<val)
00709                                 post_over+=val-pred;
00710                         else
00711                                 post_under+=pred-val;
00712                 }else{
00713                         if(pred<val)
00714                                 pre_over+=val-pred;
00715                         else
00716                                 pre_under+=pred-val;
00717                 }
00718         }
00719 
00720         if(debug&&psf&&predhist)predhist->Draw("same");
00721 
00722         //draw the predicted fit here...
00723         if(debug){
00724                 MSG("ShwFit",Msg::kDebug)<<(psf!=0)<<" "<<(c!=0)<<"\n";
00725 
00726                 if(psf && c){
00727                         c->cd();
00728                         predElec->Draw("same");
00729                         predGamma->Draw("same");
00730                         MSG("ShwFit",Msg::kDebug)<<"ASDFASDFASDF\n";
00731                 }
00732         }
00733 
00734         pp_chisq=0;
00735         pp_ndf=0;
00736         pp_igood=0;
00737         pp_p=0;
00738         if(psf){
00739                         //require nonempty histograms to prevent endless loop!
00740         
00741                 if(predhist->Integral()>0.01 && lenepl3d->Integral()>0.01)
00742                 {
00743                         MSG("ShwFit",Msg::kDebug)<<"int pred "<<(predhist->Integral()>0)<<" shw "<< (lenepl3d->Integral()>0)<<"\n";     
00744                         //pp_p=lenepl3d->Chi2TestX(predhist,pp_chisq,pp_ndf,pp_igood,"WW",(double*)0);
00745                 }
00746         }
00747         
00748         //lets directly compare the predicted histogram with the original one
00749         cmp_chisq=0;
00750         cmp_ndf = 0;
00751         peakdiff=0;
00752         if(predhist)
00753         {
00754                 cmp_ndf=predhist->GetNbinsX();
00755                 for(int k=1;k<cmp_ndf+1;k++)
00756                 {
00757                         double oc = lenepl3d->GetBinContent(k);
00758                         double pc = predhist->GetBinContent(k);
00759                         if(pc)cmp_chisq+=(oc-pc)*(oc-pc)/pc;
00760                         else if(oc)cmp_ndf--;
00761                 }
00762                 if(cmp_chisq>0)cmp_chisq=sqrt(cmp_chisq);
00763                 
00764                 int maxbin = predhist->GetMaximumBin();
00765                 if(maxbin>1)
00766                 {
00767                         peakdiff+= lenepl3d->GetBinContent(maxbin-1)-predhist->GetBinContent(maxbin-1);
00768                         peakdiff+= lenepl3d->GetBinContent(maxbin)-predhist->GetBinContent(maxbin);
00769                         peakdiff+= lenepl3d->GetBinContent(maxbin+1)-predhist->GetBinContent(maxbin+1);                                 
00770                 }else if(maxbin==1)
00771                 {
00772                         peakdiff+= lenepl3d->GetBinContent(maxbin)-predhist->GetBinContent(maxbin);
00773                         peakdiff+= lenepl3d->GetBinContent(maxbin+2)-predhist->GetBinContent(maxbin+2);
00774                         peakdiff+= lenepl3d->GetBinContent(maxbin+1)-predhist->GetBinContent(maxbin+1);                                 
00775                 }
00776         
00777         //      printf("-- %f %d %f\n",cmp_chisq,cmp_ndf, peakdiff);
00778 
00779 
00780         }
00781         
00782 
00783 
00784 
00785         MSG("ShwFit",Msg::kDebug)<<" STATUS: "<<gMinuit->fCstatu
00786                         <<" "<<efit->GetParameter(0)
00787                         <<" "<<efit->GetNDF()<<" "<<pulseheight<<endl;
00788 
00789         string fitstatus = (string)(gMinuit->fCstatu);
00790         MSG("ShwFit",Msg::kDebug)<<" STATUS: "<<fitstatus
00791                         <<", "<<efit->GetParameter(0)
00792                         <<", "<<efit->GetNDF()<<" "<<pulseheight<<endl;
00793 
00794         if(fitstatus=="CONVERGED "&&efit->GetParameter(0)<29.9&&
00795                 efit->GetNDF()>0&&pulseheight>0){
00796 
00797                 MSG("ShwFit",Msg::kDebug)<<" filling vars"<<endl;
00798 
00799                 par_a=efit->GetParameter(0);
00800                 par_b=efit->GetParameter(1);
00801                 par_e0=efit->GetParameter(2);
00802   
00803                 par_a_err=efit->GetParError(0);
00804                 par_b_err=efit->GetParError(1);
00805                 par_e0_err=efit->GetParError(2);
00806       
00807                 chisq=efit->GetChisquare();
00808      
00809                 ndf = efit->GetNDF();
00810     
00811 
00812                 //ndf should be equal to the number of bins... because we are giving 0 bins weight 1
00813                 ndf = nbins;
00814 
00815                 if(debug && psf){
00816                         c->cd();
00817                         TLatex l;
00818                         l.SetTextAlign(12);
00819                         l.SetNDC(1);
00820         
00821                         char txt[100];
00822                         sprintf(txt,"\\chi^2=%f",efit->GetChisquare());
00823                         l.DrawLatex(0.5,0.8,txt);    
00824    
00825                         sprintf(txt,"NDF=%f",ndf);
00826                         l.DrawLatex(0.5,0.75,txt);    
00827          
00828                         sprintf(txt,"\\chi^2/NDF=%f",efit->GetChisquare()/ndf);
00829                         l.DrawLatex(0.5,0.7,txt);    
00830 
00831                         sprintf(txt,"par a %f \\pm %f",par_a,par_a_err);
00832                         l.DrawLatex(0.5,0.65,txt);    
00833                 
00834                         sprintf(txt,"par b %f \\pm %f",par_b,par_b_err);
00835                         l.DrawLatex(0.5,0.6,txt);    
00836 
00837                         sprintf(txt,"e0 %f \\pm %f",par_e0,par_e0_err);
00838                         l.DrawLatex(0.5,0.55,txt);    
00839                 
00840                         sprintf(txt,"elec \\chi^2/NDF = %f ",chisq_elec);
00841                         l.DrawLatex(0.5,0.50,txt);    
00842 
00843                         sprintf(txt,"gamma \\chi^2/NDF = %f ",chisq_gamma);
00844                         l.DrawLatex(0.5,0.45,txt);  
00845 
00846                         sprintf(txt,"pre over %f under %f",pre_over,pre_under);
00847                         l.DrawLatex(0.5,0.40,txt);  
00848                         sprintf(txt,"post over %f under %f",post_over,post_under);
00849                         l.DrawLatex(0.5,0.35,txt);  
00850                 
00851                         sprintf(txt,"%f --- %f %d %f", pp_p ,pp_chisq,pp_ndf,pp_chisq/pp_ndf);
00852                         l.DrawLatex(0.5,0.3,txt);       
00853 
00854                         if(cmp_ndf)
00855                         {
00856                                 sprintf(txt,"cmp c2 %f ndf %d c2ndf %f ", cmp_chisq,cmp_ndf,cmp_chisq/(double)cmp_ndf);
00857                                 l.DrawLatex(0.5,0.25,txt);                              
00858                         }
00859 
00860                         sprintf(txt,"pd %f ",peakdiff);
00861                         l.DrawLatex(0.5,0.2,txt);       
00862                         
00863                 }
00864 
00865 
00866      
00867                 shwmax=1.*GetMaximumX(efit);
00868                 prob=efit->GetProb();
00869         
00870                 shwmaxplane=(int)(lenepl3d->
00871                         GetBinContent(lenepl3d->FindBin(shwmax)));
00872                 conv=1;
00873         }
00874 
00875 
00876         if(!psf || !debug)
00877         {
00878                 delete predElec;predElec=0;
00879                 delete predGamma;predGamma=0;
00880                 delete predhist;predhist=0;
00881         }       
00882 
00883         if(!psf)
00884         {
00885                 delete efit;efit=0;
00886 
00887         }
00888 }

Double_t ShwFit::GetMaximumX ( TF1 *  efit,
Double_t  xmin = 0,
Double_t  xmax = 0 
)

Definition at line 263 of file ShwFit.cxx.

References MuELoss::e.

Referenced by Fit3d().

00264 {                           
00265    Double_t fXmin = 0.001;
00266    Double_t fXmax = 100;
00267    Double_t fNpx = 100;                                                    
00268 
00269    Double_t x,y;
00270    if (xmin >= xmax) {xmin = fXmin; xmax = fXmax;}
00271    Double_t dx = (xmax-xmin)/fNpx;
00272    Double_t xxmax = xmin;
00273    Double_t yymax = efit->Eval(xmin+dx);
00274    for (Int_t i=0;i<fNpx;i++) {
00275       x = xmin + (i+0.5)*dx;
00276       y = efit->Eval(x);
00277       if (y > yymax) {xxmax = x; yymax = y;}
00278    }
00279    if (dx < 1.e-9*(fXmax-fXmin)) return TMath::Min(xmax,xxmax);
00280    else return GetMaximumX(efit, TMath::Max(xmin,xxmax-dx), TMath::Min(xmax,xxmax+dx));
00281 }

int ShwFit::GetPlaneFromZ ( double  z  ) 

Definition at line 171 of file ShwFit.cxx.

References planemap.

Referenced by Fit3d().

00172 {
00173         int plane=0;
00174         
00175         std::map<double,int>::iterator i = planemap->lower_bound(z);
00176         
00177         if(i!=planemap->end()) plane=i->second;
00178 
00179         
00180         return plane;
00181 }

void ShwFit::Insert ( double  ph,
double  z 
)

Definition at line 258 of file ShwFit.cxx.

References lenepl.

00259 {
00260         lenepl->Fill(z+0.5,ph);
00261 }

void ShwFit::Insert3d ( double  ph,
double  u,
double  v,
double  z 
)

Definition at line 355 of file ShwFit.cxx.

References v_ph, v_u, v_v, and v_z.

Referenced by Finder::FinalizeParticles3D(), and PrimaryShowerFinder::FindPrimaryShower().

00356 {
00357         v_ph.push_back(ph);
00358         v_u.push_back(u);
00359         v_v.push_back(v);
00360         v_z.push_back(z);
00361 }

void ShwFit::Reset (  ) 

Definition at line 198 of file ShwFit.cxx.

References chisq, cmp_chisq, cmp_ndf, conv, debug, detector, efit, lenepl, lenepl3d, ndf, par_a, par_a_err, par_b, par_b_err, par_e0, par_e0_err, peakdiff, post_over, post_under, pp_chisq, pp_igood, pp_ndf, pp_p, pre_over, pre_under, pred_b, pred_e0, pred_e_a, pred_e_chisq, pred_e_ndf, pred_g_a, pred_g_chisq, pred_g_ndf, prob, shwmax, shwmaxplane, v_ph, v_u, v_v, v_z, and zshift.

Referenced by PrimaryShowerFinder::FindPrimaryShower(), and ShwFit().

00199 {
00200 
00201         if(!debug)
00202         {
00203                 if(lenepl){delete lenepl;lenepl=0;}
00204                 if(lenepl3d){delete lenepl3d;lenepl3d=0;}
00205                 if(efit){delete efit;efit=0;}
00206         }
00207 
00208         cmp_chisq=0;
00209         cmp_ndf = 0;
00210         peakdiff=0;
00211                 
00212         pp_chisq=0;
00213         pp_ndf=0;
00214         pp_igood=0;
00215         pp_p=0;
00216         
00217         lenepl=0;
00218         par_a=0;
00219         par_b=0;
00220         par_e0=0;
00221         par_a_err=0;
00222         par_b_err=0;
00223         par_e0_err=0;   
00224         chisq=0;
00225         ndf=0;
00226         prob=0;
00227         shwmax=0;
00228         shwmaxplane=0;
00229         conv=0;
00230         
00231         
00232         pred_e_a=0;
00233         pred_g_a=0;
00234         pred_b=0;
00235         pred_e0=0;
00236         pred_e_chisq=0;
00237         pred_e_ndf=0;
00238         pred_g_chisq=0;
00239         pred_g_ndf=0;   
00240         
00241 
00242         detector=0;
00243         
00244         v_ph.clear();
00245         v_u.clear();
00246         v_v.clear();
00247         v_z.clear();
00248         
00249         zshift=0;
00250         
00251         pre_over=0;
00252         pre_under=0;
00253         post_over=0;
00254         post_under=0;
00255         
00256 }

void ShwFit::SetDetector ( int  mydet  ) 

Definition at line 350 of file ShwFit.cxx.

References detector.

00351 {
00352         detector=mydet;
00353 }


Member Data Documentation

double ShwFit::chisq
double ShwFit::conv

Definition at line 96 of file ShwFit.h.

Referenced by Fit3d(), Reset(), ShwFit(), and ~ShwFit().

Definition at line 89 of file ShwFit.h.

Referenced by Reset(), and SetDetector().

TF1 * ShwFit::efit = 0 [static]

Definition at line 77 of file ShwFit.h.

Referenced by Fit3d(), Reset(), and ~ShwFit().

TH1F * ShwFit::lenepl = 0 [static]

Definition at line 76 of file ShwFit.h.

Referenced by Insert(), Reset(), and ~ShwFit().

TH1F * ShwFit::lenepl3d = 0 [static]

Definition at line 81 of file ShwFit.h.

Referenced by Fit3d(), Reset(), and ~ShwFit().

double ShwFit::ndf
double ShwFit::par_a
double ShwFit::par_b
std::map< double, int > * ShwFit::planemap = 0 [static]

Definition at line 91 of file ShwFit.h.

Referenced by BuildPlaneMap(), GetPlaneFromZ(), and ShwFit().

double ShwFit::pp_p
double ShwFit::prob

Definition at line 32 of file ShwFit.h.

Referenced by Fit3d(), and Reset().

Definition at line 33 of file ShwFit.h.

Referenced by Fit3d(), and Reset().

std::vector<double> ShwFit::v_ph

Definition at line 84 of file ShwFit.h.

Referenced by Fit3d(), Insert3d(), and Reset().

std::vector<double> ShwFit::v_u

Definition at line 85 of file ShwFit.h.

Referenced by Fit3d(), Insert3d(), and Reset().

std::vector<double> ShwFit::v_v

Definition at line 86 of file ShwFit.h.

Referenced by Fit3d(), Insert3d(), and Reset().

std::vector<double> ShwFit::v_z

Definition at line 87 of file ShwFit.h.

Referenced by Fit3d(), Insert3d(), and Reset().

Definition at line 68 of file ShwFit.h.

Referenced by PrimaryShowerFinder::FindPrimaryShower(), and Reset().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1