00001 #include "AlgAtNuReco.h"
00002
00003 #include "Algorithm/AlgConfig.h"
00004 #include "Algorithm/AlgFactory.h"
00005 #include "Algorithm/AlgHandle.h"
00006
00007 #include "CandData/CandRecord.h"
00008 #include "Candidate/CandContext.h"
00009 #include "JobControl/JobCModuleRegistry.h"
00010 #include "JobControl/JobCommand.h"
00011 #include "MessageService/MsgService.h"
00012 #include "MinosObjectMap/MomNavigator.h"
00013 #include "UgliGeometry/UgliGeomHandle.h"
00014 #include "Validity/VldContext.h"
00015 #include "Plex/PlexPlaneId.h"
00016
00017 #include "RecoBase/CandSliceListHandle.h"
00018 #include "RecoBase/CandSliceHandle.h"
00019 #include "RecoBase/CandStripListHandle.h"
00020 #include "RecoBase/CandStripHandle.h"
00021 #include "CandDigit/CandDeMuxDigitHandle.h"
00022 #include "CandDigit/CandDeMuxDigitListHandle.h"
00023
00024 #include "CandShowerAtNuHandle.h"
00025 #include "CandShowerAtNuListHandle.h"
00026 #include "CandTrackAtNuHandle.h"
00027 #include "CandTrackAtNuListHandle.h"
00028
00029 #include "HitAtNu.h"
00030 #include "ClusterAtNu.h"
00031 #include "TrackSegmentAtNu.h"
00032 #include "ShowerSegmentAtNu.h"
00033 #include "TrackAtNu.h"
00034 #include "ShowerAtNu.h"
00035
00036 #include "CandAtNuRecoHandle.h"
00037
00038 #include <sys/time.h>
00039
00040
00041
00042
00043
00044 ClassImp(AlgAtNuReco)
00045
00046 CVSID("$Id: AlgAtNuReco.cxx,v 1.8 2006/05/22 16:44:42 rhatcher Exp $");
00047
00048 AlgAtNuReco::AlgAtNuReco()
00049 {
00050
00051 }
00052
00053 AlgAtNuReco::~AlgAtNuReco()
00054 {
00055
00056 }
00057
00058 void AlgAtNuReco::RunAlg(AlgConfig & ac, CandHandle &ch, CandContext &cx)
00059 {
00060
00061 MSG("AlgAtNuReco", Msg::kDebug) << "AlgAtNuReco::RunAlg(...)" << endl;
00062
00063
00064 CandAtNuRecoHandle& atnu_handle = dynamic_cast<CandAtNuRecoHandle&>(ch);
00065
00066
00067 Double_t fPeCut,fTimeWin;
00068 Int_t fVtxFlag,fObjMax,fObjMaxFlag;
00069 TString fListOutTrk,fListOutShw;
00070 Int_t NEWOBJCTR,DELOBJCTR;
00071
00072 fListOutTrk = ac.GetCharString("ListOutTrk");
00073 fListOutShw = ac.GetCharString("ListOutShw");
00074 fPeCut = ac.GetDouble("PeCut");
00075 fVtxFlag = ac.GetInt("VtxFlag");
00076 fTimeWin=100.0;
00077 fObjMax=100000;
00078 NEWOBJCTR=0; DELOBJCTR=0;
00079
00080
00081 CandSliceListHandle* slice_list = (CandSliceListHandle*)(cx.GetDataIn());
00082
00083
00084 Int_t itr;
00085 Int_t modtype,modnum,modpln,modstr;
00086 Int_t fEndPln,fEndStr;
00087 switch(slice_list->GetVldContext()->GetDetector()){
00088 case Detector::kFar:
00089 modtype=1; modnum=2; modpln=248; modstr=192; break;
00090 case Detector::kNear:
00091 modtype=2; modnum=1; modpln=282; modstr=96; break;
00092 case Detector::kCalib:
00093 modtype=3; modnum=1; modpln=60; modstr=24; break;
00094 default:
00095 modtype=-1; modnum=0; modpln=0; break;
00096 }
00097 fEndPln = modnum*(modpln+1); fEndStr = modstr;
00098
00099 MSG("AlgAtNuReco", Msg::kDebug)
00100 << " *** PeCut=" << fPeCut << " VtxFlag=" << fVtxFlag << endl;
00101
00102 MSG("AlgAtNuReco", Msg::kDebug)
00103 << " *** ModType=" << modtype << " ModNum=" << modnum << " "
00104 << " ModPln=" << modpln << " ModStr=" << modstr << endl;
00105
00106 MSG("AlgAtNuReco", Msg::kDebug)
00107 << " *** EndPln=" << fEndPln << " fEndStr=" << fEndStr << endl;
00108
00109
00110
00111 Double_t dtime;
00112 Double_t fTime_TrackAtNu,fTime_ShowerAtNu,fTime_AtNuReco;
00113 fTime_TrackAtNu=0.0; fTime_ShowerAtNu=0.0; fTime_AtNuReco=0.0;
00114
00115 struct timeval tpbefore;
00116 struct timeval tpafter;
00117
00118 struct timeval atnubefore;
00119 struct timeval atnuafter;
00120
00121 gettimeofday(&atnubefore,0);
00122
00123
00124 TIter sliceitr(slice_list->GetDaughterIterator());
00125 while(CandSliceHandle* slice = dynamic_cast<CandSliceHandle*>(sliceitr())){
00126 atnu_handle.AddDaughterLink(*slice);
00127 fObjMaxFlag=0;
00128
00129 MSG("AlgAtNuReco", Msg::kDebug)
00130 << " *** SLICE (" << atnu_handle.GetNDaughters() << ") strips= " << slice->GetNDaughters() << endl;
00131
00132
00133
00134
00135
00136
00137 gettimeofday(&tpbefore,0);
00138
00139
00140 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of hits *** " << endl;
00141 Int_t i,j,k,l;
00142 Bool_t cont;
00143 Double_t chg,scr=0.;
00144 Int_t id,pln,npln,mod,vuw,ctr,xtalk,assoc;
00145 Int_t Nrealhits=0,Nxtalkhits=0;
00146 TObjArray tmp;
00147 TIter stritr(slice->GetDaughterIterator());
00148 while(CandStripHandle* strip = dynamic_cast<CandStripHandle*>(stritr())){
00149 pln=strip->GetPlane(); chg=strip->GetCharge(); xtalk=0;
00150 if(pln>0 && pln<fEndPln && pln<500){
00151 TIter digitr(strip->GetDaughterIterator());
00152 while(CandDeMuxDigitHandle* digit = (CandDeMuxDigitHandle*)(digitr())){
00153 if( ( digit->GetDeMuxDigitFlagWord()<8 )
00154 && (digit->GetDeMuxDigitFlagWord() & CandDeMuxDigit::kXTalk)==(CandDeMuxDigit::kXTalk) ){
00155 xtalk=1;
00156 }
00157 }
00158 HitAtNu* hit = new HitAtNu(strip); hitbnk.Add(hit); NEWOBJCTR++;
00159 if(!xtalk && chg>fPeCut){ hitplnbnk[pln].Add(hit); Nrealhits++; }
00160 else{ hit->SetXtalkFlag(1); xtalkplnbnk[pln].Add(hit); Nxtalkhits++; }
00161 }
00162 atnu_handle.AddDaughterLink(*strip);
00163 }
00164 MSG("AlgAtNuReco", Msg::kDebug) << " *** real hits = " << Nrealhits << " xtalk hits = " << Nxtalkhits << endl;
00165
00166
00167
00168
00169
00170
00171
00172 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of clusters *** " << endl;
00173 for(i=1;i<fEndPln;i++){
00174 for(j=0;j<1+hitplnbnk[i].GetLast();j++){
00175 HitAtNu* hitj = (HitAtNu*)(hitplnbnk[i].At(j));
00176 if(hitj->GetUID()==0){
00177 ClusterAtNu* clr = new ClusterAtNu(hitj); NEWOBJCTR++;
00178 clrplnbnk[i].Add(clr); clrbnk.Add(clr);
00179 hitj->SetUID(1); ctr=1;
00180 while(ctr>0){
00181 ctr=0;
00182 for(k=0;k<1+hitplnbnk[i].GetLast();k++){
00183 HitAtNu* hitk = (HitAtNu*)(hitplnbnk[i].At(k));
00184 if(hitk->GetUID()==0 && clr->IsHitAssoc(hitk)){
00185 clr->AddHit(hitk); hitk->SetUID(1); ctr++;
00186 }
00187 }
00188 }
00189 }
00190 }
00191 }
00192
00193
00194 Int_t k0,k1,k2;
00195 Int_t nhits0,nhits1,nclrs0,nclrs1;
00196 TObjArray tmp0,tmp1,tmp2;
00197
00198 for(mod=0;mod<modnum;mod++){
00199 for(pln=0;pln<modpln+1;pln++){
00200 i=pln+mod*(modpln+1);
00201 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00202 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00203 if( clr0->GetCharge()>2.0 || clr0->GetDigits()>1.0 ){
00204 clr0->SetTrkFlag(1);
00205 }
00206 }
00207 }
00208 }
00209
00210 for(mod=0;mod<modnum;mod++){
00211 for(pln=1;pln<1+modpln;pln++){
00212 i=pln+mod*(modpln+1);
00213 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00214 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00215 nhits0=0; nhits1=0; nclrs0=0; nclrs1=0;
00216 tmp0.Clear(); tmp1.Clear();
00217 for(k=-1;k<2;k++){
00218 k0=i+2*k;
00219 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
00220 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
00221 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
00222 if(clr0->GetCharge()>1.0 && clr1->GetCharge()>1.0){
00223 assoc=clr0->IsShwAssoc(clr1);
00224 if(assoc>0){ nhits0+=1+clr1->GetHitLast(); nclrs0+=1; tmp0.Add(clr1); }
00225 if(assoc>1){ nhits1+=1+clr1->GetHitLast(); nclrs1+=1; tmp1.Add(clr1); }
00226 }
00227 }
00228 }
00229 }
00230 if(nclrs1>1 && nhits1>4){
00231 for(k=0;k<1+tmp1.GetLast();k++){
00232 ClusterAtNu* clrtmp = (ClusterAtNu*)(tmp1.At(k));
00233 clrtmp->SetShwFlag(1);
00234 }
00235 }
00236 if(nclrs0>4 && nhits0>5){
00237 for(k=0;k<1+tmp0.GetLast();k++){
00238 ClusterAtNu* clrtmp = (ClusterAtNu*)(tmp0.At(k));
00239 clrtmp->SetShwFlag(1);
00240 }
00241 }
00242 if(clr0->GetShwFlag()<1 && clr0->GetCharge()>50.0){
00243 clr0->SetShwFlag(1);
00244 }
00245 }
00246 }
00247 }
00248
00249
00250 for(mod=0;mod<modnum;mod++){
00251 for(pln=0;pln<modpln+1;pln++){
00252 i=pln+mod*(modpln+1); chg=0.0;
00253 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00254 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00255 chg+=clr0->GetCharge();
00256 }
00257 if(chg>0.0){
00258 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00259 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00260 if(1+clr1->GetHitLast()<3 && clr1->GetCharge()/chg>0.5) clr1->SetTrkPlnFlag(1);
00261 if(1+clr1->GetHitLast()>1 && clr1->GetCharge()/chg>0.5) clr1->SetShwPlnFlag(1);
00262 if(clr1->GetCharge()/chg>0.5) clr1->SetClrPlnFlag(1);
00263 }
00264 }
00265 }
00266 }
00267
00268
00269 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF HITS/CLUSTERS *** " << endl;
00270 for(i=0;i<fEndPln;i++){
00271 for(j=0;j<1+hitplnbnk[i].GetLast();j++){
00272 HitAtNu* hit = (HitAtNu*)(hitplnbnk[i].At(j));
00273 MSG("AlgAtNuReco", Msg::kVerbose)
00274 << " HIT:"
00275 << " pln=" << hit->GetPlane()
00276 << " strp=" << hit->GetStrip()
00277 << " chg=" << hit->GetCharge()
00278 << " time=" << hit->GetTime()
00279 << " dig=" << hit->GetDigits() << endl;
00280 }
00281 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00282 ClusterAtNu* clr = (ClusterAtNu*)(clrplnbnk[i].At(j));
00283 MSG("AlgAtNuReco", Msg::kVerbose)
00284 << " CLUSTER:"
00285 << " plane=" << clr->GetPlane()
00286 << " begstrip=" << clr->GetBegStrip()
00287 << " endstrip=" << clr->GetEndStrip()
00288 << " trkpln=" << clr->GetTrkPlnFlag()
00289 << " shwpln=" << clr->GetShwPlnFlag() << endl;
00290 }
00291 }
00292
00293 gettimeofday(&tpafter,0);
00294 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
00295 fTime_TrackAtNu+=0.5*dtime;
00296 fTime_ShowerAtNu+=0.5*dtime;
00297
00298
00299
00300
00301
00302
00303 gettimeofday(&tpbefore,0);
00304
00305 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of triplets *** " << endl;
00306
00307
00308 Int_t jnflag,jnflagtmp;
00309 Int_t km,kp,ktmp,ktmp1;
00310 for(mod=0;mod<modnum;mod++){
00311 for(pln=2;pln<modpln-1;pln++){
00312 i=pln+mod*(modpln+1);
00313 for(k0=0;k0<1+clrplnbnk[i].GetLast();k0++){
00314 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(k0));
00315 if(clr0->GetTrkFlag()>0){
00316 jnflag=0;
00317
00318
00319 if( 1+clrplnbnk[i-2].GetLast()>0 && 1+clrplnbnk[i+2].GetLast()>0 ){
00320 for(km=0;km<1+clrplnbnk[i-2].GetLast();km++){
00321 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-2].At(km));
00322 for(kp=0;kp<1+clrplnbnk[i+2].GetLast();kp++){
00323 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+2].At(kp));
00324 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00325 assoc=clr0->IsTrkAssoc(clrm,clrp);
00326 if(assoc>0){
00327 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00328 segplnbnk[i].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00329 if(assoc>1){
00330 clr0->SetTrkFlag(2);
00331 }
00332 }
00333 }
00334 }
00335 }
00336 }
00337
00338
00339 if( 1+clrplnbnk[i+2].GetLast()>0 ){
00340
00341
00342 if( pln>3 && 1+clrplnbnk[i-4].GetLast()>0 ){
00343 for(kp=0;kp<1+clrplnbnk[i+2].GetLast();kp++){
00344 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+2].At(kp));
00345 for(km=0;km<1+clrplnbnk[i-4].GetLast();km++){
00346 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-4].At(km));
00347 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00348 assoc = clr0->IsTrkAssoc(clrm,clrp);
00349 if(assoc>0){
00350 jnflagtmp=0;
00351 if(jnflagtmp==0 && 1+clrplnbnk[i-2].GetLast()>0){
00352 for(ktmp=0;ktmp<1+clrplnbnk[i-2].GetLast();ktmp++){
00353 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp));
00354 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00355 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrp)>0 ) jnflagtmp=1;
00356 }
00357 }
00358 if(jnflagtmp==0){
00359 for(j=0;j<1;j++){
00360 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00361 segplnbnk[i-2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00362 }
00363 if(assoc>1){
00364 clr0->SetTrkFlag(2);
00365 }
00366 }
00367 }
00368 }
00369 }
00370 }
00371 }
00372
00373
00374 if( pln>5 && 1+clrplnbnk[i-6].GetLast()>0 ){
00375 for(kp=0;kp<1+clrplnbnk[i+2].GetLast();kp++){
00376 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+2].At(kp));
00377 for(km=0;km<1+clrplnbnk[i-6].GetLast();km++){
00378 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-6].At(km));
00379 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00380 assoc = clr0->IsTrkAssoc(clrm,clrp);
00381 if(assoc>0){
00382 jnflagtmp=0;
00383 if(jnflagtmp==0 && 1+clrplnbnk[i-2].GetLast()>0){
00384 for(ktmp=0;ktmp<1+clrplnbnk[i-2].GetLast();ktmp++){
00385 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp));
00386 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00387 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrp)>0 ) jnflagtmp=1;
00388 }
00389 }
00390 if(jnflagtmp==0 && 1+clrplnbnk[i-4].GetLast()>0){
00391 for(ktmp=0;ktmp<1+clrplnbnk[i-4].GetLast();ktmp++){
00392 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-4].At(ktmp));
00393 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00394 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrp)>0 ) jnflagtmp=1;
00395 }
00396 }
00397 if(jnflagtmp==0 && 1+clrplnbnk[i-4].GetLast()>0 && 1+clrplnbnk[i-2].GetLast()>0){
00398 for(ktmp=0;ktmp<1+clrplnbnk[i-4].GetLast();ktmp++){
00399 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-4].At(ktmp));
00400 for(ktmp1=0;ktmp1<1+clrplnbnk[i-2].GetLast();ktmp1++){
00401 ClusterAtNu* clrtmp1 = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp1));
00402 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0 && clrtmp1->GetTrkFlag()>0
00403 && clrtmp->IsTrkAssoc(clrm,clrtmp1)>0 && clrtmp1->IsTrkAssoc(clrtmp,clr0)>0
00404 && clr0->IsTrkAssoc(clrtmp1,clrp)>0 ) jnflagtmp=1;
00405 }
00406 }
00407 }
00408 if(jnflagtmp==0){
00409 for(j=0;j<1;j++){
00410 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00411 segplnbnk[i-2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00412 }
00413 if(assoc>1){
00414 clr0->SetTrkFlag(2);
00415 }
00416 }
00417 }
00418 }
00419 }
00420 }
00421 }
00422 }
00423
00424
00425 if( 1+clrplnbnk[i-2].GetLast()>0 ){
00426
00427
00428 if( pln<modpln-3 && 1+clrplnbnk[i+4].GetLast()>0 ){
00429 for(km=0;km<1+clrplnbnk[i-2].GetLast();km++){
00430 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-2].At(km));
00431 for(kp=0;kp<1+clrplnbnk[i+4].GetLast();kp++){
00432 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+4].At(kp));
00433 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00434 assoc = clr0->IsTrkAssoc(clrm,clrp);
00435 if(assoc>0){
00436 jnflagtmp=0;
00437 if(jnflagtmp==0 && 1+clrplnbnk[i+2].GetLast()>0){
00438 for(ktmp=0;ktmp<1+clrplnbnk[i+2].GetLast();ktmp++){
00439 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp));
00440 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00441 && clrtmp->IsTrkAssoc(clr0,clrp)>0 && clr0->IsTrkAssoc(clrm,clrtmp)>0 ) jnflagtmp=1;
00442 }
00443 }
00444 if(jnflagtmp==0){
00445 for(j=0;j<2;j++){
00446 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00447 segplnbnk[i+2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00448 }
00449 if(assoc>1){
00450 clr0->SetTrkFlag(2);
00451 }
00452 }
00453 }
00454 }
00455 }
00456 }
00457 }
00458
00459
00460 if( pln<modpln-5 && 1+clrplnbnk[i+6].GetLast()>0 ){
00461 for(km=0;km<1+clrplnbnk[i-2].GetLast();km++){
00462 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-2].At(km));
00463 for(kp=0;kp<1+clrplnbnk[i+6].GetLast();kp++){
00464 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+6].At(kp));
00465 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00466 assoc = clr0->IsTrkAssoc(clrm,clrp);
00467 if(assoc>0){
00468 jnflagtmp=0;
00469 if(jnflagtmp==0 && 1+clrplnbnk[i+2].GetLast()>0){
00470 for(ktmp=0;ktmp<1+clrplnbnk[i+2].GetLast();ktmp++){
00471 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp));
00472 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00473 && clrtmp->IsTrkAssoc(clr0,clrp)>0 && clr0->IsTrkAssoc(clrm,clrtmp)>0 ) jnflagtmp=1;
00474 }
00475 }
00476 if(jnflagtmp==0 && 1+clrplnbnk[i+4].GetLast()>0){
00477 for(ktmp=0;ktmp<1+clrplnbnk[i+4].GetLast();ktmp++){
00478 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+4].At(ktmp));
00479 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00480 && clrtmp->IsTrkAssoc(clr0,clrp)>0 && clr0->IsTrkAssoc(clrm,clrtmp)>0 ) jnflagtmp=1;
00481 }
00482 }
00483 if(jnflagtmp==0 && 1+clrplnbnk[i+2].GetLast()>0 && 1+clrplnbnk[i+4].GetLast()>0){
00484 for(ktmp=0;ktmp<1+clrplnbnk[i+2].GetLast();ktmp++){
00485 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp));
00486 for(ktmp1=0;ktmp1<1+clrplnbnk[i+4].GetLast();ktmp1++){
00487 ClusterAtNu* clrtmp1 = (ClusterAtNu*)(clrplnbnk[i+4].At(ktmp1));
00488 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0 && clrtmp1->GetTrkFlag()>0
00489 && clr0->IsTrkAssoc(clrm,clrtmp)>0 && clrtmp->IsTrkAssoc(clr0,clrtmp1)>0
00490 && clrtmp1->IsTrkAssoc(clrtmp,clrp) ) jnflagtmp=1;
00491 }
00492 }
00493 }
00494 if(jnflagtmp==0){
00495 for(j=0;j<3;j++){
00496 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00497 segplnbnk[i+2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00498 }
00499 if(assoc>1){
00500 clr0->SetTrkFlag(2);
00501 }
00502 }
00503 }
00504 }
00505 }
00506 }
00507 }
00508 }
00509
00510
00511 if( pln>3 && 1+clrplnbnk[i-4].GetLast()>0
00512 && pln<modpln-3 && 1+clrplnbnk[i+4].GetLast()>0 ){
00513 for(km=0;km<1+clrplnbnk[i-4].GetLast();km++){
00514 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-4].At(km));
00515 for(kp=0;kp<1+clrplnbnk[i+4].GetLast();kp++){
00516 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+4].At(kp));
00517 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00518 assoc=clr0->IsTrkAssoc(clrm,clrp);
00519 if(assoc>0){
00520 jnflagtmp=0;
00521 if(jnflagtmp==0 && 1+clrplnbnk[i-2].GetLast()>0){
00522 for(ktmp=0;ktmp<1+clrplnbnk[i-2].GetLast();ktmp++){
00523 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp));
00524 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00525 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrp)>0 ) jnflagtmp=1;
00526 }
00527 }
00528 if(jnflagtmp==0 && 1+clrplnbnk[i+2].GetLast()>0){
00529 for(ktmp=0;ktmp<1+clrplnbnk[i+2].GetLast();ktmp++){
00530 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp));
00531 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00532 && clrtmp->IsTrkAssoc(clr0,clrp)>0 && clr0->IsTrkAssoc(clrm,clrtmp)>0 ) jnflagtmp=1;
00533 }
00534 }
00535 if(jnflagtmp==0 && 1+clrplnbnk[i-2].GetLast()>0 && 1+clrplnbnk[i+2].GetLast()>0){
00536 for(ktmp=0;ktmp<1+clrplnbnk[i-2].GetLast();ktmp++){
00537 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp));
00538 for(ktmp1=0;ktmp1<1+clrplnbnk[i+2].GetLast();ktmp1++){
00539 ClusterAtNu* clrtmp1 = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp1));
00540 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0 && clrtmp1->GetTrkFlag()>0
00541 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrtmp1)>0
00542 && clrtmp1->IsTrkAssoc(clr0,clrp) ) jnflagtmp=1;
00543 }
00544 }
00545 }
00546 if(jnflagtmp==0){
00547 for(j=0;j<2;j++){
00548 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00549 segplnbnk[i+2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00550 }
00551 if(assoc>1){
00552 clr0->SetTrkFlag(2);
00553 }
00554 }
00555 }
00556 }
00557 }
00558 }
00559 }
00560
00561 if( jnflag==1
00562 && 1+clr0->GetEndStrip()-clr0->GetBegStrip()<2 ){
00563 clr0->SetTrkFlag(2);
00564 }
00565
00566 }
00567 }
00568 }
00569 }
00570
00571
00572 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF TRIPLETS *** " << endl;
00573 for(i=0;i<fEndPln;i++){
00574 for(j=0;j<1+segplnbnk[i].GetLast();j++){
00575 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segplnbnk[i].At(j));
00576 ClusterAtNu* clr0 = (ClusterAtNu*)(seg->GetClrAt(0));
00577 ClusterAtNu* clr1 = (ClusterAtNu*)(seg->GetClrAt(1));
00578 ClusterAtNu* clr2 = (ClusterAtNu*)(seg->GetClrAt(2));
00579 MSG("AlgAtNuReco", Msg::kVerbose)
00580 << " (" << clr0->GetPlane() << "," << clr0->GetBegStrip() << "->" << clr0->GetEndStrip() << ") "
00581 << " (" << clr1->GetPlane() << "," << clr1->GetBegStrip() << "->" << clr1->GetEndStrip() << ") "
00582 << " (" << clr2->GetPlane() << "," << clr2->GetBegStrip() << "->" << clr2->GetEndStrip() << ") " << endl;
00583 }
00584 }
00585
00586
00587
00588
00589
00590
00591 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of track segments *** " << endl;
00592 for(mod=0;mod<modnum;mod++){
00593 for(pln=3;pln<modpln-1;pln++){
00594 i=pln+mod*(modpln+1);
00595 if(1+segplnbnk[i].GetLast()>0 && 1+segplnbnk[i+2].GetLast()>0){
00596 for(k0=0;k0<1+segplnbnk[i].GetLast();k0++){
00597 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(segplnbnk[i].At(k0));
00598 for(kp=0;kp<1+segplnbnk[i+2].GetLast();kp++){
00599 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segplnbnk[i+2].At(kp));
00600 if( seg0->IsEndAssoc(segp)>0 || segp->IsBegAssoc(seg0)>0 ){
00601 seg0->AddTriToEnd(segp); segp->AddTriToBeg(seg0);
00602 }
00603 }
00604 }
00605 }
00606 }
00607 }
00608
00609
00610
00611 Int_t jnflagm,jnflagp;
00612 TObjArray tmpm,tmpp;
00613 for(mod=0;mod<modnum;mod++){
00614 for(pln=3;pln<modpln-1;pln++){
00615 i=pln+mod*(modpln+1);
00616
00617
00618 tmpm.Clear();
00619 for(j=0;j<1+segplnbnk[i].GetLast();j++){
00620 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(segplnbnk[i].At(j));
00621 if(1+segm->GetTriEndLast()>0){
00622 tmpm.Add(segm);
00623 if(1+segm->GetTriBegLast()>0){
00624 segm->SetTmpTrkFlag(1);
00625 for(k=0;k<1+segm->GetTriBegLast();k++){
00626 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segm->GetTriBegAt(k));
00627 if( segm->GetTmpTrkFlag()<2 && segtmp->GetBegPlane()<segm->GetBegPlane() ){
00628 segm->SetTmpTrkFlag(2);
00629 }
00630 }
00631 }
00632 }
00633 }
00634
00635
00636 tmpp.Clear();
00637 for(j=0;j<1+segplnbnk[i+2].GetLast();j++){
00638 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segplnbnk[i+2].At(j));
00639 if(1+segp->GetTriBegLast()>0){
00640 tmpp.Add(segp);
00641 if(1+segp->GetTriEndLast()>0){
00642 segp->SetTmpTrkFlag(1);
00643 for(k=0;k<1+segp->GetTriEndLast();k++){
00644 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segp->GetTriEndAt(k));
00645 if( segp->GetTmpTrkFlag()<2 && segtmp->GetEndPlane()>segp->GetEndPlane() ){
00646 segp->SetTmpTrkFlag(2);
00647 }
00648 }
00649 }
00650 }
00651 }
00652
00653
00654 for(j=0;j<1+tmpp.GetLast();j++){
00655 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(tmpp.At(j));
00656 jnflag=0;
00657 for(k=0;k<1+segp->GetTriBegLast();k++){
00658 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segp->GetTriBegAt(k));
00659 if(segtmp->GetTmpTrkFlag()>1) jnflag=1;
00660 }
00661 if(jnflag>0){
00662 for(k=0;k<1+segp->GetTriBegLast();k++){
00663 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segp->GetTriBegAt(k));
00664 if(segtmp->GetTmpTrkFlag()>0 && segtmp->GetTmpTrkFlag()<2) segtmp->SetTmpTrkFlag(2);
00665 }
00666 }
00667 }
00668
00669
00670 for(j=0;j<1+tmpm.GetLast();j++){
00671 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00672 jnflag=0;
00673 for(k=0;k<1+segm->GetTriEndLast();k++){
00674 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00675 if(segtmp->GetTmpTrkFlag()>1) jnflag=1;
00676 }
00677 if(jnflag>0){
00678 for(k=0;k<1+segm->GetTriEndLast();k++){
00679 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00680 if(segtmp->GetTmpTrkFlag()>0 && segtmp->GetTmpTrkFlag()<2) segtmp->SetTmpTrkFlag(2);
00681 }
00682 }
00683 }
00684
00685
00686 if(1+tmpm.GetLast()>0 && 1+tmpp.GetLast()>0){
00687
00688
00689 for(j=0;j<1+tmpm.GetLast();j++){
00690 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00691 for(k=0;k<1+segm->GetTriEndLast();k++){
00692 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00693 if(segm->GetTmpTrkFlag()>1 && segp->GetTmpTrkFlag()>1){
00694 jnflagm=0;
00695 for(l=0;l<1+segm->GetTriBegLast();l++){
00696 TrackSegmentAtNu* segm2 = (TrackSegmentAtNu*)(segm->GetTriBegAt(l));
00697 if( jnflagm<1 && (segm2->IsEndAssoc(segp)||segp->IsBegAssoc(segm2)) ) jnflagm=1;
00698 }
00699 jnflagp=0;
00700 for(l=0;l<1+segp->GetTriEndLast();l++){
00701 TrackSegmentAtNu* segp2 = (TrackSegmentAtNu*)(segp->GetTriEndAt(l));
00702 if( jnflagp<1 && (segm->IsEndAssoc(segp2)||segp2->IsBegAssoc(segm)) ) jnflagp=1;
00703 }
00704 if(jnflagm>0 && jnflagp>0){
00705 segm->AddAssocTriToEnd(segp); segp->AddAssocTriToBeg(segm);
00706 }
00707 }
00708 }
00709 }
00710
00711
00712 for(j=0;j<1+tmpm.GetLast();j++){
00713 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00714 if(segm->GetTmpTrkFlag()>1){
00715 jnflag=0;
00716 for(k=0;k<1+segm->GetTriEndLast();k++){
00717 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00718 if(segp->GetTmpTrkFlag()>1) jnflag=1;
00719 }
00720 if(jnflag<1){
00721 for(k=0;k<1+segm->GetTriEndLast();k++){
00722 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00723 jnflagm=0;
00724 for(l=0;l<1+segm->GetTriBegLast();l++){
00725 TrackSegmentAtNu* segm2 = (TrackSegmentAtNu*)(segm->GetTriBegAt(l));
00726 if( jnflagm<1 && (segm2->IsEndAssoc(segp)||segp->IsBegAssoc(segm2)) ) jnflagm=1;
00727 }
00728 if(jnflagm>0){
00729 segm->AddAssocTriToEnd(segp); segp->AddAssocTriToBeg(segm);
00730 }
00731 }
00732 }
00733 }
00734 }
00735
00736
00737 for(j=0;j<1+tmpp.GetLast();j++){
00738 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(tmpp.At(j));
00739 if(segp->GetTmpTrkFlag()>1){
00740 jnflag=0;
00741 for(k=0;k<1+segp->GetTriBegLast();k++){
00742 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(segp->GetTriBegAt(k));
00743 if(segm->GetTmpTrkFlag()>1) jnflag=1;
00744 }
00745 if(jnflag<1){
00746 for(k=0;k<1+segp->GetTriBegLast();k++){
00747 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(segp->GetTriBegAt(k));
00748 jnflagp=0;
00749 for(l=0;l<1+segp->GetTriEndLast();l++){
00750 TrackSegmentAtNu* segp2 = (TrackSegmentAtNu*)(segp->GetTriEndAt(l));
00751 if( jnflagp<1 && (segm->IsEndAssoc(segp2)||segp2->IsBegAssoc(segm)) ) jnflagp=1;
00752 }
00753 if(jnflagp>0){
00754 segm->AddAssocTriToEnd(segp); segp->AddAssocTriToBeg(segm);
00755 }
00756 }
00757 }
00758 }
00759 }
00760
00761
00762 for(j=0;j<1+tmpm.GetLast();j++){
00763 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00764 for(k=0;k<1+segm->GetTriEndLast();k++){
00765 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00766 if(segm->GetTmpTrkFlag()<2 && segp->GetTmpTrkFlag()<2){
00767 jnflagm=0;
00768 for(l=0;l<1+segm->GetTriEndLast();l++){
00769 TrackSegmentAtNu* segmp = (TrackSegmentAtNu*)(segm->GetTriEndAt(l));
00770 if(segmp->GetTmpTrkFlag()>1) jnflagm=1;
00771 }
00772 jnflagp=0;
00773 for(l=0;l<1+segp->GetTriBegLast();l++){
00774 TrackSegmentAtNu* segpm = (TrackSegmentAtNu*)(segp->GetTriBegAt(l));
00775 if(segpm->GetTmpTrkFlag()>1) jnflagp=1;
00776 }
00777 if(jnflagm<1 && jnflagp<1){
00778 segm->AddAssocTriToEnd(segp); segp->AddAssocTriToBeg(segm);
00779 }
00780 }
00781 }
00782 }
00783
00784 }
00785
00786
00787 for(j=0;j<1+tmpm.GetLast();j++){
00788 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00789 segm->SetTmpTrkFlag(0);
00790 }
00791
00792 for(j=0;j<1+tmpp.GetLast();j++){
00793 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(tmpp.At(j));
00794 segp->SetTmpTrkFlag(0);
00795 }
00796
00797 }
00798 }
00799
00800
00801
00802 TObjArray segbnk[2];
00803 for(mod=0;mod<modnum;mod++){
00804 for(pln=3;pln<modpln-1;pln++){
00805 i=pln+mod*(modpln+1);
00806
00807 for(j=0;j<1+segplnbnk[i].GetLast();j++){
00808 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segplnbnk[i].At(j));
00809
00810
00811 if(1+seg->GetTriBegLast()==0 && 1+seg->GetTriEndLast()==0){
00812 vuw=seg->GetPlaneView(); segbnk[vuw].Add(seg); seg->SetSegment(seg);
00813 }
00814
00815
00816 if(1+seg->GetAssocTriBegLast()>0 || 1+seg->GetAssocTriEndLast()>0){
00817 jnflag=0;
00818 if(1+seg->GetAssocTriBegLast()==1){
00819 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(seg->GetAssocTriBegAt(0));
00820 if(1+segtmp->GetAssocTriEndLast()==1){
00821 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(segtmp->GetSegment());
00822 seg0->AddSegment(seg); seg->SetSegment(seg0); seg->SetUID(-1);
00823 jnflag=1;
00824 }
00825 }
00826 if(jnflag<1){
00827 vuw=seg->GetPlaneView(); segbnk[vuw].Add(seg); seg->SetSegment(seg);
00828 for(k=0;k<1+seg->GetAssocTriBegLast();k++){
00829 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(seg->GetAssocTriBegAt(k));
00830 TrackSegmentAtNu* segtmp0 = (TrackSegmentAtNu*)(segtmp->GetSegment());
00831 seg->AddSegToBeg(segtmp0); segtmp0->AddSegToEnd(seg);
00832 }
00833 }
00834 }
00835
00836 }
00837 }
00838 }
00839
00840
00841
00842 Int_t npts,npts0=4;
00843 Double_t inv_npts;
00844 for(vuw=0;vuw<2;vuw++){
00845 tmp1.Clear(); tmp2.Clear();
00846 for(i=0;i<1+segbnk[vuw].GetLast();i++){
00847 TrackSegmentAtNu* seg1 = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
00848 if(1+seg1->GetSegEndLast()==0){
00849 jnflag=0; mod=(Int_t)((seg1->GetEndPlane())/(modpln+1)); npts=npts0;
00850 if( ( (modtype==1 && seg1->GetEndStrip()>80 && seg1->GetEndStrip()<=110) )
00851 || ( 1+(seg1->GetEndPlane()-seg1->GetBegPlane())/2>5) ){
00852 inv_npts = 0.2*seg1->GetEndDir();
00853 if(inv_npts<0) inv_npts=-inv_npts; if(inv_npts<0.01) inv_npts=0.01;
00854 npts=(Int_t)(1.0/inv_npts); if(npts<13) npts=13; if(npts>25) npts=25;
00855 }
00856 for(j=1;j<npts;j++){
00857 pln=seg1->GetEndPlane()+2*j;
00858 if(jnflag==0 && pln<(modpln+1)*(mod+1)){
00859 for(k=0;k<1+segplnbnk[pln].GetLast();k++){
00860 TrackSegmentAtNu* seg2 = (TrackSegmentAtNu*)(segplnbnk[pln].At(k));
00861 if( ( seg2->GetSegment()==seg2 && 1+seg2->GetSegBegLast()==0 )
00862 && ( seg2->GetBegPlane()-seg1->GetEndPlane()>0
00863 || seg1->GetClrLast()>5 && seg2->GetClrLast()>5 )
00864 && ( j<npts0
00865 || ( modtype==1
00866 && ( ( seg2->GetBegStrip()>80 && seg2->GetBegStrip()<=100
00867 && seg1->GetEndStrip()>90 && seg1->GetEndStrip()<=110 )
00868 || ( seg2->GetBegStrip()>90 && seg2->GetBegStrip()<=110
00869 && seg1->GetEndStrip()>80 && seg1->GetEndStrip()<=100 ) ) )
00870 || ( 1+(seg1->GetEndPlane()-seg1->GetBegPlane())/2>5 ) ) ){
00871 if( seg1->IsEndAssoc(seg2) ){
00872 tmp1.Add(seg1); tmp2.Add(seg2); jnflag=1;
00873 MSG("AlgAtNuReco", Msg::kVerbose) << " NON-ADJ ASSOC: "
00874 << seg1->GetBegPlane() << " " << seg1->GetEndPlane() << " "
00875 << seg2->GetBegPlane() << " " << seg2->GetEndPlane() << endl;
00876 }
00877 }
00878 }
00879 }
00880 }
00881 }
00882 }
00883 for(i=0;i<1+tmp1.GetLast();i++){
00884 TrackSegmentAtNu* seg1 = (TrackSegmentAtNu*)(tmp1.At(i));
00885 TrackSegmentAtNu* seg2 = (TrackSegmentAtNu*)(tmp2.At(i));
00886 seg1->AddSegToEnd(seg2); seg2->AddSegToBeg(seg1);
00887 }
00888 }
00889
00890
00891 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF TRACK SEGMENTS *** " << endl;
00892 for(vuw=0;vuw<2;vuw++){
00893 MSG("AlgAtNuReco", Msg::kVerbose) << "view: " << vuw << endl;
00894 for(i=0;i<1+segbnk[vuw].GetLast();i++){
00895 TrackSegmentAtNu* segment = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
00896 MSG("AlgAtNuReco", Msg::kVerbose) << " i=" << i
00897 << " begpln: " << segment->GetBegPlane()
00898 << " begstr: " << segment->GetBegStrip()
00899 << " endpln: " << segment->GetEndPlane()
00900 << " endstr: " << segment->GetEndStrip() << endl;
00901 MSG("AlgAtNuReco", Msg::kVerbose) << " beg: " << endl;
00902 for(j=0;j<1+segment->GetSegBegLast();j++){
00903 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segment->GetSegBegAt(j));
00904 MSG("AlgAtNuReco", Msg::kVerbose) << " j=" << j
00905 << " begpln: " << segtmp->GetBegPlane()
00906 << " begstr: " << segtmp->GetBegStrip()
00907 << " endpln: " << segtmp->GetEndPlane()
00908 << " endstr: " << segtmp->GetEndStrip() << endl;
00909 }
00910 MSG("AlgAtNuReco", Msg::kVerbose) << " end: " << endl;
00911 for(j=0;j<1+segment->GetSegEndLast();j++){
00912 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segment->GetSegEndAt(j));
00913 MSG("AlgAtNuReco", Msg::kVerbose) << " j=" << j
00914 << " begpln: " << segtmp->GetBegPlane()
00915 << " begstr: " << segtmp->GetBegStrip()
00916 << " endpln: " << segtmp->GetEndPlane()
00917 << " endstr: " << segtmp->GetEndStrip() << endl;
00918 }
00919 }
00920 }
00921
00922 gettimeofday(&tpafter,0);
00923 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
00924 fTime_TrackAtNu+=dtime;
00925
00926
00927
00928
00929
00930
00931 gettimeofday(&tpbefore,0);
00932
00933 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of shower segments *** " << endl;
00934
00935
00936 Int_t pln0;
00937 TObjArray segbnk2[2];
00938 for(mod=0;mod<modnum;mod++){
00939 for(pln=1;pln<1+modpln;pln++){
00940 i=pln+mod*(modpln+1);
00941 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00942 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00943 if(clr0->GetShwFlag()>0 && clr0->GetShwFlag()<2){
00944 ShowerSegmentAtNu* seg0 = new ShowerSegmentAtNu(clr0); NEWOBJCTR++;
00945 vuw = seg0->GetPlaneView(); segbnk2[vuw].Add(seg0);
00946 clr0->SetShwFlag(2); shwsegbnk.Add(seg0); ctr=1;
00947 while(ctr>0){
00948 ctr=0;
00949 pln0 = -2 + seg0->GetBegPlane();
00950 npln = 3 + (seg0->GetEndPlane()-seg0->GetBegPlane())/2;
00951 for(k=0;k<npln;k++){
00952 k0=pln0+2*k;
00953 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
00954 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
00955 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
00956 if(clr1->GetShwFlag()>0 && clr1->GetShwFlag()<2){
00957 assoc=seg0->IsShwAssoc(clr1);
00958 if(assoc>0){
00959 seg0->AddCluster(clr1);
00960 clr1->SetShwFlag(2); ctr++;
00961 }
00962 }
00963 }
00964 }
00965 }
00966 }
00967 }
00968 }
00969 }
00970 }
00971
00972
00973 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF SHOWER SEGMENTS *** " << endl;
00974 for(vuw=0;vuw<2;vuw++){
00975 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
00976 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
00977 MSG("AlgAtNuReco", Msg::kVerbose)
00978 << "vuw=" << vuw << " i=" << i << " "
00979 << seg->GetBegPlane() << " " << seg->GetEndPlane() << " "
00980 << seg->GetBegStrip() << " " << seg->GetEndStrip() << " "
00981 << 1+seg->GetHitLast() << endl;
00982 }
00983 }
00984
00985 gettimeofday(&tpafter,0);
00986 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
00987 fTime_ShowerAtNu+=dtime;
00988
00989
00990
00991
00992
00993
00994 Int_t shwctr,trkctr,plnctr;
00995 TObjArray mysegbnk[2],mysegbnk2[2];
00996 MSG("AlgAtNuReco", Msg::kDebug) << " *** comparing views (1) *** " << endl;
00997
00998
00999 gettimeofday(&tpbefore,0);
01000
01001 for(vuw=0;vuw<2;vuw++){
01002 for(i=0;i<1+segbnk[vuw].GetLast();i++){
01003 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
01004 shwctr=0; trkctr=0; plnctr=0;
01005 for(j=0;j<1+seg->GetClrLast();j++){
01006 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01007 if(clr->GetTrkFlag()>1){ trkctr++; if(clr->GetShwFlag()<2) shwctr++; }
01008 if(clr->GetTrkPlnFlag()>0){ plnctr++; }
01009 }
01010 if( plnctr>1 ){
01011 if(shwctr>0) seg->SetUID(1); if(shwctr>2) seg->SetUID(2);
01012 }
01013 }
01014 }
01015
01016 for(i=0;i<1+segbnk[0].GetLast();i++){
01017 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(segbnk[0].At(i));
01018 for(j=0;j<1+segbnk[1].GetLast();j++){
01019 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(segbnk[1].At(j));
01020 if( segu->GetBegPlane()+1<segv->GetEndPlane() && segu->GetEndPlane()>segv->GetBegPlane()+1
01021 && segu->GetUID()>0 && segv->GetUID()>0 ){
01022 if(segv->GetUID()>1) mysegbnk[0].Add(segu); if(segu->GetUID()>1) mysegbnk[1].Add(segv);
01023 }
01024 }
01025 }
01026
01027 for(vuw=0;vuw<2;vuw++){
01028 for(i=0;i<1+mysegbnk[vuw].GetLast();i++){
01029 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(mysegbnk[vuw].At(i));
01030 seg->SetUID(2);
01031 }
01032 }
01033
01034 for(vuw=0;vuw<2;vuw++){
01035 for(i=0;i<1+segbnk[vuw].GetLast();i++){
01036 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
01037 if( seg->GetUID()>0 && seg->GetUID()<2 ) seg->SetReseedFlag(1);
01038 }
01039 }
01040
01041 gettimeofday(&tpafter,0);
01042 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01043 fTime_TrackAtNu+=dtime;
01044
01045
01046 gettimeofday(&tpbefore,0);
01047
01048 for(vuw=0;vuw<2;vuw++){
01049 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01050 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01051 trkctr=0; shwctr=0; plnctr=0;
01052 for(j=0;j<1+seg->GetClrLast();j++){
01053 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01054 if(clr->GetShwFlag()>1){ shwctr++; if(clr->GetTrkFlag()<2) trkctr++; }
01055 if(clr->GetShwPlnFlag()>0){ plnctr++; }
01056 }
01057 if( plnctr>-1 ){
01058 if(trkctr>-1) seg->SetUID(1); if(trkctr>0) seg->SetUID(2);
01059 }
01060 }
01061 }
01062
01063 for(i=0;i<1+segbnk2[0].GetLast();i++){
01064 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(segbnk2[0].At(i));
01065 for(j=0;j<1+segbnk2[1].GetLast();j++){
01066 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(segbnk2[1].At(j));
01067 if( segu->GetBegPlane()-2<segv->GetEndPlane() && segu->GetEndPlane()>segv->GetBegPlane()-2
01068 && segu->GetUID()>0 && segv->GetUID()>0 ){
01069 if(segv->GetUID()>1) mysegbnk2[0].Add(segu); if(segu->GetUID()>1) mysegbnk2[1].Add(segv);
01070 }
01071 }
01072 }
01073
01074 for(vuw=0;vuw<2;vuw++){
01075 for(i=0;i<1+mysegbnk2[vuw].GetLast();i++){
01076 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(mysegbnk2[vuw].At(i));
01077 seg->SetUID(2);
01078 }
01079 }
01080
01081 for(vuw=0;vuw<2;vuw++){
01082 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01083 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01084 if( seg->GetUID()>0 && seg->GetUID()<2 ) seg->SetReseedFlag(1);
01085 }
01086 }
01087
01088 gettimeofday(&tpafter,0);
01089 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01090 fTime_ShowerAtNu+=dtime;
01091
01092
01093
01094
01095
01096
01097 gettimeofday(&tpbefore,0);
01098
01099 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of 2D tracks *** " << endl;
01100 Int_t bpln,epln;
01101 Int_t objctr,objarrctr;
01102 Double_t dscr,topscr;
01103 Double_t strw,strdw;
01104 Int_t idm,idp;
01105 Int_t ntrks;
01106 Int_t tmpflag,trkflag;
01107 TObjArray tmpb,tmpe;
01108 TObjArray begbnk,endbnk;
01109 TObjArray tmpbegbnk,tmpendbnk;
01110 TObjArray tmpbnk,tmpbnk2;
01111
01112 if(1+segbnk[0].GetLast()>0 && 1+segbnk[1].GetLast()>0){
01113 for(vuw=0;vuw<2;vuw++){
01114
01115 cont=1; ntrks=0;
01116 while(cont){
01117 cont=0;
01118 begbnk.Clear(); endbnk.Clear();
01119
01120
01121 tmp1.Clear(); tmp2.Clear();
01122 for(i=0;i<1+segbnk[vuw].GetLast();i++){
01123 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
01124 if( seg->GetUID()>0 && seg->GetUID()<3 ){
01125 ctr=0;
01126 for(j=0;j<1+seg->GetClrLast();j++){
01127 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01128 if( clr->GetTrkPlnFlag()>0 && clr->GetTrkFlag()>1 && clr->GetTrkFlag()<3 ) ctr++;
01129 }
01130 if( ntrks>0 && ctr<3 ) seg->SetUID(0);
01131 }
01132 if(seg->GetUID()==2) tmp1.Add(seg); if(seg->GetUID()==1) tmp2.Add(seg);
01133 }
01134
01135 if(1+tmp1.GetLast()==0){
01136 for(i=0;i<1+tmp2.GetLast();i++){
01137 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp2.At(i));
01138 tmp1.Add(seg);
01139 }
01140 }
01141
01142 tmp2.Clear();
01143 if(1+tmp1.GetLast()>0){
01144 npln=-1;
01145 for(i=0;i<1+tmp1.GetLast();i++){
01146 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp1.At(i));
01147 if( 1+seg->GetEndPlane()-seg->GetBegPlane()>npln ){
01148 npln=1+seg->GetEndPlane()-seg->GetBegPlane();
01149 }
01150 }
01151 if(npln>0){
01152 for(i=0;i<1+tmp1.GetLast();i++){
01153 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp1.At(i));
01154 if( 1+seg->GetEndPlane()-seg->GetBegPlane()>12
01155 || 1+seg->GetEndPlane()-seg->GetBegPlane()>npln-5 ){
01156 tmp2.Add(seg);
01157 }
01158 }
01159 }
01160 }
01161
01162
01163 tmp1.Clear();
01164 if(1+tmp2.GetLast()>0){
01165 npln=-1;
01166 for(i=0;i<1+tmp2.GetLast();i++){
01167 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp2.At(i));
01168 tmpm.Clear(); tmpp.Clear();
01169
01170 tmpb.Clear(); tmpb.Add(seg); tmpbnk.Clear();
01171 while(1+tmpb.GetLast()>0){
01172 tmp.Clear();
01173 for(j=0;j<1+tmpb.GetLast();j++){
01174 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpb.At(j));
01175 tmp.Add(segtmp);
01176 }
01177 tmpb.Clear();
01178 for(j=0;j<1+tmp.GetLast();j++){
01179 tmpflag=0;
01180 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmp.At(j));
01181 for(k=0;k<1+segtmp->GetSegBegLast();k++){
01182 TrackSegmentAtNu* segb = (TrackSegmentAtNu*)(segtmp->GetSegBegAt(k));
01183 if(segb->GetUID()>-1 && segb->GetUID()<3){
01184 tmpflag=1;
01185 if(segb->GetTmpTrkFlag()<1){
01186 tmpb.Add(segb); segb->SetTmpTrkFlag(1); tmpbnk.Add(segb);
01187 }
01188 }
01189 }
01190 if(tmpflag==0){
01191 tmpm.Add(segtmp);
01192 }
01193 }
01194 }
01195 for(j=0;j<1+tmpbnk.GetLast();j++){
01196 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpbnk.At(j));
01197 segtmp->SetTmpTrkFlag(0);
01198 }
01199
01200 tmpe.Clear(); tmpe.Add(seg); tmpbnk.Clear();
01201 while(1+tmpe.GetLast()>0){
01202 tmp.Clear();
01203 for(j=0;j<1+tmpe.GetLast();j++){
01204 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpe.At(j));
01205 tmp.Add(segtmp);
01206 }
01207 tmpe.Clear();
01208 for(j=0;j<1+tmp.GetLast();j++){
01209 tmpflag=0;
01210 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmp.At(j));
01211 for(k=0;k<1+segtmp->GetSegEndLast();k++){
01212 TrackSegmentAtNu* sege = (TrackSegmentAtNu*)(segtmp->GetSegEndAt(k));
01213 if(sege->GetUID()>-1 && sege->GetUID()<3){
01214 tmpflag=1;
01215 if(sege->GetTmpTrkFlag()<1){
01216 tmpe.Add(sege); sege->SetTmpTrkFlag(1); tmpbnk.Add(sege);
01217 }
01218 }
01219 }
01220 if(tmpflag==0){
01221 tmpp.Add(segtmp);
01222 }
01223 }
01224 }
01225 for(j=0;j<1+tmpbnk.GetLast();j++){
01226 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpbnk.At(j));
01227 segtmp->SetTmpTrkFlag(0);
01228 }
01229
01230 npts=-1;
01231 for(j=0;j<1+tmpm.GetLast();j++){
01232 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
01233 for(k=0;k<1+tmpp.GetLast();k++){
01234 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(tmpp.At(k));
01235 if(1+segp->GetEndPlane()-segm->GetBegPlane()>npts){
01236 npts=1+segp->GetEndPlane()-segm->GetBegPlane();
01237 }
01238 }
01239 }
01240 seg->SetNPlanes(npts);
01241 if(npts>npln) npln=npts;
01242 }
01243 if(npln>0){
01244 for(i=0;i<1+tmp2.GetLast();i++){
01245 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp2.At(i));
01246 if( seg->GetNPlanes()>npln-5 ){
01247 tmp1.Add(seg);
01248 }
01249 }
01250 }
01251 }
01252
01253
01254 tmpbnk2.Clear();
01255 for(i=0;i<1+tmp1.GetLast();i++){
01256 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(tmp1.At(i));
01257
01258 MSG("AlgAtNuReco", Msg::kVerbose)
01259 << " PROPAGATE SEED SEGMENT " << seg0->GetBegPlane() << " " << seg0->GetEndPlane() << endl;
01260
01261 seg0->SetTmpTrkFlag(3); trkflag=0;
01262 tmpbnk.Clear(); tmpbnk.Add(seg0);
01263
01264
01265 tmpm.Clear(); tmpb.Clear(); tmpb.Add(seg0);
01266 while(1+tmpb.GetLast()>0){
01267 tmp.Clear();
01268 for(j=0;j<1+tmpb.GetLast();j++){
01269 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpb.At(j));
01270 tmp.Add(segtmp);
01271 }
01272 tmpb.Clear();
01273 for(j=0;j<1+tmp.GetLast();j++){
01274 tmpflag=0;
01275 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmp.At(j));
01276 for(k=0;k<1+segtmp->GetSegBegLast();k++){
01277 TrackSegmentAtNu* segb = (TrackSegmentAtNu*)(segtmp->GetSegBegAt(k));
01278 if(segb->GetUID()>-1 && segb->GetUID()<3){
01279 tmpflag=1;
01280 if(segb->GetTmpTrkFlag()<1){
01281 tmpb.Add(segb); segb->SetTmpTrkFlag(1); tmpbnk.Add(segb);
01282 }
01283 }
01284 }
01285 if(tmpflag==0){
01286 tmpm.Add(segtmp);
01287 }
01288 }
01289 }
01290
01291 npln=999;
01292 for(j=0;j<1+tmpm.GetLast();j++){
01293 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpm.At(j));
01294 if(seg->GetBegPlane()<npln) npln=seg->GetBegPlane();
01295 }
01296
01297 tmpb.Clear();
01298 for(j=0;j<1+tmpm.GetLast();j++){
01299 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpm.At(j));
01300 if(seg->GetBegPlane()<npln+5){
01301 tmpb.Add(seg);
01302 }
01303 }
01304
01305 while(1+tmpb.GetLast()>0){
01306 tmp.Clear();
01307 for(j=0;j<1+tmpb.GetLast();j++){
01308 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpb.At(j));
01309 if(seg->GetTmpTrkFlag()<2){
01310 seg->SetTmpTrkFlag(2); tmp.Add(seg);
01311 }
01312 if(seg->GetTrkFlag()<1){
01313 seg->SetTrkFlag(1); tmpbnk2.Add(seg); trkflag=1;
01314 }
01315 }
01316 tmpb.Clear();
01317 for(j=0;j<1+tmp.GetLast();j++){
01318 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp.At(j));
01319 for(k=0;k<1+seg->GetSegEndLast();k++){
01320 TrackSegmentAtNu* segb = (TrackSegmentAtNu*)(seg->GetSegEndAt(k));
01321 if(segb->GetTmpTrkFlag()>0){
01322 tmpb.Add(segb);
01323 }
01324 }
01325 }
01326 }
01327
01328
01329 tmpp.Clear(); tmpe.Clear(); tmpe.Add(seg0);
01330 while(1+tmpe.GetLast()>0){
01331 tmp.Clear();
01332 for(j=0;j<1+tmpe.GetLast();j++){
01333 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpe.At(j));
01334 tmp.Add(segtmp);
01335 }
01336 tmpe.Clear();
01337 for(j=0;j<1+tmp.GetLast();j++){
01338 tmpflag=0;
01339 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmp.At(j));
01340 for(k=0;k<1+segtmp->GetSegEndLast();k++){
01341 TrackSegmentAtNu* sege = (TrackSegmentAtNu*)(segtmp->GetSegEndAt(k));
01342 if(sege->GetUID()>-1 && sege->GetUID()<3){
01343 tmpflag=1;
01344 if(sege->GetTmpTrkFlag()<1){
01345 tmpe.Add(sege); sege->SetTmpTrkFlag(1); tmpbnk.Add(sege);
01346 }
01347 }
01348 }
01349 if(tmpflag==0){
01350 tmpp.Add(segtmp);
01351 }
01352 }
01353 }
01354
01355 npln=-999;
01356 for(j=0;j<1+tmpp.GetLast();j++){
01357 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpp.At(j));
01358 if(seg->GetEndPlane()>npln) npln=seg->GetEndPlane();
01359 }
01360
01361 tmpe.Clear();
01362 for(j=0;j<1+tmpp.GetLast();j++){
01363 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpp.At(j));
01364 if(seg->GetEndPlane()>npln-5){
01365 tmpe.Add(seg);
01366 }
01367 }
01368
01369 while(1+tmpe.GetLast()>0){
01370 tmp.Clear();
01371 for(j=0;j<1+tmpe.GetLast();j++){
01372 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpe.At(j));
01373 if(seg->GetTmpTrkFlag()<2){
01374 seg->SetTmpTrkFlag(2); tmp.Add(seg);
01375 }
01376 if(seg->GetTrkFlag()<1){
01377 seg->SetTrkFlag(1); tmpbnk2.Add(seg); trkflag=1;
01378 }
01379 }
01380 tmpe.Clear();
01381 for(j=0;j<1+tmp.GetLast();j++){
01382 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp.At(j));
01383 for(k=0;k<1+seg->GetSegBegLast();k++){
01384 TrackSegmentAtNu* sege = (TrackSegmentAtNu*)(seg->GetSegBegAt(k));
01385 if(sege->GetTmpTrkFlag()>0){
01386 tmpe.Add(sege);
01387 }
01388 }
01389 }
01390 }
01391
01392 if(trkflag){
01393
01394
01395 tmpbegbnk.Clear(); objarrctr=0;
01396 TObjArray* tmpbegsegbnk = new TObjArray(); NEWOBJCTR++;
01397 tmpbegsegbnk->Add(seg0); tmpbegbnk.Add(tmpbegsegbnk); objctr=1;
01398 while(objctr>0){
01399 objctr=0; npts=1+tmpbegbnk.GetLast();
01400 for(j=0;j<npts;j++){
01401 TObjArray* objtmp = (TObjArray*)(tmpbegbnk.At(j));
01402 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(objtmp->Last());
01403 tmpm.Clear();
01404 for(k=0;k<1+segtmp->GetSegBegLast();k++){
01405 TrackSegmentAtNu* segbeg = (TrackSegmentAtNu*)(segtmp->GetSegBegAt(k));
01406 if(segbeg->GetTmpTrkFlag()>1){
01407 tmpm.Add(segbeg);
01408 }
01409 }
01410 if(1+tmpm.GetLast()>0){
01411 for(k=1;k<1+tmpm.GetLast();k++){
01412 if(objarrctr<fObjMax){
01413 TObjArray* objnew = new TObjArray(*objtmp); NEWOBJCTR++;
01414 TrackSegmentAtNu* segbeg = (TrackSegmentAtNu*)(tmpm.At(k));
01415 objnew->Add(segbeg); tmpbegbnk.Add(objnew);
01416 objarrctr++;
01417 }
01418 else{
01419 if(!fObjMaxFlag){
01420 fObjMaxFlag=1;
01421 MSG("AlgAtNuReco", Msg::kWarning) << " *** CRAZY MEMORY GUZZLING *** " << endl;
01422 }
01423 }
01424 }
01425 TrackSegmentAtNu* segbeg = (TrackSegmentAtNu*)(tmpm.First());
01426 objtmp->Add(segbeg); objctr++;
01427 }
01428 }
01429 }
01430
01431 id=-1; topscr=-1.0;
01432 for(j=0;j<1+tmpbegbnk.GetLast();j++){
01433 TObjArray* objtmp = (TObjArray*)(tmpbegbnk.At(j));
01434 scr = seg0->GetScore(objtmp,0);
01435 if(scr>topscr){ topscr=scr; id=j; }
01436 }
01437
01438 TObjArray* objbegtmp = (TObjArray*)(tmpbegbnk.At(id));
01439 TObjArray* objbegnew = new TObjArray(*objbegtmp); begbnk.Add(objbegnew); NEWOBJCTR++;
01440 DELOBJCTR+=1+tmpbegbnk.GetLast(); tmpbegbnk.Delete();
01441
01442
01443 tmpendbnk.Clear(); objarrctr=0;
01444 TObjArray* tmpendsegbnk = new TObjArray(); NEWOBJCTR++;
01445 tmpendsegbnk->Add(seg0); tmpendbnk.Add(tmpendsegbnk); objctr=1;
01446 while(objctr>0){
01447 objctr=0; npts=1+tmpendbnk.GetLast();
01448 for(j=0;j<npts;j++){
01449 TObjArray* objtmp = (TObjArray*)(tmpendbnk.At(j));
01450 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(objtmp->Last());
01451 tmpp.Clear();
01452 for(k=0;k<1+segtmp->GetSegEndLast();k++){
01453 TrackSegmentAtNu* segend = (TrackSegmentAtNu*)(segtmp->GetSegEndAt(k));
01454 if(segend->GetTmpTrkFlag()>1){
01455 tmpp.Add(segend);
01456 }
01457 }
01458 if(1+tmpp.GetLast()>0){
01459 for(k=1;k<1+tmpp.GetLast();k++){
01460 if(objarrctr<fObjMax){
01461 TObjArray* objnew = new TObjArray(*objtmp); NEWOBJCTR++;
01462 TrackSegmentAtNu* segend = (TrackSegmentAtNu*)(tmpp.At(k));
01463 objnew->Add(segend); tmpendbnk.Add(objnew);
01464 objarrctr++;
01465 }
01466 else{
01467 if(!fObjMaxFlag){
01468 fObjMaxFlag=1;
01469 MSG("AlgAtNuReco", Msg::kWarning) << " *** CRAZY MEMORY GUZZLING *** " << endl;
01470 }
01471 }
01472 }
01473 TrackSegmentAtNu* segend = (TrackSegmentAtNu*)(tmpp.First());
01474 objtmp->Add(segend); objctr++;
01475 }
01476 }
01477 }
01478
01479 id=-1; topscr=-1.0;
01480 for(j=0;j<1+tmpendbnk.GetLast();j++){
01481 TObjArray* objtmp = (TObjArray*)(tmpendbnk.At(j));
01482 scr = seg0->GetScore(0,objtmp);
01483 if(scr>topscr){ topscr=scr; id=j; }
01484 }
01485
01486 TObjArray* objendtmp = (TObjArray*)(tmpendbnk.At(id));
01487 TObjArray* objendnew = new TObjArray(*objendtmp); endbnk.Add(objendnew); NEWOBJCTR++;
01488 DELOBJCTR+=1+tmpendbnk.GetLast(); tmpendbnk.Delete();
01489 }
01490
01491 for(j=0;j<1+tmpbnk.GetLast();j++){
01492 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpbnk.At(j));
01493 seg->SetTmpTrkFlag(0);
01494 }
01495 }
01496
01497 for(j=0;j<1+tmpbnk2.GetLast();j++){
01498 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpbnk2.At(j));
01499 seg->SetTrkFlag(0);
01500 }
01501
01502
01503 if(1+begbnk.GetLast()>0 && 1+endbnk.GetLast()>0){
01504
01505 id=-1; topscr=-1.0;
01506 for(i=0;i<1+begbnk.GetLast();i++){
01507 TObjArray* objbeg = (TObjArray*)(begbnk.At(i));
01508 TObjArray* objend = (TObjArray*)(endbnk.At(i));
01509 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(objbeg->First());
01510 scr = seg0->GetScore(objbeg,objend);
01511 if(scr>topscr){ topscr=scr; id=i; }
01512 }
01513
01514 if(1+id>0){
01515 TObjArray* objbeg = (TObjArray*)(begbnk.At(id));
01516 TObjArray* objend = (TObjArray*)(endbnk.At(id));
01517 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(objbeg->First());
01518 for(i=1;i<1+objbeg->GetLast();i++){
01519 TrackSegmentAtNu* segbeg = (TrackSegmentAtNu*)(objbeg->At(i));
01520 seg0->AddSegment(segbeg); segbeg->SetUID(-1);
01521 }
01522 for(i=1;i<1+objend->GetLast();i++){
01523 TrackSegmentAtNu* segend = (TrackSegmentAtNu*)(objend->At(i));
01524 seg0->AddSegment(segend); segend->SetUID(-1);
01525 }
01526 seg0->SetUID(3);
01527 for(i=0;i<1+seg0->GetClrLast();i++){
01528 ClusterAtNu* clr = (ClusterAtNu*)(seg0->GetClrAt(i));
01529 clr->SetTrkFlag(3);
01530 }
01531 cont=1; ntrks++;
01532 }
01533 }
01534
01535 MSG("AlgAtNuReco", Msg::kVerbose) << " ... about to delete beg/end... " << endl;
01536 DELOBJCTR+=1+begbnk.GetLast(); DELOBJCTR+=1+endbnk.GetLast();
01537 begbnk.Delete(); endbnk.Delete();
01538 MSG("AlgAtNuReco", Msg::kVerbose) << " beg/end deleted ... " << endl;
01539 }
01540 }
01541 }
01542
01543
01544 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF 2D TRACKS *** " << endl;
01545 for(vuw=0;vuw<2;vuw++){
01546 MSG("AlgAtNuReco", Msg::kVerbose) << "view: " << vuw << endl;
01547 for(i=0;i<1+segbnk[vuw].GetLast();i++){
01548 TrackSegmentAtNu* segment = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
01549 if(segment->GetUID()>2){
01550 MSG("AlgAtNuReco", Msg::kVerbose) << " i=" << i
01551 << " begpln: " << segment->GetBegPlane()
01552 << " begstr: " << segment->GetBegStrip()
01553 << " endpln: " << segment->GetEndPlane()
01554 << " endstr: " << segment->GetEndStrip() << endl;
01555 }
01556 }
01557 }
01558
01559 gettimeofday(&tpafter,0);
01560 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01561 fTime_TrackAtNu+=dtime;
01562
01563
01564
01565
01566
01567
01568 gettimeofday(&tpbefore,0);
01569
01570 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of 2D showers *** " << endl;
01571 Int_t nshws,maxshws;
01572 if(1+segbnk2[0].GetLast()>0 && 1+segbnk2[1].GetLast()>0){
01573 for(vuw=0;vuw<2;vuw++){
01574
01575 cont=1; nshws=0; maxshws=4;
01576 while(cont){
01577 cont=0;
01578
01579
01580 tmp1.Clear(); tmp2.Clear();
01581 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01582 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01583 if(seg->GetUID()>0 && seg->GetUID()<3){
01584 ctr=0;
01585 for(j=0;j<1+seg->GetClrLast();j++){
01586 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01587 if( clr->GetShwPlnFlag()>0 && clr->GetShwFlag()>1 && clr->GetShwFlag()<3 ) ctr++;
01588 }
01589 if( nshws>0 && ctr<1 ) seg->SetUID(0);
01590 }
01591 if(seg->GetUID()==2) tmp1.Add(seg); if(seg->GetUID()==1) tmp2.Add(seg);
01592 }
01593
01594 if(1+tmp1.GetLast()==0){
01595 for(i=0;i<1+tmp2.GetLast();i++){
01596 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmp2.At(i));
01597 tmp1.Add(seg);
01598 }
01599 }
01600
01601
01602 npln=0; id=-1; topscr=0.0;
01603 for(i=0;i<1+tmp1.GetLast();i++){
01604 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmp1.At(i));
01605 npts = 1+(seg->GetEndPlane()-seg->GetBegPlane())/2;
01606 if(npts>npln-3){
01607 scr=seg->GetScore();
01608 if(npts>npln || scr>topscr){
01609 npln=npts; topscr=scr;
01610 id=i;
01611 }
01612 }
01613 }
01614
01615
01616 if(1+id>0){
01617 ShowerSegmentAtNu* seg0 = (ShowerSegmentAtNu*)(tmp1.At(id));
01618
01619 mod=(Int_t)((seg0->GetEndPlane())/(modpln+1));
01620 seg0->SetUID(3);
01621 cont=1; nshws++;
01622
01623 ctr=1;
01624 while(ctr>0){
01625 ctr=0;
01626
01627
01628 pln0 = -4 + seg0->GetBegPlane();
01629 npln = 5 + (seg0->GetEndPlane()-seg0->GetBegPlane())/2;
01630 for(k=0;k<npln;k++){
01631 k0 = pln0+2*k;
01632 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
01633 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
01634 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
01635 if(clr1->GetShwFlag()<2){
01636 assoc=seg0->IsShwAssoc(clr1);
01637 if( ( assoc>0 && clr1->GetTrkFlag()<3 )
01638 || ( assoc>1 && clr1->GetPlane()-seg0->GetBegPlane()>0 && seg0->GetEndPlane()-clr1->GetPlane()>0 ) ){
01639 seg0->AddCluster(clr1);
01640 clr1->SetShwFlag(2); ctr++;
01641 }
01642 }
01643 }
01644 }
01645 }
01646
01647
01648 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01649 ShowerSegmentAtNu* segtmp = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01650 if( segtmp->GetUID()>0 && segtmp->GetUID()<3
01651 && (Int_t)((segtmp->GetEndPlane())/(modpln+1))==mod ){
01652 assoc=seg0->IsShwAssoc(segtmp);
01653 if(assoc>0){
01654 seg0->AddSegment(segtmp);
01655 segtmp->SetUID(-1); ctr++;
01656 }
01657 }
01658 }
01659
01660 }
01661
01662 for(i=0;i<1+seg0->GetClrLast();i++){
01663 ClusterAtNu* clr = (ClusterAtNu*)(seg0->GetClrAt(i));
01664 clr->SetShwFlag(3);
01665 }
01666
01667 }
01668
01669 }
01670 }
01671 }
01672
01673
01674 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF 2D SHOWERS *** " << endl;
01675 for(vuw=0;vuw<2;vuw++){
01676 MSG("AlgAtNuReco", Msg::kVerbose) << "view: " << vuw << endl;
01677 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01678 ShowerSegmentAtNu* segment = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01679 if(segment->GetUID()>2){
01680 MSG("AlgAtNuReco", Msg::kVerbose) << " i=" << i
01681 << " begpln: " << segment->GetBegPlane()
01682 << " begstr: " << segment->GetBegStrip()
01683 << " endpln: " << segment->GetEndPlane()
01684 << " endstr: " << segment->GetEndStrip() << endl;
01685 }
01686 }
01687 }
01688
01689 gettimeofday(&tpafter,0);
01690 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01691 fTime_ShowerAtNu+=dtime;
01692
01693
01694
01695
01696
01697
01698 MSG("AlgAtNuReco", Msg::kDebug) << " *** matching views (2) *** " << endl;
01699
01700 gettimeofday(&tpbefore,0);
01701
01702
01703 TObjArray tmptrk[2],tmptrkseg[2];
01704 Int_t dstr,dpln,dpln1,dpln2;
01705 Int_t bpln2,bpln1,epln1,epln2;
01706 Int_t jnflag2[2];
01707 for(i=0;i<1+segbnk[0].GetLast();i++){
01708 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(segbnk[0].At(i));
01709 if(segu->GetUID()>2){
01710 for(j=0;j<1+segbnk[1].GetLast();j++){
01711 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(segbnk[1].At(j));
01712 if(segv->GetUID()>2){
01713 if( ( (segv->GetBegPlane()-segu->GetBegPlane()>-10
01714 && segv->GetBegPlane()-segu->GetBegPlane()<10)
01715 || (segv->GetEndPlane()-segu->GetEndPlane()>-10
01716 && segv->GetEndPlane()-segu->GetEndPlane()<10) )
01717 && ( segv->GetEndPlane()-segu->GetBegPlane()>2
01718 && segu->GetEndPlane()-segv->GetBegPlane()>2 )
01719 && ( segv->GetEndTime()-segu->GetBegTime()>-fTimeWin
01720 && segu->GetEndTime()-segv->GetBegTime()>-fTimeWin ) ){
01721 tmptrkseg[0].Add(segu); tmptrkseg[1].Add(segv);
01722 }
01723 }
01724 }
01725 }
01726 }
01727
01728 if(modnum>1){
01729 TObjArray* tmpmod = new TObjArray[2*(modnum-1)]; NEWOBJCTR+=2*(modnum-1);
01730 for(i=0;i<2*(modnum-1);i++){
01731 for(j=0;j<2;j++){
01732 TObjArray* tmparr = new TObjArray(); NEWOBJCTR++;
01733 tmpmod[i].Add(tmparr);
01734 }
01735 }
01736
01737 for(i=0;i<1+tmptrkseg[0].GetLast();i++){
01738 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(tmptrkseg[0].At(i));
01739 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(tmptrkseg[1].At(i));
01740 mod=(Int_t)((segu->GetEndPlane())/(modpln+1));
01741 if( segu->GetBegPlane()<mod*(modpln+1)+20
01742 && segv->GetBegPlane()<mod*(modpln+1)+20
01743 && mod>0 ){
01744 TObjArray* tmparrum = (TObjArray*)(tmpmod[2*mod-1].At(0));
01745 TObjArray* tmparrvm = (TObjArray*)(tmpmod[2*mod-1].At(1));
01746 tmparrum->Add(segu); tmparrvm->Add(segv);
01747 }
01748 if( segu->GetEndPlane()>(mod+1)*(modpln+1)-20
01749 && segv->GetEndPlane()>(mod+1)*(modpln+1)-20
01750 && mod<modnum-1 ){
01751 TObjArray* tmparru = (TObjArray*)(tmpmod[2*mod].At(0));
01752 TObjArray* tmparrv = (TObjArray*)(tmpmod[2*mod].At(1));
01753 tmparru->Add(segu); tmparrv->Add(segv);
01754 }
01755 }
01756
01757 for(mod=0;mod<modnum-1;mod++){
01758 TObjArray* tmparru = (TObjArray*)(tmpmod[2*mod].At(0));
01759 TObjArray* tmparrv = (TObjArray*)(tmpmod[2*mod].At(1));
01760 TObjArray* tmparrup = (TObjArray*)(tmpmod[2*mod+1].At(0));
01761 TObjArray* tmparrvp = (TObjArray*)(tmpmod[2*mod+1].At(1));
01762
01763 if( 1+tmparru->GetLast()>0 && 1+tmparrup->GetLast()>0 ){
01764
01765 cont=1;
01766 while(cont){
01767 cont=0;
01768 idm=-1; idp=-1; topscr=0;
01769 for(i=0;i<1+tmparru->GetLast();i++){
01770 TrackSegmentAtNu* segum = (TrackSegmentAtNu*)(tmparru->At(i));
01771 TrackSegmentAtNu* segvm = (TrackSegmentAtNu*)(tmparrv->At(i));
01772 if( (segum->GetPartner()==segvm && segvm->GetPartner()==segum)
01773 || (segum->GetPartner()==0 && segvm->GetPartner()==0
01774 && segum->GetUID()>2 && segvm->GetUID()>2) ){
01775 for(j=0;j<1+tmparrup->GetLast();j++){
01776 TrackSegmentAtNu* segup = (TrackSegmentAtNu*)(tmparrup->At(j));
01777 TrackSegmentAtNu* segvp = (TrackSegmentAtNu*)(tmparrvp->At(j));
01778 if(segup->GetUID()>2 && segvp->GetUID()>2){
01779 if( segum->IsEndAssoc(segup) && segvm->IsEndAssoc(segvp) ){
01780 if(segum->GetBegPlane()<=segvm->GetBegPlane()){
01781 bpln2=segum->GetBegPlane(); bpln1=segvm->GetBegPlane(); }
01782 if(segvm->GetBegPlane()<=segum->GetBegPlane()){
01783 bpln2=segvm->GetBegPlane(); bpln1=segum->GetBegPlane(); }
01784 if(segup->GetEndPlane()>=segvp->GetEndPlane()){
01785 epln2=segup->GetEndPlane(); epln1=segvp->GetEndPlane(); }
01786 if(segvp->GetEndPlane()>=segup->GetEndPlane()){
01787 epln2=segvp->GetEndPlane(); epln1=segup->GetEndPlane(); }
01788 dpln1=epln2-bpln2; dpln2=epln1-bpln1;
01789 scr=dpln2*dpln2/dpln1;
01790 if(scr>topscr){
01791 idm=i; idp=j; topscr=scr;
01792 }
01793 }
01794 }
01795 }
01796 }
01797 }
01798 if(1+idm>0 && 1+idp>0){
01799 TrackSegmentAtNu* segum = (TrackSegmentAtNu*)(tmparru->At(idm));
01800 TrackSegmentAtNu* segvm = (TrackSegmentAtNu*)(tmparrv->At(idm));
01801 TrackSegmentAtNu* segup = (TrackSegmentAtNu*)(tmparrup->At(idp));
01802 TrackSegmentAtNu* segvp = (TrackSegmentAtNu*)(tmparrvp->At(idp));
01803 segup->AddSegment(segum); segup->SetPartner(segvp);
01804 segum->SetUID(-1); segum->SetPartner(0);
01805 segvp->AddSegment(segvm); segvp->SetPartner(segup);
01806 segvm->SetUID(-1); segvm->SetPartner(0);
01807 cont=1;
01808 }
01809 }
01810 }
01811 }
01812
01813 for(i=0;i<2*(modnum-1);i++){
01814 DELOBJCTR+=1+tmpmod[i].GetLast(); tmpmod[i].Delete();
01815 }
01816 delete [] tmpmod; DELOBJCTR+=2*(modnum-1);
01817 }
01818
01819
01820 cont=1;
01821 while(cont){
01822 cont=0;
01823 id=-1; topscr=0.0;
01824 for(i=0;i<1+tmptrkseg[0].GetLast();i++){
01825 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(tmptrkseg[0].At(i));
01826 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(tmptrkseg[1].At(i));
01827 if( segu->GetPartner()==segv && segv->GetPartner()==segu
01828 && segu->GetUID()<4 && segv->GetUID()<4 ){
01829 tmptrk[0].Add(segu); tmptrk[1].Add(segv);
01830 segu->SetUID(4); segv->SetUID(4);
01831 }
01832 else{
01833 if( segu->GetPartner()==0 && segv->GetPartner()==0
01834 && segu->GetUID()>2 && segv->GetUID()>2 ){
01835 if(segu->GetBegPlane()<=segv->GetBegPlane()){
01836 bpln2=segu->GetBegPlane(); bpln1=segv->GetBegPlane(); }
01837 if(segv->GetBegPlane()<=segu->GetBegPlane()){
01838 bpln2=segv->GetBegPlane(); bpln1=segu->GetBegPlane(); }
01839 if(segu->GetEndPlane()>=segv->GetEndPlane()){
01840 epln2=segu->GetEndPlane(); epln1=segv->GetEndPlane(); }
01841 if(segv->GetEndPlane()>=segu->GetEndPlane()){
01842 epln2=segv->GetEndPlane(); epln1=segu->GetEndPlane(); }
01843 dpln1=epln2-bpln2; dpln2=epln1-bpln1;
01844 scr=dpln2*dpln2/dpln1;
01845 if(scr>topscr){
01846 topscr=scr; id=i;
01847 }
01848 }
01849 }
01850 }
01851 if(1+id>0){
01852 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(tmptrkseg[0].At(id));
01853 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(tmptrkseg[1].At(id));
01854 segu->SetPartner(segv); segv->SetPartner(segu);
01855 cont=1;
01856 }
01857 }
01858
01859 ntrks=1+tmptrk[0].GetLast();
01860 for(vuw=0;vuw<2;vuw++){
01861 for(i=0;i<1+tmptrk[vuw].GetLast();i++){
01862 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
01863 for(j=0;j<1+seg->GetClrLast();j++){
01864 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01865 clr->SetTrkFlag(4);
01866 }
01867 }
01868 }
01869
01870 gettimeofday(&tpafter,0);
01871 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01872 fTime_TrackAtNu+=dtime;
01873
01874 gettimeofday(&tpbefore,0);
01875
01876
01877 TObjArray tmpshw[2],tmpshwseg[2];
01878 for(i=0;i<1+segbnk2[0].GetLast();i++){
01879 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(segbnk2[0].At(i));
01880 if(segu->GetUID()>2){
01881 for(j=0;j<1+segbnk2[1].GetLast();j++){
01882 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(segbnk2[1].At(j));
01883 if(segv->GetUID()>2){
01884 if( ( segu->GetEndPlane()-segv->GetBegPlane()>-2
01885 && segv->GetEndPlane()-segu->GetBegPlane()>-2)
01886 && ( segu->GetEndTime()-segv->GetBegTime()>-fTimeWin
01887 && segv->GetEndTime()-segu->GetBegTime()>-fTimeWin ) ){
01888 tmpshwseg[0].Add(segu); tmpshwseg[1].Add(segv);
01889 }
01890 }
01891 }
01892 }
01893 }
01894
01895
01896 TObjArray tmpclr[2],tmpseg[2];
01897 if(1+tmptrk[0].GetLast()>0){
01898 for(i=0;i<1+tmptrk[0].GetLast();i++){
01899
01900 for(vuw=0;vuw<2;vuw++){
01901 tmpclr[vuw].Clear(); tmpseg[vuw].Clear();
01902 }
01903
01904 for(vuw=0;vuw<2;vuw++){
01905 TrackSegmentAtNu* segtrk = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
01906 for(j=0;j<1+segtrk->GetClrLast();j++){
01907 ClusterAtNu* clr = (ClusterAtNu*)(segtrk->GetClrAt(j));
01908 if(clr->GetShwFlag()>2) tmpclr[vuw].Add(clr);
01909 }
01910 }
01911
01912 if(1+tmpclr[0].GetLast()>0 && 1+tmpclr[1].GetLast()>0){
01913 for(j=0;j<1+tmpshwseg[0].GetLast();j++){
01914 for(vuw=0;vuw<2;vuw++){
01915 jnflag2[vuw]=0;
01916 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmpshwseg[vuw].At(j));
01917 if(seg->GetPartner()==0){
01918 for(k=0;k<1+tmpclr[vuw].GetLast();k++){
01919 ClusterAtNu* clr = (ClusterAtNu*)(tmpclr[vuw].At(k));
01920 if(seg->ContainsClr(clr)) jnflag2[vuw]=1;
01921 }
01922 }
01923 }
01924 if(jnflag2[0]>0 && jnflag2[1]>0){
01925 for(vuw=0;vuw<2;vuw++){
01926 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmpshwseg[vuw].At(j));
01927 tmpseg[vuw].Add(seg);
01928 }
01929 }
01930 }
01931 }
01932
01933 if(1+tmpseg[0].GetLast()>0 && 1+tmpseg[1].GetLast()>0){
01934 cont=1;
01935 while(cont){
01936 cont=0;
01937 id=-1; topscr=0.0;
01938 for(j=0;j<1+tmpseg[0].GetLast();j++){
01939 ShowerSegmentAtNu* segshwu = (ShowerSegmentAtNu*)(tmpseg[0].At(j));
01940 ShowerSegmentAtNu* segshwv = (ShowerSegmentAtNu*)(tmpseg[1].At(j));
01941 if( segshwu->GetPartner()==0 && segshwv->GetPartner()==0
01942 && segshwu->GetUID()>2 && segshwv->GetUID()>2 ){
01943 dpln1=segshwu->GetEndPlane()-segshwv->GetBegPlane();
01944 dpln2=segshwv->GetEndPlane()-segshwu->GetBegPlane();
01945 if(dpln1<=dpln2) scr=dpln1; if(dpln2<=dpln1) scr=dpln2;
01946 if(scr>topscr){
01947 id=j; scr=topscr;
01948 }
01949 }
01950 }
01951 if(1+id>0){
01952 ShowerSegmentAtNu* segshwu = (ShowerSegmentAtNu*)(tmpseg[0].At(id));
01953 ShowerSegmentAtNu* segshwv = (ShowerSegmentAtNu*)(tmpseg[1].At(id));
01954 segshwu->SetPartner(segshwv); segshwv->SetPartner(segshwu);
01955 cont=1;
01956 }
01957 }
01958 }
01959 }
01960 }
01961
01962
01963 cont=1;
01964 while(cont){
01965 cont=0;
01966 id=-1; topscr=0.0;
01967 for(i=0;i<1+tmpshwseg[1].GetLast();i++){
01968 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpshwseg[0].At(i));
01969 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpshwseg[1].At(i));
01970 if( segu->GetPartner()==segv && segv->GetPartner()==segu
01971 && segu->GetUID()<4 && segv->GetUID()<4 ){
01972 tmpshw[0].Add(segu); tmpshw[1].Add(segv);
01973 segu->SetUID(4); segv->SetUID(4);
01974 }
01975 else{
01976 if( segu->GetPartner()==0 && segv->GetPartner()==0
01977 && segu->GetUID()>2 && segv->GetUID()>2 ){
01978 dpln1 = segu->GetEndPlane()-segv->GetBegPlane();
01979 dpln2 = segv->GetEndPlane()-segu->GetBegPlane();
01980 if(dpln1<=dpln2) scr=dpln1; if(dpln2<=dpln1) scr=dpln2;
01981 if(scr>topscr){
01982 topscr=scr; id=i;
01983 }
01984 }
01985 }
01986 }
01987 if(1+id>0){
01988 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpshwseg[0].At(id));
01989 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpshwseg[1].At(id));
01990 segu->SetPartner(segv); segv->SetPartner(segu);
01991 cont=1;
01992 }
01993 }
01994
01995 nshws=1+tmpshw[0].GetLast();
01996 for(vuw=0;vuw<2;vuw++){
01997 for(i=0;i<1+tmpshw[vuw].GetLast();i++){
01998 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmpshw[vuw].At(i));
01999 for(j=0;j<1+seg->GetClrLast();j++){
02000 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
02001 clr->SetShwFlag(4);
02002 }
02003 }
02004 }
02005
02006 gettimeofday(&tpafter,0);
02007 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02008 fTime_ShowerAtNu+=dtime;
02009
02010
02011
02012
02013
02014
02015 if( 1+tmpshw[0].GetLast()==0 && 1+tmpshw[1].GetLast()==0
02016 && 1+tmptrk[0].GetLast()==0 && 1+tmptrk[1].GetLast()==0 ){
02017
02018 gettimeofday(&tpbefore,0);
02019
02020
02021 MSG("AlgAtNuReco",Msg::kDebug) << " reseeding low muliplicity showers " << endl;
02022
02023 for(vuw=0;vuw<2;vuw++){
02024 tmpclr[vuw].Clear(); tmpseg[vuw].Clear(); tmpshwseg[vuw].Clear();
02025 }
02026
02027
02028 for(mod=0;mod<modnum;mod++){
02029 for(pln=1;pln<1+modpln;pln++){
02030 i=pln+mod*(modpln+1);
02031 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
02032 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
02033 clr0->SetShwFlag(0);
02034 npts=0;
02035 for(k=-4;k<5;k++){
02036 k0=i+k;
02037 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
02038 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
02039 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
02040 if(clr0->GetCharge()>1.0 && clr1->GetCharge()>1.0){
02041 npts+=1+clr1->GetHitLast();
02042 }
02043 }
02044 }
02045 }
02046 if(npts>3){
02047 vuw=clr0->GetPlaneView(); tmpclr[vuw].Add(clr0); clr0->SetShwFlag(1);
02048 }
02049 }
02050 }
02051 }
02052
02053
02054 if(1+tmpclr[0].GetLast()>0 && 1+tmpclr[1].GetLast()>0){
02055 for(vuw=0;vuw<2;vuw++){
02056 for(i=0;i<1+tmpclr[vuw].GetLast();i++){
02057 ClusterAtNu* clr0 = (ClusterAtNu*)(tmpclr[vuw].At(i));
02058 if(clr0->GetShwFlag()>0 && clr0->GetShwFlag()<2){
02059 ShowerSegmentAtNu* seg0 = new ShowerSegmentAtNu(clr0); NEWOBJCTR++;
02060 tmpseg[vuw].Add(seg0); shwsegbnk.Add(seg0);
02061 seg0->SetUID(2); seg0->SetReseedFlag(1);
02062 mod=(Int_t)((clr0->GetPlane())/(modpln+1)); clr0->SetShwFlag(2);
02063 ctr=1;
02064 while(ctr>0){
02065 ctr=0;
02066 pln0=-6+seg0->GetBegPlane();
02067 npln=7+(seg0->GetEndPlane()-seg0->GetBegPlane())/2;
02068 for(k=0;k<npln;k++){
02069 k0=pln0+2*k;
02070 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
02071 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
02072 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
02073 if( clr1->GetShwFlag()>0 && clr1->GetShwFlag()<2
02074 && seg0->IsDiffuseShwAssoc(clr1)>0 ){
02075 seg0->AddCluster(clr1);
02076 clr1->SetShwFlag(2); ctr++;
02077 }
02078 }
02079 }
02080 }
02081 }
02082 }
02083 }
02084 }
02085 }
02086
02087 MSG("AlgAtNuReco",Msg::kVerbose) << " low multiplicity 2D showers - U=" << 1+tmpseg[0].GetLast()
02088 << " V=" << 1+tmpseg[1].GetLast() << endl;
02089
02090
02091 if( 1+tmpseg[0].GetLast()>0 && 1+tmpseg[1].GetLast()>0 ){
02092
02093 for(i=0;i<1+tmpseg[0].GetLast();i++){
02094 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpseg[0].At(i));
02095 for(j=0;j<1+tmpseg[1].GetLast();j++){
02096 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpseg[1].At(j));
02097 if( 1+segu->GetHitLast()>2 || 1+segv->GetHitLast()>2
02098 && segu->GetEndPlane()-segv->GetBegPlane()>-4
02099 && segv->GetEndPlane()-segu->GetBegPlane()>-4 ){
02100 tmpshwseg[0].Add(segu); tmpshwseg[1].Add(segv);
02101 }
02102 }
02103 }
02104
02105 id=-1; topscr=0.0;
02106 for(i=0;i<1+tmpshwseg[0].GetLast();i++){
02107 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpshwseg[0].At(i));
02108 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpshwseg[1].At(i));
02109 scr=1+segu->GetHitLast()+segv->GetHitLast();
02110 if(scr>topscr){
02111 topscr=scr; id=i;
02112 }
02113 }
02114
02115 if(1+id>0){
02116 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpshwseg[0].At(id));
02117 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpshwseg[1].At(id));
02118 segu->SetPartner(segv); segv->SetPartner(segu);
02119 segu->SetUID(4); segv->SetUID(4);
02120 tmpshw[0].Add(segu); tmpshw[1].Add(segv);
02121 }
02122
02123 }
02124
02125 gettimeofday(&tpafter,0);
02126 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02127 fTime_ShowerAtNu+=dtime;
02128 }
02129
02130
02131
02132
02133
02134
02135 gettimeofday(&tpbefore,0);
02136
02137
02138 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of 3D tracks *** " << endl;
02139 Int_t km1,kp1;
02140 Double_t swx,swy,swx2,swxy,sw;
02141 Double_t w,x,y,m,c;
02142 Int_t ustrps,vstrps,trkstrps,trkplns;
02143 Int_t bplnvuw[2],eplnvuw[2];
02144 for(i=0;i<1+tmptrk[0].GetLast();i++){
02145
02146 MSG("AlgAtNuReco", Msg::kVerbose) << " making track " << i << endl;
02147
02148
02149 bpln2=-1; epln2=-1;
02150 for(vuw=0;vuw<2;vuw++){
02151 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
02152 if(bpln2<0 || seg->GetBegPlane()<bpln2) bpln2=seg->GetBegPlane();
02153 if(epln2<0 || seg->GetEndPlane()>epln2) epln2=seg->GetEndPlane();
02154 }
02155
02156 bpln1=-1; epln1=-1;
02157 for(vuw=0;vuw<2;vuw++){
02158 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
02159 if(bpln1<0 || seg->GetBegPlane()>bpln2) bpln1=seg->GetBegPlane();
02160 if(epln1<0 || seg->GetEndPlane()<epln2) epln1=seg->GetEndPlane();
02161 }
02162
02163 for(vuw=0;vuw<2;vuw++){
02164 bplnvuw[vuw]=-1; eplnvuw[vuw]=-1;
02165 }
02166
02167 MSG("AlgAtNuReco", Msg::kVerbose)
02168 << bpln2 << " (" << bpln1 << ") -> (" << epln1 << ") " << epln2 << endl;
02169
02170
02171 npln=1+epln2-bpln2;
02172 Int_t* strlst = new Int_t[npln]; NEWOBJCTR+=npln;
02173 Int_t* vuwlst = new Int_t[npln]; NEWOBJCTR+=npln;
02174 TObjArray* clrlst = new TObjArray[npln]; NEWOBJCTR+=npln;
02175 TObjArray* hitlst = new TObjArray[npln]; NEWOBJCTR+=npln;
02176 for(pln=0;pln<npln;pln++){
02177 strlst[pln]=-1; vuwlst[pln]=-1;
02178 }
02179
02180 for(vuw=0;vuw<2;vuw++){
02181 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
02182 for(j=0;j<1+seg->GetClrLast();j++){
02183 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
02184 pln = clr->GetPlane()-bpln2;
02185 if( pln>=0 && pln<npln ){
02186 clrlst[pln].Add(clr);
02187 }
02188 }
02189 }
02190
02191
02192 bpln1=bpln1-4;
02193 km=-1; pln=bpln1-bpln2; if(pln<0) pln=0;
02194 while(km<0 && pln<npln){
02195 if(1+clrlst[pln].GetLast()>0) km=pln; else pln++;
02196 }
02197 if(km>-1 && km+4<npln && 1+clrlst[km+2].GetLast()==0){
02198 if(1+clrlst[km+4].GetLast()>0){
02199 ClusterAtNu* clr = (ClusterAtNu*)(clrlst[km+4].First());
02200 if(clr->GetTrkPlnFlag()>0 && clr->GetShwPlnFlag()==0) km=bpln1-bpln2;
02201 }
02202 bpln1=bpln2+km;
02203 }
02204
02205 epln1=epln1+4;
02206 kp=-1; pln=epln1-bpln2; if(pln>npln-1) pln=npln-1;
02207 while(kp<0 && pln>-1){
02208 if(1+clrlst[pln].GetLast()>0) kp=pln; else pln--;
02209 }
02210 if(kp>-1 && kp-4>-1 && 1+clrlst[kp-2].GetLast()==0){
02211 if(1+clrlst[kp-4].GetLast()>0){
02212 ClusterAtNu* clr = (ClusterAtNu*)(clrlst[kp-4].First());
02213 if(clr->GetTrkPlnFlag()>0 && clr->GetShwPlnFlag()==0) kp=epln1-bpln2;
02214 }
02215 epln1=bpln2+kp;
02216 }
02217
02218 ustrps=0; vstrps=0; trkstrps=0;
02219 for(pln=0;pln<npln;pln++){
02220 if( bpln2+pln>bpln1 && bpln2+pln<epln1
02221 && 1+clrlst[pln].GetLast()>0 ){
02222 ClusterAtNu* clr = (ClusterAtNu*)(clrlst[pln].First());
02223 vuw = clr->GetPlaneView();
02224 vuwlst[pln]=vuw; strlst[pln]=1;
02225 if(bplnvuw[vuw]<0||bpln2+pln<bplnvuw[vuw]) bplnvuw[vuw]=bpln2+pln;
02226 if(eplnvuw[vuw]<0||bpln2+pln>eplnvuw[vuw]) eplnvuw[vuw]=bpln2+pln;
02227 if(vuw==0) ustrps++; if(vuw==1) vstrps++;
02228 if(clr->GetTrkPlnFlag()>0) trkstrps++;
02229 }
02230 }
02231
02232 MSG("AlgAtNuReco", Msg::kVerbose)
02233 << " U: (" << bplnvuw[0] << "," << eplnvuw[0] << ") "
02234 << " V: (" << bplnvuw[1] << "," << eplnvuw[1] << ") " << endl;
02235
02236 if( ustrps>2 && vstrps>2 && trkstrps>2 ){
02237
02238
02239 for(vuw=0;vuw<2;vuw++){
02240 for(pln=0;pln<npln;pln++){
02241 if(vuwlst[pln]==vuw && strlst[pln]>0){
02242
02243 ClusterAtNu* clr0 = (ClusterAtNu*)(clrlst[pln].First());
02244 k=0; km1=-1;
02245 while(pln-k>0 && km1<0){
02246 k++; if(vuwlst[pln-k]==vuw) km1=k;
02247 }
02248 k=0; kp1=-1;
02249 while(pln+k<npln-1 && kp1<0){
02250 k++; if(vuwlst[pln+k]==vuw) kp1=k;
02251 }
02252
02253 if(km1>0 && kp1>0){
02254 ClusterAtNu* clrp1 = (ClusterAtNu*)(clrlst[pln+kp1].First());
02255 ClusterAtNu* clrm1 = (ClusterAtNu*)(clrlst[pln-km1].First());
02256 if( ( clrm1->GetPlane()>bplnvuw[vuw]
02257 || ( clr0->GetPlane()-clrm1->GetPlane()<3
02258 && clr0->GetTrkPlnFlag()>0 ) )
02259 && ( clrp1->GetPlane()<eplnvuw[vuw]
02260 || ( clrp1->GetPlane()-clr0->GetPlane()<3
02261 && clr0->GetTrkPlnFlag()>0 ) ) ){
02262 if(clr0->IsTrkAssoc(clrm1,clrp1)>1){
02263 strlst[pln]=3;
02264 if( ( clrp1->GetBegStrip()-clr0->GetEndStrip()>-1
02265 && clr0->GetBegStrip()-clrm1->GetEndStrip()>-1 )
02266 || ( clrp1->GetEndStrip()-clr0->GetBegStrip()<1
02267 && clr0->GetEndStrip()-clrm1->GetBegStrip()<1 ) ){
02268 if( clr0->GetPlane()-clrm1->GetPlane()<3
02269 && clrm1->GetTrkPlnFlag()>-1 ){ strlst[pln-km1]=3; }
02270 if( clrp1->GetPlane()-clr0->GetPlane()<3
02271 && clrp1->GetTrkPlnFlag()>-1 ){ strlst[pln+kp1]=3; }
02272 }
02273 }
02274 else{
02275 if( 1+clr0->GetEndStrip()-clr0->GetBegStrip()<3
02276 && clr0->GetTrkPlnFlag()>0 ){ strlst[pln]=3; }
02277 }
02278 }
02279 }
02280
02281 }
02282 }
02283 }
02284
02285 MSG("AlgAtNuReco", Msg::kVerbose) << " *** SORT TRACK HITS *** " << endl;
02286 for(pln=0;pln<npln;pln++){
02287 if( vuwlst[pln]>-1 ){
02288 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrlst[pln].First());
02289 MSG("AlgAtNuReco", Msg::kVerbose) << " " << pln << " " << clrtmp->GetPlane()
02290 << " " << clrtmp->GetBegStrip() << " " << clrtmp->GetEndStrip()
02291 << " " << strlst[pln] << endl;
02292 }
02293 }
02294
02295
02296 for(vuw=0;vuw<2;vuw++){
02297 for(pln=0;pln<npln;pln++){
02298 if(vuwlst[pln]==vuw && strlst[pln]<2){
02299
02300 MSG("AlgAtNuReco", Msg::kVerbose) << " *** GLOBAL FIT *** pln=" << pln << endl;
02301
02302 k=0; kp=-1;
02303 while(pln+k<epln1-bpln2 && pln+k<npln-1 && kp<0){
02304 k++; if(vuwlst[pln+k]==vuw && strlst[pln+k]>2) kp=k;
02305 }
02306
02307 k=0; km=-1;
02308 while(pln-k>bpln1-bpln2 && pln-k>0 && km<0){
02309 k++; if(vuwlst[pln-k]==vuw && strlst[pln-k]>2) km=k;
02310 }
02311
02312 if(km>0){ km1=0-km-7; k1=0-km; } else{ km1=(bpln1-bpln2)-pln+1; k1=(bpln1-bpln2)-pln+1; }
02313 if(pln+km1<0){ km1=0-pln; } if(pln+k1<0){ k1=0-pln; }
02314 if(kp>0){ kp1=1+kp+7; k2=1+kp; } else{ kp1=(epln1-bpln2)-pln; k2=(epln1-bpln2)-pln; }
02315 if(pln+kp1>npln){ kp1=npln-pln; } if(pln+k2>npln){ k2=npln-pln; }
02316
02317 ClusterAtNu* clrseed = (ClusterAtNu*)(clrlst[pln].First());
02318 m=0.0; c=0.5*(clrseed->GetBegStrip()+clrseed->GetEndStrip());
02319 swx=0.0; swy=0.0; swx2=0.0; swxy=0.0; sw=0.0;
02320 for(k=km1;k<kp1;k++){
02321 if(vuwlst[pln+k]==vuw){
02322 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrlst[pln+k].First());
02323 for(l=0;l<1+clrtmp->GetHitLast();l++){
02324 HitAtNu* hittmp = (HitAtNu*)(clrtmp->GetHitAt(l));
02325 x=hittmp->GetZPos(); y=hittmp->GetStrip();
02326 w=hittmp->GetCharge()/clrtmp->GetCharge();
02327 if(strlst[pln+k]>2) w=5.0*w; else w=0.5*w;
02328 swx+=w*x; swy+=w*y; swx2+=w*x*x; swxy+=w*x*y; sw+=w;
02329 }
02330 }
02331 }
02332 if(sw>0.0){
02333 m=(sw*swxy-swx*swy)/(sw*swx2-swx*swx);
02334 c=(swy*swx2-swx*swxy)/(sw*swx2-swx*swx);
02335 }
02336
02337 for(k=k1;k<k2;k++){
02338 if(vuwlst[pln+k]==vuw && strlst[pln+k]<2){
02339 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrlst[pln+k].First());
02340 strw = m*clrtmp->GetZPos()+c;
02341 if(m>=0.0) strdw = 0.6+0.02*m; if(m<=0) strdw = 0.6-0.02*m;
02342 MSG("AlgAtNuReco", Msg::kVerbose) << " pln=" << pln+k << " strw=" << strw << " strdw=" << strdw << " (" << clrtmp->GetBegStrip() << " " << clrtmp->GetEndStrip() << ") " << endl;
02343 if( (km>0 && kp>0) || (km<0 && kp<0)
02344 || (clrtmp->GetEndStrip()>strw-1.6 && clrtmp->GetBegStrip()<strw+1.6) ){
02345 id=-1; dscr=0.0; scr=999.9;
02346 for(l=0;l<1+clrtmp->GetHitLast();l++){
02347 HitAtNu* hit = (HitAtNu*)(clrtmp->GetHitAt(l));
02348 dscr=hit->GetStrip()-strw; if(dscr<0) dscr=-dscr;
02349 if(dscr<scr){
02350 id=l; scr=dscr;
02351 }
02352 }
02353 if(1+id>0){
02354 HitAtNu* hit = (HitAtNu*)(clrtmp->GetHitAt(id));
02355 ClusterAtNu* clrnew = new ClusterAtNu(hit); clrbnk.Add(clrnew); NEWOBJCTR++;
02356 clrlst[pln+k].Add(clrnew); strlst[pln+k]=2;
02357 for(l=0;l<1+clrtmp->GetHitLast();l++){
02358 HitAtNu* hittmp = (HitAtNu*)(clrtmp->GetHitAt(l));
02359 if( ( hittmp->GetStrip()-hit->GetStrip()>0
02360 && hittmp->GetStrip()-hit->GetStrip()<strdw)
02361 || ( hittmp->GetStrip()-hit->GetStrip()<0
02362 && hittmp->GetStrip()-hit->GetStrip()>-strdw) ){
02363 clrnew->AddHit(hittmp);
02364 }
02365 }
02366 for(l=0;l<1+clrnew->GetHitLast();l++){
02367 HitAtNu* hittmp = (HitAtNu*)(clrnew->GetHitAt(l));
02368 MSG("AlgAtNuReco", Msg::kVerbose) << " ... add hit=" << hittmp->GetStrip() << " (" << pln+k << ") " << endl;
02369 }
02370 }
02371 }
02372 }
02373 }
02374 }
02375
02376 }
02377 }
02378
02379
02380 for(vuw=0;vuw<2;vuw++){
02381 for(pln=0;pln<npln;pln++){
02382 if(vuwlst[pln]==vuw && strlst[pln]>1){
02383
02384 MSG("AlgAtNuReco", Msg::kVerbose) << " *** FIT (2) *** pln=" << pln << endl;
02385
02386 ClusterAtNu* clr = (ClusterAtNu*)(clrlst[pln].Last());
02387 strw=0.5*(clr->GetBegStrip()+clr->GetEndStrip());
02388 strdw=0.6;
02389
02390 if(1+clr->GetHitLast()>1){
02391 for(k=0;k<1+clr->GetHitLast();k++){
02392 HitAtNu* hit = (HitAtNu*)(clr->GetHitAt(k));
02393 hit->SetTrkFlag(1);
02394 }
02395 swx=0.0; swy=0.0; swx2=0.0; swxy=0.0; sw=0.0;
02396 for(k=-4;k<5;k++){
02397 if( pln+k>=0 && pln+k<npln
02398 && vuwlst[pln+k]==vuw && strlst[pln+k]>1 ){
02399 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrlst[pln+k].Last());
02400 for(l=0;l<1+clrtmp->GetHitLast();l++){
02401 HitAtNu* hittmp = (HitAtNu*)(clrtmp->GetHitAt(l));
02402 x=hittmp->GetZPos(); y=hittmp->GetStrip();
02403 w=hittmp->GetCharge()/clrtmp->GetCharge();
02404 swx+=w*x; swy+=w*y; swx2+=w*x*x; swxy+=w*x*y; sw+=w;
02405 }
02406 }
02407 }
02408 if(sw>1.0){
02409 m=(sw*swxy-swx*swy)/(sw*swx2-swx*swx);
02410 c=(swy*swx2-swx*swxy)/(sw*swx2-swx*swx);
02411 strw=m*clr->GetZPos()+c;
02412 if(m>=0.0) strdw = 0.6+0.02*m; if(m<=0) strdw = 0.6-0.02*m;
02413 }
02414 }
02415
02416 MSG("AlgAtNuReco", Msg::kVerbose) << " pln=" << pln << " strw=" << strw << " strdw=" << strdw << endl;
02417
02418 id=-1; dscr=0.0; scr=999.0;
02419 for(k=0;k<1+clr->GetHitLast();k++){
02420 HitAtNu* hit = (HitAtNu*)(clr->GetHitAt(k));
02421 dscr=hit->GetStrip()-strw; if(dscr<0) dscr=-dscr;
02422 if(dscr<scr){
02423 id=k; scr=dscr;
02424 }
02425 }
02426 if(1+id>0){
02427 HitAtNu* hit = (HitAtNu*)(clr->GetHitAt(id));
02428 hitlst[pln].Add(hit); strlst[pln]=4;
02429 for(k=0;k<1+clr->GetHitLast();k++){
02430 HitAtNu* hittmp = (HitAtNu*)(clr->GetHitAt(k));
02431 if( ( hittmp->GetStrip()-hit->GetStrip()>0
02432 && hittmp->GetStrip()-hit->GetStrip()<strdw)
02433 || ( hittmp->GetStrip()-hit->GetStrip()<0
02434 && hittmp->GetStrip()-hit->GetStrip()>-strdw) ){
02435 hitlst[pln].Add(hittmp);
02436 }
02437 }
02438 }
02439 for(l=0;l<1+hitlst[pln].GetLast();l++){
02440 HitAtNu* hittmp = (HitAtNu*)(hitlst[pln].At(l));
02441 MSG("AlgAtNuReco", Msg::kVerbose) << " ... add hit=" << hittmp->GetStrip() << " (" << pln << ") " << endl;
02442 }
02443 }
02444 }
02445 }
02446
02447
02448 trkplns=0;
02449 for(pln=0;pln<npln;pln++){
02450 if( strlst[pln]>3 ){
02451 trkplns++;
02452 }
02453 }
02454
02455 if(trkplns>5){
02456 for(vuw=0;vuw<2;vuw++){
02457 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
02458 TrackAtNu* trk = new TrackAtNu(slice); trkbnk[vuw].Add(trk); NEWOBJCTR++;
02459 if(seg->GetReseedFlag()||fObjMaxFlag) jnflag2[vuw]=1; else jnflag2[vuw]=0;
02460 for(pln=0;pln<npln;pln++){
02461 if(vuwlst[pln]==vuw && strlst[pln]>3){
02462 for(k=0;k<1+hitlst[pln].GetLast();k++){
02463 HitAtNu* hit = (HitAtNu*)(hitlst[pln].At(k));
02464 hit->SetTrkFlag(2); trk->AddHit(hit);
02465 }
02466 }
02467 }
02468 }
02469 TrackAtNu* trku = (TrackAtNu*)(trkbnk[0].Last());
02470 TrackAtNu* trkv = (TrackAtNu*)(trkbnk[1].Last());
02471 trku->SetPartner(trkv); trkv->SetPartner(trku);
02472
02473 if(jnflag2[0]>0 && jnflag2[1]>0){
02474 trku->SetReseedFlag(1); trkv->SetReseedFlag(1);
02475 }
02476 }
02477
02478 }
02479
02480 delete [] strlst; DELOBJCTR+=npln;
02481 delete [] vuwlst; DELOBJCTR+=npln;
02482 delete [] clrlst; DELOBJCTR+=npln;
02483 delete [] hitlst; DELOBJCTR+=npln;
02484 }
02485
02486
02487 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF TRACKS *** " << endl;
02488 for(vuw=0;vuw<2;vuw++){
02489 for(i=0;i<1+trkbnk[vuw].GetLast();i++){
02490 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02491 MSG("AlgAtNuReco", Msg::kVerbose)
02492 << " " << vuw << " " << i
02493 << " " << trk->GetBegPlane() << " " << trk->GetBegStrip()
02494 << " " << trk->GetEndPlane() << " " << trk->GetEndStrip() << endl;
02495 }
02496 }
02497
02498 gettimeofday(&tpafter,0);
02499 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02500 fTime_TrackAtNu+=dtime;
02501
02502
02503
02504
02505
02506
02507 gettimeofday(&tpbefore,0);
02508
02509
02510 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of 3D showers *** " << endl;
02511 for(i=0;i<1+tmpshw[0].GetLast();i++){
02512
02513 MSG("AlgAtNuReco", Msg::kVerbose) << " making shower " << i << endl;
02514
02515 for(vuw=0;vuw<2;vuw++){
02516 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmpshw[vuw].At(i));
02517 ShowerAtNu* shw = new ShowerAtNu(slice); shwbnk[vuw].Add(shw); NEWOBJCTR++;
02518 if(seg->GetReseedFlag()||fObjMaxFlag) jnflag2[vuw]=1; else jnflag2[vuw]=0;
02519 for(j=0;j<1+seg->GetHitLast();j++){
02520 HitAtNu* hit = (HitAtNu*)(seg->GetHitAt(j));
02521 shw->AddHit(hit); hit->SetShwFlag(2);
02522 }
02523 }
02524 ShowerAtNu* shwu = (ShowerAtNu*)(shwbnk[0].Last());
02525 ShowerAtNu* shwv = (ShowerAtNu*)(shwbnk[1].Last());
02526 shwu->SetPartner(shwv); shwv->SetPartner(shwu);
02527
02528 if(jnflag2[0]>0 && jnflag2[1]>0){
02529 shwu->SetReseedFlag(1); shwv->SetReseedFlag(1);
02530 }
02531 }
02532
02533
02534
02535
02536
02537
02538 if(fVtxFlag){
02539 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of vertex showers *** " << endl;
02540
02541 Int_t shwflag=0;
02542 for(i=0;i<1+shwbnk[0].GetLast();i++){
02543 if(shwflag==0){
02544 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[0].At(i));
02545 if(shw->GetReseedFlag()||1) shwflag=1;
02546 }
02547 }
02548
02549 TObjArray tmphittrk[2],tmphitbeg[2],tmphitend[2];
02550 for(i=0;i<1+trkbnk[0].GetLast();i++){
02551
02552
02553 for(vuw=0;vuw<2;vuw++){
02554 tmpclr[vuw].Clear();
02555 tmphittrk[vuw].Clear(); tmphitbeg[vuw].Clear();
02556 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02557 for(j=0;j<1+trk->GetHitLast();j++){
02558 HitAtNu* hit = (HitAtNu*)(trk->GetHitAt(j));
02559 if( hit->GetPlane()<trk->GetBegPlane()+9
02560 && 2*hit->GetPlane()<trk->GetBegPlane()+trk->GetEndPlane()+1
02561 && (Int_t)((hit->GetPlane())/(modpln+1))==(Int_t)((trk->GetBegPlane())/(modpln+1)) ){
02562 tmphittrk[vuw].Add(hit);
02563 }
02564 }
02565 }
02566
02567 for(j=0;j<1+shwbnk[0].GetLast();j++){
02568 for(vuw=0;vuw<2;vuw++){
02569 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02570 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(j));
02571 jnflag2[vuw]=0;
02572 if( shw->GetBegPlane()<trk->GetBegPlane()+9
02573 && 2*shw->GetBegPlane()<trk->GetBegPlane()+trk->GetEndPlane()+1
02574 && (Int_t)((shw->GetBegPlane())/(modpln+1))==(Int_t)((trk->GetBegPlane())/(modpln+1)) ){
02575 for(k=0;k<1+tmphittrk[vuw].GetLast();k++){
02576 if(jnflag2[vuw]==0){
02577 HitAtNu* hit = (HitAtNu*)(tmphittrk[vuw].At(k));
02578 if(shw->IsDiffuseShwAssoc(hit)>0) jnflag2[vuw]=1;
02579 }
02580 }
02581 }
02582 }
02583 if(jnflag2[0]>0 || jnflag2[1]>0){
02584 for(vuw=0;vuw<2;vuw++){
02585 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(j));
02586 tmpclr[vuw].Add(shw);
02587 }
02588 }
02589 }
02590
02591 if( 1+tmpclr[0].GetLast()>0 && 1+tmpclr[1].GetLast()>0 ){
02592 id=-1; topscr=9999.9;
02593 for(j=0;j<1+tmpclr[0].GetLast();j++){
02594 ShowerAtNu* shwu = (ShowerAtNu*)(tmpclr[0].At(j));
02595 ShowerAtNu* shwv = (ShowerAtNu*)(tmpclr[1].At(j));
02596 scr = shwu->GetBegPlane()+shwv->GetBegPlane();
02597 if(scr<topscr){
02598 topscr=scr; id=j;
02599 }
02600 }
02601 if(1+id>0){
02602 MSG("AlgAtNuReco", Msg::kVerbose) << " found existing beginning vertex shower " << endl;
02603 for(vuw=0;vuw<2;vuw++){
02604 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02605 ShowerAtNu* shw = (ShowerAtNu*)(tmpclr[vuw].At(id));
02606 trk->SetBegVtxShower(shw); shw->AddAssocTrack(trk); shw->SetBegVtxFlag(1);
02607 }
02608 }
02609 }
02610
02611
02612 else{
02613 if( shwflag==0 && i==0 ){
02614 for(vuw=0;vuw<2;vuw++){
02615 tmphitbeg[vuw].Clear();
02616 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02617
02618 mod=(Int_t)((trk->GetBegPlane())/(modpln+1));
02619 for(j=-3;j<4;j++){
02620 pln=trk->GetBegPlane()+2*j;
02621 if( 2*pln<trk->GetBegPlane()+trk->GetEndPlane()+1
02622 && pln>mod*(modpln+1) && pln<(mod+1)*(modpln+1) ){
02623 for(k=0;k<1+hitplnbnk[pln].GetLast();k++){
02624 HitAtNu* hit = (HitAtNu*)(hitplnbnk[pln].At(k));
02625 if( hit->GetShwFlag()<2
02626 && ( hit->GetTrkFlag()<2 || hit->GetCharge()>50 ) ){
02627 jnflag=0;
02628 for(l=0;l<1+tmphittrk[vuw].GetLast();l++){
02629 if(jnflag==0){
02630 HitAtNu* hittrk = (HitAtNu*)(tmphittrk[vuw].At(l));
02631 if( hittrk->GetPlane()<trk->GetBegPlane()+7
02632 && hittrk->IsDiffuseShwAssoc(hit)>0 ){
02633 tmphitbeg[vuw].Add(hit); hit->SetShwFlag(1); jnflag=1;
02634 }
02635 }
02636 }
02637 }
02638 }
02639 }
02640 }
02641
02642 if(1+tmphitbeg[vuw].GetLast()>0){
02643 bpln=trk->GetBegPlane()-6; epln=trk->GetBegPlane()+6;
02644 ctr=1;
02645 while(ctr>0){
02646 ctr=0;
02647 pln0=bpln-4; npln=5+(epln-bpln)/2;
02648 for(j=0;j<npln;j++){
02649 pln=pln0+2*j;
02650 if(pln>mod*(modpln+1) && pln<(mod+1)*(modpln+1)){
02651 for(k=0;k<1+hitplnbnk[pln].GetLast();k++){
02652 HitAtNu* hit = (HitAtNu*)(hitplnbnk[pln].At(k));
02653 if(hit->GetShwFlag()<1 && hit->GetTrkFlag()<1){
02654 jnflag=0;
02655 for(l=0;l<1+tmphitbeg[vuw].GetLast();l++){
02656 if(jnflag==0){
02657 HitAtNu* hitshw = (HitAtNu*)(tmphitbeg[vuw].At(l));
02658 if( hitshw->IsDiffuseShwAssoc(hit)>0 ){
02659 tmphitbeg[vuw].Add(hit); hit->SetShwFlag(1);
02660 if(hit->GetPlane()<bpln) bpln=hit->GetPlane();
02661 if(hit->GetPlane()>epln) epln=hit->GetPlane();
02662 jnflag=1; ctr++;
02663 }
02664 }
02665 }
02666 }
02667 }
02668 }
02669 }
02670 }
02671 for(j=0;j<1+tmphitbeg[vuw].GetLast();j++){
02672 HitAtNu* hit = (HitAtNu*)(tmphitbeg[vuw].At(j));
02673 hit->SetShwFlag(0);
02674 }
02675 }
02676
02677 }
02678 }
02679 }
02680
02681
02682 for(vuw=0;vuw<2;vuw++){
02683 tmpclr[vuw].Clear();
02684 tmphittrk[vuw].Clear(); tmphitend[vuw].Clear();
02685 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02686 for(j=0;j<1+trk->GetHitLast();j++){
02687 HitAtNu* hit = (HitAtNu*)(trk->GetHitAt(j));
02688 if( hit->GetPlane()>trk->GetEndPlane()-9
02689 && 2*hit->GetPlane()>trk->GetBegPlane()+trk->GetEndPlane()-1
02690 && (Int_t)((hit->GetPlane())/(modpln+1))==(Int_t)((trk->GetEndPlane())/(modpln+1)) ){
02691 tmphittrk[vuw].Add(hit);
02692 }
02693 }
02694 }
02695
02696 for(j=0;j<1+shwbnk[0].GetLast();j++){
02697 for(vuw=0;vuw<2;vuw++){
02698 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02699 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(j));
02700 jnflag2[vuw]=0;
02701 if( shw->GetEndPlane()>trk->GetEndPlane()-9
02702 && 2*shw->GetEndPlane()>trk->GetBegPlane()+trk->GetEndPlane()-1
02703 && (Int_t)((shw->GetEndPlane())/(modpln+1))==(Int_t)((trk->GetEndPlane())/(modpln+1)) ){
02704 for(k=0;k<1+tmphittrk[vuw].GetLast();k++){
02705 if(jnflag2[vuw]==0){
02706 HitAtNu* hit = (HitAtNu*)(tmphittrk[vuw].At(k));
02707 if(shw->IsDiffuseShwAssoc(hit)>0) jnflag2[vuw]=1;
02708 }
02709 }
02710 }
02711 }
02712 if(jnflag2[0]>0 || jnflag2[1]>0){
02713 for(vuw=0;vuw<2;vuw++){
02714 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(j));
02715 tmpclr[vuw].Add(shw);
02716 }
02717 }
02718 }
02719
02720 if( 1+tmpclr[0].GetLast()>0 && 1+tmpclr[1].GetLast()>0 ){
02721 id=-1; topscr=-9999.9;
02722 for(j=0;j<1+tmpclr[0].GetLast();j++){
02723 ShowerAtNu* shwu = (ShowerAtNu*)(tmpclr[0].At(j));
02724 ShowerAtNu* shwv = (ShowerAtNu*)(tmpclr[1].At(j));
02725 scr = shwu->GetEndPlane()+shwv->GetEndPlane();
02726 if(scr>topscr){
02727 topscr=scr; id=j;
02728 }
02729 }
02730 if(1+id>0){
02731 MSG("AlgAtNuReco", Msg::kVerbose) << " found existing end vertex shower " << endl;
02732 for(vuw=0;vuw<2;vuw++){
02733 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02734 ShowerAtNu* shw = (ShowerAtNu*)(tmpclr[vuw].At(id));
02735 trk->SetEndVtxShower(shw); shw->AddAssocTrack(trk); shw->SetEndVtxFlag(1);
02736 }
02737 }
02738 }
02739
02740
02741 else{
02742 if( shwflag==0 && i==0 ){
02743 for(vuw=0;vuw<2;vuw++){
02744 tmphitend[vuw].Clear();
02745 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02746
02747 mod=(Int_t)((trk->GetEndPlane())/(modpln+1));
02748 for(j=-3;j<4;j++){
02749 pln=trk->GetEndPlane()+2*j;
02750 if( 2*pln>trk->GetBegPlane()+trk->GetEndPlane()-1
02751 && pln>mod*(modpln+1) && pln<(mod+1)*(modpln+1) ){
02752 for(k=0;k<1+hitplnbnk[pln].GetLast();k++){
02753 HitAtNu* hit = (HitAtNu*)(hitplnbnk[pln].At(k));
02754 if( hit->GetShwFlag()<2
02755 && ( hit->GetTrkFlag()<2 || hit->GetCharge()>50 ) ){
02756 jnflag=0;
02757 for(l=0;l<1+tmphittrk[vuw].GetLast();l++){
02758 if(jnflag==0){
02759 HitAtNu* hittrk = (HitAtNu*)(tmphittrk[vuw].At(l));
02760 if( hittrk->GetPlane()>trk->GetEndPlane()-7
02761 && hittrk->IsDiffuseShwAssoc(hit)>0 ){
02762 tmphitend[vuw].Add(hit); hit->SetShwFlag(1); jnflag=1;
02763 }
02764 }
02765 }
02766 }
02767 }
02768 }
02769 }
02770
02771 if(1+tmphitend[vuw].GetLast()>0){
02772 bpln=trk->GetEndPlane()-6; epln=trk->GetEndPlane()+6;
02773 ctr=1;
02774 while(ctr>0){
02775 ctr=0;
02776 pln0=bpln-4; npln=5+(epln-bpln)/2;
02777 for(j=0;j<npln;j++){
02778 pln=pln0+2*j;
02779 if(pln>mod*(modpln+1) && pln<(mod+1)*(modpln+1)){
02780 for(k=0;k<1+hitplnbnk[pln].GetLast();k++){
02781 HitAtNu* hit = (HitAtNu*)(hitplnbnk[pln].At(k));
02782 if(hit->GetShwFlag()<1 && hit->GetTrkFlag()<1){
02783 jnflag=0;
02784 for(l=0;l<1+tmphitend[vuw].GetLast();l++){
02785 if(jnflag==0){
02786 HitAtNu* hitshw = (HitAtNu*)(tmphitend[vuw].At(l));
02787 if( hitshw->IsDiffuseShwAssoc(hit)>0 ){
02788 tmphitend[vuw].Add(hit); hit->SetShwFlag(1);
02789 if(hit->GetPlane()<bpln) bpln=hit->GetPlane();
02790 if(hit->GetPlane()>epln) epln=hit->GetPlane();
02791 jnflag=1; ctr++;
02792 }
02793 }
02794 }
02795 }
02796 }
02797 }
02798 }
02799 }
02800 for(j=0;j<1+tmphitend[vuw].GetLast();j++){
02801 HitAtNu* hit = (HitAtNu*)(tmphitend[vuw].At(j));
02802 hit->SetShwFlag(0);
02803 }
02804 }
02805
02806 }
02807 }
02808 }
02809
02810 if( 1+tmphitbeg[0].GetLast()>0 && 1+tmphitbeg[1].GetLast()>0 ){
02811 MSG("AlgAtNuReco", Msg::kVerbose) << " found new vertex shower at beginning of track " << endl;
02812 for(vuw=0;vuw<2;vuw++){
02813 ShowerAtNu* shw = new ShowerAtNu(slice); NEWOBJCTR++;
02814 for(j=0;j<1+tmphitbeg[vuw].GetLast();j++){
02815 HitAtNu* hit = (HitAtNu*)(tmphitbeg[vuw].At(j));
02816 shw->AddHit(hit); hit->SetShwFlag(2);
02817 }
02818 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02819 trk->SetBegVtxShower(shw); shw->AddAssocTrack(trk);
02820 shw->SetBegVtxFlag(1); shw->SetReseedFlag(1);
02821 shwbnk[vuw].Add(shw);
02822 }
02823 }
02824
02825 if( 1+tmphitend[0].GetLast()>0 && 1+tmphitend[1].GetLast()>0 ){
02826 MSG("AlgAtNuReco", Msg::kVerbose) << " found new vertex shower at end of track " << endl;
02827 for(vuw=0;vuw<2;vuw++){
02828 ShowerAtNu* shw = new ShowerAtNu(slice); NEWOBJCTR++;
02829 for(j=0;j<1+tmphitend[vuw].GetLast();j++){
02830 HitAtNu* hit = (HitAtNu*)(tmphitend[vuw].At(j));
02831 shw->AddHit(hit); hit->SetShwFlag(2);
02832 }
02833 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02834 trk->SetEndVtxShower(shw); shw->AddAssocTrack(trk);
02835 shw->SetEndVtxFlag(1); shw->SetReseedFlag(1);
02836 shwbnk[vuw].Add(shw);
02837 }
02838 }
02839
02840 }
02841 }
02842
02843
02844
02845
02846
02847
02848
02849 MSG("AlgAtNuReco", Msg::kDebug) << " *** mop up remaining hits *** " << endl;
02850 Double_t shwwinplns,shwwinstrps;
02851 if( 1+shwbnk[0].GetLast()>0 ){
02852 shwwinplns=8; shwwinstrps=12;
02853 if(1+trkbnk[0].GetLast()==0){
02854 shwwinplns=16; shwwinstrps=25; }
02855 for(mod=0;mod<modnum;mod++){
02856 for(pln=1;pln<1+modpln;pln++){
02857 i=pln+mod*(modpln+1);
02858 for(j=0;j<1+hitplnbnk[i].GetLast();j++){
02859 HitAtNu* hit = (HitAtNu*)(hitplnbnk[i].At(j));
02860 if(hit->GetShwFlag()<1 && hit->GetTrkFlag()<1){
02861 vuw = hit->GetPlaneView();
02862 id=-1; scr=shwwinplns+shwwinstrps;
02863 for(k=0;k<1+shwbnk[vuw].GetLast();k++){
02864 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(k));
02865 if( (Int_t)((shw->GetEndPlane())/(modpln+1)) == mod ){
02866 for(l=0;l<1+shw->GetHitLast();l++){
02867 if(1+shwbnk[vuw].GetLast()>1 || id<0){
02868 HitAtNu* hittmp = (HitAtNu*)(shw->GetHitAt(l));
02869 if( (hittmp->GetShwFlag()>1)
02870 && (hit->GetTrkFlag()<1||hittmp->IsShwAssoc(hit)>0) ){
02871 dpln=hit->GetPlane()-hittmp->GetPlane(); if(dpln<0) dpln=-dpln;
02872 dstr=hit->GetStrip()-hittmp->GetStrip(); if(dstr<0) dstr=-dstr;
02873 if( dpln<shwwinplns && dstr<shwwinstrps && dpln+dstr<scr ){
02874 id=k; scr=dpln+dstr;
02875 }
02876 }
02877 }
02878 }
02879 }
02880 }
02881 if(1+id>0){
02882 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(id));
02883 shw->AddHit(hit); hit->SetShwFlag(1);
02884 }
02885 }
02886 }
02887 }
02888 }
02889 }
02890
02891
02892 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF SHOWERS *** " << endl;
02893 for(vuw=0;vuw<2;vuw++){
02894 MSG("AlgAtNuReco", Msg::kVerbose) << "view: " << vuw << endl;
02895 for(i=0;i<1+shwbnk[vuw].GetLast();i++){
02896 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(i));
02897 MSG("AlgAtNuReco", Msg::kVerbose)
02898 << " i=" << i << " "
02899 << shw->GetBegPlane() << " " << shw->GetBegStrip() << " "
02900 << shw->GetEndPlane() << " " << shw->GetEndStrip() << endl;
02901 }
02902 }
02903
02904 gettimeofday(&tpafter,0);
02905 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02906 fTime_ShowerAtNu+=dtime;
02907
02908
02909
02910
02911
02912
02913 gettimeofday(&tpbefore,0);
02914 for(i=0;i<fEndPln;i++){
02915 hitplnbnk[i].Clear();
02916 xtalkplnbnk[i].Clear();
02917 clrplnbnk[i].Clear();
02918 }
02919 DELOBJCTR+=1+clrbnk.GetLast(); clrbnk.Delete();
02920 gettimeofday(&tpafter,0);
02921 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02922 fTime_TrackAtNu+=0.5*dtime;
02923 fTime_ShowerAtNu+=0.5*dtime;
02924
02925
02926 gettimeofday(&tpbefore,0);
02927 for(i=0;i<fEndPln;i++){
02928 segplnbnk[i].Clear();
02929 }
02930 DELOBJCTR+=1+trksegbnk.GetLast(); trksegbnk.Delete();
02931 gettimeofday(&tpafter,0);
02932 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02933 fTime_TrackAtNu+=dtime;
02934
02935
02936 gettimeofday(&tpbefore,0);
02937 DELOBJCTR+=1+shwsegbnk.GetLast(); shwsegbnk.Delete();
02938 gettimeofday(&tpafter,0);
02939 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02940 fTime_ShowerAtNu+=dtime;
02941
02942 }
02943
02944
02945
02946
02947
02948 MSG("AlgAtNuReco", Msg::kDebug) << " Formation of CandXXXHandles " << endl;
02949
02950
02951 AlgFactory &af = AlgFactory::GetInstance();
02952
02953
02954 gettimeofday(&tpbefore,0);
02955
02956 TObjArray* objtrk = 0;
02957 for(itr=0;itr<1+trkbnk[0].GetLast();itr++){
02958 TrackAtNu* trku = (TrackAtNu*)(trkbnk[0].At(itr));
02959 TrackAtNu* trkv = (TrackAtNu*)(trkbnk[1].At(itr));
02960 objtrk = new TObjArray(); NEWOBJCTR++;
02961 objtrk->Add(trku);
02962 objtrk->Add(trkv);
02963 trkcxt.Add(objtrk);
02964 }
02965
02966 MSG("AlgAtNuReco", Msg::kDebug) << "*** CREATE CandTrackAtNuListHandle *** " << endl;
02967 AlgHandle ah_trk = af.GetAlgHandle("AlgTrackAtNuList", "default");
02968 CandContext cx_trk(this, cx.GetMom());
02969 cx_trk.SetCandRecord(cx.GetCandRecord());
02970 cx_trk.SetDataIn(&trkcxt);
02971 CandTrackAtNuListHandle trk_list = CandTrackAtNuList::MakeCandidate(ah_trk, cx_trk);
02972 trk_list.SetName(fListOutTrk.Data());
02973 trk_list.SetTitle(TString("Created by AtNuFindModule from ").Append(slice_list->GetName()));
02974
02975 TIter trkitr(trk_list.GetDaughterIterator());
02976 while(CandTrackAtNuHandle* trk = dynamic_cast<CandTrackAtNuHandle*>(trkitr())){
02977 atnu_handle.AddTrack(trk);
02978 }
02979
02980 gettimeofday(&tpafter,0);
02981 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02982 fTime_TrackAtNu+=dtime;
02983
02984
02985 gettimeofday(&tpbefore,0);
02986
02987 TObjArray* objshw = 0;
02988 for(itr=0;itr<1+shwbnk[0].GetLast();itr++){
02989 ShowerAtNu* shwu = (ShowerAtNu*)(shwbnk[0].At(itr));
02990 ShowerAtNu* shwv = (ShowerAtNu*)(shwbnk[1].At(itr));
02991 objshw = new TObjArray(); NEWOBJCTR++;
02992 objshw->Add(shwu);
02993 objshw->Add(shwv);
02994 shwcxt.Add(objshw);
02995 }
02996
02997 MSG("AlgAtNuReco", Msg::kDebug) << "*** CREATE CandShowerAtNuListHandle *** " << endl;
02998 AlgHandle ah_shw = af.GetAlgHandle("AlgShowerAtNuList", "default");
02999 CandContext cx_shw(this, cx.GetMom());
03000 cx_shw.SetCandRecord(cx.GetCandRecord());
03001 cx_shw.SetDataIn(&shwcxt);
03002 CandShowerAtNuListHandle shw_list = CandShowerAtNuList::MakeCandidate(ah_shw, cx_shw);
03003 shw_list.SetName(fListOutShw.Data());
03004 shw_list.SetTitle(TString("Created by AtNuFindModule from ").Append(slice_list->GetName()));
03005
03006 TIter shwitr(shw_list.GetDaughterIterator());
03007 while(CandShowerAtNuHandle* shw = dynamic_cast<CandShowerAtNuHandle*>(shwitr())){
03008 atnu_handle.AddShower(shw);
03009 }
03010
03011 gettimeofday(&tpafter,0);
03012 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
03013 fTime_ShowerAtNu+=dtime;
03014
03015
03016
03017
03018
03019
03020 gettimeofday(&tpbefore,0);
03021 DELOBJCTR+=1+hitbnk.GetLast(); hitbnk.Delete();
03022 gettimeofday(&tpafter,0);
03023 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
03024 fTime_TrackAtNu+=0.5*dtime;
03025 fTime_ShowerAtNu+=0.5*dtime;
03026
03027
03028 gettimeofday(&tpbefore,0);
03029 for(itr=0;itr<2;itr++){
03030 DELOBJCTR+=1+trkbnk[itr].GetLast(); trkbnk[itr].Delete();
03031 }
03032 DELOBJCTR+=1+trkcxt.GetLast(); trkcxt.Delete();
03033 gettimeofday(&tpafter,0);
03034 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
03035 fTime_TrackAtNu+=dtime;
03036
03037
03038 gettimeofday(&tpbefore,0);
03039 for(itr=0;itr<2;itr++){
03040 DELOBJCTR+=1+shwbnk[itr].GetLast(); shwbnk[itr].Delete();
03041 }
03042 DELOBJCTR+=1+shwcxt.GetLast(); shwcxt.Delete();
03043 gettimeofday(&tpafter,0);
03044 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
03045 fTime_ShowerAtNu+=dtime;
03046
03047
03048
03049
03050
03051
03052 MSG("AlgAtNuReco", Msg::kDebug) << " *** CANDATNURECO PARAMETERS *** " << endl;
03053
03054
03055 Double_t atnurecoscore=0.0;
03056
03057
03058 atnu_handle.SetAtNuRecoScore(atnurecoscore);
03059
03060 gettimeofday(&atnuafter,0);
03061 dtime=1000.0*(atnuafter.tv_sec-atnubefore.tv_sec)+0.001*(atnuafter.tv_usec-atnubefore.tv_usec);
03062 fTime_AtNuReco+=dtime;
03063
03064
03065
03066
03067
03068
03069 MSG("AlgAtNuReco", Msg::kDebug) << " *** MEMORY BOOK-KEEPING : "
03070 << " new=" << NEWOBJCTR
03071 << " delete=" << DELOBJCTR << endl;
03072
03073 atnu_handle.SetNBlocks(NEWOBJCTR);
03074
03075 trk_list.SetCPUTime(fTime_TrackAtNu);
03076 MSG("AlgAtNuReco", Msg::kDebug) << " *** TRACK RECO TIME = " << fTime_TrackAtNu << endl;
03077
03078 shw_list.SetCPUTime(fTime_ShowerAtNu);
03079 MSG("AlgAtNuReco", Msg::kDebug) << " *** SHOWER RECO TIME = " << fTime_ShowerAtNu << endl;
03080
03081 atnu_handle.SetCPUTime(fTime_AtNuReco);
03082 MSG("AlgAtNuReco", Msg::kDebug) << " *** TOTAL ATNURECO TIME = " << fTime_AtNuReco << endl;
03083
03084
03085 CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
03086 candrec->SecureCandHandle(trk_list);
03087 candrec->SecureCandHandle(shw_list);
03088
03089 return;
03090
03091 }
03092
03093 void AlgAtNuReco::Trace(const char * ) const
03094 {
03095
03096 }
03097
03098
