[807] | 1 | // %%%%%%%%%%%%% |
---|
| 2 | // gemc headers |
---|
| 3 | // %%%%%%%%%%%%% |
---|
| 4 | #include "DC_hitprocess.h" |
---|
| 5 | #include "CLHEP/Random/RandGaussT.h" |
---|
| 6 | |
---|
| 7 | PH_output DC_HitProcess :: ProcessHit(MHit* aHit, gemc_opts Opt) |
---|
| 8 | { |
---|
| 9 | string hd_msg = Opt.args["LOG_MSG"].args + " DC Hit Process " ; |
---|
| 10 | double HIT_VERBOSITY = Opt.args["HIT_VERBOSITY"].arg; |
---|
| 11 | |
---|
| 12 | PH_output out; |
---|
| 13 | out.identity = aHit->GetId(); |
---|
| 14 | |
---|
| 15 | HCname = "DC Hit Process"; |
---|
| 16 | // %%%%%%%%%%%%%%%%%%% |
---|
| 17 | // Raw hit information |
---|
| 18 | // %%%%%%%%%%%%%%%%%%% |
---|
| 19 | int nsteps = aHit->GetPos().size(); |
---|
| 20 | |
---|
| 21 | |
---|
| 22 | |
---|
| 23 | // Identifying the fastest track in the hit |
---|
| 24 | int trackId = -1; |
---|
| 25 | double minTime = 10000; |
---|
| 26 | vector<int> stepTrackId = aHit->GetTrackId(); |
---|
| 27 | vector<double> stepTime = aHit->GetTime(); |
---|
| 28 | vector<G4double> Edep = aHit->GetEdep(); |
---|
| 29 | |
---|
| 30 | for(int s=0; s<nsteps; s++) |
---|
| 31 | { |
---|
| 32 | if(stepTrackId[s] != trackId && stepTime[s]/ns < minTime && Edep[s] >= aHit->GetThreshold()) |
---|
| 33 | trackId = stepTrackId[s]; |
---|
| 34 | } |
---|
| 35 | |
---|
| 36 | // If no step pass the threshold, getting the fastest track |
---|
| 37 | if(trackId == -1) |
---|
| 38 | { |
---|
| 39 | for(int s=0; s<nsteps; s++) |
---|
| 40 | { |
---|
| 41 | if(stepTrackId[s] != trackId && stepTime[s]/ns < minTime) |
---|
| 42 | trackId = stepTrackId[s]; |
---|
| 43 | } |
---|
| 44 | } |
---|
| 45 | |
---|
| 46 | // Get Total Energy deposited of the fastest track |
---|
| 47 | double Etot = 0; |
---|
| 48 | int count_track_steps = 0; |
---|
| 49 | for(int s=0; s<nsteps; s++) |
---|
| 50 | if(stepTrackId[s] == trackId) |
---|
| 51 | { |
---|
| 52 | Etot = Etot + Edep[s]; |
---|
| 53 | count_track_steps++; |
---|
| 54 | } |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | // average global, local positions of fastest track |
---|
| 58 | double x, y, z; |
---|
| 59 | double lx, ly, lz; |
---|
| 60 | x = y = z = lx = ly = lz = 0; |
---|
| 61 | vector<G4ThreeVector> pos = aHit->GetPos(); |
---|
| 62 | vector<G4ThreeVector> Lpos = aHit->GetLPos(); |
---|
| 63 | |
---|
| 64 | /* if(Etot>0)*/ |
---|
| 65 | for(int s=0; s<nsteps; s++) |
---|
| 66 | { |
---|
| 67 | if(stepTrackId[s] == trackId) |
---|
| 68 | { |
---|
| 69 | /* x = x + pos[s].x()*Edep[s]/Etot; |
---|
| 70 | y = y + pos[s].y()*Edep[s]/Etot; |
---|
| 71 | z = z + pos[s].z()*Edep[s]/Etot; |
---|
| 72 | lx = lx + Lpos[s].x()*Edep[s]/Etot; |
---|
| 73 | ly = ly + Lpos[s].y()*Edep[s]/Etot; |
---|
| 74 | lz = lz + Lpos[s].z()*Edep[s]/Etot;*/ |
---|
| 75 | x = x + pos[s].x()/count_track_steps; |
---|
| 76 | y = y + pos[s].y()/count_track_steps; |
---|
| 77 | z = z + pos[s].z()/count_track_steps; |
---|
| 78 | lx = lx + Lpos[s].x()/count_track_steps; |
---|
| 79 | ly = ly + Lpos[s].y()/count_track_steps; |
---|
| 80 | lz = lz + Lpos[s].z()/count_track_steps; |
---|
| 81 | } |
---|
| 82 | } |
---|
| 83 | // else |
---|
| 84 | // { |
---|
| 85 | // x = pos[0].x(); |
---|
| 86 | // y = pos[0].y(); |
---|
| 87 | // z = pos[0].z(); |
---|
| 88 | // lx = Lpos[0].x(); |
---|
| 89 | // ly = Lpos[0].y(); |
---|
| 90 | // lz = Lpos[0].z(); |
---|
| 91 | // } |
---|
| 92 | |
---|
| 93 | // average time of fastest track |
---|
| 94 | double time = 0; |
---|
| 95 | vector<G4double> times = aHit->GetTime(); |
---|
| 96 | for(int s=0; s<nsteps; s++) if(stepTrackId[s] == trackId) time = time + times[s]/nsteps; |
---|
| 97 | |
---|
| 98 | // Energy of the track |
---|
| 99 | double Ene = aHit->GetE(); |
---|
| 100 | |
---|
| 101 | out.raws.push_back(Etot); |
---|
| 102 | out.raws.push_back(x); |
---|
| 103 | out.raws.push_back(y); |
---|
| 104 | out.raws.push_back(z); |
---|
| 105 | out.raws.push_back(lx); |
---|
| 106 | out.raws.push_back(ly); |
---|
| 107 | out.raws.push_back(lz); |
---|
| 108 | out.raws.push_back(time); |
---|
| 109 | out.raws.push_back((double) aHit->GetPID()); |
---|
| 110 | out.raws.push_back(aHit->GetVert().getX()); |
---|
| 111 | out.raws.push_back(aHit->GetVert().getY()); |
---|
| 112 | out.raws.push_back(aHit->GetVert().getZ()); |
---|
| 113 | out.raws.push_back(Ene); |
---|
| 114 | out.raws.push_back((double) aHit->GetmPID()); |
---|
| 115 | out.raws.push_back(aHit->GetmVert().getX()); |
---|
| 116 | out.raws.push_back(aHit->GetmVert().getY()); |
---|
| 117 | out.raws.push_back(aHit->GetmVert().getZ()); |
---|
| 118 | |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 122 | // Distance 1: smear by 300 um |
---|
| 123 | // Drift Time 1: using 50 um/ns |
---|
| 124 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 125 | |
---|
| 126 | // Finding Wire position |
---|
| 127 | double NWIRES = 112; |
---|
| 128 | double ylength = aHit->GetDetector().dimensions[3]; |
---|
| 129 | double deltay = 2.0*ylength/NWIRES; |
---|
| 130 | vector<identifier> iden = aHit->GetId(); |
---|
| 131 | int nwire = iden[3].id; |
---|
| 132 | double WIRE_Y = nwire*deltay; |
---|
| 133 | |
---|
| 134 | // Finding DOCA |
---|
| 135 | double doca = 10000; |
---|
| 136 | double LR = 0; |
---|
| 137 | |
---|
| 138 | for(int s=0; s<nsteps; s++) |
---|
| 139 | { |
---|
| 140 | G4ThreeVector DOCA(0, Lpos[s].y() + ylength - WIRE_Y, Lpos[s].z()); |
---|
| 141 | if(DOCA.mag() <= doca && stepTrackId[s] == trackId ) |
---|
| 142 | { |
---|
| 143 | doca = DOCA.mag(); |
---|
| 144 | if(DOCA.y() >=0 ) LR = 1; |
---|
| 145 | else LR = -1; |
---|
| 146 | |
---|
| 147 | } |
---|
| 148 | } |
---|
| 149 | // if(aHit->GetPID() == 2112) |
---|
| 150 | // cout << " wire=" << nwire << " layer=" << iden[2].id |
---|
| 151 | // << " doca=" << doca << " nsteps=" << count_track_steps << endl; |
---|
| 152 | |
---|
| 153 | |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 157 | // Digitization |
---|
| 158 | // DC ID: |
---|
| 159 | // sector, SL, Layer, Wire |
---|
| 160 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 161 | int sector = out.identity[0].id; |
---|
| 162 | int SL = out.identity[1].id; |
---|
| 163 | int Layer = out.identity[2].id; |
---|
| 164 | int Wire = out.identity[3].id; |
---|
| 165 | |
---|
| 166 | double drift_velocity = 0; |
---|
| 167 | if(SL == 1 || SL == 2) drift_velocity = 0.053; |
---|
| 168 | if(SL == 3 || SL == 4) drift_velocity = 0.026; |
---|
| 169 | if(SL == 5 || SL == 6) drift_velocity = 0.036; |
---|
| 170 | |
---|
| 171 | double sdoca = -1; |
---|
| 172 | while(sdoca < 0) sdoca = CLHEP::RandGauss::shoot(doca, 0.3); |
---|
| 173 | double time1 = doca/drift_velocity; |
---|
| 174 | double stime1 = sdoca/drift_velocity; |
---|
| 175 | |
---|
| 176 | out.raws.push_back(LR); |
---|
| 177 | out.raws.push_back(doca); |
---|
| 178 | out.raws.push_back(sdoca); |
---|
| 179 | out.raws.push_back(time1); |
---|
| 180 | out.raws.push_back(stime1); |
---|
| 181 | |
---|
| 182 | |
---|
| 183 | out.dgtz.push_back(sector); |
---|
| 184 | out.dgtz.push_back(SL); |
---|
| 185 | out.dgtz.push_back(Layer); |
---|
| 186 | out.dgtz.push_back(Wire); |
---|
| 187 | |
---|
| 188 | if(HIT_VERBOSITY>4) |
---|
| 189 | cout << hd_msg << " Sector: " << sector |
---|
| 190 | << " Superlayer: " << SL |
---|
| 191 | << " Layer: " << Layer |
---|
| 192 | << " Wire: " << Wire |
---|
| 193 | << " L/R: " << LR |
---|
| 194 | << " doca: " << doca |
---|
| 195 | << " sdoca: " << sdoca |
---|
| 196 | << " time1" << time1 |
---|
| 197 | << " stime1=" << stime1 << endl; |
---|
| 198 | |
---|
| 199 | return out; |
---|
| 200 | } |
---|
| 201 | |
---|
| 202 | vector<identifier> DC_HitProcess :: ProcessID(vector<identifier> id, G4Step* aStep, detector Detector) |
---|
| 203 | { |
---|
| 204 | double NWIRES = 112; |
---|
| 205 | |
---|
| 206 | G4StepPoint *prestep = aStep->GetPreStepPoint(); |
---|
| 207 | G4StepPoint *poststep = aStep->GetPostStepPoint(); |
---|
| 208 | G4VTouchable* TH = (G4VTouchable*) aStep->GetPreStepPoint()->GetTouchable(); |
---|
| 209 | |
---|
| 210 | string name = TH->GetVolume(0)->GetName(); |
---|
| 211 | G4ThreeVector xyz = poststep->GetPosition(); |
---|
| 212 | G4ThreeVector Lxyz = prestep->GetTouchableHandle()->GetHistory() |
---|
| 213 | ->GetTopTransform().TransformPoint(xyz); |
---|
| 214 | |
---|
| 215 | double ylength = Detector.dimensions[3]; |
---|
| 216 | double deltay = 2.0*ylength/NWIRES; |
---|
| 217 | double loc_y = Lxyz.y() + ylength; |
---|
| 218 | |
---|
| 219 | int nwire = (int) floor(loc_y/deltay); |
---|
| 220 | vector<identifier> yid = id; |
---|
| 221 | |
---|
| 222 | if(nwire > NWIRES) cout << " SuperLayer: " << yid[1].id << " layer: " << yid[2].id |
---|
| 223 | << " wire: " << nwire << " length: " << ylength |
---|
| 224 | << " ly: " << Lxyz.y() << " deltay: " << deltay |
---|
| 225 | << " loc_y: " << loc_y << " nwire*deltay: " << fabs( nwire*deltay - loc_y ) << endl; |
---|
| 226 | |
---|
| 227 | |
---|
| 228 | yid[3].id = nwire; |
---|
| 229 | |
---|
| 230 | if(fabs( (nwire+1)*deltay - loc_y ) < fabs( nwire*deltay - loc_y ) ) |
---|
| 231 | yid[3].id = nwire + 1; |
---|
| 232 | |
---|
| 233 | return yid; |
---|
| 234 | } |
---|
| 235 | |
---|
| 236 | |
---|
| 237 | |
---|
| 238 | |
---|
| 239 | |
---|