00001 #include "TMath.h" 00002 00003 #include "MessageService/MsgService.h" 00004 00005 #include "AtNuEvent/AtmosEvent.h" 00006 #include "AtNuEvent/AtmosScintHit.h" 00007 #include "AtNuEvent/AtmosStrip.h" 00008 #include "AtNuEvent/AtmosMC.h" 00009 00010 #include "AtNuUtils/UtilStrip.h" 00011 00012 CVSID("$Id: UtilStrip.cxx,v 1.2 2006/07/19 22:41:29 gmieg Exp $"); 00013 00014 const double C45 = TMath::Sqrt(2.) / 2.; //Cosine of 45 degrees 00015 00016 std::vector<int> UtilStrip::GangedStripsNegative(int strip) { 00017 std::vector<int> Gang; 00018 if(strip<0 || strip>191) return Gang; 00019 Gang.push_back(strip); 00020 00021 int AddStrip = strip; 00022 //Climb Up 00023 while (true) { 00024 if(AddStrip%24 == 0) AddStrip += 47; 00025 else AddStrip += 23; 00026 if(AddStrip > 191) break; 00027 Gang.push_back(AddStrip); 00028 } 00029 00030 //Climb Down 00031 AddStrip = strip; 00032 while (true) { 00033 if(AddStrip<47) AddStrip -= 23; 00034 else if((AddStrip-47)%24 == 0) AddStrip -= 47; 00035 else AddStrip -= 23; 00036 if(AddStrip < 0) break; 00037 Gang.push_back(AddStrip); 00038 } 00039 00040 return Gang; 00041 } 00042 00043 std::vector<int> UtilStrip::GangedStripsPositive(int strip) { 00044 std::vector<int> Gang; 00045 if(strip<0 || strip>191) return Gang; 00046 Gang.push_back(strip); 00047 00048 int AddStrip = strip; 00049 while (true) { 00050 AddStrip += 24; 00051 if(AddStrip > 191) AddStrip -= 192; 00052 if(AddStrip == strip) break; 00053 Gang.push_back(AddStrip); 00054 } 00055 00056 return Gang; 00057 } 00058 00059 double UtilStrip::SumEnergy(std::set<SetStrip> strips) { 00060 double sumenergy = 0.0; 00061 std::set<SetStrip>::iterator itr = strips.begin(); 00062 for(; itr!=strips.end(); itr++) sumenergy += itr->Energy; 00063 return sumenergy; 00064 } 00065 00066 int UtilStrip::NStrips(std::set<SetStrip> strips) { 00067 int nstrips = 0; 00068 std::set<SetStrip>::iterator itr = strips.begin(); 00069 for(; itr!=strips.end(); itr++) if(itr->Energy>0) nstrips++; 00070 return nstrips; 00071 } 00072 00073 int UtilStrip::NDigits(std::set<SetStrip> strips) { 00074 int ndigits = 0; 00075 std::set<SetStrip>::iterator itr = strips.begin(); 00076 for(; itr!=strips.end(); itr++) 00077 if(itr->Energy>0 && itr->DoubleEnded) ndigits++; 00078 return ndigits; 00079 } 00080 00081 /* 00082 * SetStrip class 00083 */ 00084 SetStrip::SetStrip(int i, const AtmosStrip *strip) { 00085 TCIndex = i; 00086 00087 #ifdef PECORR 00088 Energy = strip->QPEcorr[0] + strip->QPEcorr[1]; 00089 #else 00090 Energy = strip->QPE[0] + strip->QPE[1]; 00091 #endif 00092 00093 DoubleEnded = (bool)(strip->Ndigits - 1); 00094 } 00095 00096 /* 00097 * StripCluster class 00098 */ 00099 void StripCluster::Clear() { 00100 NStrips = 0; 00101 SumU = SumV = SumZ = SumSqU = SumSqV = SumSqZ = 0.0; 00102 00103 SumPE = 0.0; 00104 WSumU = WSumV = WSumZ = WSumSqU = WSumSqV = WSumSqZ = 0.0; 00105 00106 Strips.clear(); 00107 } 00108 00109 void StripCluster::AddStripDList(const StripDList &sdl) { 00110 NStrips++; 00111 SumPE += sdl.PE; 00112 00113 SumU += sdl.U; 00114 SumV += sdl.V; 00115 SumZ += sdl.Z; 00116 00117 SumSqU += sdl.U * sdl.U; 00118 SumSqV += sdl.V * sdl.V; 00119 SumSqZ += sdl.Z * sdl.Z; 00120 00121 WSumU += sdl.U * sdl.PE; 00122 WSumV += sdl.V * sdl.PE; 00123 WSumZ += sdl.Z * sdl.PE; 00124 00125 WSumSqU += sdl.U * sdl.U * sdl.PE; 00126 WSumSqV += sdl.V * sdl.V * sdl.PE; 00127 WSumSqZ += sdl.Z * sdl.Z * sdl.PE; 00128 00129 Strips.insert(SetStrip(sdl.TCIndex,false,sdl.PE)); 00130 } 00131 00132 double StripCluster::RMSU() { 00133 return TMath::Sqrt(SumSqU/(double)NStrips + (CentU()*CentU())); 00134 } 00135 00136 double StripCluster::RMSV() { 00137 return TMath::Sqrt(SumSqV/(double)NStrips + (CentV()*CentV())); 00138 } 00139 00140 double StripCluster::RMSZ() { 00141 return TMath::Sqrt(SumSqZ/(double)NStrips + (CentZ()*CentZ())); 00142 } 00143 00144 StripDList::StripDList(int i) { 00145 TCIndex = i; 00146 InCluster = false; 00147 ClusterCount = 0; 00148 U = V = Z = 0.0; 00149 } 00150 00151 bool StripDList::SetInCluster() { 00152 ClusterCount++; 00153 if(InCluster) return(true); 00154 InCluster = true; 00155 return(false); 00156 } 00157 00158 bool operator<(const StripDList &sdl1, const StripDList &sdl2) { 00159 if(sdl1.VecDStrip.size() == 0 || 00160 sdl2.VecDStrip.size() == 0) return false; 00161 return(sdl1.VecDStrip[0].Disp < sdl2.VecDStrip[0].Disp); 00162 } 00163 00164 bool operator==(const StripDList &sdl1, const StripDList &sdl2) { 00165 return(sdl1.TCIndex == sdl2.TCIndex); 00166 } 00167 00168 bool operator!=(const StripDList &sdl1, const StripDList &sdl2) { 00169 return(sdl1.TCIndex != sdl2.TCIndex); 00170 }