#include <DmxUtilities.h>
Public Member Functions | |
| DmxUtilities () | |
| virtual | ~DmxUtilities () |
| void | FillPlaneArray (DmxStatus *status, CandDeMuxDigitListHandle &cdlh, AlgConfig &acd) |
| Float_t | CheckForXTalk (Int_t pixel, Float_t *planeSet) |
| void | FillHitPixels (CandDeMuxDigitHandleItr crdhi, Int_t chip, Float_t *eastSet, Float_t *westSet) |
| Int_t | GetDigitPixel (Int_t channel) |
| Int_t | FindVertexPlane (DmxPlaneItr &planeItr) |
| Int_t | FindEndPlane (DmxPlaneItr &planeItr) |
| Int_t | FindMuonStartPlane (DmxPlaneItr &planeItr, Float_t peGainConversion) |
| Int_t | FindBeamMuonStartPlane (DmxPlaneItr &planeItr, Float_t peGainConversion) |
| Bool_t | IsValidFit (DmxPlaneItr &planeItr, Double_t a1, Double_t a2, Double_t a3, Double_t a4, Float_t offset) |
| Bool_t | IsValidFit (DmxPlaneItr &planeItr, Double_t a1, Double_t a2, Double_t a3, Double_t a4) |
| Int_t | FindPlanesOffFit (DmxPlaneItr &planeItr, Int_t strayCut) |
| Bool_t | IsOverlappingMultiple (DmxPlaneItr &planeItr, Float_t vertexZ, Float_t slopeCutOff, Float_t interceptCutOff, Float_t peakCutOff) |
| Int_t | CheckFit (DmxPlaneItr &planeItr) |
| Float_t | CheckFitWithTiming (DmxPlaneItr &planeItr) |
| void | FindCubicFit (Double_t *x, Double_t *y, Double_t *weight, Int_t nPoints, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Double_t &chiSq) |
| void | FindQuadraticFit (Double_t *x, Double_t *y, Double_t *weight, Int_t nPoints, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &chiSq) |
| void | FindLinearFit (Double_t *x, Double_t *y, Double_t *weight, Int_t nPoints, Double_t &a1, Double_t &a2, Double_t &chiSq) |
|
|
Definition at line 42 of file DmxUtilities.cxx. 00043 {
00044 }
|
|
|
Definition at line 235 of file DmxUtilities.cxx. 00236 {
00237 //MSG("Dmx", Msg::kDebug) << "deleting DmxUtilities object " << endl;
00238 }
|
|
|
Definition at line 764 of file DmxUtilities.cxx. References DmxPlane::GetCoG(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlane::IsStray(), DmxPlane::SetStray(), and DmxPlane::SetStrips(). Referenced by AlgDeMuxGolden::RunAlg(). 00765 {
00766 DmxPlane *plane = 0;
00767
00768 //compare the best hyps in valid shower planes to the fit hyp for the amount of signal
00769 //mated in the plane. pick the one with the most mated signal, as long as the difference
00770 //is significant - like 0.2 or more. also only use the best hyps if they have more than
00771 //85% of the signal mated in them
00772 Float_t firstBestMated = 0.;
00773 Float_t firstSMated = 0.;
00774 Float_t firstTMated = 0.;
00775 Float_t firstMated = 0.;
00776 Float_t lastBestMated = 0.;
00777 Float_t lastSMated = 0.;
00778 Float_t lastTMated = 0.;
00779 Float_t lastMated = 0.;
00780
00781 Int_t ctr = 0;
00782 Int_t totPlanes = planeItr.SizeSelect();
00783 Int_t fixedPlanes = 0;
00784 planeItr.ResetFirst();
00785
00786 //make an array to hold the CoG positions of the planes
00787 Float_t firstPlaneCoG = -1.;
00788 Float_t firstPlaneBCoG = -1.;
00789 Float_t firstPlaneSCoG = -1.;
00790 Float_t firstPlaneTCoG = -1.;
00791 Float_t lastPlaneCoG = -1.;
00792 Float_t lastPlaneBCoG = -1.;
00793 Float_t lastPlaneSCoG = -1.;
00794 Float_t lastPlaneTCoG = -1.;
00795 Float_t firstPlaneSpace = 0.;
00796 Float_t lastPlaneSpace = 0.;
00797
00798 while(planeItr.IsValid()){
00799
00800 plane = planeItr.Ptr();
00801 if(ctr == 0){
00802 firstPlaneCoG = plane->GetCoG();
00803 firstPlaneSpace = plane->GetPlaneNumber();
00804 if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
00805
00806 firstPlaneBCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
00807 firstBestMated = dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
00808
00809 firstPlaneSCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
00810 firstSMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
00811
00812 firstPlaneTCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
00813 firstTMated = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
00814
00815 firstMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio();
00816 }
00817 else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
00818 firstPlaneBCoG = firstPlaneCoG;
00819 firstPlaneSCoG = firstPlaneCoG;
00820 firstPlaneTCoG = firstPlaneCoG;
00821 firstMated = 1.;
00822 firstBestMated = 1.;
00823 firstSMated = 1.;
00824 firstTMated = 1.;
00825 }
00826 }
00827 else if(ctr == 1) firstPlaneSpace -= plane->GetPlaneNumber();
00828 else if(ctr == totPlanes-2) lastPlaneSpace = plane->GetPlaneNumber();
00829 else if(ctr == totPlanes-1){
00830 lastPlaneCoG = plane->GetCoG();
00831 lastPlaneSpace -= plane->GetPlaneNumber();
00832 if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
00833
00834 lastPlaneBCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
00835 lastBestMated = dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
00836
00837 lastPlaneSCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
00838 lastSMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
00839
00840 lastPlaneTCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
00841 lastTMated = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
00842
00843 lastMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio();
00844 }
00845 else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
00846 lastPlaneBCoG = lastPlaneCoG;
00847 lastPlaneSCoG = lastPlaneCoG;
00848 lastPlaneTCoG = lastPlaneCoG;
00849 lastMated = 1.;
00850 lastBestMated = 1.;
00851 lastSMated = 1.;
00852 lastTMated = 1.;
00853 }
00854 }
00855 ++ctr;
00856 planeItr.Next();
00857 }
00858 planeItr.ResetFirst();
00859
00860 //which is the better hypothesis among the 3 best?
00861 bool useFirstB = false;
00862 bool useFirstS = false;
00863 bool useFirstT = false;
00864 bool useLastB = false;
00865 bool useLastS = false;
00866 bool useLastT = false;
00867
00868 if(firstBestMated>=firstSMated && firstBestMated>=firstTMated
00869 && TMath::Abs(firstBestMated-firstMated)>0.2) useFirstB = true;
00870 else if(firstSMated>firstBestMated && firstSMated>=firstTMated
00871 && TMath::Abs(firstSMated-firstMated)>0.2) useFirstS = true;
00872 else if(firstTMated>firstSMated && firstTMated>firstBestMated
00873 && TMath::Abs(firstTMated-firstMated)>0.2) useFirstT = true;
00874
00875 if(lastBestMated>=lastSMated && lastBestMated>=lastTMated
00876 && TMath::Abs(lastBestMated-lastMated)>0.2) useLastB = true;
00877 else if(lastSMated>lastBestMated && lastSMated>=lastTMated
00878 && TMath::Abs(lastSMated-lastMated)>0.2) useLastS = true;
00879 else if(lastTMated>lastSMated && lastTMated>lastBestMated
00880 && TMath::Abs(lastTMated-lastMated)>0.2) useLastT = true;
00881
00882
00883 // cout << "first plane: " << firstPlaneCoG << " "
00884 // << firstPlaneBCoG << " " << firstPlaneSCoG
00885 // << " " << firstPlaneTCoG << " " << firstMated
00886 // << " " << firstBestMated << " "
00887 // << firstSMated << " " << firstTMated << " "
00888 // << useFirstB << " " << useFirstS << " "
00889 // << useFirstT << " " << firstPlaneSpace << endl;
00890
00891 // cout << "last plane: " << lastPlaneCoG << " " << lastPlaneBCoG
00892 // << " " << lastPlaneSCoG << " " << lastPlaneTCoG << " " << lastMated
00893 // << " " << lastBestMated << " " << lastSMated << " " << lastTMated
00894 // << " " << useLastB << " " << useLastS << " " << useLastT
00895 // << " " << lastPlaneSpace << endl;
00896
00897 //the first and last planes in the view should be reset to one of the
00898 //best hypothesis
00899 //
00900 //1. the difference between the best hypothesis and the nearest neighboring
00901 // valid plane hypothesis is less than 0.25m/plane
00902 //
00903 // and
00904 //
00905 //2. the best hypothesis has a lot more signal mated
00906
00907 //plane iterator was previously set to the first plane in the set
00908
00909 if(useFirstB
00910 && TMath::Abs(firstPlaneSpace)*0.25>TMath::Abs(firstPlaneCoG-firstPlaneBCoG)){
00911 plane = planeItr.Ptr();
00912 if(plane->IsStray()){
00913 if(plane->GetPlaneType() == DmxPlaneTypes::kShower)
00914 dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("best");
00915 else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon) plane->SetStrips();
00916 plane->SetStray(false);
00917 ++fixedPlanes;
00918 }
00919 }
00920 else if(useFirstS
00921 && TMath::Abs(firstPlaneSpace)*0.25>TMath::Abs(firstPlaneCoG-firstPlaneSCoG)){
00922 plane = planeItr.Ptr();
00923 if(plane->IsStray()){
00924 if(plane->GetPlaneType() == DmxPlaneTypes::kShower)
00925 dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("second");
00926 plane->SetStray(false);
00927 ++fixedPlanes;
00928 }
00929 }
00930 else if(useFirstT
00931 && TMath::Abs(firstPlaneSpace)*0.25>TMath::Abs(firstPlaneCoG-firstPlaneTCoG)){
00932 plane = planeItr.Ptr();
00933 if(plane->IsStray()){
00934 if(plane->GetPlaneType() == DmxPlaneTypes::kShower)
00935 dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("third");
00936 plane->SetStray(false);
00937 ++fixedPlanes;
00938 }
00939 }
00940
00941 //set the iterator to the last plane
00942 planeItr.ResetLast();
00943
00944 if(useLastB && TMath::Abs(lastPlaneSpace)*0.25>TMath::Abs(lastPlaneCoG-lastPlaneBCoG)){
00945 plane = planeItr.Ptr();
00946 if(plane->IsStray()){
00947 if(plane->GetPlaneType() == DmxPlaneTypes::kShower)
00948 dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("best");
00949 else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon) plane->SetStrips();
00950 plane->SetStray(false);
00951 ++fixedPlanes;
00952 }
00953 }
00954 else if(useLastS
00955 && TMath::Abs(lastPlaneSpace)*0.25>TMath::Abs(lastPlaneCoG-lastPlaneSCoG)){
00956 plane = planeItr.Ptr();
00957 if(plane->IsStray()){
00958 if(plane->GetPlaneType() == DmxPlaneTypes::kShower)
00959 dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("second");
00960 plane->SetStray(false);
00961 ++fixedPlanes;
00962 }
00963 }
00964 else if(useLastT
00965 && TMath::Abs(lastPlaneSpace)*0.25>TMath::Abs(lastPlaneCoG-lastPlaneTCoG)){
00966 plane = planeItr.Ptr();
00967 if(plane->IsStray()){
00968 if(plane->GetPlaneType() == DmxPlaneTypes::kShower)
00969 dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("third");
00970 plane->SetStray(false);
00971 ++fixedPlanes;
00972 }
00973 }
00974
00975 return fixedPlanes;
00976 }
|
|
|
Definition at line 1160 of file DmxUtilities.cxx. References DmxPlane::GetCoG(), DmxPlane::GetPlaneView(), and DmxPlane::GetTimingOffset(). Referenced by AlgDeMuxCosmics::RunAlg(), and AlgDeMuxBeam::RunAlg(). 01161 {
01162
01163 //the idea here is to get the offset from the center line from timing for a plane
01164 //in one view and compare that to the center of gravity from the demuxer of the neighboring
01165 //plane in the other view.
01166
01167 //take the average of the centers of gravity for the planes on either side of the plane you
01168 //are interested in
01169
01170 planeItr.ResetFirst();
01171
01172 Int_t planeCtr = 0;
01173 Int_t numPlanes = planeItr.SizeSelect();
01174 Int_t forwardCtr = 0;
01175 Int_t backwardCtr = 0;
01176 bool opposite = false;
01177 Float_t offset = 0.;
01178 Float_t avOffsetDiff = 0.;
01179 Float_t avCoG = 0.;
01180 Float_t avCoGCnt = 0.;
01181
01182 DmxPlane *plane = 0;
01183 DmxPlane *nextPlane = 0;
01184 DmxPlane *prevPlane = 0;
01185
01186 while( planeItr.IsValid() ){
01187
01188 plane = planeItr.Ptr();
01189 //MSG("DmxUtil", Msg::kDebug) << "plane of offset = " << plane->GetPlaneNumber() << endl;
01190
01191 offset = plane->GetTimingOffset();
01192 //MSG("DmxUtil", Msg::kDebug) << " offset = " << offset << " is valid = " << (Int_t)plane->IsValid() << endl;
01193 forwardCtr = 0;
01194 backwardCtr = 0;
01195 avCoG = 0.;
01196 avCoGCnt = 0.;
01197 //this is the first plane, so you can only look at the one after it
01198 if(planeCtr == 0 && offset != -10.){
01199
01200 //advance to the next plane in the opposite view
01201 while(planeItr.IsValid() && opposite == false){
01202 planeItr.Next();
01203 ++forwardCtr;
01204 nextPlane = planeItr.Ptr();
01205 if(nextPlane->GetPlaneView() != plane->GetPlaneView()){
01206 //MSG("DmxUtil", Msg::kDebug) << "next opposite plane to vertex = " << nextPlane->GetPlaneNumber()
01207 // << " cog = " << nextPlane->GetCoG() << endl;
01208 opposite = true;
01209 avOffsetDiff += TMath::Abs(TMath::Abs(offset)-TMath::Abs(nextPlane->GetCoG()));
01210 }
01211 }//end loop to find next plane in opposite view
01212
01213 while(forwardCtr > 0){
01214 planeItr.Prev();
01215 --forwardCtr;
01216 }//end loop to return planeItr to original plane
01217
01218 }//end if first plane in set
01219 else if(planeCtr > 0 && planeCtr<numPlanes-1 && offset != -10.){
01220
01221 opposite = false;
01222 //advance to the next plane in the opposite view
01223 while(planeItr.IsValid() && !opposite){
01224 planeItr.Next();
01225 ++forwardCtr;
01226 if(planeItr.IsValid()){
01227 nextPlane = planeItr.Ptr();
01228 if(nextPlane->GetPlaneView() != plane->GetPlaneView()){
01229 //MSG("DmxUtil", Msg::kDebug) << "next opposite plane = " << nextPlane->GetPlaneNumber()
01230 // << " cog = " << nextPlane->GetCoG() << endl;
01231 opposite = true;
01232 avCoG += nextPlane->GetCoG();
01233 avCoGCnt += 1.;
01234 }
01235 }
01236 }//end loop to find next plane in opposite view
01237
01238 while(forwardCtr > 0){
01239 planeItr.Prev();
01240 --forwardCtr;
01241 }//end loop to return planeItr to original plane
01242
01243 opposite = false;
01244
01245 //if(planeItr.IsValid() )MSG("DmxUtil", Msg::kDebug) << "start backward look from plane "
01246 // << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber()
01247 // << endl;
01248 while(planeItr.IsValid() && !opposite){
01249 planeItr.Prev();
01250 ++backwardCtr;
01251 if(planeItr.IsValid()){
01252 prevPlane = planeItr.Ptr();
01253 if(prevPlane->GetPlaneView() != plane->GetPlaneView()){
01254 //MSG("DmxUtil", Msg::kDebug) << "prev opposite plane = " << prevPlane->GetPlaneNumber()
01255 // << " cog = " << prevPlane->GetCoG() << endl;
01256 opposite = true;
01257 avCoG += prevPlane->GetCoG();
01258 avCoGCnt += 1.;
01259 }
01260 }
01261 }//end loop to find prev plane in opposite view
01262
01263 while(backwardCtr > 0){
01264 planeItr.Next();
01265 --backwardCtr;
01266 }//end loop to return planeItr to original plane
01267
01268 //get the average CoG for the planes on either side and subtract the timing offset
01269 if(avCoGCnt > 0.) avOffsetDiff += TMath::Abs(TMath::Abs(offset)-TMath::Abs(avCoG/avCoGCnt));
01270
01271 }//end if plane between end planes
01272 else if(planeCtr == numPlanes-1 && offset != -10.){
01273 opposite = false;
01274 while(planeItr.IsValid() && !opposite){
01275 planeItr.Prev();
01276 ++backwardCtr;
01277 if(planeItr.IsValid()){
01278 prevPlane = planeItr.Ptr();
01279 if(prevPlane->GetPlaneView() != plane->GetPlaneView()){
01280 //MSG("DmxUtil", Msg::kDebug) << "prev opposite plane to end = " <<prevPlane->GetPlaneNumber()
01281 // << " cog = " << prevPlane->GetCoG() << endl;
01282 opposite = true;
01283 avOffsetDiff += TMath::Abs(TMath::Abs(offset)-TMath::Abs(prevPlane->GetCoG()));
01284 }
01285 }
01286 }//end loop to find prev plane in opposite view
01287
01288 while(backwardCtr > 0){
01289 planeItr.Next();
01290 --backwardCtr;
01291 }//end loop to return planeItr to original plane
01292
01293 }//end if plane is end plane
01294
01295 ++planeCtr;
01296 planeItr.Next();
01297 }//end loop over planes
01298
01299 if(numPlanes>0) avOffsetDiff /= (1.0*numPlanes);
01300
01301 planeItr.ResetFirst();
01302
01303 return avOffsetDiff;
01304 }
|
|
||||||||||||
|
Definition at line 241 of file DmxUtilities.cxx. Referenced by DmxDeMuxModule::Ana(), DmxDeMuxCosmicsModule::Ana(), and FillPlaneArray(). 00242 {
00243
00244 //go through and look to see if the digit is from cross talk
00245 Float_t crossTalk = -1.;
00246
00247 if(pixel == 0){
00248 crossTalk = planeSet[1] + planeSet[4];
00249 }
00250 if(pixel == 1) {
00251 crossTalk = planeSet[0] + planeSet[2] + planeSet[5];
00252 }
00253 if(pixel == 2){
00254 crossTalk = planeSet[1] + planeSet[3] + planeSet[6];
00255 }
00256 if(pixel == 3){
00257 crossTalk = planeSet[2] + planeSet[7];
00258 }
00259 if(pixel == 4){
00260 crossTalk = planeSet[0] + planeSet[5] + planeSet[8];
00261 }
00262 if(pixel == 5){
00263 crossTalk = planeSet[1] + planeSet[4] + planeSet[6] + planeSet[9];
00264 }
00265 if(pixel == 6){
00266 crossTalk = planeSet[2] + planeSet[5] + planeSet[7] + planeSet[10];
00267 }
00268 if(pixel == 7){
00269 crossTalk = planeSet[3] + planeSet[6] + planeSet[11];
00270 }
00271 if(pixel == 8){
00272 crossTalk = planeSet[4] + planeSet[9] + planeSet[12];
00273 }
00274 if(pixel == 9){
00275 crossTalk = planeSet[5] + planeSet[8] + planeSet[10] + planeSet[13];
00276 }
00277 if(pixel == 10){
00278 crossTalk = planeSet[6] + planeSet[9] + planeSet[11] + planeSet[14];
00279 }
00280 if(pixel == 11){
00281 crossTalk = planeSet[7] + planeSet[10] + planeSet[15];
00282 //MSG("DmxUtil", Msg::kDebug)<< fUtilities.GetEventNumber() << "\t"
00283 // << pixel << "\t" << "\t" << planeSet[11] << "\t"
00284 // <<planeSet[7]/planeSet[11] << "\t"
00285 // << planeSet[10]/planeSet[11] << "\t"
00286 // << planeSet[15]/planeSet[11] << endl;
00287 }
00288 if(pixel == 12){
00289 crossTalk = planeSet[8] + planeSet[13];
00290 }
00291 if(pixel == 13){
00292 crossTalk = planeSet[9] + planeSet[12] + planeSet[14];
00293 }
00294 if(pixel == 14){
00295 crossTalk = planeSet[10] + planeSet[13] + planeSet[15];
00296 //MSG("DmxUtil", Msg::kDebug)<< fUtilities.GetEventNumber() << "\t"
00297 // << pixel << "\t" << "\t" << planeSet[14] << "\t"
00298 // <<planeSet[13]/planeSet[14] << "\t"
00299 // << planeSet[10]/planeSet[14] << "\t"
00300 // << planeSet[15]/planeSet[14] << endl;
00301 }
00302 if(pixel == 15){
00303 crossTalk = planeSet[11] + planeSet[14];
00304 }
00305
00306 return crossTalk;
00307 }
|
|
||||||||||||||||||||
|
Definition at line 333 of file DmxUtilities.cxx. References digit(), CandDigitHandle::GetChannelId(), CandDigitHandle::GetCharge(), GetDigitPixel(), PlexSEIdAltL::GetEnd(), CandDigitHandle::GetPlexSEIdAltL(), RawChannelId::GetVaChannel(), and RawChannelId::GetVaChip(). Referenced by DmxDeMuxModule::Ana(), DmxDeMuxCosmicsModule::Ana(), and FillPlaneArray(). 00335 {
00336 Int_t pixel = -1;
00337 //loop over the digits and see if they have been set
00338 while( crdhi.IsValid() ){
00339
00340 //get the digit and it's strip, plane, and side
00341 CandDeMuxDigitHandle *digit = crdhi.Ptr();
00342 Int_t digitChip = digit->GetChannelId().GetVaChip();
00343 Int_t channel = digit->GetChannelId().GetVaChannel();
00344 pixel = GetDigitPixel(channel);
00345 Float_t charge = digit->GetCharge();
00346
00347 if( digitChip == chip){
00348 if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast){
00349 eastSet[pixel] += charge;
00350 }
00351 if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest){
00352 westSet[pixel] += charge;
00353 }
00354 }
00355 //MSG("Dmx2", Msg::kDebug)<< pixel << "\t" << channel<< "\t" << chip << "\t"
00356 // << planeSet[pixel] << endl;
00357
00358 crdhi.Next();
00359 }
00360
00361 crdhi.ResetFirst();
00362 return;
00363 }
|
|
||||||||||||||||
|
Definition at line 47 of file DmxUtilities.cxx. References CheckForXTalk(), digit(), FillHitPixels(), CandDigitHandle::GetChannelId(), CandDigitHandle::GetCharge(), CandHandle::GetDaughterIterator(), CandDeMuxDigitHandle::GetDeMuxDigitFlagWord(), GetDigitPixel(), Registry::GetDouble(), PlexSEIdAltL::GetEnd(), Registry::GetInt(), CandHandle::GetNDaughters(), PlexSEIdAltL::GetPlane(), CandDigitHandle::GetPlexSEIdAltL(), CandDigitHandle::GetPlexSEIdAltLWritable(), RawChannelId::GetVaChannel(), RawChannelId::GetVaChip(), PlexSEIdAltL::IsVetoShield(), MSG, CandDeMuxDigitHandle::SetDeMuxDigitFlagBit(), PlexSEIdAltL::SetDemuxVetoFlag(), and DmxStatus::SetPlaneArray(). Referenced by AlgDeMuxGolden::RunAlg(), AlgDeMuxCosmics::RunAlg(), and AlgDeMuxBeam::RunAlg(). 00049 {
00050
00051 TObjArray *planeArray = new TObjArray();
00052 MSG("DmxUtil", Msg::kDebug) << cdlh.GetNDaughters() << " digits in DmxUtil" << endl;
00053
00054 //make two iterators - one over the whole set and one to iterate over each plane
00055 //plane in turn
00056 CandDeMuxDigitHandleItr fCdhit = cdlh.GetDaughterIterator();
00057
00058 //create the KeyFunc
00059 CandDeMuxDigitHandleKeyFunc *planeKF = fCdhit.CreateKeyFunc();
00060
00061 //program the KeyFunc with the sort function
00062 planeKF->SetFun(KeyFromPlane);
00063
00064 //get the NavSet from the iterator and pass the KeyFunc to it
00065 fCdhit.GetSet()->AdoptSortKeyFunc(planeKF);
00066
00067 //clear the KF pointer because we no longer own the KeyFunc
00068 planeKF = 0;
00069
00070 //initialize variable to keep track of the current plane number
00071 Int_t planeNum = 0;
00072
00073 //MSG("DmxUtil", Msg::kDebug) << "current plane = " << planeNum << endl;
00074
00075 CandDeMuxDigitHandle *currentDigit = 0;
00076 CandDeMuxDigitHandle *digit = 0;
00077
00078 //go through the sorted Digits
00079 while( (currentDigit = fCdhit()) ){
00080
00081 Int_t plane = 0;
00082
00083 //find the plane number corresponding to that Digit
00084 plane = currentDigit->GetPlexSEIdAltL().GetPlane();
00085
00086 //MSG("DmxUtil", Msg::kDebug) << "plane number of current digit = " << plane
00087 // << endl;
00088
00089 //if the plane number of the current digit is different from the
00090 //current plane number, then add another DmxPlaneDebug object to
00091 //the TObjArray container for that plane number and increment the
00092 //fNumberOfPlanes variable
00093
00094 if(plane != planeNum && !currentDigit->GetPlexSEIdAltL().IsVetoShield()){
00095
00096 planeNum = plane;
00097
00098 //make a new iterator for this plane's digit. dont move this outside the
00099 //loop because then you have one iterator for all planes versus 1 iterator
00100 //per plane. makes a big difference.
00101 CandDeMuxDigitHandleItr planeCdhit(cdlh.GetDaughterIterator());
00102 CandDeMuxDigitHandleKeyFunc *plane1KF = planeCdhit.CreateKeyFunc();
00103 plane1KF->SetFun(KeyFromPlane);
00104 planeCdhit.GetSet()->AdoptSortKeyFunc(plane1KF,true,false);
00105 plane1KF = 0;
00106
00107 //get the slice for this plane
00108 planeCdhit.GetSet()->Slice(planeNum);
00109 planeCdhit.GetSet()->Update();
00110
00111 Int_t numDigits = 0;
00112 if( acd.GetInt("UseCandDigitMask")==1 ){
00113 //loop over the digits in this plane and figure out which ones might be
00114 //from cross talk. set the cutoff value for signal fraction to be in
00115 //AlgConfigDeMux
00116 Float_t eastSet0[] = {0., 0., 0., 0., 0., 0., 0., 0.,
00117 0., 0., 0., 0., 0., 0., 0., 0.};
00118 Float_t westSet0[] = {0., 0., 0., 0., 0., 0., 0., 0.,
00119 0., 0., 0., 0., 0., 0., 0., 0.};
00120 Float_t eastSet1[] = {0., 0., 0., 0., 0., 0., 0., 0.,
00121 0., 0., 0., 0., 0., 0., 0., 0.};
00122 Float_t westSet1[] = {0., 0., 0., 0., 0., 0., 0., 0.,
00123 0., 0., 0., 0., 0., 0., 0., 0.};
00124 Float_t eastSet2[] = {0., 0., 0., 0., 0., 0., 0., 0.,
00125 0., 0., 0., 0., 0., 0., 0., 0.};
00126 Float_t westSet2[] = {0., 0., 0., 0., 0., 0., 0., 0.,
00127 0., 0., 0., 0., 0., 0., 0., 0.};
00128
00129 FillHitPixels(planeCdhit, 0, eastSet0, westSet0);
00130 FillHitPixels(planeCdhit, 1, eastSet1, westSet1);
00131 FillHitPixels(planeCdhit, 2, eastSet2, westSet2);
00132 Int_t pixel = -1;
00133 Float_t crossTalk = 0.;
00134
00135 numDigits = planeCdhit.SizeSelect();
00136
00137 //loop over the digits in this plane to identify low signal and xtalk
00138 //digits. also set the DeMuxVetoFlag for all digits in this plane - if the
00139 //digits have a strip end in the demuxed hypothesis, the DeMuxVetoFlag will
00140 //be unset when the digit is demuxed.
00141 while( (digit = planeCdhit()) ){
00142
00143 pixel = GetDigitPixel(digit->GetChannelId().GetVaChannel());
00144 crossTalk = 0.;
00145
00146 currentDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(1);
00147
00148 if(digit->GetCharge()
00149 <acd.GetDouble("AveragePEGainConversion")*acd.GetDouble("IgnoreSignalLimit")){
00150 digit->SetDeMuxDigitFlagBit(CandDeMuxDigit::kLowSignal);
00151 //MSG("DmxUtil", Msg::kDebug) << "digit with charge " << digit->GetCharge()
00152 // << " has too little charge" << endl;
00153 --numDigits;
00154 }
00155
00156 if( digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast
00157 && digit->GetDeMuxDigitFlagWord() == 0){
00158
00159 if(digit->GetChannelId().GetVaChip() == 0)
00160 crossTalk = CheckForXTalk(pixel, eastSet0);
00161 else if(digit->GetChannelId().GetVaChip() == 1)
00162 crossTalk = CheckForXTalk(pixel, eastSet1);
00163 else if(digit->GetChannelId().GetVaChip() == 2)
00164 crossTalk = CheckForXTalk(pixel, eastSet2);
00165
00166 //mark it if its cross talk and subtract its charge to the total for the side
00167 if(digit->GetCharge()<(acd.GetDouble("XTalkFractionLimit")*crossTalk)){
00168 digit->SetDeMuxDigitFlagBit(CandDeMuxDigit::kXTalk);
00169 --numDigits;
00170 //MSG("DmxUtil", Msg::kDebug) << "digit with charge = " << digit->GetCharge()
00171 // << " is xtalk" << endl;
00172 }
00173 }//end if digit is from the east side and not already flagged
00174 else if( digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest &&
00175 digit->GetDeMuxDigitFlagWord() == 0){
00176
00177 if(digit->GetChannelId().GetVaChip() == 0)
00178 crossTalk = CheckForXTalk(pixel, westSet0);
00179 else if(digit->GetChannelId().GetVaChip() == 1)
00180 crossTalk = CheckForXTalk(pixel, westSet1);
00181 else if(digit->GetChannelId().GetVaChip() == 2)
00182 crossTalk = CheckForXTalk(pixel, westSet2);
00183
00184 //mark it if its not cross talk and add its charge to the total for the side
00185 if(digit->GetCharge()<(acd.GetDouble("XTalkFractionLimit")*crossTalk)){
00186 digit->SetDeMuxDigitFlagBit(CandDeMuxDigit::kXTalk);
00187 --numDigits;
00188 //MSG("DmxUtil", Msg::kDebug) << "digit with charge = " << digit->GetCharge()
00189 // << " is xtalk" << endl;
00190 }
00191 }//end if digit is from the west side and not already marked
00192
00193 }//end loop over digits to mark xtalk digits
00194 }
00195 //if not marking xtalk digits,
00196 //then set numDigits to size of the whole set
00197 else numDigits = planeCdhit.SizeSelect();
00198
00199 planeCdhit.ResetFirst();
00200
00201 //determine the type of plane object by the number of digits it
00202 //has. if > 2, its a shower plane, if <= 2 its
00203 //a muon plane. dont do a thing if the plane is a veto shield plane
00204
00205 if(numDigits > 2 && planeNum>-1){
00206 planeArray->AddLast(new DmxShowerPlane(acd, planeCdhit, planeNum));
00207 //MSG("DmxUtil", Msg::kDebug) << "Plane " << planeNum << " is a Shower Plane"
00208 // << " with " << numDigits << " digits" << endl;
00209 }
00210 else if(numDigits<=2 && planeNum>-1){
00211 planeArray->AddLast(new DmxMuonPlane(acd, planeCdhit, planeNum));
00212 //MSG("DmxUtil", Msg::kDebug) << "Plane " << planeNum << " is a Muon Plane"
00213 // << " with " << numDigits << " digits" << endl;
00214 }
00215
00216 //reset the iterator over each plane's digits
00217 planeCdhit.GetSet()->ClearSlice(false);
00218
00219 }//end if on a new plane
00220
00221 }//end loop over all digits
00222
00223 //cout << planeArray->GetLast() << endl;
00224 status->SetPlaneArray(planeArray);
00225
00226 //zero the pointer to the planeArray because the DmxStatus object now owns the
00227 //array
00228 planeArray = 0;
00229
00230 //MSG("DmxUtil", Msg::kDebug) << "done filling plane array" << endl;
00231 return;
00232 }
|
|
||||||||||||
|
Definition at line 575 of file DmxUtilities.cxx. References DmxPlane::GetNumberOfStrips(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlaneTypes::GetPlaneTypeString(), DmxPlane::IsValid(), and MSG. Referenced by AlgDeMuxBeam::RunAlg(). 00576 {
00577 //now locate where the muon starts, if anywhere - define the start
00578 //plane to be the first muon plane in a set of 4/5 muon planes with
00579 //3 of the 4 having at least 3 pe's each
00580
00581 //this method assumes that the planeItr was sliced to be from the vertex to the
00582 //end planes of the event
00583
00584 Int_t fMuonStartPlaneNumber = -1;
00585
00586 Int_t ctr = 0;
00587 Int_t muonCtr = 0;
00588 Int_t startPlane = 0;
00589 Int_t endPlane = 0;
00590 Int_t peCtr = 0;
00591 Int_t lastEventPlane = 0;
00592
00593 DmxPlane *currentPlane = 0;
00594
00595 while( (currentPlane = planeItr()) ){
00596
00597 MSG("DmxUtilities", Msg::kDebug) << "plane " << currentPlane->GetPlaneNumber()
00598 << " type " << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType())
00599 << " charge " << currentPlane->GetPlaneCharge()
00600 << " strips " << currentPlane->GetNumberOfStrips()
00601 << " is valid " << (Int_t)currentPlane->IsValid() << endl;
00602
00603 lastEventPlane = currentPlane->GetPlaneNumber();
00604 }
00605
00606 planeItr.ResetFirst();
00607
00608 while(planeItr.IsValid() && fMuonStartPlaneNumber == -1){
00609 ctr = 0;
00610
00611 muonCtr = 0;
00612 startPlane = planeItr.Ptr()->GetPlaneNumber();
00613 endPlane = startPlane+10;
00614 peCtr = 0;
00615
00616 MSG("DmxUtilities", Msg::kDebug) << "plane number " << startPlane
00617 << " charge " << planeItr.Ptr()->GetPlaneCharge()
00618 << " type "
00619 << DmxPlaneTypes::GetPlaneTypeString(planeItr.Ptr()->GetPlaneType())
00620 << endl;
00621 planeItr.Next();
00622
00623 //check 4 planes at a time
00624 while(ctr < 4 && planeItr.IsValid() ){
00625 currentPlane = planeItr.Ptr();
00626
00627 MSG("DmxUtilities", Msg::kDebug) << "plane number " << currentPlane->GetPlaneNumber()
00628 << " charge " << currentPlane->GetPlaneCharge()
00629 << " type "
00630 << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType())
00631 << " num strips " << currentPlane->GetNumberOfStrips()
00632 << endl;
00633 endPlane = currentPlane->GetPlaneNumber();
00634
00635 //check to see that its a muon plane or only has 2 strips hit
00636 if( currentPlane->GetNumberOfStrips() < 3 ) ++muonCtr;
00637
00638 if(currentPlane->GetPlaneCharge() > (peGainConversion * 3.)) ++peCtr;
00639
00640 MSG("DmxUtilities", Msg::kDebug) << " muon planes = " << muonCtr
00641 << " peCtr = " << peCtr << endl;
00642 planeItr.Next();
00643 ++ctr;
00644 }
00645
00646 if(muonCtr == 4 && peCtr >= 3 && endPlane-startPlane<=5)
00647 fMuonStartPlaneNumber = startPlane;
00648
00649 //if there are still planes left to examine, the iterator will be valid
00650 //otherwise, you ran out of planes, so break the cycle
00651 if(planeItr.IsValid() ){
00652 planeItr.Prev();
00653 planeItr.Prev();
00654 planeItr.Prev();
00655 planeItr.Prev();
00656 }
00657 }//end loop to find the muon track starting plane
00658
00659 //found the most likely start plane for the muon. but lets make sure
00660 //we arent being fooled into thinking a shower event is really a muon event.
00661 //loop over the planes and see the number of planes that are shower like and if
00662 //there are lots more shower planes than muon planes
00663 planeItr.GetSet()->ClearSlice(false);
00664
00665 planeItr.GetSet()->Slice(fMuonStartPlaneNumber, lastEventPlane);
00666 planeItr.GetSet()->Update();
00667 planeItr.ResetFirst();
00668
00669 Int_t numShowerPlanes = 0;
00670 Int_t numMuonPlanes = 0;
00671 while( (currentPlane = planeItr()) ){
00672
00673 if(currentPlane->GetPlaneType() == DmxPlaneTypes::kShower
00674 && currentPlane->GetNumberOfStrips()>2) ++numShowerPlanes;
00675 else ++numMuonPlanes;
00676
00677 }
00678
00679 if(1.*numShowerPlanes>0.25*numMuonPlanes) fMuonStartPlaneNumber = -1;
00680
00681 planeItr.GetSet()->ClearSlice(false);
00682
00683 return fMuonStartPlaneNumber;
00684 }
|
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 1309 of file DmxUtilities.cxx. 01312 {
01313
01314 //use method of determinants to do the least squares fit
01315 //this is from bevington chapter 7
01316 //
01317 //f_1(x) = 1
01318 //f_2(x) = x
01319 //f_3(x) = x^2
01320 //f_4(x) = x^3
01321
01322 //need to find the sums
01323 //SIGMA y_i*f_1(x_i)/weight_i
01324 //SIGMA y_i*f_2(x_i)/weight_i
01325 //SIGMA y_i*f_3(x_i)/weight_i
01326 //SIGMA y_i*f_4(x_i)/weight_i
01327 //SIGMA f_1(x_i)*f_1(x_i)/weight_i
01328 //SIGMA f_1(x_i)*f_2(x_i)/weight_i
01329 //SIGMA f_1(x_i)*f_3(x_i)/weight_i
01330 //SIGMA f_1(x_i)*f_4(x_i)/weight_i
01331 //SIGMA f_2(x_i)*f_1(x_i)/weight_i
01332 //SIGMA f_2(x_i)*f_2(x_i)/weight_i
01333 //SIGMA f_2(x_i)*f_3(x_i)/weight_i
01334 //SIGMA f_2(x_i)*f_4(x_i)/weight_i
01335 //SIGMA f_3(x_i)*f_1(x_i)/weight_i
01336 //SIGMA f_3(x_i)*f_2(x_i)/weight_i
01337 //SIGMA f_3(x_i)*f_3(x_i)/weight_i
01338 //SIGMA f_3(x_i)*f_4(x_i)/weight_i
01339 //SIGMA f_4(x_i)*f_1(x_i)/weight_i
01340 //SIGMA f_4(x_i)*f_2(x_i)/weight_i
01341 //SIGMA f_4(x_i)*f_3(x_i)/weight_i
01342 //SIGMA f_4(x_i)*f_4(x_i)/weight_i
01343 //
01344 //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
01345 //10 instead of 16 variables
01346
01347 Double_t yf1divWeight = 0.;
01348 Double_t yf2divWeight = 0.;
01349 Double_t yf3divWeight = 0.;
01350 Double_t yf4divWeight = 0.;
01351 Double_t f1f1divWeight = 0.;
01352 Double_t f1f2divWeight = 0.;
01353 Double_t f1f3divWeight = 0.;
01354 Double_t f1f4divWeight = 0.;
01355 Double_t f2f2divWeight = 0.;
01356 Double_t f2f3divWeight = 0.;
01357 Double_t f2f4divWeight = 0.;
01358 Double_t f3f3divWeight = 0.;
01359 Double_t f3f4divWeight = 0.;
01360 Double_t f4f4divWeight = 0.;
01361
01362 for(Int_t i=0; i<nPoints; i++){
01363 yf1divWeight += y[i]/(weight[i]*weight[i]);
01364 yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
01365 yf3divWeight += (y[i]*x[i]*x[i])/(weight[i]*weight[i]);
01366 yf4divWeight += (y[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01367 f1f1divWeight += 1./(weight[i]*weight[i]);
01368 f1f2divWeight += x[i]/(weight[i]*weight[i]);
01369 f1f3divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01370 f1f4divWeight += (x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01371 f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01372 f2f3divWeight += (x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01373 f2f4divWeight += (x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01374 f3f3divWeight += (x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01375 f3f4divWeight += (x[i]*x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01376 f4f4divWeight += (x[i]*x[i]*x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01377 }
01378
01379 //make a bunch of matricies to hold these values
01380 TMatrix deltaMatrix(4,4);
01381 TMatrix a1Matrix(4,4);
01382 TMatrix a2Matrix(4,4);
01383 TMatrix a3Matrix(4,4);
01384 TMatrix a4Matrix(4,4);
01385
01386 //fill the matricies
01387 deltaMatrix(0,0) = f1f1divWeight;
01388 deltaMatrix(0,1) = f1f2divWeight;
01389 deltaMatrix(0,2) = f1f3divWeight;
01390 deltaMatrix(0,3) = f1f4divWeight;
01391 deltaMatrix(1,0) = f1f2divWeight;
01392 deltaMatrix(1,1) = f2f2divWeight;
01393 deltaMatrix(1,2) = f2f3divWeight;
01394 deltaMatrix(1,3) = f2f4divWeight;
01395 deltaMatrix(2,0) = f1f3divWeight;
01396 deltaMatrix(2,1) = f2f3divWeight;
01397 deltaMatrix(2,2) = f3f3divWeight;
01398 deltaMatrix(2,3) = f3f4divWeight;
01399 deltaMatrix(3,0) = f1f4divWeight;
01400 deltaMatrix(3,1) = f2f4divWeight;
01401 deltaMatrix(3,2) = f3f4divWeight;
01402 deltaMatrix(3,3) = f4f4divWeight;
01403
01404 a1Matrix(0,0) = yf1divWeight;
01405 a1Matrix(0,1) = f1f2divWeight;
01406 a1Matrix(0,2) = f1f3divWeight;
01407 a1Matrix(0,3) = f1f4divWeight;
01408 a1Matrix(1,0) = yf2divWeight;
01409 a1Matrix(1,1) = f2f2divWeight;
01410 a1Matrix(1,2) = f2f3divWeight;
01411 a1Matrix(1,3) = f2f4divWeight;
01412 a1Matrix(2,0) = yf3divWeight;
01413 a1Matrix(2,1) = f2f3divWeight;
01414 a1Matrix(2,2) = f3f3divWeight;
01415 a1Matrix(2,3) = f3f4divWeight;
01416 a1Matrix(3,0) = yf4divWeight;
01417 a1Matrix(3,1) = f2f4divWeight;
01418 a1Matrix(3,2) = f3f4divWeight;
01419 a1Matrix(3,3) = f4f4divWeight;
01420
01421 a2Matrix(0,0) = f1f1divWeight;
01422 a2Matrix(0,1) = yf1divWeight;
01423 a2Matrix(0,2) = f1f3divWeight;
01424 a2Matrix(0,3) = f1f4divWeight;
01425 a2Matrix(1,0) = f1f2divWeight;
01426 a2Matrix(1,1) = yf2divWeight;
01427 a2Matrix(1,2) = f2f3divWeight;
01428 a2Matrix(1,3) = f2f4divWeight;
01429 a2Matrix(2,0) = f1f3divWeight;
01430 a2Matrix(2,1) = yf3divWeight;
01431 a2Matrix(2,2) = f3f3divWeight;
01432 a2Matrix(2,3) = f3f4divWeight;
01433 a2Matrix(3,0) = f1f4divWeight;
01434 a2Matrix(3,1) = yf4divWeight;
01435 a2Matrix(3,2) = f3f4divWeight;
01436 a2Matrix(3,3) = f4f4divWeight;
01437
01438 a3Matrix(0,0) = f1f1divWeight;
01439 a3Matrix(0,1) = f1f2divWeight;
01440 a3Matrix(0,2) = yf1divWeight;
01441 a3Matrix(0,3) = f1f4divWeight;
01442 a3Matrix(1,0) = f1f2divWeight;
01443 a3Matrix(1,1) = f2f2divWeight;
01444 a3Matrix(1,2) = yf2divWeight;
01445 a3Matrix(1,3) = f2f4divWeight;
01446 a3Matrix(2,0) = f1f3divWeight;
01447 a3Matrix(2,1) = f2f3divWeight;
01448 a3Matrix(2,2) = yf3divWeight;
01449 a3Matrix(2,3) = f3f4divWeight;
01450 a3Matrix(3,0) = f1f4divWeight;
01451 a3Matrix(3,1) = f2f4divWeight;
01452 a3Matrix(3,2) = yf4divWeight;
01453 a3Matrix(3,3) = f4f4divWeight;
01454
01455 a4Matrix(0,0) = f1f1divWeight;
01456 a4Matrix(0,1) = f1f2divWeight;
01457 a4Matrix(0,2) = f1f3divWeight;
01458 a4Matrix(0,3) = yf1divWeight;
01459 a4Matrix(1,0) = f1f2divWeight;
01460 a4Matrix(1,1) = f2f2divWeight;
01461 a4Matrix(1,2) = f2f3divWeight;
01462 a4Matrix(1,3) = yf2divWeight;
01463 a4Matrix(2,0) = f1f3divWeight;
01464 a4Matrix(2,1) = f2f3divWeight;
01465 a4Matrix(2,2) = f3f3divWeight;
01466 a4Matrix(2,3) = yf3divWeight;
01467 a4Matrix(3,0) = f1f4divWeight;
01468 a4Matrix(3,1) = f2f4divWeight;
01469 a4Matrix(3,2) = f3f4divWeight;
01470 a4Matrix(3,3) = yf4divWeight;
01471
01472 a1 = -10000.;
01473 a2 = -10000.;
01474 a3 = -10000.;
01475 a4 = -10000.;
01476 Double_t delta = deltaMatrix.Determinant();
01477
01478 if(delta != 0.){
01479 a1 = a1Matrix.Determinant()/delta;
01480 a2 = a2Matrix.Determinant()/delta;
01481 a3 = a3Matrix.Determinant()/delta;
01482 a4 = a4Matrix.Determinant()/delta;
01483 }
01484
01485 //find the chi^2 value for the fit
01486 chiSq = 0.;
01487 Double_t chi = 0.;
01488 for(Int_t j = 0; j<nPoints; j++){
01489
01490 chi = (y[j] - (a1 + a2*x[j] + a3*(x[j]*x[j]) + a4*(x[j]*x[j]*x[j])))/weight[j];
01491 chiSq += chi*chi;
01492
01493 }
01494
01495 //get the chiSq per degree of freedom
01496
01497 chiSq /= (1.*nPoints);
01498
01499 return;
01500 }
|
|
|
Definition at line 463 of file DmxUtilities.cxx. Referenced by AlgDeMuxGolden::RunAlg(), AlgDeMuxCosmics::RunAlg(), and AlgDeMuxBeam::RunAlg(). 00464 {
00465 Int_t endPlaneNumber = 500;
00466
00467 while( planeItr.IsValid() && endPlaneNumber == 500){
00468 Int_t planeNumber = planeItr.Ptr()->GetPlaneNumber();
00469
00470 planeItr.Next();
00471 if( planeItr.IsValid() ){
00472 //if the difference between the current plane number and the previous is >= 10
00473 //call the previous plane the last one in the event
00474 if( TMath::Abs(planeItr.Ptr()->GetPlaneNumber()
00475 - planeNumber) >= 5 ){ endPlaneNumber = planeNumber; }
00476 }
00477 else{
00478 //didn't find a cutoff plane so call the last plane it
00479 endPlaneNumber = planeNumber;
00480 }
00481 }
00482
00483 return endPlaneNumber;
00484 }
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 1634 of file DmxUtilities.cxx. References MSG. Referenced by AlgDeMuxCosmics::FindFit(), AlgDeMuxBeam::FindFitForHypothesisSet(), and AlgDeMuxCosmics::FindWindowCosmicSolution(). 01637 {
01638
01639 //use method of determinants to do the least squares fit
01640 //this is from bevington chapter 6
01641 //
01642 //f_1(x) = 1
01643 //f_2(x) = x
01644
01645 //need to find the sums
01646 //SIGMA y_i*f_1(x_i)/weight_i
01647 //SIGMA y_i*f_2(x_i)/weight_i
01648 //SIGMA f_1(x_i)*f_1(x_i)/weight_i
01649 //SIGMA f_1(x_i)*f_2(x_i)/weight_i
01650 //SIGMA f_2(x_i)*f_1(x_i)/weight_i
01651 //SIGMA f_2(x_i)*f_2(x_i)/weight_i
01652 //
01653 //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
01654 //3 instead of 4 variables
01655
01656 Double_t yf1divWeight = 0.;
01657 Double_t yf2divWeight = 0.;
01658 Double_t f1f1divWeight = 0.;
01659 Double_t f1f2divWeight = 0.;
01660 Double_t f2f2divWeight = 0.;
01661
01662 for(Int_t i=0; i<nPoints; ++i){
01663 MSG("DmxUtilities", Msg::kDebug) << " x = " << x[i]
01664 << " y = " << y[i]
01665 << " weight = " << weight[i] << endl;
01666 if(weight[i]>0.){
01667 yf1divWeight += y[i]/(weight[i]*weight[i]);
01668 yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
01669 f1f1divWeight += 1./(weight[i]*weight[i]);
01670 f1f2divWeight += x[i]/(weight[i]*weight[i]);
01671 f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01672 }
01673 }
01674
01675 //make a bunch of matricies to hold these values
01676 TMatrix deltaMatrix(2,2);
01677 TMatrix a1Matrix(2,2);
01678 TMatrix a2Matrix(2,2);
01679
01680 //fill the matricies
01681 deltaMatrix(0,0) = f1f1divWeight;
01682 deltaMatrix(0,1) = f1f2divWeight;
01683 deltaMatrix(1,0) = f1f2divWeight;
01684 deltaMatrix(1,1) = f2f2divWeight;
01685
01686 a1Matrix(0,0) = yf1divWeight;
01687 a1Matrix(0,1) = f1f2divWeight;
01688 a1Matrix(1,0) = yf2divWeight;
01689 a1Matrix(1,1) = f2f2divWeight;
01690
01691 a2Matrix(0,0) = f1f1divWeight;
01692 a2Matrix(0,1) = yf1divWeight;
01693 a2Matrix(1,0) = f1f2divWeight;
01694 a2Matrix(1,1) = yf2divWeight;
01695
01696
01697 a1 = -10000.;
01698 a2 = -10000.;
01699 Double_t delta = deltaMatrix.Determinant();
01700
01701 if(delta != 0.){
01702 a1 = a1Matrix.Determinant()/delta;
01703 a2 = a2Matrix.Determinant()/delta;
01704 }
01705
01706 //find the chi^2 value for the fit
01707 chiSq = 0.;
01708 Double_t chi = 0.;
01709 for(Int_t j = 0; j<nPoints; ++j){
01710
01711 chi = (y[j] - (a1 + a2*x[j]))/weight[j];
01712 chiSq += chi*chi;
01713
01714 }
01715
01716 //get the chiSq per degree of freedom
01717
01718 chiSq /= (1.*nPoints);
01719
01720 return;
01721 }
|
|
||||||||||||
|
Definition at line 487 of file DmxUtilities.cxx. References DmxPlane::GetNumberOfStrips(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlaneTypes::GetPlaneTypeString(), DmxPlane::IsValid(), and MSG. 00488 {
00489 //now locate where the muon starts, if anywhere - define the start
00490 //plane to be the first muon plane in a set of 4/5 muon planes with
00491 //3 of the 4 having at least 3 pe's each
00492
00493 //this method assumes that the planeItr was sliced to be from the vertex to the
00494 //end planes of the event
00495
00496 Int_t fMuonStartPlaneNumber = -1;
00497
00498 Int_t ctr = 0;
00499 Int_t muonCtr = 0;
00500 Int_t startPlane = 0;
00501 Int_t endPlane = 0;
00502 Int_t peCtr = 0;
00503
00504 DmxPlane *currentPlane = 0;
00505
00506 while( (currentPlane = planeItr()) ){
00507
00508 MSG("DmxUtilities", Msg::kDebug) << "plane " << currentPlane->GetPlaneNumber()
00509 << " type " << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType())
00510 << " charge " << currentPlane->GetPlaneCharge()
00511 << " strips " << currentPlane->GetNumberOfStrips()
00512 << " is valid " << (Int_t)currentPlane->IsValid() << endl;
00513 }
00514
00515 planeItr.ResetFirst();
00516
00517 while(planeItr.IsValid() && fMuonStartPlaneNumber == -1){
00518 ctr = 0;
00519
00520 muonCtr = 0;
00521 startPlane = planeItr.Ptr()->GetPlaneNumber();
00522 endPlane = startPlane+10;
00523 peCtr = 0;
00524
00525 MSG("DmxUtilities", Msg::kDebug) << "plane number " << startPlane
00526 << " charge " << planeItr.Ptr()->GetPlaneCharge()
00527 << " type "
00528 << DmxPlaneTypes::GetPlaneTypeString(planeItr.Ptr()->GetPlaneType())
00529 << endl;
00530 planeItr.Next();
00531
00532 //check 4 planes at a time
00533 while(ctr < 4 && planeItr.IsValid() ){
00534 currentPlane = planeItr.Ptr();
00535
00536 MSG("DmxUtilities", Msg::kDebug) << "plane number " << currentPlane->GetPlaneNumber()
00537 << " charge " << currentPlane->GetPlaneCharge()
00538 << " type "
00539 << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType())
00540 << " num strips " << currentPlane->GetNumberOfStrips()
00541 << endl;
00542 endPlane = currentPlane->GetPlaneNumber();
00543
00544 //check to see that its a muon plane or only has 2 strips hit
00545 if( currentPlane->GetNumberOfStrips() < 3 ) ++muonCtr;
00546
00547 if(currentPlane->GetPlaneCharge() > (peGainConversion * 3.)) ++peCtr;
00548
00549 MSG("DmxUtilities", Msg::kDebug) << " muon planes = " << muonCtr
00550 << " peCtr = " << peCtr << endl;
00551 planeItr.Next();
00552 ++ctr;
00553 }
00554
00555 if(muonCtr == 4 && peCtr >= 3 && endPlane-startPlane<=5)
00556 fMuonStartPlaneNumber = startPlane;
00557
00558 //if there are still planes left to examine, the iterator will be valid
00559 //otherwise, you ran out of planes, so break the cycle
00560 if(planeItr.IsValid() ){
00561 planeItr.Prev();
00562 planeItr.Prev();
00563 planeItr.Prev();
00564 planeItr.Prev();
00565 }
00566 }//end loop to find the muon track starting plane
00567
00568 planeItr.GetSet()->ClearSlice(false);
00569 planeItr.Reset();
00570
00571 return fMuonStartPlaneNumber;
00572 }
|
|
||||||||||||
|
Definition at line 1089 of file DmxUtilities.cxx. References DmxPlane::GetCoG(), DmxPlane::GetPlaneType(), DmxPlane::IsValid(), and DmxPlane::SetStray(). Referenced by AlgDeMuxCosmics::RunAlg(), AlgDeMuxBeam::RunAlg(), AlgDeMuxGolden::UseGoldenFit(), AlgDeMuxCosmics::UseSingleFit(), and AlgDeMuxCosmics::UseSlidingWindow(). 01090 {
01091
01092 planeItr.Reset();
01093
01094 Int_t strayPlanes = 0;
01095
01096 DmxPlane *plane = 0;
01097
01098 //MSG("DmxUtil", Msg::kDebug) << "checking for stray planes in event " << endl;
01099 while( planeItr.IsValid() ){
01100
01101 plane = planeItr.Ptr();
01102
01103 if( plane->IsValid() ){
01104 //&& plane->GetPlaneNumber()>= firstPlane && plane->GetPlaneNumber()<= lastPlane){
01105 Double_t fitCoG = -1.;
01106
01107
01108 if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
01109 Double_t diff1 = 0.;
01110 Double_t diff2 = 0.;
01111 Double_t diff3 = 0.;
01112 fitCoG = plane->GetCoG();
01113
01114 if( dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis() ){
01115 diff1 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG());
01116 //MSG("DmxUtil", Msg::kDebug)<<"best " << dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG()<<endl;
01117 }
01118 if( dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis() ){
01119 diff2 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG());
01120 //MSG("DmxUtil", Msg::kDebug)<<"2nd " <<dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG()<<endl;
01121 }
01122 if( dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis() ){
01123 diff3 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG());
01124 //MSG("DmxUtil", Msg::kDebug)<<"3rd " << dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG()<<endl;
01125 }
01126
01127 //fStrayCut is in strip numbers, the 0.042 factor is the width of a strip
01128 if(diff1>strayCut*.042 && diff2>strayCut*.042 && diff3>strayCut*.042){
01129 //MSG("DmxUtil", Msg::kDebug) << "plane = " << plane->GetPlaneNumber() << " is stray "
01130 // << diff1 << " " << diff2 << " " << diff3 << endl;
01131
01132 ++strayPlanes;
01133 plane->SetStray(true);
01134 // MSG("DmxUtil", Msg::kDebug) << "\tstray plane = " << plane->GetPlaneNumber()
01135 // << "\tfit cog = " << fitCoG << endl;
01136 }
01137 }//end if shower plane
01138 // else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
01139 // fitCoG = plane->GetCoG();
01140 // if( TMath::Abs(fitCoG-dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG()) > strayCut*0.042){
01141 // //MSG("DmxUtil", Msg::kDebug)<<"muon " << plane->GetCoG()<<endl;
01142 // ++strayPlanes;
01143 // plane->SetStray(true);
01144 // MSG("DmxUtil", Msg::kDebug) << "\tstray plane = " << plane->GetPlaneNumber()
01145 // << "\tfit cog = " << fitCoG << " initial cog = "
01146 // << dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG()
01147 // <<endl;
01148 // }
01149 // }//end if muon plane
01150 }//end if it is a valid plane
01151 planeItr.Next();
01152 }//end loop over planes
01153
01154 planeItr.Reset();
01155
01156 return strayPlanes;
01157 }
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 1505 of file DmxUtilities.cxx. 01508 {
01509
01510 //use method of determinants to do the least squares fit
01511 //this is from bevington chapter 7
01512 //
01513 //f_1(x) = 1
01514 //f_2(x) = x
01515 //f_3(x) = x^2
01516
01517 //need to find the sums
01518 //SIGMA y_i*f_1(x_i)/weight_i
01519 //SIGMA y_i*f_2(x_i)/weight_i
01520 //SIGMA y_i*f_3(x_i)/weight_i
01521 //SIGMA f_1(x_i)*f_1(x_i)/weight_i
01522 //SIGMA f_1(x_i)*f_2(x_i)/weight_i
01523 //SIGMA f_1(x_i)*f_3(x_i)/weight_i
01524 //SIGMA f_2(x_i)*f_1(x_i)/weight_i
01525 //SIGMA f_2(x_i)*f_2(x_i)/weight_i
01526 //SIGMA f_2(x_i)*f_3(x_i)/weight_i
01527 //SIGMA f_3(x_i)*f_1(x_i)/weight_i
01528 //SIGMA f_3(x_i)*f_2(x_i)/weight_i
01529 //SIGMA f_3(x_i)*f_3(x_i)/weight_i
01530
01531 //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
01532 //6 instead of 9 variables
01533
01534 Double_t yf1divWeight = 0.;
01535 Double_t yf2divWeight = 0.;
01536 Double_t yf3divWeight = 0.;
01537 Double_t f1f1divWeight = 0.;
01538 Double_t f1f2divWeight = 0.;
01539 Double_t f1f3divWeight = 0.;
01540 Double_t f2f2divWeight = 0.;
01541 Double_t f2f3divWeight = 0.;
01542 Double_t f3f3divWeight = 0.;
01543
01544 for(Int_t i=0; i<nPoints; i++){
01545 yf1divWeight += y[i]/(weight[i]*weight[i]);
01546 yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
01547 yf3divWeight += (y[i]*x[i]*x[i])/(weight[i]*weight[i]);
01548 f1f1divWeight += 1./(weight[i]*weight[i]);
01549 f1f2divWeight += x[i]/(weight[i]*weight[i]);
01550 f1f3divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01551 f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01552 f2f3divWeight += (x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01553 f3f3divWeight += (x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01554 }
01555
01556 //make a bunch of matricies to hold these values
01557 TMatrix deltaMatrix(3,3);
01558 TMatrix a1Matrix(3,3);
01559 TMatrix a2Matrix(3,3);
01560 TMatrix a3Matrix(3,3);
01561
01562 //fill the matricies
01563 deltaMatrix(0,0) = f1f1divWeight;
01564 deltaMatrix(0,1) = f1f2divWeight;
01565 deltaMatrix(0,2) = f1f3divWeight;
01566 deltaMatrix(1,0) = f1f2divWeight;
01567 deltaMatrix(1,1) = f2f2divWeight;
01568 deltaMatrix(1,2) = f2f3divWeight;
01569 deltaMatrix(2,0) = f1f3divWeight;
01570 deltaMatrix(2,1) = f2f3divWeight;
01571 deltaMatrix(2,2) = f3f3divWeight;
01572
01573 a1Matrix(0,0) = yf1divWeight;
01574 a1Matrix(0,1) = f1f2divWeight;
01575 a1Matrix(0,2) = f1f3divWeight;
01576 a1Matrix(1,0) = yf2divWeight;
01577 a1Matrix(1,1) = f2f2divWeight;
01578 a1Matrix(1,2) = f2f3divWeight;
01579 a1Matrix(2,0) = yf3divWeight;
01580 a1Matrix(2,1) = f2f3divWeight;
01581 a1Matrix(2,2) = f3f3divWeight;
01582
01583 a2Matrix(0,0) = f1f1divWeight;
01584 a2Matrix(0,1) = yf1divWeight;
01585 a2Matrix(0,2) = f1f3divWeight;
01586 a2Matrix(1,0) = f1f2divWeight;
01587 a2Matrix(1,1) = yf2divWeight;
01588 a2Matrix(1,2) = f2f3divWeight;
01589 a2Matrix(2,0) = f1f3divWeight;
01590 a2Matrix(2,1) = yf3divWeight;
01591 a2Matrix(2,2) = f3f3divWeight;
01592
01593 a3Matrix(0,0) = f1f1divWeight;
01594 a3Matrix(0,1) = f1f2divWeight;
01595 a3Matrix(0,2) = yf1divWeight;
01596 a3Matrix(1,0) = f1f2divWeight;
01597 a3Matrix(1,1) = f2f2divWeight;
01598 a3Matrix(1,2) = yf2divWeight;
01599 a3Matrix(2,0) = f1f3divWeight;
01600 a3Matrix(2,1) = f2f3divWeight;
01601 a3Matrix(2,2) = yf3divWeight;
01602
01603 a1 = -10000.;
01604 a2 = -10000.;
01605 a3 = -10000.;
01606 Double_t delta = deltaMatrix.Determinant();
01607
01608 if(delta != 0.){
01609 a1 = a1Matrix.Determinant()/delta;
01610 a2 = a2Matrix.Determinant()/delta;
01611 a3 = a3Matrix.Determinant()/delta;
01612 }
01613
01614 //find the chi^2 value for the fit
01615 chiSq = 0.;
01616 Double_t chi = 0.;
01617 for(Int_t j = 0; j<nPoints; j++){
01618
01619 chi = (y[j] - (a1 + a2*x[j] + a3*(x[j]*x[j])))/weight[j];
01620 chiSq += chi*chi;
01621
01622 }
01623
01624 //get the chiSq per degree of freedom
01625
01626 chiSq /= (1.*nPoints);
01627
01628 return;
01629 }
|
|
|
Definition at line 367 of file DmxUtilities.cxx. References MSG. Referenced by AlgDeMuxGolden::RunAlg(), AlgDeMuxCosmics::RunAlg(), and AlgDeMuxBeam::RunAlg(). 00368 {
00369 MSG("DmxUtil", Msg::kDebug) << "In FindVertexPlane" << endl;
00370
00371 planeItr.ResetFirst();
00372
00373 Int_t vertexPlaneNumber = -1;
00374 Int_t firstPlaneNumber = 0;
00375 if( planeItr.IsValid() )
00376 firstPlaneNumber = planeItr.Ptr()->GetPlaneNumber();
00377
00378 planeItr.ResetFirst();
00379
00380 //find the vertex. make it the first plane in a set of 5 planes where at least 3
00381 //of the planes are valid. use the number 3 since it is the least number of planes
00382 //you can have and still drop one to refine a fit. use 5 because the trigger is 4/5
00383 Int_t beginPlaneNumber = 0;
00384 Int_t endPlaneNumber = 0;
00385 Int_t validCtr = 0;
00386 Int_t cnt = 0;
00387
00388 while( planeItr.IsValid() ){
00389
00390 beginPlaneNumber = planeItr.Ptr()->GetPlaneNumber();
00391 endPlaneNumber = 0;
00392 validCtr = 0;
00393
00394 cnt = 0;
00395 while( planeItr.IsValid() && cnt < 5 && validCtr < 3){
00396 if( planeItr.Ptr()->IsValid() ){ ++validCtr;}
00397 ++cnt;
00398 endPlaneNumber = planeItr.Ptr()->GetPlaneNumber();
00399 MSG("DmxUtil", Msg::kDebug)<<"Start Plane = " << beginPlaneNumber
00400 << "\tThis Plane = "
00401 << endPlaneNumber << "\tValid = "
00402 << (Int_t)planeItr.Ptr()->IsValid()
00403 << "\tValidCtr = " << validCtr << endl;
00404 planeItr.Next();
00405 }
00406
00407 if( validCtr >= 3 && (endPlaneNumber-beginPlaneNumber) <= 5 ){
00408 vertexPlaneNumber = beginPlaneNumber;
00409 MSG("DmxUtil", Msg::kDebug)<<"\tVertex Plane = " << vertexPlaneNumber << endl;
00410 break;
00411 }
00412 else{
00413 for(Int_t a = 0; a < cnt-1; a++){
00414 planeItr.Prev();
00415 }
00416 }
00417 } // end loop to find vertex
00418
00419 //see if you can find any contiguos planes before the vertex, if so, make them the vertex
00420 //slice from 0 to the vertex and then iterate backwards. if there is a gap in signal
00421 //of more than 2 planes, stop
00422
00423 planeItr.Reset();
00424 planeItr.GetSet()->ClearSlice(false);
00425
00426 Int_t nextPlane = 0;
00427
00428 if( vertexPlaneNumber > -1 && firstPlaneNumber < vertexPlaneNumber){
00429 planeItr.GetSet()->Slice(firstPlaneNumber, vertexPlaneNumber);
00430 planeItr.GetSet()->Update();
00431
00432 planeItr.ResetLast();
00433
00434 while( planeItr.IsValid() ){
00435
00436 Int_t currentPlane = planeItr.Ptr()->GetPlaneNumber();
00437 MSG("DmxUtil", Msg::kDebug) << "\tcurrentPlaneNumber = " << currentPlane;
00438 planeItr.Prev();
00439 nextPlane = 0;
00440 if(planeItr.IsValid()){
00441 nextPlane = planeItr.Ptr()->GetPlaneNumber();
00442 MSG("DmxUtil", Msg::kDebug) << "\tnextPlaneNumber = " << nextPlane << endl;
00443
00444 if( vertexPlaneNumber - nextPlane < 3)
00445 vertexPlaneNumber = nextPlane;
00446 else break;
00447 }
00448
00449 }//end loop over planes
00450
00451 planeItr.GetSet()->ClearSlice(false);
00452 planeItr.ResetFirst();
00453 }
00454
00455 planeItr.GetSet()->ClearSlice(false);
00456 planeItr.ResetFirst();
00457
00458 MSG("DmxUtil", Msg::kDebug) << "\tvertexPlaneNumber = " << vertexPlaneNumber << endl;
00459 return vertexPlaneNumber;
00460 }
|
|
|
Definition at line 310 of file DmxUtilities.cxx. Referenced by DmxDeMuxModule::Ana(), DmxDeMuxCosmicsModule::Ana(), FillHitPixels(), and FillPlaneArray(). 00311 {
00312
00313 //find the pixel that the digit came from. make stupid arrays
00314 //to hold the pixel to va channel info, as for some reason it
00315 //isnt available from the framework
00316 Int_t vaChip0[] = {-1, -1, 14, 0, 15, 1, 10, 4, 11,
00317 5, 6, 8, 7, 9, 2, 13, 3, 12};
00318 //Int_t vaChip1[] = {-1, -1, 22, 8, 23, 9, 18, 12, 19,
00319 // 13, 14, 16, 15, 17, 10, 21, 11, 20};
00320 //Int_t vaChip2[] = {-1, -1, 6, 16, 7, 17, 2, 20, 3, 21,
00321 // 22, 0, 23, 1, 18, 5, 19, 4};
00322
00323 Int_t pixel = vaChip0[channel];
00324 //if(chip == 0){ pixel = vaChip0[channel];}
00325 //if(chip == 1){ pixel = vaChip1[channel];}
00326 //if(chip == 2){ pixel = vaChip2[channel];}
00327
00328 //MSG("Dmx2", Msg::kDebug)<< "pixel = " << pixel << endl;
00329 return pixel;
00330 }
|
|
||||||||||||||||||||||||
|
Definition at line 980 of file DmxUtilities.cxx. References DmxPlane::GetPlaneNumber(), and DmxPlane::GetPlaneType(). Referenced by AlgDeMuxGolden::RunAlg(), and AlgDeMuxCosmics::RunAlg(). 00983 {
00984
00985 planeItr.ResetFirst();
00986 DmxPlane *plane = 0;
00987
00988 //make a histogram for each supermodule. this is because the method is trickier in
00989 //tpos-z space than strip-plane space, so just deal with the gap by pretending its not
00990 //there
00991 TH2F houghSpaceSM1 = TH2F("houghSpaceSM1", "Hough Space for Multiples",
00992 30, -30, 30, 500, -500, 500);
00993 TH2F houghSpaceSM2 = TH2F("houghSpaceSM2", "Hough Space for Multiples",
00994 30, -30, 30, 500, -500, 500);
00995 TH2F houghSpaceCutSM1 = TH2F("hS75SM1", "Hough Space for U with cut",
00996 30, -30, 30, 500, -500, 500);
00997 TH2F houghSpaceCutSM2 = TH2F("hS75SM2", "Hough Space for U with cut",
00998 30, -30, 30, 500, -500, 500);
00999
01000 Float_t planeZ = 0.;
01001 Float_t cog = 0.;
01002 while( (plane = planeItr()) ){
01003
01004 planeZ = 1.*plane->GetPlaneNumber();
01005
01006 if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
01007
01008 for(Int_t j = 0; j <3; j++){
01009 if( j==0 && dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis() ){
01010 cog = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
01011 }
01012 else if( j==1 && dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis() ){
01013 cog = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
01014 }
01015 else if( j==2 && dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis() ){
01016 cog = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
01017 }
01018 //loop over possible slopes
01019 for(Int_t m = -29; m<30; m+=2){
01020
01021 Float_t c = cog - m*1.*(planeZ - vertexZ);
01022
01023 if(planeZ<249) houghSpaceSM1.Fill(m,c);
01024 else if(planeZ>249) houghSpaceSM2.Fill(m,c);
01025
01026 }//end loop over slopes
01027 }//end loop over hypotheses
01028
01029 }//end if shower plane
01030 else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
01031 cog = dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG();
01032
01033 //loop over possible slopes
01034 for(Int_t m = -29; m<30; m+=2){
01035
01036 Float_t c = cog - m*1.*(planeZ - vertexZ);
01037
01038 //fill the space with this entry 3 times so it is weighted the same as a shower plane
01039 if(planeZ<249){
01040 houghSpaceSM1.Fill(m,c);
01041 houghSpaceSM1.Fill(m,c);
01042 houghSpaceSM1.Fill(m,c);
01043 }
01044 else if(planeZ>249){
01045 houghSpaceSM2.Fill(m,c);
01046 houghSpaceSM2.Fill(m,c);
01047 houghSpaceSM2.Fill(m,c);
01048 }
01049 }//end loop over slopes
01050
01051 }//end if muon plane
01052 }//end loop over planes
01053
01054 planeItr.ResetFirst();
01055
01056 //now fill a second histogram with the bins that are at least peakCutOff of the maximum
01057 //bin in the previous histogram
01058 Stat_t binContent = 0;
01059 Double_t maximumSM1 = houghSpaceSM1.GetMaximum();
01060 Double_t maximumSM2 = houghSpaceSM2.GetMaximum();
01061
01062 for(Int_t x = 1; x < 30+1; x++){
01063 //cout << x << endl;
01064 for(Int_t y = 1; y < 500+1; y++){
01065 binContent = houghSpaceSM1.GetBinContent(x,y);
01066 if(binContent>= peakCutOff*maximumSM1) houghSpaceCutSM1.SetBinContent(x,y,binContent);
01067 binContent = houghSpaceSM2.GetBinContent(x,y);
01068 if(binContent>= peakCutOff*maximumSM2) houghSpaceCutSM2.SetBinContent(x,y,binContent);
01069
01070 }
01071 }
01072
01073 //if the rms in the slope is less than the value determined by tests
01074 //and the rms in the intercept is greater than the cutoff, mark the event
01075 //as a possible overlapping multiple
01076 bool multiple = false;
01077 if(houghSpaceCutSM1.GetRMS(1)<=slopeCutOff
01078 && houghSpaceCutSM1.GetRMS(2)>=interceptCutOff)
01079 multiple = true;
01080 else if(houghSpaceCutSM2.GetRMS(1)<=slopeCutOff
01081 && houghSpaceCutSM2.GetRMS(2)>=interceptCutOff)
01082 multiple = true;
01083
01084 return multiple;
01085 }
|
|
||||||||||||||||||||||||
|
Definition at line 726 of file DmxUtilities.cxx. References DmxPlane::GetZPosition(). 00728 {
00729 Int_t nonPhysicalCtr = 0;
00730 bool validFit = true;
00731 Int_t planesInEvent = planeItr.SizeSelect();
00732 //MSG("DmxUtil", Msg::kDebug)<<"\tin IsValidFit " << a << "\t" << b << endl;
00733 planeItr.Reset();
00734
00735 DmxPlane *plane = 0;
00736 Double_t fitCoG = 0.;
00737
00738 while( planeItr.IsValid() ){
00739
00740 plane = planeItr.Ptr();
00741
00742 fitCoG = (a1 + a2*(plane->GetZPosition())
00743 + a3*TMath::Power(plane->GetZPosition(),2)
00744 + a4*TMath::Power(plane->GetZPosition(),3));
00745 //MSG("DmxUtil", Msg::kDebug)<<"\t" << plane->GetZPosition() << endl;
00746
00747 //if the fit CoG for the plane is outside of the detector, it is non Physical
00748 if( TMath::Abs(fitCoG) >= 4.5) ++nonPhysicalCtr;
00749
00750 planeItr.Next();
00751 }
00752
00753 planeItr.Reset();
00754
00755 if(planesInEvent<8 && (1.*nonPhysicalCtr)/(1.*planesInEvent)>0.33) validFit = false;
00756 else if(planesInEvent>=8 && nonPhysicalCtr>=4) validFit = false;
00757
00758 return validFit;
00759 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 689 of file DmxUtilities.cxx. References DmxPlane::GetZPosition(). Referenced by AlgDeMuxCosmics::FindCosmicSolution(), and AlgDeMuxBeam::FindFit(). 00691 {
00692 Int_t nonPhysicalCtr = 0;
00693 bool validFit = true;
00694 Int_t planesInEvent = planeItr.SizeSelect();
00695 //MSG("DmxUtil", Msg::kDebug)<<"\tin IsValidFit " << a << "\t" << b << endl;
00696 planeItr.Reset();
00697
00698 DmxPlane *plane = 0;
00699 Double_t fitCoG = 0.;
00700
00701 while( planeItr.IsValid() ){
00702
00703 plane = planeItr.Ptr();
00704
00705 fitCoG = (a1 + a2*(plane->GetZPosition()-offset)
00706 + a3*TMath::Power(plane->GetZPosition()-offset,2)
00707 + a4*TMath::Power(plane->GetZPosition()-offset,3));
00708 //MSG("DmxUtil", Msg::kDebug)<<"\t" << plane->GetZPosition() << endl;
00709
00710 //if the fit CoG for the plane is outside of the detector, it is non Physical
00711 if( TMath::Abs(fitCoG) >= 4.5) ++nonPhysicalCtr;
00712
00713 planeItr.Next();
00714 }
00715
00716 planeItr.Reset();
00717
00718 if(planesInEvent<8 && (1.*nonPhysicalCtr)/(1.*planesInEvent)>0.33) validFit = false;
00719 else if(planesInEvent>=8 && nonPhysicalCtr>=4) validFit = false;
00720
00721 return validFit;
00722 }
|
1.3.9.1