AtmosCalculator Class Reference

#include <AtmosCalculator.h>

List of all members.

Public Member Functions

 AtmosCalculator ()
virtual ~AtmosCalculator ()
virtual void EventProperties (AtmosEvent *myevent, TClonesArray *StripList)
virtual void TrackProperties (AtmosTrack *mytrack, TClonesArray *StripList)
virtual void ShowerProperties (AtmosShower *myshower, TClonesArray *StripList)

Private Member Functions

bool CalculateTrace (double *m, double *c, double *coord, double *trace)
void DetectorSides (double *m, double *c, double *position, int side)

Detailed Description

Definition at line 9 of file AtmosCalculator.h.


Constructor & Destructor Documentation

AtmosCalculator::AtmosCalculator (  ) 

Definition at line 17 of file AtmosCalculator.cxx.

References Msg::kDebug, and MSG.

00018 {
00019   MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::AtmosCalculator()  " << endl;
00020 }

AtmosCalculator::~AtmosCalculator (  )  [virtual]

Definition at line 22 of file AtmosCalculator.cxx.

References Msg::kDebug, and MSG.

00023 {
00024   MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::~AtmosCalculator()  " << endl;
00025 }


Member Function Documentation

bool AtmosCalculator::CalculateTrace ( double *  m,
double *  c,
double *  coord,
double *  trace 
) [private]

Definition at line 1135 of file AtmosCalculator.cxx.

References DetectorSides(), Msg::kDebug, and MSG.

Referenced by TrackProperties().

01137 { 
01138   MSG("AtmosCalculator",Msg::kDebug) << "CalculateTrace" << endl;
01139 
01140   bool TraceSuccess=false;
01141   double TempTrace[3]={0.,0.,0.}; double TempSum=0.;
01142   double TraceLength(1.e99);
01143 
01144   for(int i=0; i<8; ++i) { // Consider the eight sides in turn
01145     TempSum=0.;
01146 
01147     DetectorSides(m,c,TempTrace, i);
01148 
01149     for(int j=0; j<3; ++j) {TempSum+=pow(Coord[j]-TempTrace[j],2);}
01150     TempSum=pow(TempSum,0.5);
01151 
01152     if(TempSum<TraceLength) {
01153       for(int j=0; j<3; ++j) {Trace[j]=TempTrace[j]-Coord[j];}
01154       TraceLength=TempSum;
01155       TraceSuccess=true;
01156     }
01157   }
01158 
01159   return TraceSuccess;
01160 }

void AtmosCalculator::DetectorSides ( double *  m,
double *  c,
double *  position,
int  side 
) [private]

Definition at line 1166 of file AtmosCalculator.cxx.

Referenced by CalculateTrace().

01167 {
01168   position[2]=-1.e99;
01169   
01170   switch(side) {
01171   case 0: if(m[0]!=-m[1]) {position[2] = (pow(32.,0.5)-c[0]-c[1])/(m[0]+m[1]);}  break;
01172   case 1: if(m[0]!=0.)    {position[2] = (4.-c[0])/m[0];}                        break;
01173   case 2: if(m[0]!=m[1])  {position[2] = (pow(32.,0.5)-c[0]+c[1])/(m[0]-m[1]);}  break;
01174   case 3: if(m[1]!=0.)    {position[2] = (-4.-c[1])/m[1];}                       break;
01175   case 4: if(m[0]!=-m[1]) {position[2] = (-pow(32.,0.5)-c[0]-c[1])/(m[0]+m[1]);} break;
01176   case 5: if(m[0]!=0.)    {position[2] = (-4.-c[0])/m[0];}                       break;
01177   case 6: if(m[0]!=m[1])  {position[2] = (-pow(32.,0.5)-c[0]+c[1])/(m[0]-m[1]);} break;
01178   case 7: if(m[1]!=0.)    {position[2] = (4.-c[1])/m[1];}                        break;
01179   default: position[2] = -1.e99; break;
01180   }
01181   
01182   if(position[2]!=-1.e99) {
01183     position[0]=m[0]*position[2]+c[0];
01184     position[1]=m[1]*position[2]+c[1];  
01185   }
01186   else {position[0]=-1.e99; position[1]=-1.e99;}
01187 }

void AtmosCalculator::EventProperties ( AtmosEvent myevent,
TClonesArray *  StripList 
) [virtual]

Definition at line 27 of file AtmosCalculator.cxx.

References AtmosReco::DoubleEndedCharge, AtmosReco::DoubleEndedChargeUview, AtmosReco::DoubleEndedChargeVview, AtmosReco::EastSideCharge, Msg::kDebug, MSG, AtmosReco::MultiMuDirCosU, AtmosReco::MultiMuDirCosV, AtmosReco::MultiMuDirCosX, AtmosReco::MultiMuDirCosY, AtmosReco::MultiMuDirCosZ, AtmosReco::MultiMuReco, AtmosReco::MultiMuTracks, AtmosStrip::Ndigits, AtmosStrip::Plane, AtmosEvent::RecoInfo, AtmosStrip::Sigcorr, AtmosReco::SingleEndedCharge, AtmosReco::SingleEndedChargeUview, AtmosReco::SingleEndedChargeVview, AtmosReco::TotalCharge, AtmosStrip::View, and AtmosReco::WestSideCharge.

Referenced by NtpMaker::FillCandInfo().

00028 {
00029   MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::EventProperties(...)  " << endl;
00030   
00031   // pulse height stuff
00032 
00033   double dQ(0.0),dQe(0.0),dQw(0.0);
00034   double totalQ(0.0),eastQ(0.0),westQ(0.0);
00035   double singleQ(0.0),singleQU(0.0),singleQV(0.0);
00036   double doubleQ(0.0),doubleQU(0.0),doubleQV(0.0);
00037 
00038   int ListEnd = StripList->GetLast()+1;
00039   for(int strp=0; strp<ListEnd; ++strp)
00040     {
00041       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00042       if(strip.Plane>-1 && strip.Plane<486)
00043         {
00044           dQe=strip.Sigcorr[0]/65.0;
00045           dQw=strip.Sigcorr[1]/65.0;
00046           dQ=(strip.Sigcorr[0]+strip.Sigcorr[1])/65.0; 
00047             // hack a rough conversion from SigCorr -> PEcorr
00048 
00049           totalQ += dQ;
00050           eastQ += dQe; westQ += dQw;
00051 
00052           if( strip.Ndigits==1 )
00053             {
00054               singleQ += dQ;
00055               if( strip.View==0 ) singleQU += dQ;
00056               if( strip.View==1 ) singleQV += dQ;
00057             }
00058           if( strip.Ndigits==2 )
00059             {
00060               doubleQ += dQ;
00061               if( strip.View==0 ) doubleQU += dQ;
00062               if( strip.View==1 ) doubleQV += dQ;
00063             }
00064 
00065         }
00066     }
00067 
00068   MyEvent->RecoInfo.TotalCharge = totalQ;
00069   MyEvent->RecoInfo.EastSideCharge = eastQ;
00070   MyEvent->RecoInfo.WestSideCharge = westQ;
00071   MyEvent->RecoInfo.SingleEndedCharge = singleQ;
00072   MyEvent->RecoInfo.SingleEndedChargeUview = singleQU;
00073   MyEvent->RecoInfo.SingleEndedChargeVview = singleQV;
00074   MyEvent->RecoInfo.DoubleEndedCharge = doubleQ;
00075   MyEvent->RecoInfo.DoubleEndedChargeUview = doubleQU;
00076   MyEvent->RecoInfo.DoubleEndedChargeVview = doubleQV;
00077 
00078   // multiple muon reconstruction goes here
00079   MyEvent->RecoInfo.MultiMuReco = 0;
00080   MyEvent->RecoInfo.MultiMuTracks = 0;
00081   MyEvent->RecoInfo.MultiMuDirCosU = 0.0;
00082   MyEvent->RecoInfo.MultiMuDirCosV = 0.0;
00083   MyEvent->RecoInfo.MultiMuDirCosX = 0.0;
00084   MyEvent->RecoInfo.MultiMuDirCosY = 0.0;
00085   MyEvent->RecoInfo.MultiMuDirCosZ = 0.0;
00086 
00087   return;
00088 }

void AtmosCalculator::ShowerProperties ( AtmosShower myshower,
TClonesArray *  StripList 
) [virtual]

Definition at line 874 of file AtmosCalculator.cxx.

References AtmosShower::AssocShwPH, AtmosShower::AssocShwPHfrac, ShowerMOI::EigenValues, ShowerMOI::EigenVectors, AtmosShower::Index, Msg::kDebug, Msg::kVerbose, AtmosStrip::L, AtmosShower::MaxChargeInPlane, AtmosShower::MaxPlaneNumber, AtmosShower::MaxStripsInPlane, AtmosShower::MeanChargePerPlane, AtmosShower::MeanDCoreHit, AtmosShower::MeanDCosHit, AtmosShower::MeanStripsPerPlane, AtmosShower::MinPlaneNumber, AtmosShower::MOIUVZEigenValue0, AtmosShower::MOIUVZEigenValue1, AtmosShower::MOIUVZEigenValue2, AtmosShower::MOIUVZEigenVector0, AtmosShower::MOIUVZEigenVector1, AtmosShower::MOIUVZEigenVector2, AtmosShower::MOIUZEigenValue0, AtmosShower::MOIUZEigenValue1, AtmosShower::MOIUZEigenVector0, AtmosShower::MOIUZEigenVector1, AtmosShower::MOIVZEigenValue0, AtmosShower::MOIVZEigenValue1, AtmosShower::MOIVZEigenVector0, AtmosShower::MOIVZEigenVector1, MSG, AtmosShower::Ndigits, AtmosStrip::Ndigits, AtmosShower::NonFidFrac, nPlanes, AtmosShower::NplanesUview, AtmosShower::NplanesVview, AtmosShower::NstripsDoubleEnded, AtmosShower::NstripsSingleEnded, AtmosStrip::Plane, ShowerTrace::ReTrace, ShowerTrace::ReTraceZ, AtmosShower::RMSChargePerPlane, AtmosShower::RMSDCoreHit, AtmosShower::RMSDCosHit, AtmosShower::RMSStripsPerPlane, AtmosStrip::Shw, AtmosStrip::Sigcorr, AtmosStrip::T, ShowerTrace::Trace, ShowerTrace::TraceZ, AtmosShower::TrkLikePlanes, ShowerTrace::UpTrace, ShowerTrace::UpTraceZ, AtmosShower::UVassymetry, AtmosStrip::View, AtmosShower::VtxDirCosU, AtmosShower::VtxDirCosV, AtmosShower::VtxDirCosZ, AtmosShower::VtxDistToEdge, AtmosShower::VtxForTrace, AtmosShower::VtxForTraceZ, AtmosShower::VtxRevTrace, AtmosShower::VtxRevTraceZ, AtmosShower::VtxU, AtmosShower::VtxUpTrace, AtmosShower::VtxUpTraceZ, AtmosShower::VtxV, AtmosShower::VtxX, AtmosShower::VtxY, AtmosShower::VtxZ, and AtmosStrip::Z.

Referenced by NtpMaker::FillCandInfo().

00875 {
00876   MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::ShowerProperties(...)  " << endl;
00877 
00878   ShowerTrace ShwTrace(MyShower);
00879   MyShower->VtxForTrace = ShwTrace.Trace;
00880   MyShower->VtxForTraceZ = ShwTrace.TraceZ;
00881   MyShower->VtxRevTrace = ShwTrace.ReTrace;
00882   MyShower->VtxRevTraceZ = ShwTrace.ReTraceZ;
00883   MyShower->VtxUpTrace = ShwTrace.UpTrace;
00884   MyShower->VtxUpTraceZ = ShwTrace.UpTraceZ;
00885 
00886   // Get Shower Index
00887   Int_t index = MyShower->Index;
00888 
00889   // Closest distance to edge of detector
00890   MSG("AtmosCalculator",Msg::kVerbose) << "    Distance from Edge of Detector" << endl;
00891   double vtxdiru=MyShower->VtxDirCosU;
00892   double vtxdirv=MyShower->VtxDirCosV;
00893   double vtxdirz=MyShower->VtxDirCosZ;
00894   double vtxu=MyShower->VtxU;  
00895   double vtxv=MyShower->VtxV;  
00896   double vtxx=MyShower->VtxX;  
00897   double vtxy=MyShower->VtxY;
00898   double vtxz=MyShower->VtxZ;
00899 
00900   Double_t rmin,temp;
00901 
00902   rmin=4.0;
00903   temp=4.0-vtxu; if(temp<rmin) rmin=temp;
00904   temp=4.0+vtxu; if(temp<rmin) rmin=temp;
00905   temp=4.0-vtxv; if(temp<rmin) rmin=temp;
00906   temp=4.0+vtxv; if(temp<rmin) rmin=temp;
00907   temp=4.0-vtxx; if(temp<rmin) rmin=temp;
00908   temp=4.0+vtxx; if(temp<rmin) rmin=temp;
00909   temp=4.0-vtxy; if(temp<rmin) rmin=temp;
00910   temp=4.0+vtxy; if(temp<rmin) rmin=temp;
00911   MyShower->VtxDistToEdge = rmin;//*************
00912 
00913   // Highest and lowest plane number in shower
00914   MSG("AtmosCalculator",Msg::kVerbose) << "    Highest and Lowest Plane Number" << endl;
00915   int thispln;
00916   int minplane=-1,maxplane=-1;
00917   Int_t nSingleEnds(0), nDoubleEnds(0), nDigits(0);
00918   Double_t dQ;
00919   Double_t nonfidShwPH(0.), totShwPH(0.), totPlnPH(0.);
00920   Int_t shwplnstp[500];
00921   Double_t shwplnchg[500];
00922   Double_t plnchg[500];
00923 
00924   memset(shwplnstp, 0, 500*sizeof(Int_t));
00925   memset(shwplnchg, 0, 500*sizeof(Double_t));
00926   memset(plnchg, 0, 500*sizeof(Double_t));
00927 
00928   vector<double> UStripsU, UStripsZ, UStripsE;
00929   vector<double> VStripsV, VStripsZ, VStripsE;
00930   vector<double> AllStripsU, AllStripsV, AllStripsZ, AllStripsE;
00931 
00932   double UStpVtxShift, VStpVtxShift, ZStpVtxShift, RStpVtxShift;
00933   double DotProd, DCosHit, DCoreHit;
00934   double SumDCosHit(0.0), SumSqDCosHit(0.0), SumDCoreHit(0.), SumSqDCoreHit(0.0);
00935   
00936   int ListEnd = StripList->GetLast()+1;
00937 
00938   double ucShwStrips(0.0), ShwStrips(0.0);
00939   bool UCsides(false);
00940   bool UCends(false);
00941   double Xdis(0.3);
00942   int Xpln(5);
00943   double stripX, stripY, stripMaxUVXY(-1.0);
00944 
00945   for(int strp=0; strp<ListEnd; ++strp)
00946     {
00947       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00948       thispln=strip.Plane;
00949       if(thispln<0 || thispln>485) continue;
00950 
00951       dQ=(strip.Sigcorr[0]+strip.Sigcorr[1])/65.0; // hack a rough conversion 
00952       plnchg[thispln] += dQ;
00953 
00954       totPlnPH += dQ;
00955 
00956       //Everything past this point is only for strips associated with
00957       //the shower
00958       if((strip.Shw&index)!=index) continue;
00959 
00960       if (strip.View == 0) {
00961         AllStripsU.push_back(strip.T);
00962         AllStripsV.push_back(strip.L);
00963 
00964         UStripsU.push_back(strip.T);
00965         UStripsZ.push_back(strip.Z);
00966         UStripsE.push_back(dQ);
00967       }
00968       else {
00969         AllStripsU.push_back(strip.L);
00970         AllStripsV.push_back(strip.T);
00971 
00972         VStripsV.push_back(strip.T);
00973         VStripsZ.push_back(strip.Z);
00974         VStripsE.push_back(dQ);
00975       }
00976 
00977       AllStripsE.push_back(dQ);
00978       AllStripsZ.push_back(strip.Z);
00979 
00980       shwplnstp[thispln]++;
00981       shwplnchg[thispln] += dQ;
00982 
00983       totShwPH += dQ;
00984       if(strip.Ndigits==1) ++nSingleEnds;
00985       if(strip.Ndigits==2) ++nDoubleEnds;
00986       nDigits += strip.Ndigits;
00987 
00988       if(minplane<0||thispln<minplane) minplane=thispln;
00989       if(maxplane<0||thispln>maxplane) maxplane=thispln;
00990 
00991       if(fabs(strip.T)>stripMaxUVXY) stripMaxUVXY = fabs(strip.T);
00992       if(fabs(strip.L)>stripMaxUVXY) stripMaxUVXY = fabs(strip.L);
00993       stripX = fabs(strip.T + strip.L) / sqrt(2.0);
00994       if(fabs(stripX)>stripMaxUVXY) stripMaxUVXY = fabs(stripX);
00995       stripY = fabs(strip.T - strip.L) / sqrt(2.0);
00996       if(fabs(stripY)>stripMaxUVXY) stripMaxUVXY = fabs(stripY);
00997       UCsides = ( stripMaxUVXY>(4-Xdis) );
00998       UCends  = ( (strip.Plane<=Xpln) || ((strip.Plane>=(249-Xpln)) &&
00999                 (strip.Plane<=(249+Xpln))) || (strip.Plane>=(486-Xpln)) );
01000       if( UCsides || UCends ) nonfidShwPH += dQ;
01001       if( UCsides || UCends ) ucShwStrips += 1.0;
01002       ShwStrips += 1.0;
01003 
01004       //Shift to vertex frame
01005       UStpVtxShift = (strip.View==0?strip.T:strip.L) - vtxu;
01006       VStpVtxShift = (strip.View==1?strip.T:strip.L) - vtxv;
01007       ZStpVtxShift = strip.Z - vtxz;
01008       RStpVtxShift = sqrt(UStpVtxShift*UStpVtxShift +
01009                           VStpVtxShift*VStpVtxShift +
01010                           ZStpVtxShift*ZStpVtxShift);
01011 
01012       if (RStpVtxShift != 0.0) {
01013         //Find the angle between and Shower Direction
01014         DotProd = fabs(UStpVtxShift*vtxdiru +
01015                        VStpVtxShift*vtxdirv + ZStpVtxShift*vtxdirz);
01016         DCosHit = DotProd / RStpVtxShift;
01017         if(RStpVtxShift <= DotProd) DCoreHit = 0.0;
01018         else DCoreHit = sqrt(RStpVtxShift*RStpVtxShift - DotProd*DotProd);
01019 
01020         SumDCosHit +=  DCosHit * dQ;
01021         SumSqDCosHit += DCosHit * DCosHit * dQ;
01022         SumDCoreHit +=  DCoreHit * dQ;
01023         SumSqDCoreHit += DCoreHit * DCoreHit * dQ;
01024       }
01025     }
01026 
01027   if( totShwPH>0.0 ) MyShower->NonFidFrac = nonfidShwPH/totShwPH;
01028 
01029   double ShwMeanDCosHit = SumDCosHit / totPlnPH;
01030   double ShwRMSDCosHit = (SumSqDCosHit/totPlnPH) - (ShwMeanDCosHit*ShwMeanDCosHit);
01031   if(ShwRMSDCosHit<=0) ShwRMSDCosHit = 0.0;
01032   else ShwRMSDCosHit = sqrt(ShwRMSDCosHit);
01033 
01034   MyShower->MeanDCosHit = ShwMeanDCosHit;
01035   MyShower->RMSDCosHit = ShwRMSDCosHit;
01036 
01037   double ShwMeanDCoreHit = SumDCoreHit / totPlnPH;
01038   double ShwRMSDCoreHit = (SumSqDCoreHit/totPlnPH) - (ShwMeanDCoreHit*ShwMeanDCoreHit);
01039   if(ShwRMSDCoreHit<=0) ShwRMSDCoreHit = 0.0;
01040   else ShwRMSDCoreHit = sqrt(ShwRMSDCoreHit);
01041 
01042   MyShower->MeanDCoreHit = ShwMeanDCoreHit;
01043   MyShower->RMSDCoreHit = ShwRMSDCoreHit;
01044 
01045   //Solve the 3D UVZ MOI
01046   ShowerMOI MOIUVZ(AllStripsU, AllStripsV, AllStripsZ, AllStripsE);
01047   MyShower->MOIUVZEigenValue0 = MOIUVZ.EigenValues[0];
01048   MyShower->MOIUVZEigenValue1 = MOIUVZ.EigenValues[1];
01049   MyShower->MOIUVZEigenValue2 = MOIUVZ.EigenValues[2];
01050   for (int i=0; i<3; i++) {
01051     MyShower->MOIUVZEigenVector0[i] = MOIUVZ.EigenVectors[0][i];
01052     MyShower->MOIUVZEigenVector1[i] = MOIUVZ.EigenVectors[1][i];
01053     MyShower->MOIUVZEigenVector2[i] = MOIUVZ.EigenVectors[2][i];
01054   }
01055 
01056   ShowerMOI MOIUZ(UStripsU, UStripsZ, UStripsE);
01057   MyShower->MOIUZEigenValue0 = MOIUZ.EigenValues[0];
01058   MyShower->MOIUZEigenValue1 = MOIUZ.EigenValues[1];
01059   for (int i=0; i<2; i++) {
01060     MyShower->MOIUZEigenVector0[i] = MOIUZ.EigenVectors[0][i];
01061     MyShower->MOIUZEigenVector1[i] = MOIUZ.EigenVectors[1][i];
01062   }
01063 
01064   ShowerMOI MOIVZ(VStripsV, VStripsZ, VStripsE);
01065   MyShower->MOIVZEigenValue0 = MOIVZ.EigenValues[0];
01066   MyShower->MOIVZEigenValue1 = MOIVZ.EigenValues[1];
01067   for (int i=0; i<2; i++) {
01068     MyShower->MOIVZEigenVector0[i] = MOIVZ.EigenVectors[0][i];
01069     MyShower->MOIVZEigenVector1[i] = MOIVZ.EigenVectors[1][i];
01070   }
01071 
01072   MyShower->MinPlaneNumber = minplane;
01073   MyShower->MaxPlaneNumber = maxplane;
01074             
01075   MyShower->Ndigits = nDigits;
01076   MyShower->NstripsSingleEnded = nSingleEnds;
01077   MyShower->NstripsDoubleEnded = nDoubleEnds;
01078 
01079   MyShower->AssocShwPH = totShwPH;
01080   MyShower->AssocShwPHfrac = totShwPH / totPlnPH;
01081   
01082   int uPlanes(0),vPlanes(0);
01083   double nPlanes(0.);
01084   int maxplnstp(0); double maxplnchg(0.0); 
01085   int sumstp(0); double sumsqstp(0.0); 
01086   double sumchg(0); double sumsqchg(0.0); 
01087   int ntrkplns(0);
01088   for(int ipln=minplane; ipln<=maxplane; ipln++)
01089     {
01090       if(shwplnstp[ipln] == 0) continue;
01091       if((ipln/249)^(ipln%2)) uPlanes++;
01092       else vPlanes++;
01093       nPlanes += 1.0;
01094 
01095       //Count up the tracklike planes
01096       if(plnchg[ipln]>0.0 && plnchg[ipln]<80.0) {
01097         if(shwplnstp[ipln] <=2) ntrkplns++;
01098       }
01099 
01100       if(shwplnstp[ipln] > maxplnstp) maxplnstp = shwplnstp[ipln];
01101       if(shwplnchg[ipln] > maxplnchg) maxplnchg = shwplnchg[ipln];
01102 
01103       sumstp += shwplnstp[ipln];
01104       sumsqstp += shwplnstp[ipln]*shwplnstp[ipln];
01105 
01106       sumchg += shwplnchg[ipln];
01107       sumsqchg += shwplnchg[ipln]*shwplnchg[ipln];
01108     }
01109 
01110   MyShower->NplanesUview = uPlanes;
01111   MyShower->NplanesVview = vPlanes;
01112   if(uPlanes+vPlanes>0) MyShower->UVassymetry = abs(uPlanes-vPlanes)/(0.5*(uPlanes+vPlanes));
01113   MyShower->TrkLikePlanes = ntrkplns;
01114 
01115   MyShower->MaxStripsInPlane = maxplnstp;
01116   MyShower->MaxChargeInPlane = maxplnchg;
01117 
01118   double ShwMeanStripsPerPlane = sumstp/nPlanes;
01119   double ShwRMSStripsPerPlane = (sumsqstp/nPlanes) - ShwMeanStripsPerPlane*ShwMeanStripsPerPlane;
01120   if(ShwRMSStripsPerPlane <= 0.0) ShwRMSStripsPerPlane = 0.0;
01121   else ShwRMSStripsPerPlane = sqrt(ShwRMSStripsPerPlane);
01122   MyShower->MeanStripsPerPlane = ShwMeanStripsPerPlane;
01123   MyShower->RMSStripsPerPlane = ShwRMSStripsPerPlane;
01124     
01125   double ShwMeanChargePerPlane = sumchg/nPlanes;
01126   double ShwRMSChargePerPlane = (sumsqchg/nPlanes) - ShwMeanChargePerPlane*ShwMeanChargePerPlane;
01127   if(ShwRMSChargePerPlane <= 0.0) ShwRMSChargePerPlane = 0.0;
01128   else ShwRMSChargePerPlane = sqrt(ShwRMSChargePerPlane);
01129   MyShower->MeanChargePerPlane = ShwMeanChargePerPlane;
01130   MyShower->RMSChargePerPlane = ShwRMSChargePerPlane;
01131     
01132   return;
01133 }

void AtmosCalculator::TrackProperties ( AtmosTrack mytrack,
TClonesArray *  StripList 
) [virtual]

Definition at line 90 of file AtmosCalculator.cxx.

References AtmosTrack::AssocTrkPH, AtmosTrack::AssocTrkPHfrac, CalculateTrace(), AtmosTrack::EndDirCosU, AtmosTrack::EndDirCosV, AtmosTrack::EndDirCosZ, AtmosTrack::EndDistToEdge, AtmosTrack::EndDistToEdgeDigits, AtmosTrack::EndLinearDirFitChisqU, AtmosTrack::EndLinearDirFitChisqV, AtmosTrack::EndLinearDirFitNdfU, AtmosTrack::EndLinearDirFitNdfV, AtmosTrack::EndPlane, AtmosTrack::EndPlaneDigits, AtmosTrack::EndQmax, AtmosTrack::EndRmax, AtmosTrack::EndTrace, AtmosTrack::EndTraceZ, AtmosTrack::EndU, AtmosTrack::EndUmean, AtmosTrack::EndUwidth, AtmosTrack::EndV, AtmosTrack::EndVmean, AtmosTrack::EndVwidth, AtmosTrack::EndX, AtmosTrack::EndY, AtmosTrack::EndZ, AtmosTrack::Index, Msg::kDebug, Msg::kVerbose, AtmosStrip::L, AtmosTrack::LinearDirCosU, AtmosTrack::LinearDirCosV, AtmosTrack::LinearDirCosZ, AtmosTrack::LinearDirFitChisq, AtmosTrack::LinearDirFitChisqU, AtmosTrack::LinearDirFitChisqV, AtmosTrack::LinearDirFitNdf, AtmosTrack::LinearDirFitNdfU, AtmosTrack::LinearDirFitNdfV, Munits::m, AtmosTrack::MaxPlaneNumber, AtmosTrack::MinPlaneNumber, MSG, AtmosStrip::Ndigits, AtmosTrack::Ndigits, AtmosTrack::NonFidFrac, AtmosTrack::NplanesTrackGaps, AtmosTrack::NplanesTrackOnly, AtmosTrack::NplanesUview, AtmosTrack::NplanesVview, AtmosTrack::NstripsDoubleEnded, AtmosTrack::NstripsSingleEnded, AtmosStrip::Plane, AtmosStrip::Sigcorr, Particle::Sign(), AtmosStrip::Strip, AtmosStrip::T, AtmosStrip::Trk, AtmosTrack::TrkLikePlanes, AtmosTrack::TrkPH, AtmosTrack::UVassymetry, AtmosStrip::View, AtmosTrack::VtxDirCosU, AtmosTrack::VtxDirCosV, AtmosTrack::VtxDirCosZ, AtmosTrack::VtxDistToEdge, AtmosTrack::VtxDistToEdgeDigits, AtmosTrack::VtxLinearDirFitChisqU, AtmosTrack::VtxLinearDirFitChisqV, AtmosTrack::VtxLinearDirFitNdfU, AtmosTrack::VtxLinearDirFitNdfV, AtmosTrack::VtxPlane, AtmosTrack::VtxPlaneDigits, AtmosTrack::VtxQmax, AtmosTrack::VtxRmax, AtmosTrack::VtxTrace, AtmosTrack::VtxTraceZ, AtmosTrack::VtxU, AtmosTrack::VtxUmean, AtmosTrack::VtxUwidth, AtmosTrack::VtxV, AtmosTrack::VtxVmean, AtmosTrack::VtxVwidth, AtmosTrack::VtxX, AtmosTrack::VtxY, AtmosTrack::VtxZ, AtmosStrip::Z, AtmosTrack::ZeroCurveChi2, and AtmosTrack::ZeroCurveNdf.

Referenced by NtpMaker::FillCandInfo().

00091 {
00092   MSG("AtmosCalculator",Msg::kDebug) << "  AtmosCalculator::TrackProperties(...)  " << endl;
00093 
00094   // Get Track Index
00095   int index=MyTrack->Index;
00096 
00097   // minimum/maximum plane number
00098   int vtxpln=MyTrack->VtxPlane; 
00099   int endpln=MyTrack->EndPlane;
00100   int minpln=(MyTrack->VtxPlane<MyTrack->EndPlane)?MyTrack->VtxPlane:MyTrack->EndPlane;
00101   int maxpln=(MyTrack->VtxPlane<MyTrack->EndPlane)?MyTrack->EndPlane:MyTrack->VtxPlane;
00102   double dir=(MyTrack->VtxPlane>MyTrack->EndPlane)?-1.0:1.0;
00103   MyTrack->MinPlaneNumber = minpln;
00104   MyTrack->MaxPlaneNumber = maxpln;
00105 
00106   // Set vertex/end position/direction
00107   double vtxdiru=MyTrack->VtxDirCosU; double enddiru=MyTrack->EndDirCosU;
00108   double vtxdirv=MyTrack->VtxDirCosV; double enddirv=MyTrack->EndDirCosV;
00109   double vtxdirz=MyTrack->VtxDirCosZ; double enddirz=MyTrack->EndDirCosZ;
00110   double vtxu=MyTrack->VtxU; double endu=MyTrack->EndU;
00111   double vtxv=MyTrack->VtxV; double endv=MyTrack->EndV;
00112   double vtxx=MyTrack->VtxX; double endx=MyTrack->EndX;
00113   double vtxy=MyTrack->VtxY; double endy=MyTrack->EndY;
00114   double vtxz=MyTrack->VtxZ; double endz=MyTrack->EndZ;
00115   
00116   // Closest distance to edge of detector
00117   MSG("AtmosCalculator",Msg::kVerbose) << "    Distance from Edge of Detector" << endl;
00118   double rmin,temp;
00119 
00120   rmin=4.0; 
00121   temp=4.0-vtxu; if(temp<rmin) rmin=temp;
00122   temp=4.0+vtxu; if(temp<rmin) rmin=temp;
00123   temp=4.0-vtxv; if(temp<rmin) rmin=temp;
00124   temp=4.0+vtxv; if(temp<rmin) rmin=temp;
00125   temp=4.0-vtxx; if(temp<rmin) rmin=temp;
00126   temp=4.0+vtxx; if(temp<rmin) rmin=temp;
00127   temp=4.0-vtxy; if(temp<rmin) rmin=temp;
00128   temp=4.0+vtxy; if(temp<rmin) rmin=temp;
00129   MyTrack->VtxDistToEdge = rmin; // VtxDistToEdge
00130 
00131   rmin=4.0; 
00132   temp=4.0-endu; if(temp<rmin) rmin=temp;
00133   temp=4.0+endu; if(temp<rmin) rmin=temp;
00134   temp=4.0-endv; if(temp<rmin) rmin=temp;
00135   temp=4.0+endv; if(temp<rmin) rmin=temp;
00136   temp=4.0-endx; if(temp<rmin) rmin=temp;
00137   temp=4.0+endx; if(temp<rmin) rmin=temp;
00138   temp=4.0-endy; if(temp<rmin) rmin=temp;
00139   temp=4.0+endy; if(temp<rmin) rmin=temp; 
00140   MyTrack->EndDistToEdge = rmin; // EndDistToEdge
00141  
00142   // Calculate Traces
00143   // ================
00144 
00145   double invSqrt2 = pow(1./2.,0.5);
00146   double m[2]; double c[2]; double Sign;
00147   double Coord[3]; double Trace[3];
00148   bool TraceSuccess;
00149 
00150   // For vertex
00151   if(vtxdirz!=0.) {
00152     m[0]=vtxdiru/vtxdirz;
00153     m[1]=vtxdirv/vtxdirz;
00154     c[0]=vtxu-(m[0]*vtxz);
00155     c[1]=vtxv-(m[1]*vtxz);
00156     Coord[0] = vtxu; Coord[1] = vtxv; Coord[2] = vtxz;
00157     Trace[0] = 1.e99; Trace[1] = 1.e99; Trace[2] = 1.e99;
00158 
00159     TraceSuccess = CalculateTrace(m,c,Coord,Trace);
00160     if(TraceSuccess==true) {
00161       Sign=1.;
00162       if(fabs(Coord[0])>4. || fabs(Coord[1])>4. || fabs(invSqrt2*(Coord[0]+Coord[1]))>4. 
00163          || fabs(invSqrt2*(Coord[0]-Coord[1]))>4.) {Sign=-1.;}
00164         
00165       MyTrack->VtxTrace=(Sign*pow(pow(Trace[0],2) + pow(Trace[1],2) + pow(Trace[2],2), 0.5));
00166       MyTrack->VtxTraceZ=(Sign*fabs(Trace[2]));
00167     }
00168   }
00169     
00170   // For end
00171   if(enddirz!=0.) {
00172     m[0] = enddiru/enddirz;
00173     m[1] = enddirv/enddirz;
00174     c[0] = endu-(m[0]*endz);
00175     c[1] = endv-(m[1]*endz);
00176     Coord[0] = endu; Coord[1] = endv; Coord[2] = endz;
00177     Trace[0] = 1.e99; Trace[1]=1.e99; Trace[2]=1.e99;
00178     
00179     TraceSuccess = CalculateTrace(m,c,Coord,Trace);
00180     if(TraceSuccess==true) {
00181       Sign=1.;
00182       if(fabs(Coord[0])>4. || fabs(Coord[1])>4. || fabs(invSqrt2*(Coord[0]+Coord[1]))>4. 
00183          || fabs(invSqrt2*(Coord[0]-Coord[1]))>4.) {Sign=-1.;}
00184 
00185       MyTrack->EndTrace=(Sign*pow(pow(Trace[0],2) + pow(Trace[1],2) + pow(Trace[2],2), 0.5));
00186       MyTrack->EndTraceZ=(Sign*fabs(Trace[2]));
00187     }
00188   }
00189 
00191 
00192 
00193   // Track containment variables
00194   // ===========================
00195 
00196   MSG("AtmosCalculator",Msg::kVerbose) << "    Track Containment Variables" << endl;
00197   Int_t nSingleEnds(0),nDoubleEnds(0),nDigits(0);
00198   Double_t dT,dU,dV,dQ,vtxdUmax(0.),vtxdVmax(0.),enddUmax(0.),enddVmax(0.);
00199   Double_t vtxUwt2(0.),vtxUwt(0.),vtxUw(0.);
00200   Double_t vtxVwt2(0.),vtxVwt(0.),vtxVw(0.);
00201   Double_t endUwt2(0.),endUwt(0.),endUw(0.);
00202   Double_t endVwt2(0.),endVwt(0.),endVw(0.);
00203   Double_t totTrkPH(0.),totTrkAssocPH(0.),totPlnPH(0.);
00204   Double_t trkAssocPH[500],trkPH[500],plnPH[500];
00205   Int_t bstrp[500],estrp[500];
00206   Int_t nstrps[500],ntrkstrps[500];
00207   Int_t planeview[500];
00208   memset(trkAssocPH, 0, 500*sizeof(Double_t));
00209   memset(trkPH, 0, 500*sizeof(Double_t));
00210   memset(plnPH, 0, 500*sizeof(Double_t));
00211   memset(bstrp, -1, 500*sizeof(Int_t));
00212   memset(estrp, -1, 500*sizeof(Int_t));
00213   memset(nstrps, 0, 500*sizeof(Int_t));
00214   memset(ntrkstrps, 0, 500*sizeof(Int_t));
00215   memset(planeview, -1, 500*sizeof(Int_t));
00216   int thispln;
00217 
00218   int ListEnd = StripList->GetLast()+1;
00219   for(int strp=0; strp<ListEnd; ++strp)
00220     {
00221       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00222       thispln=strip.Plane;
00223       if(thispln>-1 && thispln<486)
00224         {
00225 
00226           ++nstrps[thispln];
00227  
00228           if(strip.View==0)
00229             {
00230               if(planeview[thispln]<0) planeview[thispln]=0;
00231             }
00232 
00233           if(strip.View==1)
00234             {
00235               if(planeview[thispln]<0) planeview[thispln]=1;
00236             }
00237 
00238           if((strip.Trk&index)==index) 
00239             {
00240               if(bstrp[thispln]<0||strip.Strip<bstrp[thispln]) bstrp[thispln]=strip.Strip;
00241               if(estrp[thispln]<0||strip.Strip>estrp[thispln]) estrp[thispln]=strip.Strip;
00242               if(strip.Ndigits==1) ++nSingleEnds; if(strip.Ndigits==2) ++nDoubleEnds;
00243               nDigits += strip.Ndigits;
00244               ++ntrkstrps[thispln];
00245             }
00246         }
00247     }
00248 
00249   for(int strp=0; strp<ListEnd; ++strp)
00250     {
00251       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00252       thispln=strip.Plane;
00253       dQ=(strip.Sigcorr[0]+strip.Sigcorr[1])/65.0; // hack a rough conversion 
00254                                                    // from SigCorr -> PEcorr
00255 
00256       // Record pulse height information
00257       if(thispln>-1 && thispln<486)
00258         {
00259           // strips associated with track (+/-1 strip window...)
00260           if( bstrp[thispln]>-1 && estrp[thispln]>-1 && strip.Strip>bstrp[thispln]-2 && strip.Strip<estrp[thispln]+2)
00261             {
00262               trkAssocPH[thispln]+=dQ;
00263               totTrkAssocPH+=dQ;
00264             }
00265 
00266           // strips on track (+/-0 strip window...)
00267           if( bstrp[thispln]>-1 && estrp[thispln]>-1 && strip.Strip>bstrp[thispln]-1 && strip.Strip<estrp[thispln]+1)
00268             {
00269               trkPH[thispln]+=dQ;
00270               totTrkPH+=dQ;
00271             }
00272 
00273           // total pulse height
00274           plnPH[thispln]+=dQ;
00275           totPlnPH+=dQ;
00276         }
00277       
00278       // Handle the Track containment information
00279       // (could require the strips to NOT be on the track here...)
00280       if( abs(vtxpln-thispln)<5 && !(vtxpln<249 && thispln>249) && !(vtxpln>249 && thispln<249) )
00281         {//strip must be within four of the vtx and not on the other side of the SM gap - could do this more neatly with Z
00282           if(strip.View==0)
00283             { 
00284               dU=strip.T-vtxu; vtxdUmax=(dU>vtxdUmax)?dU:vtxdUmax;
00285               dT=(strip.T)-(vtxu+dir*vtxdiru*(strip.Z-vtxz));
00286               vtxUwt2+=dQ*dT*dT; vtxUwt+=dQ*dT; vtxUw+=dQ;
00287             }
00288           if(strip.View==1)
00289             { 
00290               dV=strip.T-vtxv; vtxdVmax=(dV>vtxdVmax)?dV:vtxdVmax;
00291               dT=(strip.T)-(vtxv+dir*vtxdirv*(strip.Z-vtxz));
00292               vtxVwt2+=dQ*dT*dT; vtxVwt+=dQ*dT; vtxVw+=dQ;
00293             }
00294         }
00295 
00296       if( abs(endpln-thispln)<5 && !(endpln<249 && thispln>249) && !(endpln>249 && thispln<249) )
00297         {//strip must be within four of the end and not on the other side of the SM gap - could do this more neatly with Z
00298           if(strip.View==0)
00299             { 
00300               dU=strip.T-endu; enddUmax=(dU>enddUmax)?dU:enddUmax;
00301               dT=(strip.T)-(endu+dir*enddiru*(strip.Z-endz));
00302               endUwt2+=dQ*dT*dT; endUwt+=dQ*dT; endUw+=dQ;
00303             }
00304           if(strip.View==1)
00305             { 
00306               dV=strip.T-endv; enddVmax=(dV>enddVmax)?dV:enddVmax;
00307               dT=(strip.T)-(endv+dir*enddirv*(strip.Z-endz));
00308               endVwt2+=dQ*dT*dT; endVwt+=dQ*dT; endVw+=dQ;
00309             }
00310         }
00311     }
00312 
00313   // VtxRmax, EndRmax
00314   MyTrack->VtxRmax = pow( vtxdUmax*vtxdUmax + vtxdVmax*vtxdVmax, 0.5 );
00315   MyTrack->EndRmax = pow( enddUmax*enddUmax + enddVmax*enddVmax, 0.5 );
00316 
00317   // VtxUmean, VtxUwidth
00318   if(vtxUw>0.0)
00319     {
00320       MyTrack->VtxUmean=vtxUwt/vtxUw;
00321       if((vtxUwt2/vtxUw)>(vtxUwt/vtxUw)*(vtxUwt/vtxUw))
00322         {
00323           MyTrack->VtxUwidth=sqrt((vtxUwt2/vtxUw)-(vtxUwt/vtxUw)*(vtxUwt/vtxUw));
00324         }
00325     }
00326 
00327   // VtxVmean, VtxVwidth
00328   if(vtxVw>0.0)
00329     {
00330       MyTrack->VtxVmean=vtxVwt/vtxVw;
00331       if((vtxVwt2/vtxVw)>(vtxVwt/vtxVw)*(vtxVwt/vtxVw))
00332         {
00333           MyTrack->VtxVwidth = sqrt((vtxVwt2/vtxVw)-(vtxVwt/vtxVw)*(vtxVwt/vtxVw));
00334         }
00335     }
00336 
00337 
00338   // EndUmean, EndUwidth
00339   if(endUw>0.0)
00340     {
00341       MyTrack->EndUmean=endUwt/endUw;
00342       if((endUwt2/endUw)>(endUwt/endUw)*(endUwt/endUw))
00343         {
00344           MyTrack->EndUwidth=sqrt((endUwt2/endUw)-(endUwt/endUw)*(endUwt/endUw));
00345         }
00346     }
00347 
00348   // EndVmean, EndVwidth
00349   if(endVw>0.0)
00350     {
00351       MyTrack->EndVmean=endVwt/endVw;  
00352       if((endVwt2/endVw)>(endVwt/endVw)*(endVwt/endVw))
00353         {
00354           MyTrack->EndVwidth=sqrt((endVwt2/endVw)-(endVwt/endVw)*(endVwt/endVw));
00355         }
00356     }
00357 
00358   // Pulse height variables
00359   MSG("AtmosCalculator",Msg::kVerbose) << "    Pulse Height Variables" << endl;
00360   double vtxqmax=0.0,endqmax=0.0;
00361   for(int ipln=0; ipln<500; ++ipln)
00362     {
00363       if( abs(vtxpln-ipln)<5 && !(vtxpln<249 && ipln>249) && !(vtxpln>249 && ipln<249) )
00364         {
00365           if( plnPH[ipln]>vtxqmax ) vtxqmax=plnPH[ipln];
00366         }
00367       if( abs(endpln-ipln)<5 && !(endpln<249 && ipln>249) && !(endpln>249 && ipln<249) )
00368         {
00369           if( plnPH[ipln]>endqmax ) endqmax=plnPH[ipln];
00370         }
00371     }
00372 
00373   MyTrack->VtxQmax = vtxqmax;          // VtxQmax
00374   MyTrack->EndQmax = endqmax;          // EndQmax
00375   MyTrack->TrkPH = totTrkPH;           // TrkPH
00376   MyTrack->AssocTrkPH = totTrkAssocPH; // AssocTrkPH
00377   if(totPlnPH>0.0) MyTrack->AssocTrkPHfrac = totTrkAssocPH/totPlnPH;
00378 
00379   // Additional plane/strip variables;
00380   MSG("AtmosCalculator",Msg::kVerbose) << "    Additional Plane/Strip Variables" << endl;
00381   int uPlanes(0),vPlanes(0);
00382   int nPlanesTrkOnly(0),nPlanesTrkGaps(0);
00383   for(int ipln=0; ipln<500; ++ipln)
00384     {
00385       if(ipln>=minpln && ipln<=maxpln)
00386         { 
00387           if(ntrkstrps[ipln]>0)
00388             {
00389               if(planeview[ipln]==0) ++uPlanes;
00390               if(planeview[ipln]==1) ++vPlanes;
00391               if(ntrkstrps[ipln]==nstrps[ipln]) ++nPlanesTrkOnly;
00392             }
00393           else
00394             {
00395               ++nPlanesTrkGaps;
00396             }
00397         }
00398     }
00399 
00400   MyTrack->Ndigits = nDigits;                 // Ndigits
00401   MyTrack->NstripsSingleEnded = nSingleEnds;  // NstripsSingleEnded
00402   MyTrack->NstripsDoubleEnded = nDoubleEnds;  // NstripsDoubleEnded
00403   MyTrack->NplanesTrackOnly = nPlanesTrkOnly; // NplanesTrackOnly
00404   MyTrack->NplanesTrackGaps = nPlanesTrkGaps; // NplanesTrackGaps
00405   MyTrack->NplanesUview = uPlanes;            // NplanesUview
00406   MyTrack->NplanesVview = vPlanes;            // NplanesVview
00407   if(uPlanes+vPlanes>0) MyTrack->UVassymetry = abs(uPlanes-vPlanes)/(0.5*(uPlanes+vPlanes));                                 // UVassymetry
00408   
00409   // Number of track-like planes
00410   Int_t ntrkplns(0);
00411   for(int ipln=0; ipln<500; ++ipln)
00412     {
00413       if( plnPH[ipln]>0.0 && plnPH[ipln]<80.0 ) 
00414         {
00415           if(trkAssocPH[ipln]/plnPH[ipln]>0.8) ++ntrkplns;
00416         }
00417     }
00418 
00419   MyTrack->TrkLikePlanes = ntrkplns; // TrkLikePlanes
00420 
00422   
00423   // Linear Fits
00424   // ===========
00425 
00426   //double dstrp = 0.041/sqrt(12.0);
00427   double dstrpSQ = 0.041*0.041/12.0;
00428 
00429   double u,v,x,y,z,w,du,dv;
00430   double mu,mv,cu,cv,denom;
00431   double Uchi2,Vchi2;
00432   int Undf,Vndf;
00433   int Upts,Vpts;
00434 
00435   // Linear Fit 
00436   // (Separate Linear Fit In Each View)
00437   MSG("AtmosCalculator",Msg::kVerbose) << "    Linear Fit " << endl;
00438   double Uw(0.0), Uwz(0.0), Uwzz(0.0), Uwu(0.0), Uwuu(0.0), Uwzu(0.0);
00439   double Vw(0.0), Vwz(0.0), Vwzz(0.0), Vwv(0.0), Vwvv(0.0), Vwzv(0.0);
00440   
00441   Upts=0; Vpts=0; 
00442 
00443   for(int strp=0; strp<ListEnd; ++strp)
00444     {
00445       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00446       thispln=strip.Plane;
00447       if(thispln>-1 && thispln<486)
00448         {
00449           if((strip.Trk&index)==index) 
00450             {
00451               z=strip.Z; 
00452               w=1.0;
00453 
00454               if(strip.View==0)
00455                 {
00456                   u=strip.T;
00457                   Uw+=w;
00458                   Uwu+=w*u; Uwz+=w*z; Uwuu+=w*u*u;
00459                   Uwzu+=w*u*z; Uwzz+=w*z*z;
00460                   ++Upts;
00461                 }
00462               
00463               if(strip.View==1)
00464                 {
00465                   v=strip.T;
00466                   Vw+=w;
00467                   Vwv+=w*v; Vwz+=w*z; Vwvv+=w*v*v;
00468                   Vwzv+=w*v*z; Vwzz+=w*z*z;
00469                   ++Vpts;
00470                 }
00471               
00472             }
00473         }
00474     }
00475 
00476   mu=0.0; mv=0.0; cu=0; cv=0;
00477 
00478   if(Upts>1 && Vpts>1)
00479     {
00480       denom = Uw*Uwzz-Uwz*Uwz;
00481       if(fabs(denom)>0.0) { mu=(Uw*Uwzu-Uwz*Uwu)/denom; cu=(Uwu*Uwzz-Uwz*Uwzu)/denom; }
00482       denom = Vw*Vwzz-Vwz*Vwz;
00483       if(fabs(denom)>0.0) { mv=(Vw*Vwzv-Vwz*Vwv)/denom; cv=(Vwv*Vwzz-Vwz*Vwzv)/denom; }
00484       // (calculate error properly now)
00485       // err=(Uwuu+mu*mu*Uwzz+Uw*cu*cu+2.0*(mu*cu*Uwz-mu*Uwzu-cu*Uwu))+(Vwvv+mv*mv*Vwzz+Vw*cv*cv+2.0*(mv*cv*Vwz-mv*Vwzv-cv*Vwv));
00486       // chi2=err;
00487       // ndf=Upts+Vpts-4;
00488     }
00489   
00490   double precou=mu; double precov=mv; double precoz=1.0;
00491   double precos=pow(precou*precou+precov*precov+1.0, -0.5);
00492   precou*=precos; precov*=precos; precoz*=precos;
00493   MyTrack->LinearDirCosU=precou*dir; // LinearDirCosU
00494   MyTrack->LinearDirCosV=precov*dir; // LinearDirCosV
00495   MyTrack->LinearDirCosZ=precoz*dir; // LinearDirCosZ
00496 
00497   // Unconstrained Linear Fit 
00498   //  (chi2 about linear fit to track)
00499   MSG("AtmosCalculator",Msg::kVerbose) << "    Linear Fit (unconstrained) " << endl;
00500 
00501   mu = 0.0;
00502   if( MyTrack->LinearDirCosZ!=0 )
00503     mu = MyTrack->LinearDirCosU/MyTrack->LinearDirCosZ; 
00504   cu = vtxu-mu*vtxz;
00505 
00506   if( MyTrack->LinearDirCosZ!=0 )
00507     mv = MyTrack->LinearDirCosV/MyTrack->LinearDirCosZ; 
00508   cv = vtxv-mv*vtxz;
00509 
00510   Uchi2=0.0; Vchi2=0.0; 
00511   Undf=0; Vndf=0;
00512 
00513   for(int strp=0; strp<ListEnd; ++strp)
00514     {
00515       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00516       thispln=strip.Plane;
00517       if(thispln>-1 && thispln<486)
00518         {
00519           if((strip.Trk&index)==index) 
00520             {
00521               if( strip.View==0 )
00522                 {
00523                   du=strip.T-(mu*strip.Z+cu);
00524                   Uchi2+=(du*du)/dstrpSQ; //(dstrp*dstrp); 
00525                   ++Undf;
00526                 }
00527 
00528               if( strip.View==1 )
00529                 {
00530                   dv=strip.T-(mv*strip.Z+cv);
00531                   Vchi2+=dv*dv/dstrpSQ; //(dstrp*dstrp); 
00532                   ++Vndf;
00533                 }
00534 
00535             }
00536         }
00537     }
00538   
00539   MyTrack->LinearDirFitChisqU = Uchi2;       // LinearDirFitChisqU
00540   MyTrack->LinearDirFitNdfU   = Undf;        // LinearDirFitNdfU
00541   MyTrack->LinearDirFitChisqV = Vchi2;       // LinearDirFitChisqV
00542   MyTrack->LinearDirFitNdfV   = Vndf;        // LinearDirFitNdfV
00543   MyTrack->LinearDirFitChisq  = Uchi2+Vchi2; // LinearDirFitChisq
00544   MyTrack->LinearDirFitNdf    = Undf+Vndf;   // LinearDirFitNdf
00545 
00546   // Constrained Linear Fits 
00547   //  (chi2 about linear fit using reconstructed track direction)
00548   MSG("AtmosCalculator",Msg::kVerbose) << "    Linear Fit (vertex) " << endl;
00549 
00550   mu = 0.0;
00551   if( MyTrack->VtxDirCosZ!=0 )
00552     mu = MyTrack->VtxDirCosU/MyTrack->VtxDirCosZ;
00553   cu = vtxu-mu*vtxz;
00554 
00555   mv = 0.0;
00556   if( MyTrack->VtxDirCosZ!=0 ) 
00557     mv = MyTrack->VtxDirCosV/MyTrack->VtxDirCosZ; 
00558   cv = vtxv-mv*vtxz;
00559 
00560   Uchi2=0.0; Vchi2=0.0; 
00561   Undf=0; Vndf=0;
00562 
00563   for(int strp=0; strp<ListEnd; ++strp)
00564     {
00565       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00566       thispln=strip.Plane;
00567       if(thispln>-1 && thispln<486)
00568         {
00569           if((strip.Trk&index)==index) 
00570             {
00571               if( strip.View==0 )
00572                 {
00573                   du=strip.T-(mu*strip.Z+cu);
00574                   Uchi2+=(du*du)/dstrpSQ; //(dstrp*dstrp); 
00575                   ++Undf;
00576                 }
00577 
00578               if( strip.View==1 )
00579                 {
00580                   dv=strip.T-(mv*strip.Z+cv);
00581                   Vchi2+=dv*dv/dstrpSQ; //(dstrp*dstrp); 
00582                   ++Vndf;
00583                 }
00584 
00585             }
00586         }
00587     }
00588 
00589   MyTrack->VtxLinearDirFitChisqU = Uchi2; // VtxLinearDirFitChisqU
00590   MyTrack->VtxLinearDirFitNdfU   = Undf;  // VtxLinearDirFitNdfU
00591   MyTrack->VtxLinearDirFitChisqV = Vchi2; // VtxLinearDirFitChisqV
00592   MyTrack->VtxLinearDirFitNdfV   = Vndf;  // VtxLinearDirFitNdfV
00593 
00594   mu = 0.0;
00595   if( MyTrack->EndDirCosZ!=0 )
00596     mu = MyTrack->EndDirCosU/MyTrack->EndDirCosZ;
00597   cu = endu-mu*endz;
00598 
00599   mv = 0.0;
00600   if( MyTrack->EndDirCosZ!=0 ) 
00601     mv = MyTrack->EndDirCosV/MyTrack->EndDirCosZ; 
00602   cv = endv-mv*endz;
00603 
00604   Uchi2=0.0; Vchi2=0.0; 
00605   Undf=0; Vndf=0;
00606 
00607   for(int strp=0; strp<ListEnd; ++strp)
00608     {
00609       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00610       thispln=strip.Plane;
00611       if(thispln>-1 && thispln<486)
00612         {
00613           if((strip.Trk&index)==index) 
00614             {
00615               if( strip.View==0 )
00616                 {
00617                   du=strip.T-(mu*strip.Z+cu);
00618                   Uchi2+=(du*du)/dstrpSQ; // (dstrp*dstrp); 
00619                   ++Undf;
00620                 }
00621 
00622               if( strip.View==1 )
00623                 {
00624                   dv=strip.T-(mv*strip.Z+cv);
00625                   Vchi2+=dv*dv/dstrpSQ;   // (dstrp*dstrp); 
00626                   ++Vndf;
00627                 }
00628 
00629             }
00630         }
00631     }
00632 
00633   MyTrack->EndLinearDirFitChisqU = Uchi2; // EndLinearDirFitChisqU
00634   MyTrack->EndLinearDirFitNdfU   = Undf;  // EndLinearDirFitNdfU
00635   MyTrack->EndLinearDirFitChisqV = Vchi2; // EndLinearDirFitChisqV
00636   MyTrack->EndLinearDirFitNdfV   = Vndf;  // EndLinearDirFitNdfV
00637   
00638   // Zero Curvature Fit
00639   //  (join up the track vertex and end positions - 
00640   //   calculate the RMS of the strips about this line)
00641   MSG("AtmosCalculator",Msg::kVerbose) << "    Zero Curvature Fit" << endl;
00642 
00643   mu = (endu-vtxu)/(endz-vtxz); 
00644   cu = vtxu-mu*vtxz;
00645   mv = (endv-vtxv)/(endz-vtxz);
00646   cv = vtxv-mv*vtxz;
00647  
00648   Uchi2=0.0; Vchi2=0.0; 
00649   Undf=0; Vndf=0;
00650 
00651   for(int strp=0; strp<ListEnd; ++strp)
00652     {
00653       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00654       thispln=strip.Plane;
00655       if(thispln>-1 && thispln<486)
00656         {
00657           if((strip.Trk&index)==index) 
00658             {
00659               if( strip.View==0 )
00660                 {
00661                   du=strip.T-(mu*strip.Z+cu);
00662                   Uchi2+=(du*du)/dstrpSQ; //(dstrp*dstrp); 
00663                   ++Undf;
00664                 }
00665 
00666               if( strip.View==1 )
00667                 {
00668                   dv=strip.T-(mv*strip.Z+cv);
00669                   Vchi2+=dv*dv/dstrpSQ; //(dstrp*dstrp); 
00670                   ++Vndf;
00671                 }
00672 
00673             }
00674         }
00675     }
00676 
00677   MyTrack->ZeroCurveChi2 = Uchi2 + Vchi2;
00678   MyTrack->ZeroCurveNdf = Undf + Vndf;
00679 
00681 
00682   // Digit Containment Variables
00683   // ===========================
00684 
00685   MSG("AtmosCalculator",Msg::kVerbose) << "    Digit Containment Variables" << endl;
00686   double wPlaneChargeSigCorr[500],wPlaneMeanTPos[500],wPlaneZPos[500];
00687   memset(wPlaneChargeSigCorr, 0, 500*sizeof(double));
00688   memset(wPlaneMeanTPos, 0, 500*sizeof(double));
00689   memset(wPlaneZPos, 0, 500*sizeof(double));
00690   double stripU(0.0),stripV(0.0);
00691   double nonfidTrkQ(0.0),totTrkQ(0.0);
00692   bool UCsides(false),UCends(false);
00693   double Xdis(0.3);
00694   int Xpln(5);
00695   double VtxZdigit[3]={vtxx,vtxy,MyTrack->VtxPlane};
00696   double EndZdigit[3]={endx,endy,MyTrack->EndPlane};
00697   double VtxRdigit[3]={vtxx,vtxy,MyTrack->VtxPlane};
00698   double EndRdigit[3]={endx,endy,MyTrack->EndPlane};
00699   int fPlaneMin = 500,fPlaneMax = 0;
00700  
00701   for(int strp=0; strp<ListEnd; ++strp)
00702     {
00703       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00704       wPlaneZPos[strip.Plane] = strip.Z;
00705       wPlaneChargeSigCorr[strip.Plane] += (strip.Sigcorr[0]+strip.Sigcorr[1]);
00706       wPlaneMeanTPos[strip.Plane] += (strip.T*(strip.Sigcorr[0]+strip.Sigcorr[1]));
00707       if(strip.Plane>fPlaneMax) fPlaneMax = strip.Plane;
00708       if(strip.Plane<fPlaneMin) fPlaneMin = strip.Plane;
00709     }
00710 
00711   for(int ipln=fPlaneMin; ipln<=fPlaneMax; ++ipln)
00712     {
00713       if(wPlaneChargeSigCorr[ipln]>0.0)
00714         {
00715           wPlaneMeanTPos[ipln]/=wPlaneChargeSigCorr[ipln];
00716         }
00717     }
00718  
00719   for(int strp=0; strp<ListEnd; ++strp)
00720     {
00721       AtmosStrip& strip = *((AtmosStrip*)(StripList->At(strp)));
00722 
00723       const int plane = strip.Plane;
00724       const double tpos = strip.T;
00725       
00726       double nz[2]={-1.0, -1.0};
00727       double ntpos[2]={-1.0, -1.0};
00728       double stp_guess_z=wPlaneZPos[plane];
00729       double stp_guess_tpos=-9999.0;
00730 
00731       if(plane>0 && plane<500 && wPlaneChargeSigCorr[plane-1]>0.0 && wPlaneChargeSigCorr[plane+1]>0.0)
00732         {
00733           nz[0] = wPlaneZPos[plane-1];
00734           nz[1] = wPlaneZPos[plane+1];
00735           ntpos[0] = wPlaneMeanTPos[plane-1];
00736           ntpos[1] = wPlaneMeanTPos[plane+1];
00737         }
00738       else if(plane>2  && plane<500 && wPlaneChargeSigCorr[plane-1]>0.0 && wPlaneChargeSigCorr[plane-3]>0.0)
00739         {
00740           nz[0] = wPlaneZPos[plane-1];
00741           nz[1] = wPlaneZPos[plane-3];
00742           ntpos[0] = wPlaneMeanTPos[plane-1];
00743           ntpos[1] = wPlaneMeanTPos[plane-3];
00744         }
00745       else if(plane<497  && plane>0 && wPlaneChargeSigCorr[plane+1]>0.0 && wPlaneChargeSigCorr[plane+3]>0.0)
00746         {
00747           nz[0] = wPlaneZPos[plane+1];
00748           nz[1] = wPlaneZPos[plane+3];
00749           ntpos[0] = wPlaneMeanTPos[plane+1];
00750           ntpos[1] = wPlaneMeanTPos[plane+3];      
00751         }
00752       else if(plane>1 && wPlaneChargeSigCorr[plane-1]>0.0)
00753         {stp_guess_tpos = wPlaneMeanTPos[plane-1];}
00754       else if(plane<497 && wPlaneChargeSigCorr[plane+1]>0.0)
00755         {stp_guess_tpos = wPlaneMeanTPos[plane+1];}
00756       else{/*lone hit*/ } 
00757        
00758       //
00759       //Extrapolate position
00760       //
00761       if(stp_guess_tpos<-1000.0 &&  nz[0]>0.0)
00762         {
00763           stp_guess_tpos =((ntpos[0]-ntpos[1])/(nz[0]-nz[1]))*stp_guess_z + 
00764             ((ntpos[1]*nz[0] - ntpos[0]*nz[1]) / (nz[0]-nz[1]));
00765         }
00766     
00767       //calc XY coordinates
00768       if(stp_guess_tpos>-1000.0)
00769         {
00770           u = (strip.View==0)?tpos:stp_guess_tpos;
00771           v = (strip.View==1)?tpos:stp_guess_tpos;
00772           y = (u+v)*0.707106781;
00773           x = (u-v)*0.707106781;
00774           
00775           /*
00776             Want to record pretrackPH, posttrackPH, first hit before vtx, last hit after end
00777           */
00778           if((strip.Trk&index)==index)
00779             {
00780               stripU = (strip.View==0)?strip.T:strip.L;
00781               stripV = (strip.View==1)?strip.T:strip.L;
00782               UCsides = ( (fabs(stripU)>(4-Xdis)) || (fabs(stripV)>(4-Xdis)) || (((fabs(stripU)+fabs(stripV))/sqrt(2.0))>(4-Xdis)) );
00783               UCends  = ( (strip.Plane<=Xpln) || ((strip.Plane>=(249-Xpln)) && (strip.Plane<=(249+Xpln))) || (strip.Plane>=(486-Xpln)) );
00784               if( UCsides || UCends ) nonfidTrkQ+=strip.Sigcorr[0]+strip.Sigcorr[1];
00785               totTrkQ+=strip.Sigcorr[0]+strip.Sigcorr[1];
00786             }
00787 
00788           // don't bother using very low PH hits
00789           if( (strip.Sigcorr[0]+strip.Sigcorr[1])<100.0 ) continue; 
00790 
00791           if(MyTrack->VtxPlane<MyTrack->EndPlane)
00792             {
00793               if(strip.Plane<=MyTrack->VtxPlane)
00794                 {
00795                   if(VtxZdigit[2]>double(strip.Plane))
00796                     {
00797                       VtxZdigit[0]=x;
00798                       VtxZdigit[1]=y;
00799                       VtxZdigit[2]=double(strip.Plane);
00800                     }
00801                   if( ( pow(VtxRdigit[0]*VtxRdigit[0] + VtxRdigit[1]*VtxRdigit[1] ,0.5)<pow(u*u + v*v ,0.5)) )
00802                     {
00803                       VtxRdigit[0]=x;
00804                       VtxRdigit[1]=y;
00805                       VtxRdigit[2]=double(strip.Plane);
00806                     }
00807                 }
00808               
00809               if(strip.Plane>=MyTrack->EndPlane)
00810                 {
00811                   if(EndZdigit[2]<double(strip.Plane))
00812                     {
00813                       EndZdigit[0]=x;
00814                       EndZdigit[1]=y;
00815                       EndZdigit[2]=double(strip.Plane);
00816                     }
00817                   if(pow(EndRdigit[0]*EndRdigit[0] + EndRdigit[1]*EndRdigit[1] ,0.5)<pow(u*u + v*v ,0.5))
00818                     {
00819                       EndRdigit[0]=x;
00820                       EndRdigit[1]=y;
00821                       EndRdigit[2]=double(strip.Plane);
00822                     }
00823                 }
00824               
00825             }
00826           else//neglect the equality case
00827             {
00828               if(strip.Plane>=MyTrack->VtxPlane)
00829                 {
00830                   if(VtxZdigit[2]<double(strip.Plane))
00831                     {
00832                       VtxZdigit[0]=x;
00833                       VtxZdigit[1]=y;
00834                       VtxZdigit[2]=double(strip.Plane);
00835                     }
00836                   if(pow(VtxRdigit[0]*VtxRdigit[0] + VtxRdigit[1]*VtxRdigit[1] ,0.5)<pow(u*u + v*v ,0.5))
00837                     {
00838                       VtxRdigit[0]=x;
00839                       VtxRdigit[1]=y;
00840                       VtxRdigit[2]=double(strip.Plane);
00841                     }
00842                 }
00843               if(strip.Plane<=MyTrack->EndPlane)
00844                 {
00845                   //Trace stuff in here
00846                   if(EndZdigit[2]>double(strip.Plane))
00847                     {
00848                       EndZdigit[0]=x;
00849                       EndZdigit[1]=y;
00850                       EndZdigit[2]=double(strip.Plane);
00851                     }
00852                   if(pow(EndRdigit[0]*EndRdigit[0] + EndRdigit[1]*EndRdigit[1] ,0.5)<pow(u*u + v*v ,0.5))
00853                     {
00854                       EndRdigit[0]=x;
00855                       EndRdigit[1]=y;
00856                       EndRdigit[2]=double(strip.Plane);
00857                     }
00858                 }
00859             }
00860         }           
00861     }
00862   //here decide on the values to use...
00863   MyTrack->VtxDistToEdgeDigits = 4.0 - pow(VtxRdigit[0]*VtxRdigit[0] + VtxRdigit[1]*VtxRdigit[1] ,0.5); // VtxDistToEdgeDigits
00864   MyTrack->VtxPlaneDigits = int(VtxZdigit[2]); // VtxPlaneDigits
00865 
00866   MyTrack->EndDistToEdgeDigits = 4.0 - pow(EndRdigit[0]*EndRdigit[0] + EndRdigit[1]*EndRdigit[1] ,0.5); // EndDistToEdgeDigits
00867   MyTrack->EndPlaneDigits = int(EndZdigit[2]); // EndPlaneDigits
00868 
00869   //calc the amount of track outside the fiducial volume
00870   if( totTrkQ>0.0 ) MyTrack->NonFidFrac=nonfidTrkQ/totTrkQ; // NonFidFrac
00871   return;
00872 }


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1