#include <AtmosCalculator.h>
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) |
|
|
Definition at line 17 of file AtmosCalculator.cxx. References MSG. 00018 {
00019 MSG("AtmosCalculator",Msg::kDebug) << " AtmosCalculator::AtmosCalculator() " << endl;
00020 }
|
|
|
Definition at line 22 of file AtmosCalculator.cxx. References MSG. 00023 {
00024 MSG("AtmosCalculator",Msg::kDebug) << " AtmosCalculator::~AtmosCalculator() " << endl;
00025 }
|
|
||||||||||||||||||||
|
Definition at line 1135 of file AtmosCalculator.cxx. References DetectorSides(), 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 }
|
|
||||||||||||||||||||
|
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 }
|
|
||||||||||||
|
Definition at line 27 of file AtmosCalculator.cxx. References AtmosReco::DoubleEndedCharge, AtmosReco::DoubleEndedChargeUview, AtmosReco::DoubleEndedChargeVview, AtmosReco::EastSideCharge, 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 }
|
|
||||||||||||
|
Definition at line 874 of file AtmosCalculator.cxx. References abs(), AtmosShower::AssocShwPH, AtmosShower::AssocShwPHfrac, ShowerMOI::EigenValues, ShowerMOI::EigenVectors, AtmosShower::Index, 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 }
|
|
||||||||||||
|
Definition at line 90 of file AtmosCalculator.cxx. References abs(), 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, AtmosStrip::L, AtmosTrack::LinearDirCosU, AtmosTrack::LinearDirCosV, AtmosTrack::LinearDirCosZ, AtmosTrack::LinearDirFitChisq, AtmosTrack::LinearDirFitChisqU, AtmosTrack::LinearDirFitChisqV, AtmosTrack::LinearDirFitNdf, AtmosTrack::LinearDirFitNdfU, AtmosTrack::LinearDirFitNdfV, AtmosTrack::MaxPlaneNumber, AtmosTrack::MinPlaneNumber, MSG, AtmosTrack::Ndigits, AtmosStrip::Ndigits, AtmosTrack::NonFidFrac, AtmosTrack::NplanesTrackGaps, AtmosTrack::NplanesTrackOnly, AtmosTrack::NplanesUview, AtmosTrack::NplanesVview, AtmosTrack::NstripsDoubleEnded, AtmosTrack::NstripsSingleEnded, AtmosStrip::Plane, AtmosStrip::Sigcorr, 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 }
|
1.3.9.1