#include <FracVarAna.h>
Inheritance diagram for FracVarAna:
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 | |
| FracVar & | fFracVar |
| Int_t | fDisplay |
| TNtuple * | display |
Static Private Attributes | |
| TMultiLayerPerceptron * | fneuralNet = 0 |
|
|
Definition at line 76 of file FracVarAna.cxx. References display, 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 }
|
|
|
Definition at line 111 of file FracVarAna.cxx. References display. 00112 {
00113 if(display){
00114 delete display;
00115 display=0;
00116 }
00117
00118 }
|
|
||||||||||||
|
Implements NueAnaBase. Definition at line 120 of file FracVarAna.cxx. References abs(), NtpSRPlane::beg, FracVar::dis2stp, display, NtpSRPlane::end, 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(), SntpHelpers::GetStrip(), SntpHelpers::GetTrack(), MAXMSG, NtpSRPlane::n, NtpSREvent::nstrip, NtpSREvent::ntrack, FracVar::passcuts, NtpSREvent::ph, NtpSRStrip::ph0, NtpSRStrip::ph1, FracVar::pid, FracVar::pid1, NtpSRStrip::plane, NtpSRVertex::plane, NtpSREvent::plane, NtpSRTrack::plane, NtpSRStrip::planeview, FracVar::Reset(), FracVar::shw_max, FracVar::shw_npl, FracVar::shw_nstp, FracVar::shw_slp, NtpSRPulseHeight::sigcor, NtpSREvent::stp, NtpSRStrip::strip, NtpSRStrip::tpos, mlpANN::value(), NtpSREvent::vtx, 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 }
|
|
|
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 }
|
|
||||||||||||||||
|
Definition at line 722 of file FracVarAna.cxx. References SntpHelpers::GetShower(), NtpSRStripPulseHeight::gev, NtpSREvent::nshower, 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 }
|
|
||||||||||||||||
|
Definition at line 657 of file FracVarAna.cxx. References abs(), NtpSRPlane::beg, NtpSRPlane::end, SntpHelpers::GetTrack(), NtpSREvent::ntrack, 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 }
|
|
||||||||||||||||||||
|
Definition at line 680 of file FracVarAna.cxx. References SntpHelpers::GetShower(), SntpHelpers::GetTrack(), NtpSRStripPulseHeight::gev, NtpSRShowerPulseHeight::linCCgev, NtpSREvent::nshower, NtpSRShower::ph, NtpSREvent::shw, NtpSRShower::shwph, NtpSRTrack::vtx, NtpSRShower::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 }
|
|
|
Definition at line 756 of file FracVarAna.cxx. 00756 {
00757 }
|
|
|
Definition at line 32 of file FracVarAna.h. References fDisplay. Referenced by NueDisplayModule::NueDisplayModule(). 00032 {fDisplay = fd;};
|
|
|
Definition at line 39 of file FracVarAna.h. Referenced by Analyze(), Draw(), FracVarAna(), and ~FracVarAna(). |
|
|
Definition at line 38 of file FracVarAna.h. Referenced by SetDisplay(). |
|
|
Definition at line 37 of file FracVarAna.h. Referenced by Analyze(). |
|
|
Definition at line 27 of file FracVarAna.cxx. Referenced by Analyze(), and FracVarAna(). |
1.3.9.1