FracVarAna Class Reference

#include <FracVarAna.h>

Inheritance diagram for FracVarAna:
NueAnaBase

List of all members.

Public Member Functions

 FracVarAna (FracVar &fv)
virtual ~FracVarAna ()
void Analyze (int evtn, RecRecordImp< RecCandHeader > *srobj)
bool LoadLargestTrackFromEvent (NtpSREvent *event, RecRecordImp< RecCandHeader > *record, int &trkidx)
bool LoadShowerAtTrackVertex (NtpSREvent *event, RecRecordImp< RecCandHeader > *record, int trkidx, int &shwidx)
bool LoadLargestShowerFromEvent (NtpSREvent *event, RecRecordImp< RecCandHeader > *record, int &shwidx)
void SetDisplay (Int_t fd)
void Draw (TPad *pad)
void Print (TPad *pad)

Private Attributes

FracVarfFracVar
Int_t fDisplay
TNtuple * display

Static Private Attributes

static TMultiLayerPerceptron * fneuralNet = 0

Detailed Description

Definition at line 13 of file FracVarAna.h.


Constructor & Destructor Documentation

FracVarAna::FracVarAna ( FracVar fv  ) 

Definition at line 76 of file FracVarAna.cxx.

References display, exit(), and fneuralNet.

00076                                  :
00077   fFracVar(fv),
00078   fDisplay(0)
00079 {
00080   display = new TNtuple("display","display","tpos:z:ph:uv:core");
00081   static bool first = true;
00082                                                                                 
00083   if(first){
00084     char *srt_dir = getenv("SRT_PRIVATE_CONTEXT");
00085     char annfile[10000];
00086     sprintf(annfile,"%s/NueAna/data/fvann032007.root",srt_dir);
00087     ifstream Test(annfile);
00088     if (!Test){
00089       srt_dir = getenv("SRT_PUBLIC_CONTEXT");
00090       sprintf(annfile,"%s/NueAna/data/fvann032007.root",srt_dir);
00091       ifstream Test_again(annfile);
00092       if (!Test_again){
00093         cout<<"Couldn't find ANN object for FracVarAna, blame Tingjun"<<endl;
00094         exit(0);
00095       }
00096     }
00097                                                                                 
00098     static TFile *f = TFile::Open(annfile);
00099     fneuralNet = (TMultiLayerPerceptron*) f->Get("mlp");
00100     cout<<"Reading FracVar ann from : "<<annfile<<endl;
00101     first = false;
00102   }
00103                                                                                 
00104   if (!fneuralNet) {
00105     cout<<"Couldn't find ANN object for FracVarAna, blame Tingjun"<<endl;
00106     exit(0);
00107   }
00108   
00109 }

FracVarAna::~FracVarAna (  )  [virtual]

Definition at line 111 of file FracVarAna.cxx.

References display.

00112 {
00113   if(display){
00114     delete display;
00115     display=0;
00116   }
00117 
00118 }


Member Function Documentation

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

Implements NueAnaBase.

Definition at line 120 of file FracVarAna.cxx.

References FracVar::dis2stp, display, NueAnaBase::evtstp0mip, NueAnaBase::evtstp1mip, fDisplay, fFracVar, fneuralNet, FracVar::fract_10_counters, FracVar::fract_12_counters, FracVar::fract_14_counters, FracVar::fract_16_counters, FracVar::fract_18_counters, FracVar::fract_1_plane, FracVar::fract_20_counters, FracVar::fract_2_counters, FracVar::fract_2_planes, FracVar::fract_3_planes, FracVar::fract_4_counters, FracVar::fract_4_planes, FracVar::fract_5_planes, FracVar::fract_6_counters, FracVar::fract_6_planes, FracVar::fract_8_counters, FracVar::fract_asym, FracVar::fract_road, SntpHelpers::GetEvent(), VHS::GetStrip(), SntpHelpers::GetTrack(), Msg::kWarning, MAXMSG, NtpSRPlane::n, FracVar::passcuts, NtpSRStrip::ph0, NtpSRStrip::ph1, FracVar::pid, FracVar::pid1, NtpSRStrip::plane, NtpSRTrack::plane, NtpSRStrip::planeview, FracVar::Reset(), FracVar::shw_max, FracVar::shw_npl, FracVar::shw_nstp, FracVar::shw_slp, NtpSRPulseHeight::sigcor, NueAnaBase::sigcormeu, NtpSRStrip::strip, NtpSRStrip::tpos, mlpANN::value(), FracVar::vtxph, WeightedFit(), and NtpSRStrip::z.

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

00121 {
00122   if(srobj==0){
00123     return;
00124   }
00125 
00126   if(((dynamic_cast<NtpStRecord *>(srobj))==0)&&
00127      ((dynamic_cast<NtpSRRecord *>(srobj))==0)){
00128     return;
00129   }
00130 
00131   fFracVar.Reset();
00132   if (fDisplay) display->Reset();
00133 
00134   fFracVar.passcuts = 1;
00135   NtpSREvent *event = SntpHelpers::GetEvent(evtn,srobj);
00136   
00137   if ((event->ph.sigcor+event->ph.sigcor)<2e4)  fFracVar.passcuts = 0;
00138   if ((event->ph.sigcor+event->ph.sigcor)>1e5)  fFracVar.passcuts = 0;
00139 
00140   Int_t ntrks = event->ntrack;
00141   for (int itrk = 0; itrk<ntrks; itrk++){
00142     NtpSRTrack *track = SntpHelpers::GetTrack(itrk,srobj);
00143     if (track->plane.n>13) fFracVar.passcuts = 0;
00144   }
00145 
00146   //Int_t nshws = event->nshower;
00147   if (true){ //nshws){//we should have at lease one shower in the event for the rest of the analysis
00148 
00149     if (fDisplay){
00150       for (int istp = 0; istp<event->nstrip; istp++){
00151         Int_t index = event->stp[istp];
00152         NtpSRStrip*strip = SntpHelpers::GetStrip(index,srobj);
00153         if(!strip) continue;
00154 //        if(!evtstp0mip){
00155 //          MSG("FracVarAna",Msg::kError)<<"No mip strip information"<<endl;
00156 //          continue;
00157 //        }
00158 
00159         float x1 = strip->tpos;
00160         float x2 = strip->z;
00161         float x3 = (strip->ph0.sigcor+strip->ph1.sigcor)/sigcormeu;
00162                    //evtstp0mip[index] + evtstp1mip[index];;
00163         int x4 = int(strip->planeview);
00164         display->Fill(x1,x2,x3,x4,0);
00165       }
00166     }
00167 
00168 //    //find the primary shower, now only require the maximum ph
00169 //    Int_t pshw = 0; //primary shower index in the shower array
00170 //    Float_t maxshwph = 0;
00171 //    for (int ishw = 0; ishw<nshws; ishw++){
00172 //      Int_t shwindex = event->shw[ishw];
00173 //      NtpSRShower *shower = SntpHelpers::GetShower(shwindex,srobj);
00174 //      if (shower->ph.sigcor>maxshwph){
00175 //      pshw = shwindex;
00176 //      maxshwph = shower->ph.sigcor;
00177 //      }
00178 //    }
00179 
00180     //Get the shower close to the track vertex or the biggest shower if
00181     //there is no track. Pull out from Mad.
00182 /*
00183     int track_index   = -1;
00184     int shower_index  = -1;
00185     if(ReleaseType::IsBirch(release)){
00186       if (LoadLargestTrackFromEvent(event,srobj,track_index)){
00187         if(!LoadShowerAtTrackVertex(event,srobj,track_index,shower_index)){
00188           LoadLargestShowerFromEvent(event,srobj,shower_index);
00189         }
00190       }
00191       else LoadLargestShowerFromEvent(event,srobj,shower_index);
00192     }
00193     else {//cedar ...
00194       if (ntrks) track_index = event->trk[0];
00195       //not using event->primshw, the concern is the 2GeV cut off
00196       if (nshws) shower_index = event->shw[0];
00197     }
00198     if (shower_index == -1) return; //no decent shower was found
00199 
00200     NtpSRShower *shower = SntpHelpers::GetShower(shower_index,srobj);
00201     if (!shower) return;     //no decent shower was found
00202 */
00203 
00204     if (event->plane.n<5) fFracVar.passcuts = 0;
00205     if (event->plane.n>15) fFracVar.passcuts = 0;
00206     //get rid of cosmics
00207     if (event->plane.beg>event->plane.end) return;
00208     //get rid of events with shw.ph>evt.ph
00209 //    if (event->ph.mip>event->ph.mip){
00210 //      MAXMSG("FracVarAna",Msg::kWarning,5)
00211 //      <<"event ph("<<event->ph.mip<<"mip)>event ph("<<event->ph.mip<<"mip)"<<endl;
00212 //      return;
00213 //    }
00214       
00215     float vtx_ph = 0;
00216     //Int_t shwnstp = event->nstrip;
00217     Float_t plane[14];
00218     for (int i = 0; i<14; i++){
00219       plane[i] = 0;
00220     }
00221 
00222     vector<Float_t> stpph;
00223     vector<Float_t> stpph2; //for sorting purpose
00224     vector<Float_t>::iterator p;
00225     vector<Int_t> stppl;
00226     vector<Int_t> stpstp;
00227     vector<Float_t> stpt;
00228     vector<Float_t> stpz;
00229     vector<Int_t> planev;
00230 
00231     Float_t total_ph = 0;
00232     Int_t vtxpl = event->vtx.plane;
00233 //    if(ReleaseType::IsCedar(release)){
00234 //      NtpVtxFinder vtxf;
00235 //      NtpStRecord *st = dynamic_cast<NtpStRecord *>(srobj);
00236 //      vtxf.SetTargetEvent(evtn, st);
00237 //      if(vtxf.FindVertex() > 0){
00238 //      vtxpl = vtxf.VtxPlane();
00239 //      }
00240 //    }
00241 
00242     for (Int_t istp = 0; istp<event->nstrip; istp++){
00243       Int_t index = event->stp[istp];
00244       NtpSRStrip *strip = SntpHelpers::GetStrip(index,srobj);
00245       if(!strip) continue;
00246       //float charge = strip->ph0.sigcor+strip->ph1.sigcor;
00247  
00248       float charge = 0;
00249       //if (event->stpph1mip[istp]>0) charge += event->stpph1mip[istp];
00250       //if (event->stpph0mip[istp]>0) charge += event->stpph0mip[istp];
00251 
00252       if (evtstp0mip[index] > 0) charge += evtstp0mip[index];
00253       if (evtstp1mip[index] > 0) charge += evtstp1mip[index];
00254 
00255       if (abs(strip->plane-vtxpl)<3){
00256         vtx_ph+=charge;
00257       }
00258       stpph.push_back(charge);
00259       stpph2.push_back(charge);
00260       stppl.push_back(strip->plane);
00261       stpstp.push_back(strip->strip);
00262       stpt.push_back(strip->tpos);
00263       stpz.push_back(strip->z);
00264       planev.push_back(int(strip->planeview));
00265       total_ph += charge;
00266     }
00267 
00268     Float_t maxph = 0;
00269 
00270     Float_t maxph1 = -1;
00271     Float_t maxph2 = -1;
00272     int maxpl1 = -1;
00273     int maxpl2 = -1;
00274 
00275     for(unsigned int istp=0; istp<stpph.size(); istp++){//Loop over all strips
00276       if (stppl[istp]>=vtxpl-2&&stppl[istp]<vtxpl+12)
00277         plane[stppl[istp]-vtxpl+2] += stpph[istp];//fill plane energy information
00278       if (stpph[istp]>maxph){
00279         maxph = stpph[istp];
00280         fFracVar.shw_max = stppl[istp] - vtxpl;
00281       }
00282       if (stpph[istp]>maxph1){
00283         maxph1 = stpph[istp];
00284         maxpl1 = stppl[istp];
00285       }
00286       else if (stpph[istp]>maxph2){
00287         maxph2 = stpph[istp];
00288         maxpl2 = stppl[istp];
00289       }
00290     }
00291     if (maxpl1!=-1&&maxpl2!=-1) fFracVar.dis2stp = abs(maxpl1-maxpl2);
00292     if (fFracVar.shw_max>100) fFracVar.shw_max = -2;
00293     //I don't understand why this is happening. This should be a temporary solution.
00294     if (fFracVar.shw_max<-100) fFracVar.shw_max = -2;
00295 
00296     if (fFracVar.shw_max<0){
00297       MAXMSG("FracVarAna",Msg::kWarning,5)
00298         <<"fFracVar.shw_max<0"<<endl;
00299       fFracVar.shw_max = 0;
00300     }
00301 
00302     sort(stpph2.begin(), stpph2.end());
00303     
00304     Float_t sum_1_plane = 0,sum_2_planes = 0,sum_3_planes = 0;
00305     Float_t sum_4_planes = 0,sum_5_planes = 0,sum_6_planes = 0;
00306     Float_t sum_2_counts = 0, sum_4_counts = 0, sum_6_counts = 0;
00307     Float_t sum_8_counts = 0, sum_10_counts = 0, sum_12_counts = 0;
00308     Float_t sum_14_counts = 0, sum_16_counts = 0, sum_18_counts = 0;
00309     Float_t sum_20_counts = 0;
00310 
00311     //calculate fFracVar.fract_1_plane ...
00312     
00313     if (total_ph>0.){
00314       for (Int_t i = 0; i<14; i++){
00315         if (i<14){
00316           sum_1_plane = plane[i];
00317           if (sum_1_plane/total_ph > fFracVar.fract_1_plane) {
00318             fFracVar.fract_1_plane = sum_1_plane/total_ph;
00319           }
00320         }
00321         if (i<13){
00322           sum_2_planes = plane[i]+plane[i+1];
00323           if (sum_2_planes/total_ph > fFracVar.fract_2_planes) 
00324             fFracVar.fract_2_planes = sum_2_planes/total_ph;
00325         }
00326         if (i<12){
00327           sum_3_planes = plane[i]+plane[i+1]+plane[i+2];
00328           if (sum_3_planes/total_ph > fFracVar.fract_3_planes) 
00329             fFracVar.fract_3_planes = sum_3_planes/total_ph;
00330         }
00331         if (i<11){
00332           sum_4_planes = plane[i]+plane[i+1]+plane[i+2]+plane[i+3];
00333           if (sum_4_planes/total_ph > fFracVar.fract_4_planes) 
00334             fFracVar.fract_4_planes = sum_4_planes/total_ph;
00335         }
00336         if (i<10) {
00337           sum_5_planes = plane[i]+plane[i+1]+plane[i+2]+plane[i+3]+plane[i+4];
00338           if (sum_5_planes/total_ph > fFracVar.fract_5_planes) 
00339             fFracVar.fract_5_planes = sum_5_planes/total_ph;
00340         }
00341         if (i<9) {
00342           sum_6_planes = plane[i]+plane[i+1]+plane[i+2]+plane[i+3]+plane[i+4]+plane[i+5];
00343           if (sum_6_planes/total_ph > fFracVar.fract_6_planes) 
00344             fFracVar.fract_6_planes = sum_6_planes/total_ph;
00345         }
00346       }
00347     }//if(total_ph>0)
00348     
00349     p = stpph2.end();
00350     while (p!=stpph2.begin()&&p>stpph2.end()-20){
00351       p--;
00352       if (p>=stpph2.end()-2) sum_2_counts += *p;
00353       if (p>=stpph2.end()-4) sum_4_counts += *p;
00354       if (p>=stpph2.end()-6) sum_6_counts += *p;
00355       if (p>=stpph2.end()-8) sum_8_counts += *p;
00356       if (p>=stpph2.end()-10) sum_10_counts += *p;
00357       if (p>=stpph2.end()-12) sum_12_counts += *p;
00358       if (p>=stpph2.end()-14) sum_14_counts += *p;
00359       if (p>=stpph2.end()-16) sum_16_counts += *p;
00360       if (p>=stpph2.end()-18) sum_18_counts += *p;
00361       if (p>=stpph2.end()-20) sum_20_counts += *p;
00362       
00363       //cout<<*p<<endl;
00364     }
00365     
00366     if (total_ph>0) {
00367       fFracVar.fract_2_counters = sum_2_counts/total_ph;
00368       fFracVar.fract_4_counters = sum_4_counts/total_ph;
00369       fFracVar.fract_6_counters = sum_6_counts/total_ph;
00370       fFracVar.fract_8_counters = sum_8_counts/total_ph;
00371       fFracVar.fract_10_counters = sum_10_counts/total_ph;
00372       fFracVar.fract_12_counters = sum_12_counts/total_ph;
00373       fFracVar.fract_14_counters = sum_14_counts/total_ph;
00374       fFracVar.fract_16_counters = sum_16_counts/total_ph;
00375       fFracVar.fract_18_counters = sum_18_counts/total_ph;
00376       fFracVar.fract_20_counters = sum_20_counts/total_ph;
00377     }
00378 
00379     vector<Double_t> upos;
00380     vector<Double_t> zupos;
00381     vector<Double_t> vpos;
00382     vector<Double_t> zvpos;
00383     vector<Double_t> ucharge;
00384     vector<Double_t> vcharge;
00385     vector<Int_t> ustrip;
00386     vector<Int_t> uplane;
00387     vector<Int_t> vstrip;
00388     vector<Int_t> vplane;
00389     
00390     vector<Double_t> ufit;
00391     vector<Double_t> zufit;
00392     vector<Double_t> vfit;
00393     vector<Double_t> zvfit;
00394     vector<Double_t> ufitcharge;
00395     vector<Double_t> vfitcharge;
00396 
00397     vector<Int_t> ifitu;
00398     vector<Int_t> ifitv;
00399     
00400     vector<Double_t> ushwstrip;
00401     vector<Double_t> ushwplane;
00402     vector<Double_t> vshwstrip;
00403     vector<Double_t> vshwplane;
00404     vector<Int_t> shwplane;
00405     
00406     Double_t total_charge = 0.;
00407     
00408     Float_t toler1 = 3*0.0412;
00409     Float_t toler2 = 1.5*0.0412;
00410     
00411     Double_t thr1 = 0.05; //threshold to determine shw_beg and shw_end
00412     
00413     Int_t begplane = event->plane.beg;
00414     Int_t endplane = event->plane.end;
00415 
00416     if (begplane>endplane) {
00417       MAXMSG("FracVarAna",Msg::kWarning,5)
00418         <<"event->plane.beg= "<<begplane
00419         <<"event->plane.end= "<<endplane<<endl;
00420       return;
00421     }
00422     vector<Double_t> planeph(endplane-begplane+1);
00423     for (int i = 0; i<endplane-begplane+1; i++){
00424       planeph[i] = 0;
00425     }
00426     
00427     for (unsigned int i = 0; i<stpph.size(); i++){
00428       if (planev[i] == 2){//u view
00429         upos.push_back(stpt[i]);
00430         zupos.push_back(stpz[i]);
00431         ustrip.push_back(stpstp[i]);
00432         uplane.push_back(stppl[i]);
00433         ucharge.push_back(stpph[i]);
00434         ifitu.push_back(0);
00435       }
00436       
00437       if (planev[i] == 3){//v view
00438         vpos.push_back(stpt[i]);
00439         zvpos.push_back(stpz[i]);
00440         vstrip.push_back(stpstp[i]);
00441         vplane.push_back(stppl[i]);
00442         vcharge.push_back(stpph[i]);
00443         ifitv.push_back(0);
00444       }
00445       
00446       if (stppl[i]>=begplane&&stppl[i]<=endplane){
00447         planeph[stppl[i]-begplane] += stpph[i];
00448       }
00449       
00450       total_charge += stpph[i];
00451       
00452     }
00453    
00454     // calculate the u and v charge, then asymmetry
00455     Float_t ucharge_tot = 0, vcharge_tot = 0;
00456   
00457     for (unsigned int i = 0; i<ucharge.size(); i++) {
00458       ucharge_tot += ucharge[i];
00459     }   
00460 
00461     for (unsigned int i = 0; i<vcharge.size(); i++) {
00462       vcharge_tot += vcharge[i];
00463     }
00464 
00465     if (ucharge_tot + vcharge_tot > 0) fFracVar.fract_asym = TMath::Abs(ucharge_tot-vcharge_tot)/(ucharge_tot+vcharge_tot); 
00466  
00467     //find event range
00468     Int_t shw_beg = 485;
00469     Int_t shw_end = 1;
00470     
00471     for (Int_t ipl = 0; ipl<endplane-begplane+1; ipl++){
00472       if (planeph[ipl]>total_charge*thr1/2 && ipl+begplane<shw_beg
00473           && ipl<endplane-begplane && planeph[ipl+1]>total_charge*thr1)
00474         shw_beg = ipl + begplane;
00475       if (planeph[ipl]>total_charge*thr1/2&&ipl+begplane>shw_end)
00476         shw_end = ipl + begplane;
00477     }
00478     
00479     if (shw_beg == 485) shw_beg = begplane;
00480     if (shw_end == 1) shw_end = endplane;
00481     //cout<<begplane<<" "<<endplane<<endl;
00482     //cout<<shw_beg<<" "<<shw_end<<endl;
00483     
00484     for (unsigned int i = 0; i<ifitu.size(); i++){
00485       if (uplane[i]>=shw_beg&&uplane[i]<=shw_end){
00486         ufit.push_back(upos[i]);
00487         zufit.push_back(zupos[i]);
00488         ufitcharge.push_back(ucharge[i]);
00489       }
00490     }
00491     
00492     for (unsigned int i = 0; i<ifitv.size(); i++){
00493       if (vplane[i]>=shw_beg&&vplane[i]<=shw_end){
00494         vfit.push_back(vpos[i]);
00495         zvfit.push_back(zvpos[i]);
00496         vfitcharge.push_back(vcharge[i]);
00497       }
00498     }
00499     
00500     //perform weighted fit of u, v vs z positions
00501     Double_t uparm[2];
00502     Double_t vparm[2];
00503     Int_t fituok = 0;
00504     Int_t fitvok = 0;
00505     
00506     if (ufit.size()) {
00507       fituok = WeightedFit(ufit.size(),&zufit[0],&ufit[0],&ufitcharge[0],&uparm[0]);
00508       //cout<<"u "<<fituok<<" "<<uparm[0]<<" "<<uparm[1]<<endl;
00509     }
00510     if (vfit.size()) {
00511       fitvok = WeightedFit(vfit.size(),&zvfit[0],&vfit[0],&vfitcharge[0],&vparm[0]);
00512       //cout<<"v "<<fitvok<<" "<<vparm[0]<<" "<<vparm[1]<<endl;
00513     }
00514     
00515     ufit.clear();
00516     zufit.clear();
00517     ufitcharge.clear();
00518     vfit.clear();
00519     zvfit.clear();
00520     vfitcharge.clear();
00521     
00522     for (unsigned int i = 0; i<ifitu.size(); i++){
00523         if (!fituok&&TMath::Abs((upos[i]-(uparm[0]+zupos[i]*uparm[1]))*cos(atan(uparm[1])))<toler1&&uplane[i]>=shw_beg&&uplane[i]<=shw_end){
00524         ifitu[i] = 1;
00525         ushwstrip.push_back(ustrip[i]+0.5);
00526         ushwplane.push_back(uplane[i]+0.5);
00527         ufit.push_back(upos[i]);
00528         zufit.push_back(zupos[i]);
00529         ufitcharge.push_back(ucharge[i]);
00530       }
00531     }
00532     
00533     for (unsigned int i = 0; i<ifitv.size(); i++){
00534         if (!fitvok&&TMath::Abs((vpos[i]-(vparm[0]+zvpos[i]*vparm[1]))*cos(atan(vparm[1])))<toler1&&vplane[i]>=shw_beg&&vplane[i]<=shw_end){
00535         ifitv[i] = 1;
00536         vshwstrip.push_back(vstrip[i]+0.5);
00537         vshwplane.push_back(vplane[i]+0.5);
00538         vfit.push_back(vpos[i]);
00539         zvfit.push_back(zvpos[i]);
00540         vfitcharge.push_back(vcharge[i]);
00541       }
00542     }
00543     
00544     for(int itr = 0; itr<2; itr++){
00545       
00546       if (ufit.size()) {
00547         fituok = WeightedFit(ufit.size(),&zufit[0],&ufit[0],&ufitcharge[0],&uparm[0]);
00548         //cout<<"u "<<fituok<<" "<<uparm[0]<<" "<<uparm[1]<<endl;
00549       }
00550       if (vfit.size()) {
00551         fitvok = WeightedFit(vfit.size(),&zvfit[0],&vfit[0],&vfitcharge[0],&vparm[0]);
00552         //cout<<"v "<<fitvok<<" "<<vparm[0]<<" "<<vparm[1]<<endl;
00553       }
00554       
00555       ushwstrip.clear();
00556       ushwplane.clear();
00557       vshwstrip.clear();
00558       vshwplane.clear();
00559       ufit.clear();
00560       zufit.clear();
00561       ufitcharge.clear();
00562       vfit.clear();
00563       zvfit.clear();
00564       vfitcharge.clear();
00565       
00566       
00567       for (unsigned int i = 0; i<ifitu.size(); i++){
00568         ifitu[i] = 0;
00569           if (!fituok&&TMath::Abs((upos[i]-(uparm[0]+zupos[i]*uparm[1]))*cos(atan(uparm[1])))<toler2&&uplane[i]>=shw_beg&&uplane[i]<=shw_end){
00570           ifitu[i] = 1;
00571           ushwstrip.push_back(ustrip[i]+0.5);
00572           ushwplane.push_back(uplane[i]+0.5);
00573           ufit.push_back(upos[i]);
00574           zufit.push_back(zupos[i]);
00575           ufitcharge.push_back(ucharge[i]);
00576           if (itr==1) {
00577             shwplane.push_back(uplane[i]);
00578             display->Fill(upos[i],zupos[i],100000,2,1);
00579           }
00580         }
00581       }
00582       
00583       for (unsigned int i = 0; i<ifitv.size(); i++){
00584         ifitv[i] = 0;
00585           if (!fitvok&&TMath::Abs((vpos[i]-(vparm[0]+zvpos[i]*vparm[1]))*cos(atan(vparm[1])))<toler2&&vplane[i]>=shw_beg&&vplane[i]<=shw_end){
00586           ifitv[i] = 1;
00587           vshwstrip.push_back(vstrip[i]+0.5);
00588           vshwplane.push_back(vplane[i]+0.5);
00589           vfit.push_back(vpos[i]);
00590           zvfit.push_back(zvpos[i]);
00591           vfitcharge.push_back(vcharge[i]);
00592           if (itr==1) {
00593             shwplane.push_back(vplane[i]);
00594             display->Fill(vpos[i],zvpos[i],100000,3,1);
00595           }
00596         }
00597       }
00598     }
00599     
00600     Double_t shwcoreph = 0;
00601     fFracVar.shw_nstp = 0;
00602     for (unsigned int i = 0; i<ufitcharge.size(); i++){
00603       shwcoreph += ufitcharge[i];
00604       fFracVar.shw_nstp ++;
00605     }
00606     for (unsigned int i = 0; i<vfitcharge.size(); i++){
00607       shwcoreph += vfitcharge[i];
00608       fFracVar.shw_nstp ++;
00609     }
00610     if (total_ph) {
00611       fFracVar.fract_road = shwcoreph/total_ph;
00612       fFracVar.vtxph = vtx_ph/total_ph;
00613     }
00614     Int_t shw_begpl = 1000;
00615     Int_t shw_endpl = 0;
00616     for (unsigned int i = 0; i<shwplane.size(); i++){
00617       if (shw_begpl>shwplane[i]) shw_begpl = shwplane[i];
00618       if (shw_endpl<shwplane[i]) shw_endpl = shwplane[i];
00619    }
00620     if (shw_endpl>=shw_begpl) fFracVar.shw_npl = shw_endpl - shw_begpl + 1;
00621 
00622     double slope = -1;
00623     if (!fituok&&!fitvok&&TMath::Abs(uparm[1])<100&&TMath::Abs(vparm[1])<100){
00624       slope = sqrt(pow(uparm[1],2)+pow(vparm[1],2));
00625     }
00626     fFracVar.shw_slp = slope;
00627     //calculate ANN pid
00628     mlpANN ann;
00629     Double_t params[14];
00630     params[0] = fFracVar.fract_1_plane;
00631     params[1] = fFracVar.fract_2_planes;
00632     params[2] = fFracVar.fract_3_planes;
00633     params[3] = fFracVar.fract_4_planes;
00634     params[4] = fFracVar.fract_5_planes;
00635     params[5] = fFracVar.fract_6_planes;
00636     params[6] = fFracVar.fract_2_counters;
00637     params[7] = fFracVar.fract_4_counters;
00638     params[8] = fFracVar.fract_6_counters;
00639     params[9] = fFracVar.fract_8_counters;
00640     params[10] = fFracVar.fract_10_counters;
00641     params[11] = fFracVar.fract_12_counters;
00642     params[12] = fFracVar.fract_road;
00643     params[13] = fFracVar.shw_max;
00644     
00645     fFracVar.pid = ann.value(0,params[0],params[1],params[2],params[3],params[4],params[5],params[6],params[7],params[8],params[9],params[10],params[11],params[12],params[13]);
00646     fFracVar.pid1 = fneuralNet->Evaluate(0,params);
00647     
00648   }//if (nshws)
00649   else {
00650     fFracVar.passcuts = 0;
00651   }
00652 }

void FracVarAna::Draw ( TPad *  pad  ) 

Definition at line 744 of file FracVarAna.cxx.

References display.

Referenced by NueDisplayModule::Analyze().

00744                               {
00745   pad->cd(1);
00746   //TLatex *t1 = new TLatex(0.5,0.5,Form("%.3f",fFracVar.fract_road));
00747   //t1->Draw();
00748   display->Draw("tpos:z","ph*(uv==2&&!core)","colz");
00749   display->Draw("tpos:z","ph*(uv==2&&core)","box same");
00750   pad->cd(2);
00751   display->Draw("tpos:z","ph*(uv==3&&!core)","colz");
00752   display->Draw("tpos:z","ph*(uv==3&&core)","box same");
00753   pad->Modified();
00754 }

bool FracVarAna::LoadLargestShowerFromEvent ( NtpSREvent event,
RecRecordImp< RecCandHeader > *  record,
int &  shwidx 
)

Definition at line 722 of file FracVarAna.cxx.

References SntpHelpers::GetShower(), NtpSRStripPulseHeight::gev, NtpSRShower::ph, and NtpSREvent::shw.

00725 {
00726   bool status = false;
00727   float largest_ph = 0;
00728   for (int i=0; i<event->nshower; i++){
00729     NtpSRShower *shower = SntpHelpers::GetShower(event->shw[i],record);
00730     if (!shower) continue;
00731     if (shower->ph.gev>largest_ph){
00732       shwidx = event->shw[i];
00733       largest_ph = shower->ph.gev;
00734     }
00735   }
00736   if (largest_ph>0){
00737     status = true;
00738   }
00739   return status;
00740 }

bool FracVarAna::LoadLargestTrackFromEvent ( NtpSREvent event,
RecRecordImp< RecCandHeader > *  record,
int &  trkidx 
)

Definition at line 657 of file FracVarAna.cxx.

References NtpSRPlane::beg, NtpSRPlane::end, SntpHelpers::GetTrack(), NtpSRTrack::plane, and NtpSREvent::trk.

00660 {
00661   bool status = false;
00662   float longest_trk = 0;
00663   for (int i=0; i<event->ntrack; i++){
00664     NtpSRTrack *track = SntpHelpers::GetTrack(event->trk[i],record);
00665     if (!track) continue;
00666     int trklen = abs(track->plane.end-track->plane.beg)+1;
00667     if (trklen>longest_trk){
00668       trkidx = event->trk[i];
00669       longest_trk = trklen;
00670     }
00671   }
00672   if (longest_trk>0){
00673     status = true;
00674   }
00675   return status;
00676 }

bool FracVarAna::LoadShowerAtTrackVertex ( NtpSREvent event,
RecRecordImp< RecCandHeader > *  record,
int  trkidx,
int &  shwidx 
)

Definition at line 680 of file FracVarAna.cxx.

References SntpHelpers::GetShower(), SntpHelpers::GetTrack(), NtpSRStripPulseHeight::gev, NtpSRShowerPulseHeight::linCCgev, NtpSRShower::ph, NtpSREvent::shw, NtpSRShower::shwph, NtpSRShower::vtx, NtpSRTrack::vtx, and NtpSRVertex::z.

00684 {
00685   bool status = false;
00686   NtpSRTrack *track = SntpHelpers::GetTrack(trkidx,record);
00687   if (!track) return false;
00688   float closest_ph = 0.;
00689   float closest_vtx = 1000.;
00690   const float close_enough = 7*0.06; //about 1 interaction length
00691   for (int i=0; i<event->nshower; i++){
00692     NtpSRShower *shower = SntpHelpers::GetShower(event->shw[i],record);
00693     if (!shower) continue;
00694     float vtx_sep = TMath::Abs(shower->vtx.z-track->vtx.z);
00695 
00696     //if this is the closest shower
00697     if (vtx_sep<closest_vtx){
00698       shwidx = event->shw[i];
00699       closest_vtx = vtx_sep;
00700       if (shower->shwph.linCCgev>0) closest_ph = shower->shwph.linCCgev;
00701       else closest_ph = shower->ph.gev;
00702     }
00703     // if this isn't the closest shower but it's within close_enough m of the
00704     // track vertex and has a larger pulseheight
00705     else if(vtx_sep<close_enough && shower->ph.gev>closest_ph){
00706       shwidx = event->shw[i];
00707       closest_vtx = vtx_sep;
00708 
00709       if (shower->shwph.linCCgev>0) closest_ph = shower->shwph.linCCgev;
00710       else closest_ph = shower->ph.gev;
00711     }
00712   }
00713   if (closest_vtx<close_enough){
00714     status = true;
00715   }
00716   else shwidx = -1;
00717   return status;
00718 }

void FracVarAna::Print ( TPad *  pad  ) 

Definition at line 756 of file FracVarAna.cxx.

00756                                    {
00757 }

void FracVarAna::SetDisplay ( Int_t  fd  )  [inline]

Definition at line 32 of file FracVarAna.h.

References fDisplay.

Referenced by NueDisplayModule::NueDisplayModule().

00032 {fDisplay = fd;};


Member Data Documentation

TNtuple* FracVarAna::display [private]

Definition at line 39 of file FracVarAna.h.

Referenced by Analyze(), Draw(), FracVarAna(), and ~FracVarAna().

Int_t FracVarAna::fDisplay [private]

Definition at line 38 of file FracVarAna.h.

Referenced by Analyze(), and SetDisplay().

Definition at line 37 of file FracVarAna.h.

Referenced by Analyze().

TMultiLayerPerceptron * FracVarAna::fneuralNet = 0 [static, private]

Definition at line 40 of file FracVarAna.h.

Referenced by Analyze(), and FracVarAna().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1