loon: AlgShowerCam.cxx Source File

00001 
00002 // Package: AtNuReco
00003 //
00004 // AlgShowerCam.cxx
00005 //
00006 // chapman@hep.phy.cam.ac.uk
00007 // blake@hep.phy.cam.ac.uk
00009 #include "Algorithm/AlgConfig.h"
00010 #include "Algorithm/AlgFactory.h"
00011 #include "Algorithm/AlgHandle.h"
00012 
00013 #include "MessageService/MsgService.h" 
00014 #include "JobControl/JobCModuleRegistry.h"
00015 #include "Calibrator/Calibrator.h"
00016 #include "RecoBase/CandStripHandle.h"
00017 #include "CandDigit/CandDigitHandle.h"
00018 #include "CandDigit/CandDigitListHandle.h"
00019 #include "Candidate/CandContext.h"
00020 
00021 #include "UgliGeometry/UgliGeomHandle.h"
00022 #include "Validity/VldContext.h"
00023 #include "Plex/PlexPlaneId.h"
00024 
00025 #include "AtNuReco/AlgShowerCam.h"
00026 #include "AtNuReco/CandShowerAtNuHandle.h"
00027 #include "AtNuReco/HitCamAtNu.h"
00028 #include "AtNuReco/ShowerCamAtNu.h"
00029 #include "AtNuReco/TrackCamAtNu.h"
00030 
00031 //
00032 // $Log $
00033 //
00034 
00035 ClassImp(AlgShowerCam)
00036 
00037 CVSID("$Id: AlgShowerCam.cxx,v 1.5 2009/01/25 17:00:03 blake Exp $");
00039 AlgShowerCam::AlgShowerCam()
00040 {
00041 }
00043 AlgShowerCam::~AlgShowerCam()
00044 {
00045 }
00047 void AlgShowerCam::RunAlg(AlgConfig& ac, CandHandle &ch, CandContext &cx)
00048 {
00049 
00050   MSG("AlgShowerCam", Msg::kDebug) << "AlgShowerCam::RunAlg(...)" << endl;
00051 
00052   // Get CandShowerAtNuHandle
00053   CandShowerAtNuHandle& shower = dynamic_cast<CandShowerAtNuHandle&>(ch);
00054 
00056   // Unpack AlgConfig
00057   double fFibreIndex=1.77;
00058   int fCamAnaOnOff=0;
00059   if( ac.KeyExists("FibreIndex") )  fFibreIndex  = ac.GetDouble("FibreIndex");
00060   if( ac.KeyExists("CamAnaOnOff") ) fCamAnaOnOff = ac.GetInt("CamAnaOnOff");
00061   MSG("AlgShowerCam", Msg::kDebug) << " FibreIndex=" << fFibreIndex << endl;
00062   MSG("AlgShowerCam", Msg::kDebug) << " CamAnaOnOff=" << fCamAnaOnOff << endl;
00063 
00065   // Unpack CandContext
00066   const TObjArray *arr = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00067   ShowerCamAtNu* shwu = (ShowerCamAtNu*)(arr->At(0));
00068   ShowerCamAtNu* shwv = (ShowerCamAtNu*)(arr->At(1)); 
00069   shower.SetCandSlice(shwu->GetCandSliceHandle());
00070 
00071   CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
00072   VldContext *vldc = (VldContext*)(candrec->GetVldContext());
00073   UgliGeomHandle ugh(*vldc);
00074 
00076   // Reset Calibrator
00077   Calibrator& cal = Calibrator::Instance();
00078   cal.ReInitialise(*vldc);
00079 
00081   // Get Geometry
00082   int atmosflag,modtype,modnum,modpln,modstr;
00083   switch(candrec->GetVldContext()->GetDetector()){
00084     case Detector::kFar:
00085       atmosflag=1; modtype=1; modnum=2; modpln=248; modstr=192; break;
00086     case Detector::kNear:
00087       atmosflag=0; modtype=2; modnum=1; modpln=282; modstr=96; break;
00088     case Detector::kCalib:
00089       atmosflag=0; modtype=3; modnum=1; modpln=60; modstr=24; break;
00090     default:
00091       atmosflag=0; modtype=-1; modnum=0; modpln=0; modstr=0; break;
00092   }
00093 
00094   int i,j;//,k;
00095   //int vuw;
00096 
00097   double begZsm1=-9999.9,endZsm1=-9999.9;
00098   double begZsm2=-9999.9,endZsm2=-9999.9;
00099   int begsm1=-1,endsm1=-1,begsm2=-1,endsm2=-1;
00100 
00101   for(j=0;j<486;j++){
00102     PlexPlaneId myplaneid(vldc->GetDetector(),j,1);
00103     UgliSteelPlnHandle myplanehandle = ugh.GetSteelPlnHandle(myplaneid);
00104     if(myplanehandle.IsValid()){
00105       if(j<=modpln && begZsm1<0.0){ begsm1=j; begZsm1=myplanehandle.GetZ0(); }
00106       if(j<=modpln && myplanehandle.GetZ0()>endZsm1){ endsm1=j; endZsm1=myplanehandle.GetZ0(); }
00107       if(j>modpln && begZsm2<0.0){ begsm2=j; begZsm2=myplanehandle.GetZ0(); }
00108       if(j>modpln && myplanehandle.GetZ0()>endZsm2){ endsm2=j; endZsm2=myplanehandle.GetZ0(); }
00109     }
00110   }
00111 
00112   MSG("AlgShowerCam", Msg::kDebug) << " *** detector geometry *** " << endl
00113                                    << "  SM1 :  pln " << begsm1 << " -> " << endsm1
00114                                    << "    Z " << begZsm1 << " -> " << endZsm1 << endl
00115                                    << "  SM2 :  pln " << begsm2 << " -> " << endsm2 
00116                                    << "    Z " << begZsm2 << " -> " << endZsm2 << endl;
00117 
00119   // Find the extent of the shower in planes
00120   int bpln,epln;
00121   if(shwu->GetBegPlane()<shwv->GetBegPlane()) bpln=shwu->GetBegPlane(); else bpln=shwv->GetBegPlane();
00122   if(shwu->GetEndPlane()>shwv->GetEndPlane()) epln=shwu->GetEndPlane(); else epln=shwv->GetEndPlane();
00123   const int npln = 1+epln-bpln;  
00124 
00125 
00127   // Fill the vector with an empty TrkPlaneInfo object for each plane in the range spanned by the track
00128   PlaneInfo tmpplane;
00129   tmpplane.plnvuw=-1; tmpplane.plnnum=0;
00130   tmpplane.U=0.0; tmpplane.V=0.0; tmpplane.Z=0.0; tmpplane.Q=0.0; 
00131   tmpplane.Qm=0.0; tmpplane.Qp=0.0; tmpplane.CTm=0.0; tmpplane.CTp=0.0; 
00132   tmpplane.Wm=0.0; tmpplane.Wp=0.0;
00133   for(i=0;i<npln;i++) { fPlaneInfo.push_back(tmpplane); }
00134 
00136   // Calculate Track Co-ordinates
00137 
00138   MSG("AlgShowerCam", Msg::kDebug) << " *** calculate co-ordinates *** " << endl;
00139   ExtractHitProperties(shwu, bpln, shower);
00140   ExtractHitProperties(shwv, bpln, shower);
00141   int nplanes=0,nstrips=0;
00142   double totph=0.0;
00143   for(i=0;i<npln;i++)
00144     {
00145       if(fPlaneInfo[i].plnvuw>-1 && fPlaneInfo[i].Q>0.0)
00146         {
00147           totph+=fPlaneInfo[i].Q;
00148           fPlaneInfo[i].U/=fPlaneInfo[i].Q; fPlaneInfo[i].V/=fPlaneInfo[i].Q;
00149           nstrips+=fPlaneInfo[i].plnnum; ++nplanes;
00150         }
00151     }
00152 
00153   // Interpolate to get V for U planes and vice versa
00154   FindTransversePositions();
00155   // Calculate The Timing Info For Each Shower StripEnd
00156   SetupTimingInfo(&ugh, fFibreIndex);
00157   // Set U, V for each plane/stripend in the shower
00158   SetShowerCoordinates(bpln, shower);
00159 
00160 //   MSG("AlgShowerCam", Msg::kDebug) << " *** shower hits *** " << endl;
00161 //   for(pln=0;pln<npln;pln++){
00162 //     if( plnvuw[pln]>-1 ){
00163 //       MSG("AlgShowerCam",Msg::kDebug)
00164 //         << " pln=" << pln << " vuw=" << plnvuw[pln]
00165 //         << " z=" << Z[pln] << " u=" << U[pln] << " v=" << V[pln] << endl;
00166 //     }
00167 //   }
00168 
00169 
00170   // ***********************************
00171   // * D E T E R M I N E   V E R T E X * 
00172   // ***********************************
00173   vtxu=0.0,vtxv=0.0,vtxx=0.0,vtxy=0.0,vtxz=0.0; 
00174   vtxpln=-1; vtxshw=false;
00175   double trkscore=0.0,showerscore=0.0,dirscore=0.0;
00176   FindShowerVertex(shwu, shwv, trkscore, bpln, epln);
00177 
00178   //*****************************************
00179   //* D E T E R M I N E   D I R E C T I O N *
00180   //*****************************************
00181 
00182   diru=0.0; dirv=0.0; dirz=0.0; dirs=0.0; 
00183   direrru=0.0; direrrv=0.0; 
00184   double offset=0.0;
00185   DetermineDirection( shwu, shwv, offset, trkscore, showerscore, dirscore);
00186 
00187   // ***********************************
00188   // * C A L C U L A T E   E N E R G Y *
00189   // ***********************************
00190   double energy(0.0);
00191   CalculateShowerEnergy( shower, totph, energy );
00192 
00193   // determine overall direction 
00194   double dir=1.0;
00195   //if(dirscore>0.0) dir=1.0; 
00196   if(dirscore<0.0) dir=-1.0;  
00197 
00198   // *******************************************
00199   // * S E T   S H O W E R   V A R I A B L E S *
00200   // *******************************************
00201 
00202   // set shower variables
00203   shower.SetTimeSlope((dir)/(3.0e8)); 
00204   shower.SetTimeOffset((offset)/(3.0e8)); 
00205   shower.SetVtxU(vtxu);
00206   shower.SetVtxV(vtxv);
00207   shower.SetVtxZ(vtxz);
00208   //shower.SetVtxR(vtxr);
00209   shower.SetVtxT((offset)/(3.0e8));
00210   shower.SetVtxPlane(vtxpln);
00211   shower.SetDirCosU(dir*diru);
00212   shower.SetDirCosV(dir*dirv);
00213   shower.SetDirCosZ(dir*dirz);
00214   shower.SetEnergy(energy);
00215 
00216   MSG("AlgShowerCam", Msg::kDebug) 
00217     << " *** CREATING NEW SHOWER *** " << endl
00218     << "  NStrips = " << shower.GetNDaughters() << endl
00219     << "  TimeSlope = " << shower.GetTimeSlope() << endl
00220     << "  TimeOffset = " << shower.GetTimeOffset() << endl
00221     << "  VtxU = " << shower.GetVtxU() << endl
00222     << "  VtxV = " << shower.GetVtxV() << endl
00223     << "  VtxZ = " << shower.GetVtxZ() << endl
00224     << "  VtxT = " << shower.GetVtxT() << endl
00225     << "  VtxPln = " << shower.GetVtxPlane() << endl
00226     << "  DirCosU = " << shower.GetDirCosU() << endl
00227     << "  DirCosV = " << shower.GetDirCosV() << endl
00228     << "  DirCosZ = " << shower.GetDirCosZ() << endl
00229     << "  Energy = " << shower.GetEnergy() << endl;
00230 
00231   for(i=0;i<npln;i++)
00232     {
00233       fHitArr[i].clear();
00234     }
00235   fPlaneInfo.clear();
00236   
00237   return;
00238 }
00240 
00241 
00243 void AlgShowerCam::Trace(const char * /* c */) const
00244 {
00245 
00246 }
00248 
00249 
00251 void AlgShowerCam::ExtractHitProperties(const ShowerCamAtNu* shw, const int& bpln, CandShowerHandle& shower)
00252 {
00253   //Loop over all the HitCamAtNu objects in this ShowerCamAtNu object
00254   //Extract U/V/Z postion info as will as charge from each one
00255   //Add daughter links for the strips they represent to the CandShowerHandle.
00256   unsigned int nhits = shw-> GetNHits();
00257   int pln;
00258   for(unsigned int i=0;i<nhits;++i)
00259     {
00260       HitCamAtNu* hit = shw->GetHit(i);
00261       CandStripHandle* strip = hit->GetCandStripHandle();
00262       pln = hit->GetPlane()-bpln;
00263       fHitArr[pln].push_back(hit);
00264       if(fPlaneInfo[pln].plnvuw<0)
00265         {
00266           fPlaneInfo[pln].plnvuw = hit->GetPlaneView();
00267           fPlaneInfo[pln].Z = hit->GetZPos();
00268         }
00269       if(fPlaneInfo[pln].plnvuw==0) fPlaneInfo[pln].U += hit->GetCharge()*hit->GetTPos();
00270       if(fPlaneInfo[pln].plnvuw==1) fPlaneInfo[pln].V += hit->GetCharge()*hit->GetTPos();
00271       fPlaneInfo[pln].Q += hit->GetCharge();
00272       shower.AddDaughterLink(*strip);
00273       //shower.SetInShower(strip,hit->GetShwFlag());
00274       (fPlaneInfo[pln].plnnum)++; 
00275     }
00276 }
00278 
00279 
00281 void AlgShowerCam::FindTransversePositions()
00282 {
00283   const int npln = (int)fPlaneInfo.size();
00284   // interpolate between views for U,V
00285   int flag,vuw;
00286   int pln;
00287   int k,km1,kp1,km2,kp2;
00288   for(pln=0;pln<npln;pln++)
00289     {
00290       if(fPlaneInfo[pln].plnvuw>-1)
00291         {
00292           flag=0;
00293           vuw=fPlaneInfo[pln].plnvuw;
00294           k=0; km1=-1; km2=-1;
00295           while(pln-k>0 && km2<0)
00296             {
00297               k++; 
00298               if(fPlaneInfo[pln-k].plnvuw>-1 && fPlaneInfo[pln-k].plnvuw!=vuw)
00299                 {
00300                   if(km1<0) km1=k; else km2=k;
00301                 }
00302             }
00303           k=0; kp1=-1; kp2=-1;
00304           while(pln+k<npln-1 && kp2<0)
00305             {
00306               k++;
00307               if(fPlaneInfo[pln+k].plnvuw>-1 && fPlaneInfo[pln+k].plnvuw!=vuw)
00308                 {
00309                   if(kp1<0) kp1=k; else kp2=k;
00310                 }
00311             }
00312           if(km1>0 && kp1>0)
00313             { 
00314               if(vuw==0) fPlaneInfo[pln].V=(double(kp1)*fPlaneInfo[pln-km1].V+double(km1)*fPlaneInfo[pln+kp1].V)/double(km1+kp1);
00315               if(vuw==1) fPlaneInfo[pln].U=(double(kp1)*fPlaneInfo[pln-km1].U+double(km1)*fPlaneInfo[pln+kp1].U)/double(km1+kp1);
00316               flag=1;
00317             }
00318           else
00319             {
00320               if(km1>0)
00321                 {
00322                   if(vuw==0) 
00323                     {
00324                       if(km2>0) fPlaneInfo[pln].V=fPlaneInfo[pln-km1].V+(fPlaneInfo[pln].Z-fPlaneInfo[pln-km1].Z)*(fPlaneInfo[pln-km1].V-fPlaneInfo[pln-km2].V)/(fPlaneInfo[pln-km1].Z-fPlaneInfo[pln-km2].Z); 
00325                       else fPlaneInfo[pln].V=fPlaneInfo[pln-km1].V;
00326                     }
00327                   if(vuw==1) 
00328                     {
00329                       if(km2>0) fPlaneInfo[pln].U=fPlaneInfo[pln-km1].U+(fPlaneInfo[pln].Z-fPlaneInfo[pln-km1].Z)*(fPlaneInfo[pln-km1].U-fPlaneInfo[pln-km2].U)/(fPlaneInfo[pln-km1].Z-fPlaneInfo[pln-km2].Z); 
00330                       else fPlaneInfo[pln].U=fPlaneInfo[pln-km1].U;
00331                     }
00332                   flag=1;
00333                 }
00334               if(kp1>0)
00335                 {
00336                   if(vuw==0) 
00337                     {
00338                       if(kp2>0) fPlaneInfo[pln].V=fPlaneInfo[pln+kp1].V+(fPlaneInfo[pln].Z-fPlaneInfo[pln+kp1].Z)*(fPlaneInfo[pln+kp1].V-fPlaneInfo[pln+kp2].V)/(fPlaneInfo[pln+kp1].Z-fPlaneInfo[pln+kp2].Z); 
00339                       else fPlaneInfo[pln].V=fPlaneInfo[pln+kp1].V;
00340                     }
00341                   if(vuw==1) 
00342                     {
00343                       if(kp2>0) fPlaneInfo[pln].U=fPlaneInfo[pln+kp1].U+(fPlaneInfo[pln].Z-fPlaneInfo[pln+kp1].Z)*(fPlaneInfo[pln+kp1].U-fPlaneInfo[pln+kp2].U)/(fPlaneInfo[pln+kp1].Z-fPlaneInfo[pln+kp2].Z); 
00344                       else fPlaneInfo[pln].U=fPlaneInfo[pln+kp1].U;
00345                     }
00346                   flag=1;
00347                 }
00348             }
00349           if(flag==0){ fPlaneInfo[pln].plnvuw=-1; }
00350         }
00351     }
00352   return;
00353 //   // estimate transverse co-ordinates
00354 //   int km,kp;
00355 //   double q,sq,sqt;
00356 //   for(pln=0;pln<npln;pln++){
00357 //     if(plnvuw[pln]>-1){
00358 //       vuw=plnvuw[pln];
00359 //       sq=0.0; sqt=0.0; 
00360 //       k=0; km=-1;
00361 //       while(pln-k>0 && km<0){
00362 //         k++; if(plnvuw[pln-k]>-1 && plnvuw[pln-k]!=vuw) km=k;
00363 //       }
00364 //       k=0; kp=-1;
00365 //       while(pln+k<npln-1 && kp<0){
00366 //         k++; if(plnvuw[pln+k]>-1 && plnvuw[pln+k]!=vuw) kp=k; 
00367 //       }
00368 
00369 //       if(km>0){
00370 //         if(vuw==0){ q=1.0/km; sq+=q*Q[pln-km]; sqt+=q*Q[pln-km]*V[pln-km]; }
00371 //         if(vuw==1){ q=1.0/km; sq+=q*Q[pln-km]; sqt+=q*Q[pln-km]*U[pln-km]; }
00372 //       }
00373 //       if(kp>0){
00374 //         if(vuw==0){ q=1.0/kp; sq+=q*Q[pln+kp]; sqt+=q*Q[pln+kp]*V[pln+kp]; }
00375 //         if(vuw==1){ q=1.0/kp; sq+=q*Q[pln+kp]; sqt+=q*Q[pln+kp]*U[pln+kp]; }
00376 //       }
00377 
00378 //       if(sq>0.0){
00379 //         if(vuw==0){ V[pln]=sqt/sq; } if(vuw==1){ U[pln]=sqt/sq; }
00380 //       }
00381 //       else{
00382 //         if(vuw==0){ V[pln]=0.0; } if(vuw==1){ U[pln]=0.0; }
00383 //       }
00384 //     }
00385 //   }
00386 }
00388 
00389 
00391 void AlgShowerCam::SetupTimingInfo(UgliGeomHandle* ugh, const double& fFibreIndex)
00392 {
00393   double tpos = 0., fibre = 0., ctime =0., digchg = 0.;
00394   HitCamAtNu* hit = 0; CandStripHandle* strip = 0;
00395   const unsigned int npln = fPlaneInfo.size();
00396   for(unsigned int pln=0;pln<npln;pln++)
00397     {
00398       for(unsigned int g=0;g<fHitArr[pln].size();++g)
00399         {
00400           hit = fHitArr[pln][g];
00401           strip = hit->GetCandStripHandle();  
00402           if( strip->GetPlaneView()==PlaneView::kU
00403               || strip->GetPlaneView()==PlaneView::kX ) tpos = -fPlaneInfo[pln].V;
00404           if( strip->GetPlaneView()==PlaneView::kV
00405               || strip->GetPlaneView()==PlaneView::kY ) tpos = fPlaneInfo[pln].U;
00406           PlexStripEndId stripid = strip->GetStripEndId();
00407           UgliStripHandle striphandle = ugh->GetStripHandle(stripid);
00408 
00409           if(strip->GetNDigit(StripEnd::kPositive)>0)
00410             { 
00411               fibre = striphandle.ClearFiber(StripEnd::kPositive)+striphandle.WlsPigtail(StripEnd::kPositive)+striphandle.GetHalfLength()-tpos; 
00412               ctime = 3.0e8*strip->GetTime(StripEnd::kPositive)-fFibreIndex*fibre; 
00413               digchg=strip->GetCharge(StripEnd::kPositive); 
00414               fPlaneInfo[pln].Qp+=digchg; fPlaneInfo[pln].CTp+=digchg*ctime;
00415             }
00416         
00417           if(strip->GetNDigit(StripEnd::kNegative)>0)
00418             { 
00419               fibre = striphandle.ClearFiber(StripEnd::kNegative)+striphandle.WlsPigtail(StripEnd::kNegative)+striphandle.GetHalfLength()+tpos; 
00420               ctime = 3.0e8*strip->GetTime(StripEnd::kNegative)-fFibreIndex*fibre; 
00421               digchg=strip->GetCharge(StripEnd::kNegative); 
00422               fPlaneInfo[pln].Qm+=digchg; fPlaneInfo[pln].CTm+=digchg*ctime;
00423             }   
00424         }
00425     }
00426 }
00428 
00429 
00431 void AlgShowerCam::SetShowerCoordinates(const int& bpln, CandShowerAtNuHandle& shower)
00432 {//int& plnm, int& plnp, 
00433   int pln;
00434   const int npln = (int)fPlaneInfo.size();
00435   for(int i=0;i<npln;i++)
00436     {
00437       pln=bpln+i;
00438       if(fPlaneInfo[i].plnvuw>-1)
00439         {
00440           shower.SetU(pln,fPlaneInfo[i].U); shower.SetV(pln,fPlaneInfo[i].V);
00441           if(fPlaneInfo[i].Qm>0.0)
00442             {
00443               fPlaneInfo[i].CTm/=fPlaneInfo[i].Qm;
00444             }
00445           if(fPlaneInfo[i].Qp>0.0)
00446             {
00447               fPlaneInfo[i].CTp/=fPlaneInfo[i].Qp;
00448             }
00449           //if(plnm<0||i<plnm) plnm=i; if(plnp<0||i>plnp) plnp=i;
00450         }
00451     }
00452   return;
00453 
00454 //   // set map values
00455 //   for(i=0;i<npln;i++){
00456 //     pln=bpln+i;
00457 //     if(plnvuw[i]>-1){
00458 //       shower.SetU(pln,U[i]);
00459 //       shower.SetV(pln,V[i]);    
00460 //       if(Qm[i]>0.0){
00461 //         CTm[i]=CTm[i]/Qm[i];
00462 //       }
00463 //       if(Qp[i]>0.0){
00464 //         CTp[i]=CTp[i]/Qp[i];
00465 //       }
00466 //       if(Qm[i]+Qp[i]>0){
00467 //         CT[i]=(Qm[i]*CTm[i]+Qp[i]*CTp[i])/(Qm[i]+Qp[i]);
00468 //       }
00469 //     }
00470 //   }
00471 }
00473 
00474 
00476 void AlgShowerCam::FindShowerVertex( ShowerCamAtNu* shwu, ShowerCamAtNu* shwv, double& trkscr, const int& bpln, const int& epln)
00477 {
00478 
00479   MSG("AlgShowerCam", Msg::kDebug) << " *** determine vertex  *** " << endl;
00480   int trkbpln=-1,trkepln=-1,shwbpln=-1,shwepln=-1;
00481   unsigned int npln =  fPlaneInfo.size();
00482   unsigned int NHitsInPlane;
00483   // determine vertex plane
00484   double totchgpln=0.0,totchg=0.0;
00485   int pln,k,km,kp;
00486   for(unsigned int i=0;i<npln;i++)
00487     {
00488       NHitsInPlane= fHitArr[i].size();
00489       for(unsigned int j=0;j<NHitsInPlane;j++)
00490         {
00491           HitCamAtNu* hit = fHitArr[i][j];
00492           if(hit->GetShwFlag()>1)
00493             {
00494               totchg+=hit->GetCharge();
00495               totchgpln+=hit->GetCharge()*hit->GetPlane(); 
00496               if(shwbpln<0||hit->GetPlane()<shwbpln) shwbpln=hit->GetPlane();
00497               if(shwepln<0||hit->GetPlane()>shwepln) shwepln=hit->GetPlane(); 
00498             } 
00499         }
00500     }
00501   if(totchg>0.0)
00502     {
00503       vtxpln=(int)(totchgpln/totchg);
00504     }
00505   else
00506     {
00507       vtxpln=(int)(0.5*(bpln+epln));
00508       shwbpln=bpln; shwepln=epln;
00509     }
00510 
00511   const TrackCamAtNu* vtxtrkU = shwu->GetAssocTrackAt(0);
00512   const TrackCamAtNu* vtxtrkV = shwv->GetAssocTrackAt(0);
00513   if(vtxtrkU && vtxtrkV)
00514     {
00515       vtxshw=1;
00516       if(vtxtrkU->GetBegPlane()<vtxtrkV->GetBegPlane()) trkbpln=vtxtrkU->GetBegPlane(); else trkbpln=vtxtrkV->GetBegPlane();
00517       if(vtxtrkU->GetEndPlane()>vtxtrkV->GetEndPlane()) trkepln=vtxtrkU->GetEndPlane(); else trkepln=vtxtrkV->GetEndPlane();
00518       trkscr = (trkbpln+trkepln-2.0*vtxpln)/(1.0+trkepln-trkbpln);
00519       if(trkscr>0.0){ trkscr = (shwbpln+shwepln-2.0*trkbpln)/(1.0+shwepln-shwbpln); }
00520       if(trkscr<0.0){ trkscr = (shwbpln+shwepln-2.0*trkepln)/(1.0+shwepln-shwbpln); }
00521       if( ( shwu->GetBegVtxFlag() && shwv->GetBegVtxFlag())
00522           && !( shwu->GetEndVtxFlag() && shwv->GetEndVtxFlag()) ) vtxpln=trkbpln;
00523       if( ( shwu->GetEndVtxFlag() && shwv->GetEndVtxFlag())
00524           && !( shwu->GetBegVtxFlag() && shwv->GetBegVtxFlag()) ) vtxpln=trkepln;
00525     }
00526   
00527   // determine vertex position
00528   pln=vtxpln-bpln;
00529   if( pln>-1 && pln<(int)npln && fPlaneInfo[pln].plnvuw>-1 )
00530     {
00531       vtxu=fPlaneInfo[pln].U; vtxv=fPlaneInfo[pln].V; vtxz=fPlaneInfo[pln].Z;
00532     }
00533   else
00534     {
00535       k=0; km=-1;
00536       while(pln-k>0 && km<0)
00537         {
00538           k++; if(pln-k<(int)npln && fPlaneInfo[pln-k].plnvuw>-1) km=k;
00539         }
00540       k=0; kp=-1;
00541       while(pln+k<(int)npln-1 && kp<0)
00542         {
00543           k++; if(pln+k>-1 && fPlaneInfo[pln+k].plnvuw>-1) kp=k; 
00544     }
00545       if(km>0 && kp>0)
00546         {
00547           vtxu=(km*fPlaneInfo[pln+kp].U+kp*fPlaneInfo[pln-km].U)/(km+kp); 
00548           vtxv=(km*fPlaneInfo[pln+kp].V+kp*fPlaneInfo[pln-km].V)/(km+kp); 
00549           vtxz=(km*fPlaneInfo[pln+kp].Z+kp*fPlaneInfo[pln-km].Z)/(km+kp); 
00550         }
00551       else
00552         {
00553           if(km>0){ vtxu=fPlaneInfo[pln-km].U; vtxv=fPlaneInfo[pln-km].V; vtxz=fPlaneInfo[pln-km].Z; } 
00554           if(kp>0){ vtxu=fPlaneInfo[pln+kp].U; vtxv=fPlaneInfo[pln+kp].V; vtxz=fPlaneInfo[pln+kp].Z; }
00555         }
00556     }
00557   
00558   vtxx=0.7071*(vtxu-vtxv); vtxy=0.7071*(vtxu+vtxv); 
00559   
00560   MSG("AlgShowerCam", Msg::kVerbose) 
00561     << "   ... vertex: "
00562     << " plane = " << vtxpln
00563     << " position = (" << vtxu << "," << vtxv << "," << vtxz << ")" << endl;
00564   
00565   return;
00566 }
00568 
00569 
00570 
00572 void AlgShowerCam::DetermineDirection( ShowerCamAtNu* shwu, ShowerCamAtNu* shwv, double& offset, double& trkscr, double& shwscr, double& dirscr)
00573 {
00574   MSG("AlgShowerCam", Msg::kDebug) << " *** determine direction *** " << endl;
00575 
00576   // measure overall timeslope
00577   double ctimeslope=0.0,ctimeoffset=0.0,ctimeaverage=0.0,ctimescatter=0.0;
00578   double sw,swz,swt,swzz,swzt,swtt;
00579   double du,dv,dz;
00580   double z,t,w;
00581 
00582   int flag,npts;
00583 
00584   npts=0;
00585   sw=0.0; swz=0.0; swzz=0.0; swt=0.0; swtt=0.0; swzt=0.0;
00586   unsigned int npln = fPlaneInfo.size();
00587   for(unsigned int i=0;i<npln;i++)
00588     {
00589       if( fPlaneInfo[i].plnvuw>-1 ){
00590         z=fPlaneInfo[i].Z; t=fPlaneInfo[i].CTm; w=fPlaneInfo[i].Qm;
00591         swz+=w*z; swzz+=w*z*z; swt+=w*t; swtt+=w*t*t; swzt+=w*z*t;
00592         sw+=w; 
00593         t=fPlaneInfo[i].CTp; w=fPlaneInfo[i].Qp;
00594         swz+=w*z; swzz+=w*z*z; swt+=w*t; swtt+=w*t*t; swzt+=w*z*t;
00595         sw+=w; 
00596         npts++;
00597       }
00598     }
00599   if(npts>1)
00600     {
00601       ctimeslope=(sw*swzt-swz*swt)/(sw*swzz-swz*swz); 
00602       ctimeoffset=(swt*swzz-swz*swzt)/(sw*swzz-swz*swz);
00603       ctimeaverage=(swt/sw);
00604       if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 )
00605         {
00606           ctimescatter=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00607         }
00608     }
00609   offset=ctimeaverage;
00610   shwscr=ctimescatter;
00611 
00612   if(vtxshw) dirscr=trkscr; else dirscr=shwscr;
00613 
00614   npts=0;
00615   diru=0.0; direrru=0.0;
00616   if(1+shwu->GetEndPlane()-shwu->GetBegPlane()>1)
00617     {
00618       sw=0.0; swz=0.0; swt=0.0; 
00619       swzz=0.0; swzt=0.0; swtt=0.0;
00620       for(unsigned int i=0;i<npln;i++)
00621         {
00622           flag=0;
00623           for(unsigned int j=0;j<fHitArr[i].size();j++)
00624             {
00625               HitCamAtNu* hit = fHitArr[i][j];
00626               if(hit->GetShwFlag()>1 && hit->GetPlaneView()==0)
00627                 {
00628                   du=hit->GetTPos()-vtxu; dz=hit->GetZPos()-vtxz; w=hit->GetCharge();
00629                   swz+=w*dz; swt+=w*du; 
00630                   swzz+=w*dz*dz; swzt+=w*dz*du; swtt+=w*du*du;
00631                   sw+=w; flag=1;
00632                 }           
00633             }
00634           if(flag) npts++;
00635         }
00636       if(npts>1)
00637         {
00638           diru=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00639           if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 )
00640             {
00641               direrru=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00642             }
00643         }
00644     }
00645   
00646   npts=0;
00647   dirv=0.0; direrrv=0.0;
00648   if(1+shwv->GetEndPlane()-shwv->GetBegPlane()>1)
00649     {
00650       sw=0.0; swz=0.0; swt=0.0; 
00651       swzz=0.0; swzt=0.0; swtt=0.0;
00652       for(unsigned int i=0;i<npln;i++)
00653         {
00654           flag=0;
00655           for(unsigned int j=0;j<fHitArr[i].size();j++)
00656             {
00657               HitCamAtNu* hit = fHitArr[i][j];
00658               if(hit->GetShwFlag()>1 && hit->GetPlaneView()==1)
00659                 {
00660                   dv=hit->GetTPos()-vtxv; dz=hit->GetZPos()-vtxz; w=hit->GetCharge();
00661                   swz+=w*dz; swt+=w*dv; 
00662                   swzz+=w*dz*dz; swzt+=w*dz*dv; swtt+=w*dv*dv;
00663                   sw+=w; flag=1;
00664                 }
00665             }
00666           if(flag) npts++;
00667         }
00668       if(npts>1)
00669         {
00670           dirv=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00671           if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 )
00672             {
00673               direrrv=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00674             }
00675         }
00676     }
00677   
00678   dirs=sqrt(diru*diru+dirv*dirv+1.0);
00679   diru=diru/dirs; dirv=dirv/dirs; dirz=1.0/dirs;
00680   
00681   MSG("AlgShowerCam", Msg::kVerbose) 
00682     << "   ... direction: "
00683     << "  (" << diru << "," << dirv << "," << dirz << ")" << endl;
00684   return;
00685 }
00687 
00688 
00689 
00691 void AlgShowerCam::CalculateShowerEnergy( CandShowerHandle& shower, const double& totph, double& energy )
00692 {
00693   Calibrator& cal = Calibrator::Instance();
00694   MSG("AlgShowerCam", Msg::kDebug) << " *** calculate energy *** " << endl;
00695   // do the SIGMAP/SIGMIP calibration
00696   PlexStripEndId plexseid;
00697   FloatErr sigcorr,sigmip,sigmap;
00698   FloatErr lpos;
00699   unsigned int npln = fPlaneInfo.size();
00700   for(unsigned int pln=0;pln<npln;pln++)
00701     {
00702       for(unsigned int i=0;i<fHitArr[pln].size();i++)
00703         {
00704           HitCamAtNu* hit = fHitArr[pln][i];
00705           CandStripHandle* strip = hit->GetCandStripHandle();
00706           
00707           lpos = 0.0;
00708           if( strip->GetPlaneView()==PlaneView::kU || strip->GetPlaneView()==PlaneView::kX ) lpos = fPlaneInfo[pln].V;
00709           if( strip->GetPlaneView()==PlaneView::kV || strip->GetPlaneView()==PlaneView::kY ) lpos = fPlaneInfo[pln].U;
00710           lpos.SetError(1.0);
00711           
00712           if(strip->GetNDigit(StripEnd::kPositive)>0)
00713             {
00714               plexseid=strip->GetStripEndId(StripEnd::kPositive);
00715               sigcorr=strip->GetCharge(StripEnd::kPositive,CalDigitType::kSigCorr);
00716               sigmap=cal.GetAttenCorrectedTpos(sigcorr,lpos,plexseid);
00717               sigmip=cal.GetMIP(sigmap,plexseid);      
00718               shower.CalibrateSigMapped(plexseid.GetEncoded(),sigmap);
00719               shower.CalibrateMIP(plexseid.GetEncoded(),sigmip);
00720             }
00721           
00722           if(strip->GetNDigit(StripEnd::kNegative)>0)
00723             { 
00724               plexseid=strip->GetStripEndId(StripEnd::kNegative);
00725               sigcorr=strip->GetCharge(StripEnd::kNegative,CalDigitType::kSigCorr);
00726               sigmap=cal.GetAttenCorrectedTpos(sigcorr,lpos,plexseid);
00727               sigmip=cal.GetMIP(sigmap,plexseid);      
00728               shower.CalibrateSigMapped(plexseid.GetEncoded(),sigmap);
00729               shower.CalibrateMIP(plexseid.GetEncoded(),sigmip);
00730             }
00731         }
00732     }
00733   // calculate energy from pulse height
00734   if(totph>0.0)
00735     { 
00736       energy = totph; // NEEDS TUNING
00737       energy = 0.3 + totph/150.0;
00738     }
00739   return;
00740 }