#include <AlgShowerAtNu.h>
Inheritance diagram for AlgShowerAtNu:

Public Member Functions | |
| AlgShowerAtNu () | |
| ~AlgShowerAtNu () | |
| void | RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| void | Trace (const char *c) const |
Private Attributes | |
| TObjArray | fHitArr [500] |
|
|
Definition at line 33 of file AlgShowerAtNu.cxx. 00034 {
00035
00036 }
|
|
|
Definition at line 38 of file AlgShowerAtNu.cxx. 00039 {
00040
00041 }
|
|
||||||||||||||||
|
Implements AlgBase. Definition at line 43 of file AlgShowerAtNu.cxx. References CandHandle::AddDaughterLink(), UgliStripHandle::ClearFiber(), fHitArr, ShowerAtNu::GetAssocTrackAt(), ObjTrackAtNu::GetBegPlane(), ObjShowerAtNu::GetBegPlane(), ShowerAtNu::GetBegVtxFlag(), CandContext::GetCandRecord(), ShowerAtNu::GetCandSliceHandle(), HitAtNu::GetCandStripHandle(), CandStripHandle::GetCharge(), HitAtNu::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), CandShowerAtNuHandle::GetDirCosErrU(), CandShowerAtNuHandle::GetDirCosErrV(), CandRecoHandle::GetDirCosU(), CandRecoHandle::GetDirCosV(), CandRecoHandle::GetDirCosZ(), CandShowerAtNuHandle::GetDirTimeScore(), CandShowerAtNuHandle::GetDirVtxShwScore(), Registry::GetDouble(), PlexSEIdAltL::GetEnd(), ObjTrackAtNu::GetEndPlane(), ObjShowerAtNu::GetEndPlane(), ShowerAtNu::GetEndVtxFlag(), CandShowerHandle::GetEnergy(), UgliStripHandle::GetHalfLength(), ObjAtNu::GetHitAt(), ObjAtNu::GetHitLast(), CandShowerAtNuHandle::GetMaxPlane(), CandShowerAtNuHandle::GetMinPlane(), CandShowerAtNuHandle::GetNPlanes(), CandShowerAtNuHandle::GetNStrips(), HitAtNu::GetPlane(), CandStripHandle::GetPlaneView(), HitAtNu::GetPlaneView(), CandDigitHandle::GetPlexSEIdAltL(), CandShowerAtNuHandle::GetReseedFlag(), ObjShowerAtNu::GetReseedFlag(), HitAtNu::GetShwFlag(), UgliGeomHandle::GetSteelPlnHandle(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandStripHandle::GetTime(), CandRecoHandle::GetTimeOffset(), CandRecoHandle::GetTimeSlope(), HitAtNu::GetTPos(), RecMinos::GetVldContext(), CandRecoHandle::GetVtxPlane(), CandShowerAtNuHandle::GetVtxR(), CandShowerAtNuHandle::GetVtxShw(), CandRecoHandle::GetVtxT(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), UgliSteelPlnHandle::GetZ0(), HitAtNu::GetZPos(), UgliSteelPlnHandle::IsValid(), MSG, CandRecoHandle::SetCandSlice(), CandShowerAtNuHandle::SetDirCosErrU(), CandShowerAtNuHandle::SetDirCosErrV(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandShowerAtNuHandle::SetDirTimeScore(), CandShowerAtNuHandle::SetDirVtxShwScore(), CandShowerHandle::SetEnergy(), CandShowerAtNuHandle::SetMaxPlane(), CandShowerAtNuHandle::SetMinPlane(), CandShowerAtNuHandle::SetNPlanes(), CandShowerAtNuHandle::SetNStrips(), CandShowerAtNuHandle::SetReseedFlag(), CandRecoHandle::SetTimeOffset(), CandRecoHandle::SetTimeSlope(), CandShowerHandle::SetU(), CandShowerHandle::SetV(), CandRecoHandle::SetVtxPlane(), CandShowerAtNuHandle::SetVtxR(), CandShowerAtNuHandle::SetVtxShw(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), and UgliStripHandle::WlsPigtail(). 00044 {
00045
00046 MSG("AlgShowerAtNu", Msg::kDebug) << "AlgShowerAtNu::RunAlg(...)" << endl;
00047
00048 // Get CandShowerAtNuHandle
00049 CandShowerAtNuHandle& shower = dynamic_cast<CandShowerAtNuHandle&>(ch);
00050
00051 // Unpack AlgConfig
00052 Double_t fFibreIndex;
00053 fFibreIndex = ac.GetDouble("FibreIndex");
00054 MSG("AlgShowerAtNu", Msg::kDebug) << " FibreIndex=" << fFibreIndex << endl;
00055
00056 // Unpack CandContext
00057 const TObjArray *arr = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00058 ShowerAtNu* shwu = (ShowerAtNu*)(arr->At(0));
00059 ShowerAtNu* shwv = (ShowerAtNu*)(arr->At(1));
00060 shower.SetCandSlice(shwu->GetCandSliceHandle());
00061
00062 CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
00063 VldContext *vldc = (VldContext*)(candrec->GetVldContext());
00064 UgliGeomHandle ugh(*vldc);
00065
00066 // Get Geometry
00067 Int_t atmosflag,modtype,modnum,modpln,modstr;
00068 switch(candrec->GetVldContext()->GetDetector()){
00069 case DetectorType::kFar:
00070 atmosflag=1; modtype=1; modnum=2; modpln=248; modstr=192; break;
00071 case DetectorType::kNear:
00072 atmosflag=0; modtype=2; modnum=1; modpln=282; modstr=96; break;
00073 case DetectorType::kCalib:
00074 atmosflag=0; modtype=3; modnum=1; modpln=60; modstr=24; break;
00075 default:
00076 atmosflag=0; modtype=-1; modnum=0; modpln=0; modstr=0; break;
00077 }
00078
00079 Int_t i,j,k;
00080 Int_t vuw;
00081 Int_t pln,bpln,epln,npln;
00082 Double_t begZsm1=-9999.9,endZsm1=-9999.9;
00083 Double_t begZsm2=-9999.9,endZsm2=-9999.9;
00084 Int_t begsm1=-1,endsm1=-1,begsm2=-1,endsm2=-1;
00085
00086 for(j=0;j<486;j++){
00087 PlexPlaneId myplaneid(vldc->GetDetector(),j,1);
00088 UgliSteelPlnHandle myplanehandle = ugh.GetSteelPlnHandle(myplaneid);
00089 if(myplanehandle.IsValid()){
00090 if(j<=modpln && begZsm1<0.0){ begsm1=j; begZsm1=myplanehandle.GetZ0(); }
00091 if(j<=modpln && myplanehandle.GetZ0()>endZsm1){ endsm1=j; endZsm1=myplanehandle.GetZ0(); }
00092 if(j>modpln && begZsm2<0.0){ begsm2=j; begZsm2=myplanehandle.GetZ0(); }
00093 if(j>modpln && myplanehandle.GetZ0()>endZsm2){ endsm2=j; endZsm2=myplanehandle.GetZ0(); }
00094 }
00095 }
00096
00097 MSG("AlgShowerAtNu", Msg::kDebug) << " *** detector geometry *** " << endl
00098 << " SM1 : pln " << begsm1 << " -> " << endsm1
00099 << " Z " << begZsm1 << " -> " << endZsm1 << endl
00100 << " SM2 : pln " << begsm2 << " -> " << endsm2
00101 << " Z " << begZsm2 << " -> " << endZsm2 << endl;
00102
00103
00104 if(shwu->GetBegPlane()<shwv->GetBegPlane()) bpln=shwu->GetBegPlane(); else bpln=shwv->GetBegPlane();
00105 if(shwu->GetEndPlane()>shwv->GetEndPlane()) epln=shwu->GetEndPlane(); else epln=shwv->GetEndPlane();
00106
00107 // ***********************************************
00108 // * C A L C U L A T E C O - O R D I N A T E S *
00109 // ***********************************************
00110
00111 MSG("AlgShowerAtNu", Msg::kDebug) << " *** calculate co-ordinates *** " << endl;
00112
00113 npln = 1+epln-bpln;
00114 Double_t* Z = new Double_t[npln];
00115 Double_t* U = new Double_t[npln];
00116 Double_t* V = new Double_t[npln];
00117 Double_t* CT = new Double_t[npln];
00118 Double_t* Q = new Double_t[npln];
00119 Int_t* plnvuw = new Int_t[npln];
00120 Int_t* plnnum = new Int_t[npln];
00121 for(i=0;i<npln;i++){
00122 plnvuw[i]=-1; plnnum[i]=0;
00123 U[i]=0.0; V[i]=0.0; Z[i]=0.0;
00124 Q[i]=0.0; CT[i]=0.0;
00125 }
00126
00127 // loop over hits and fill arrays
00128 for(i=0;i<1+shwu->GetHitLast();i++){
00129 HitAtNu* hit = (HitAtNu*)(shwu->GetHitAt(i));
00130 CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
00131 pln = hit->GetPlane()-bpln;
00132 fHitArr[pln].Add(hit);
00133 if(plnvuw[pln]<0){
00134 plnvuw[pln]=hit->GetPlaneView();
00135 Z[pln]=hit->GetZPos();
00136 }
00137 U[pln]+=hit->GetTPos()*hit->GetCharge();
00138 Q[pln]+=hit->GetCharge();
00139 shower.AddDaughterLink(*strip);
00140 plnnum[pln]++;
00141 }
00142
00143 for(i=0;i<1+shwv->GetHitLast();i++){
00144 HitAtNu* hit = (HitAtNu*)(shwv->GetHitAt(i));
00145 CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
00146 pln = hit->GetPlane()-bpln;
00147 fHitArr[pln].Add(hit);
00148 if(plnvuw[pln]<0){
00149 plnvuw[pln]=hit->GetPlaneView();
00150 Z[pln]=hit->GetZPos();
00151 }
00152 V[pln]+=hit->GetTPos()*hit->GetCharge();
00153 Q[pln]+=hit->GetCharge();
00154 shower.AddDaughterLink(*strip);
00155 plnnum[pln]++;
00156 }
00157
00158 Int_t nplanes=0,nstrips=0;
00159 for(i=0;i<npln;i++){
00160 if(plnvuw[i]>-1){
00161 if(Q[i]>0.0){
00162 U[i]=U[i]/Q[i];
00163 V[i]=V[i]/Q[i];
00164 nstrips+=plnnum[i];
00165 nplanes++;
00166 }
00167 }
00168 }
00169
00170 // estimate transverse co-ordinates
00171 Int_t km,kp;
00172 Double_t q,sq,sqt;
00173 for(pln=0;pln<npln;pln++){
00174 if(plnvuw[pln]>-1){
00175 vuw=plnvuw[pln];
00176 sq=0.0; sqt=0.0;
00177 k=0; km=-1;
00178 while(pln-k>0 && km<0){
00179 k++; if(plnvuw[pln-k]>-1 && plnvuw[pln-k]!=vuw) km=k;
00180 }
00181 k=0; kp=-1;
00182 while(pln+k<npln-1 && kp<0){
00183 k++; if(plnvuw[pln+k]>-1 && plnvuw[pln+k]!=vuw) kp=k;
00184 }
00185
00186 if(km>0){
00187 if(vuw==0){ q=1.0/km; sq+=q*Q[pln-km]; sqt+=q*Q[pln-km]*V[pln-km]; }
00188 if(vuw==1){ q=1.0/km; sq+=q*Q[pln-km]; sqt+=q*Q[pln-km]*U[pln-km]; }
00189 }
00190 if(kp>0){
00191 if(vuw==0){ q=1.0/kp; sq+=q*Q[pln+kp]; sqt+=q*Q[pln+kp]*V[pln+kp]; }
00192 if(vuw==1){ q=1.0/kp; sq+=q*Q[pln+kp]; sqt+=q*Q[pln+kp]*U[pln+kp]; }
00193 }
00194
00195 if(sq>0.0){
00196 if(vuw==0){ V[pln]=sqt/sq; } if(vuw==1){ U[pln]=sqt/sq; }
00197 }
00198 else{
00199 if(vuw==0){ V[pln]=0.0; } if(vuw==1){ U[pln]=0.0; }
00200 }
00201 }
00202 }
00203
00204 // get timing information
00205 Double_t tpos,fibre,corr;
00206 Double_t tot_evtchg=0.0,tot_evtchg_fibrecorr=0.0;
00207 Double_t digchg,tot_digchg;
00208 Double_t ctime,tot_ctime;
00209 Double_t indx=fFibreIndex;
00210 for(pln=0;pln<npln;pln++){
00211 tot_ctime=0.0; tot_digchg=0.0;
00212 for(i=0;i<1+fHitArr[pln].GetLast();i++){
00213 HitAtNu* hit = (HitAtNu*)(fHitArr[pln].At(i));
00214
00215 // timing/charge
00216 CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
00217 if(strip->GetPlaneView()==PlaneView::kU) tpos = -V[pln];
00218 if(strip->GetPlaneView()==PlaneView::kV) tpos = U[pln];
00219 PlexStripEndId stripid = strip->GetStripEndId();
00220 UgliStripHandle striphandle = ugh.GetStripHandle(stripid);
00221 TIter digitr(strip->GetDaughterIterator());
00222 while(CandDigitHandle* digit = dynamic_cast<CandDigitHandle*>(digitr())){
00223 ctime=0.0; fibre=0.0; digchg=0.0;
00224 if(digit->GetPlexSEIdAltL().GetEnd()==StripEnd::kPositive){ fibre = striphandle.ClearFiber(StripEnd::kPositive)+striphandle.WlsPigtail(StripEnd::kPositive)+striphandle.GetHalfLength()-tpos; ctime = 3.0e8*strip->GetTime(StripEnd::kPositive); digchg=strip->GetCharge(StripEnd::kPositive); }
00225 if(digit->GetPlexSEIdAltL().GetEnd()==StripEnd::kNegative){ fibre = striphandle.ClearFiber(StripEnd::kNegative)+striphandle.WlsPigtail(StripEnd::kNegative)+striphandle.GetHalfLength()+tpos; ctime = 3.0e8*strip->GetTime(StripEnd::kNegative); digchg=strip->GetCharge(StripEnd::kNegative); }
00226 ctime = ctime-indx*fibre; corr = exp(-(8.0-fibre)/7.05);
00227 tot_ctime += digchg*ctime; tot_digchg += digchg;
00228 tot_evtchg += digchg; tot_evtchg_fibrecorr += digchg*corr;
00229 }
00230 }
00231 if(tot_digchg>0.0){
00232 CT[pln] = tot_ctime/tot_digchg;
00233 }
00234 }
00235
00236 // set map values
00237 for(i=0;i<npln;i++){
00238 pln=bpln+i;
00239 if(plnvuw[i]>-1){
00240 shower.SetU(pln,U[i]);
00241 shower.SetV(pln,V[i]);
00242 }
00243 }
00244
00245 MSG("AlgShowerAtNu", Msg::kDebug) << " *** shower hits *** " << endl;
00246 for(pln=0;pln<npln;pln++){
00247 if( plnvuw[pln]>-1 ){
00248 MSG("AlgShowerAtNu",Msg::kDebug)
00249 << " pln=" << pln << " vuw=" << plnvuw[pln]
00250 << " z=" << Z[pln] << " u=" << U[pln] << " v=" << V[pln] << endl;
00251 }
00252 }
00253
00254
00255 // ***********************************
00256 // * D E T E R M I N E V E R T E X *
00257 // ***********************************
00258
00259 MSG("AlgShowerAtNu", Msg::kDebug) << " *** determine vertex *** " << endl;
00260
00261 Bool_t vtxshw=0;
00262 Double_t trkscr=0.0,shwscr=0.0,dirscr=0.0;
00263 Double_t vtxu=0.0,vtxv=0.0,vtxx=0.0,vtxy=0.0,vtxz=0.0;
00264 Int_t vtxpln=-1,trkbpln=-1,trkepln=-1,shwbpln=-1,shwepln=-1;
00265
00266 // determine vertex plane
00267 Double_t totchgpln=0.0,totchg=0.0;
00268 for(i=0;i<npln;i++){
00269 for(j=0;j<1+fHitArr[i].GetLast();j++){
00270 HitAtNu* hit = (HitAtNu*)(fHitArr[i].At(j));
00271 if(hit->GetShwFlag()>1){
00272 totchg+=hit->GetCharge();
00273 totchgpln+=hit->GetCharge()*hit->GetPlane();
00274 if(shwbpln<0||hit->GetPlane()<shwbpln) shwbpln=hit->GetPlane();
00275 if(shwepln<0||hit->GetPlane()>shwepln) shwepln=hit->GetPlane();
00276 }
00277 }
00278 }
00279 if(totchg>0.0){
00280 vtxpln=(Int_t)(totchgpln/totchg);
00281 }
00282 else{
00283 vtxpln=(Int_t)(0.5*(bpln+epln));
00284 shwbpln=bpln; shwepln=epln;
00285 }
00286
00287 TrackAtNu* vtxtrkU = (TrackAtNu*)(shwu->GetAssocTrackAt(0));
00288 TrackAtNu* vtxtrkV = (TrackAtNu*)(shwv->GetAssocTrackAt(0));
00289 if(vtxtrkU && vtxtrkV){
00290 vtxshw=1;
00291 if(vtxtrkU->GetBegPlane()<vtxtrkV->GetBegPlane()) trkbpln=vtxtrkU->GetBegPlane(); else trkbpln=vtxtrkV->GetBegPlane();
00292 if(vtxtrkU->GetEndPlane()>vtxtrkV->GetEndPlane()) trkepln=vtxtrkU->GetEndPlane(); else trkepln=vtxtrkV->GetEndPlane();
00293 trkscr = (trkbpln+trkepln-2.0*vtxpln)/(1.0+trkepln-trkbpln);
00294 if(trkscr>0.0){ trkscr = (shwbpln+shwepln-2.0*trkbpln)/(1.0+shwepln-shwbpln); }
00295 if(trkscr<0.0){ trkscr = (shwbpln+shwepln-2.0*trkepln)/(1.0+shwepln-shwbpln); }
00296 if( ( shwu->GetBegVtxFlag() && shwv->GetBegVtxFlag())
00297 && !( shwu->GetEndVtxFlag() && shwv->GetEndVtxFlag()) ) vtxpln=trkbpln;
00298 if( ( shwu->GetEndVtxFlag() && shwv->GetEndVtxFlag())
00299 && !( shwu->GetBegVtxFlag() && shwv->GetBegVtxFlag()) ) vtxpln=trkepln;
00300 }
00301
00302 // determine vertex position
00303 pln=vtxpln-bpln;
00304 if( pln>-1 && pln<npln
00305 && plnvuw[pln]>-1 ){
00306 vtxu=U[pln]; vtxv=V[pln]; vtxz=Z[pln];
00307 }
00308 else{
00309 k=0; km=-1;
00310 while(pln-k>0 && km<0){
00311 k++; if(pln-k<npln && plnvuw[pln-k]>-1) km=k;
00312 }
00313 k=0; kp=-1;
00314 while(pln+k<npln-1 && kp<0){
00315 k++; if(pln+k>-1 && plnvuw[pln+k]>-1) kp=k;
00316 }
00317 if(km>0 && kp>0){
00318 vtxu=(km*U[pln+kp]+kp*U[pln-km])/(km+kp);
00319 vtxv=(km*V[pln+kp]+kp*V[pln-km])/(km+kp);
00320 vtxz=(km*Z[pln+kp]+kp*Z[pln-km])/(km+kp);
00321 }
00322 else{
00323 if(km>0){ vtxu=U[pln-km]; vtxv=V[pln-km]; vtxz=Z[pln-km]; }
00324 if(kp>0){ vtxu=U[pln+kp]; vtxv=V[pln+kp]; vtxz=Z[pln+kp]; }
00325 }
00326 }
00327
00328 vtxx=0.7071*(vtxu-vtxv); vtxy=0.7071*(vtxu+vtxv);
00329
00330 MSG("AlgShowerAtNu", Msg::kVerbose)
00331 << " ... vertex: "
00332 << " plane = " << vtxpln
00333 << " position = (" << vtxu << "," << vtxv << "," << vtxz << ")" << endl;
00334
00335
00336 // *************************************
00337 // * C O N T A I N M E N T S T U F F *
00338 // *************************************
00339
00340 MSG("AlgShowerAtNu", Msg::kDebug) << " *** containment stuff *** " << endl;
00341
00342 Double_t vtxr=0.0;
00343
00344 if(atmosflag){
00345
00346 // closest distance to edge
00347 Double_t xm,xp,ym,yp,um,up,vm,vp;
00348 Double_t rmin;
00349
00350 rmin=4.0;
00351 up=4.0-vtxu; if(up<rmin) rmin=up;
00352 um=4.0+vtxu; if(um<rmin) rmin=um;
00353 vp=4.0-vtxv; if(vp<rmin) rmin=vp;
00354 vm=4.0+vtxv; if(vm<rmin) rmin=vm;
00355 xp=4.0-vtxx; if(xp<rmin) rmin=xp;
00356 xm=4.0+vtxx; if(xm<rmin) rmin=xm;
00357 yp=4.0-vtxy; if(yp<rmin) rmin=yp;
00358 ym=4.0+vtxy; if(ym<rmin) rmin=ym;
00359 vtxr=rmin;
00360
00361 }
00362
00363 //*****************************************
00364 //* D E T E R M I N E D I R E C T I O N *
00365 //*****************************************
00366
00367 MSG("AlgShowerAtNu", Msg::kDebug) << " *** determine direction *** " << endl;
00368
00369 // measure overall timeslope
00370 Double_t ctimeslope=0.0,ctimeoffset=0.0,ctimeaverage=0.0,ctimescatter=0.0;
00371 Double_t sw,swz,swt,swzz,swzt,swtt;
00372 Double_t du,dv,dz;
00373 Double_t z,t,w;
00374 Double_t diru,dirv,dirz,dirs;
00375 Double_t direrru,direrrv;
00376 Double_t offset;
00377 Int_t flag,npts;
00378
00379 npts=0;
00380 sw=0.0; swz=0.0; swzz=0.0; swt=0.0; swtt=0.0; swzt=0.0;
00381 for(i=0;i<npln;i++){
00382 if( plnvuw[i]>-1 ){
00383 z=Z[i]; t=CT[i]; w=Q[i];
00384 swz+=w*z; swzz+=w*z*z; swt+=w*t; swtt+=w*t*t; swzt+=w*z*t;
00385 sw+=w; npts++;
00386 }
00387 }
00388 if(npts>1){
00389 ctimeslope=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00390 ctimeoffset=(swt*swzz-swz*swzt)/(sw*swzz-swz*swz);
00391 ctimeaverage=(swt/sw);
00392 if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 ){
00393 ctimescatter=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00394 }
00395 }
00396 offset=ctimeaverage;
00397 shwscr=ctimescatter;
00398
00399 if(vtxshw) dirscr=trkscr; else dirscr=shwscr;
00400
00401 npts=0;
00402 diru=0.0; direrru=0.0;
00403 if(1+shwu->GetEndPlane()-shwu->GetBegPlane()>1){
00404 sw=0.0; swz=0.0; swt=0.0;
00405 swzz=0.0; swzt=0.0; swtt=0.0;
00406 for(i=0;i<npln;i++){
00407 flag=0;
00408 for(j=0;j<1+fHitArr[i].GetLast();j++){
00409 HitAtNu* hit = (HitAtNu*)(fHitArr[i].At(j));
00410 if(hit->GetShwFlag()>1 && hit->GetPlaneView()==0){
00411 du=hit->GetTPos()-vtxu; dz=hit->GetZPos()-vtxz; w=hit->GetCharge();
00412 swz+=w*dz; swt+=w*du;
00413 swzz+=w*dz*dz; swzt+=w*dz*du; swtt+=w*du*du;
00414 sw+=w; flag=1;
00415 }
00416 }
00417 if(flag) npts++;
00418 }
00419 if(npts>1){
00420 diru=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00421 if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 ){
00422 direrru=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00423 }
00424 }
00425 }
00426
00427 npts=0;
00428 dirv=0.0; direrrv=0.0;
00429 if(1+shwv->GetEndPlane()-shwv->GetBegPlane()>1){
00430 sw=0.0; swz=0.0; swt=0.0;
00431 swzz=0.0; swzt=0.0; swtt=0.0;
00432 for(i=0;i<npln;i++){
00433 flag=0;
00434 for(j=0;j<1+fHitArr[i].GetLast();j++){
00435 HitAtNu* hit = (HitAtNu*)(fHitArr[i].At(j));
00436 if(hit->GetShwFlag()>1 && hit->GetPlaneView()==1){
00437 dv=hit->GetTPos()-vtxv; dz=hit->GetZPos()-vtxz; w=hit->GetCharge();
00438 swz+=w*dz; swt+=w*dv;
00439 swzz+=w*dz*dz; swzt+=w*dz*dv; swtt+=w*dv*dv;
00440 sw+=w; flag=1;
00441 }
00442 }
00443 if(flag) npts++;
00444 }
00445 if(npts>1){
00446 dirv=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00447 if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 ){
00448 direrrv=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00449 }
00450 }
00451 }
00452
00453 dirs=sqrt(diru*diru+dirv*dirv+1.0);
00454 diru=diru/dirs; dirv=dirv/dirs; dirz=1.0/dirs;
00455
00456 MSG("AlgShowerAtNu", Msg::kVerbose)
00457 << " ... direction: "
00458 << " (" << diru << "," << dirv << "," << dirz << ")" << endl;
00459
00460
00461 // ***********************************
00462 // * C A L C U L A T E E N E R G Y *
00463 // ***********************************
00464
00465 MSG("AlgShowerAtNu", Msg::kDebug) << " *** calculate energy *** " << endl;
00466
00467 // calculate energy from pulse height
00468 Double_t energy;
00469 Double_t tot_ph;
00470 if(tot_evtchg>0.0){
00471 tot_ph = tot_evtchg;
00472 energy = tot_ph; // NEEDS TUNING
00473 energy = 0.3 + tot_ph/150.0;
00474 }
00475
00476 // *******************************************
00477 // * S E T S H O W E R V A R I A B L E S *
00478 // *******************************************
00479
00480 // set reseed flag
00481 Bool_t reseedflag=0;
00482 if(shwu->GetReseedFlag()||shwv->GetReseedFlag()) reseedflag=1;
00483
00484 // determine overall direction
00485 Double_t dir=1.0;
00486 if(dirscr>0.0) dir=1.0; if(dirscr<0.0) dir=-1.0;
00487
00488 // set shower variables
00489 shower.SetMinPlane(bpln);
00490 shower.SetMaxPlane(epln);
00491 shower.SetNPlanes(nplanes);
00492 shower.SetNStrips(nstrips);
00493 shower.SetTimeSlope((dir)/(3.0e8));
00494 shower.SetTimeOffset((offset)/(3.0e8));
00495 shower.SetDirTimeScore(shwscr);
00496 shower.SetDirVtxShwScore(dirscr);
00497 shower.SetVtxShw(vtxshw);
00498 shower.SetVtxU(vtxu);
00499 shower.SetVtxV(vtxv);
00500 shower.SetVtxZ(vtxz);
00501 shower.SetVtxR(vtxr);
00502 shower.SetVtxT((offset)/(3.0e8));
00503 shower.SetVtxPlane(vtxpln);
00504 shower.SetDirCosU(dir*diru);
00505 shower.SetDirCosV(dir*dirv);
00506 shower.SetDirCosZ(dir*dirz);
00507 shower.SetDirCosErrU(dir*direrru);
00508 shower.SetDirCosErrV(dir*direrrv);
00509 shower.SetEnergy(energy);
00510 shower.SetReseedFlag(reseedflag);
00511
00512 MSG("AlgShowerAtNu", Msg::kDebug)
00513 << " *** CREATING NEW SHOWER *** " << endl
00514 << " MinPlane = " << shower.GetMinPlane() << endl
00515 << " MaxPlane = " << shower.GetMaxPlane() << endl
00516 << " NPlanes = " << shower.GetNPlanes() << endl
00517 << " NStrips = " << shower.GetNStrips() << endl
00518 << " TimeSlope = " << shower.GetTimeSlope() << endl
00519 << " TimeOffset = " << shower.GetTimeOffset() << endl
00520 << " DirTimeScore = " << shower.GetDirTimeScore() << endl
00521 << " DirVtxShwScore= " << shower.GetDirVtxShwScore() << endl
00522 << " VtxShw = " << shower.GetVtxShw() << endl
00523 << " VtxU = " << shower.GetVtxU() << endl
00524 << " VtxV = " << shower.GetVtxV() << endl
00525 << " VtxZ = " << shower.GetVtxZ() << endl
00526 << " VtxR = " << shower.GetVtxR() << endl
00527 << " VtxT = " << shower.GetVtxT() << endl
00528 << " VtxPln = " << shower.GetVtxPlane() << endl
00529 << " DirCosU = " << shower.GetDirCosU() << endl
00530 << " DirCosV = " << shower.GetDirCosV() << endl
00531 << " DirCosZ = " << shower.GetDirCosZ() << endl
00532 << " DirCosErrU = " << shower.GetDirCosErrU() << endl
00533 << " DirCosErrV = " << shower.GetDirCosErrV() << endl
00534 << " Energy = " << shower.GetEnergy() << endl
00535 << " ReseedFlag = " << shower.GetReseedFlag() << endl;
00536
00537
00538 for(i=0;i<npln;i++){
00539 fHitArr[i].Clear();
00540 }
00541
00542 delete [] Z;
00543 delete [] U;
00544 delete [] V;
00545 delete [] CT;
00546 delete [] Q;
00547 delete [] plnvuw;
00548 delete [] plnnum;
00549
00550 return;
00551 }
|
|
|
Reimplemented from AlgBase. Definition at line 553 of file AlgShowerAtNu.cxx. 00554 {
00555
00556 }
|
|
|
Definition at line 18 of file AlgShowerAtNu.h. Referenced by RunAlg(). |
1.3.9.1