AlgShowerAtNu Class Reference

#include <AlgShowerAtNu.h>

Inheritance diagram for AlgShowerAtNu:
AlgBase

List of all members.

Public Member Functions

 AlgShowerAtNu ()
 ~AlgShowerAtNu ()
void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
void Trace (const char *c) const

Private Attributes

TObjArray fHitArr [500]

Detailed Description

Definition at line 8 of file AlgShowerAtNu.h.


Constructor & Destructor Documentation

AlgShowerAtNu::AlgShowerAtNu (  ) 

Definition at line 34 of file AlgShowerAtNu.cxx.

00035 {
00036 
00037 }

AlgShowerAtNu::~AlgShowerAtNu (  ) 

Definition at line 39 of file AlgShowerAtNu.cxx.

00040 {
00041 
00042 }


Member Function Documentation

void AlgShowerAtNu::RunAlg ( AlgConfig ac,
CandHandle ch,
CandContext cx 
) [virtual]

Implements AlgBase.

Definition at line 44 of file AlgShowerAtNu.cxx.

References CandHandle::AddDaughterLink(), CandRecoHandle::CalibrateMIP(), CandRecoHandle::CalibrateSigMapped(), UgliStripHandle::ClearFiber(), fHitArr, ShowerAtNu::GetAssocTrackAt(), Calibrator::GetAttenCorrectedTpos(), ObjShowerAtNu::GetBegPlane(), ObjTrackAtNu::GetBegPlane(), ShowerAtNu::GetBegVtxFlag(), CandContext::GetCandRecord(), ShowerAtNu::GetCandSliceHandle(), HitAtNu::GetCandStripHandle(), HitAtNu::GetCharge(), CandStripHandle::GetCharge(), CandContext::GetDataIn(), VldContext::GetDetector(), CandShowerAtNuHandle::GetDirCosErrU(), CandShowerAtNuHandle::GetDirCosErrV(), CandRecoHandle::GetDirCosU(), CandRecoHandle::GetDirCosV(), CandRecoHandle::GetDirCosZ(), CandShowerAtNuHandle::GetDirTimeScore(), CandShowerAtNuHandle::GetDirVtxShwScore(), Registry::GetDouble(), PlexStripEndId::GetEncoded(), ObjShowerAtNu::GetEndPlane(), ObjTrackAtNu::GetEndPlane(), ShowerAtNu::GetEndVtxFlag(), CandShowerHandle::GetEnergy(), UgliStripHandle::GetHalfLength(), ObjAtNu::GetHitAt(), ObjAtNu::GetHitLast(), Registry::GetInt(), CandShowerAtNuHandle::GetMaxPlane(), CandShowerAtNuHandle::GetMinPlane(), Calibrator::GetMIP(), CandStripHandle::GetNDigit(), CandShowerAtNuHandle::GetNPlanes(), CandShowerAtNuHandle::GetNStrips(), VHS::GetPlane(), HitAtNu::GetPlane(), HitAtNu::GetPlaneView(), CandStripHandle::GetPlaneView(), ObjShowerAtNu::GetReseedFlag(), CandShowerAtNuHandle::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(), Calibrator::Instance(), Detector::kCalib, Msg::kDebug, Detector::kFar, Munits::km, Detector::kNear, StripEnd::kNegative, StripEnd::kPositive, CalDigitType::kSigCorr, PlaneView::kU, PlaneView::kV, Msg::kVerbose, PlaneView::kX, PlaneView::kY, MSG, Calibrator::ReInitialise(), CandRecoHandle::SetCandSlice(), CandShowerAtNuHandle::SetDirCosErrU(), CandShowerAtNuHandle::SetDirCosErrV(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandShowerAtNuHandle::SetDirTimeScore(), CandShowerAtNuHandle::SetDirVtxShwScore(), CandShowerHandle::SetEnergy(), ValueErr< T >::SetError(), 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(), UgliStripHandle::WlsPigtail(), and MinosMaterial::Z().

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

void AlgShowerAtNu::Trace ( const char *  c  )  const [virtual]

Reimplemented from AlgBase.

Definition at line 617 of file AlgShowerAtNu.cxx.

00618 {
00619 
00620 }


Member Data Documentation

TObjArray AlgShowerAtNu::fHitArr[500] [private]

Definition at line 18 of file AlgShowerAtNu.h.

Referenced by RunAlg().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1