| File: | /home/sdobbs/work/clang/halld_recon/src/libraries/BCAL/DBCALClump_factory.cc |
| Warning: | line 310, column 7 Value stored to 'Mt' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | /* |
| 2 | * DBCALClump.h |
| 3 | * |
| 4 | * Created by Beni Zihlmann Tue Mar 12 2013 |
| 5 | * |
| 6 | */ |
| 7 | |
| 8 | |
| 9 | #include <JANA/JApplication.h> |
| 10 | using namespace jana; |
| 11 | |
| 12 | #include <DBCALClump_factory.h> |
| 13 | |
| 14 | |
| 15 | //---------------- |
| 16 | // init |
| 17 | //---------------- |
| 18 | jerror_t DBCALClump_factory::init(void) |
| 19 | { |
| 20 | return NOERROR; |
| 21 | } |
| 22 | |
| 23 | //---------------- |
| 24 | // init |
| 25 | //---------------- |
| 26 | jerror_t DBCALClump_factory::brun(JEventLoop *loop, int32_t runnumber) |
| 27 | { |
| 28 | map<string, double> bcalparms; |
| 29 | |
| 30 | if ( !loop->GetCalib("BCAL/mc_parms", bcalparms)){ |
| 31 | cout<<"DBCALClump_factory: loading values from TOF data base"<<endl; |
| 32 | } else { |
| 33 | cout << "DBCALClumpo_factory: Error loading values from BCAL MC data base" <<endl; |
| 34 | |
| 35 | VELOCITY = 15.; // set to some reasonable value |
| 36 | |
| 37 | return NOERROR; |
| 38 | } |
| 39 | |
| 40 | VELOCITY = bcalparms["C_EFFECTIVE"]; |
| 41 | |
| 42 | |
| 43 | return NOERROR; |
| 44 | } |
| 45 | |
| 46 | |
| 47 | |
| 48 | //---------------- |
| 49 | // evnt |
| 50 | //---------------- |
| 51 | jerror_t DBCALClump_factory::evnt(JEventLoop *loop, uint64_t eventnumber) { |
| 52 | |
| 53 | const DBCALHit *BcalMatrixU[48*4][4]; // Up stream |
| 54 | const DBCALHit *BcalMatrixD[48*4][4]; // Down stream |
| 55 | |
| 56 | int MatrixUD[48*4][4]; |
| 57 | |
| 58 | float EU[48*4]; |
| 59 | float ED[48*4]; |
| 60 | int EUi[48*4]; |
| 61 | int EDi[48*4]; |
| 62 | |
| 63 | vector<const DBCALHit*> AllBcalHits; |
| 64 | loop->Get(AllBcalHits); |
| 65 | |
| 66 | for (int i=0;i<48*4;i++){ |
| 67 | EU[i] = 0.; |
| 68 | ED[i] = 0.; |
| 69 | EUi[i] = 0; |
| 70 | EDi[i] = 0; |
| 71 | for (int m=0;m<4;m++){ |
| 72 | BcalMatrixU[i][m] = NULL__null; |
| 73 | BcalMatrixD[i][m] = NULL__null; |
| 74 | MatrixUD[i][m] = 0; |
| 75 | } |
| 76 | } |
| 77 | |
| 78 | // setup BCAL Matrix |
| 79 | |
| 80 | for (int i=0; i< (int)AllBcalHits.size();i++){ |
| 81 | const DBCALHit* hit = AllBcalHits[i]; |
| 82 | //cout<<"Module: "<<hit->module<<" Sector: "<<hit->sector<<" Layer: "<<hit->layer<<endl; |
| 83 | int chan = hit->sector-1 + (hit->module-1)*4; |
| 84 | if ((hit->t<200.) && (hit->t>0.)) { // HARD CODED |
| 85 | if (hit->end == DBCALGeometry::kUpstream){ |
| 86 | BcalMatrixU[chan][hit->layer-1] = hit; |
| 87 | } else{ |
| 88 | BcalMatrixD[chan][hit->layer-1] = hit; |
| 89 | } |
| 90 | } |
| 91 | } |
| 92 | |
| 93 | // Build the up/down stream coincidence Matrix |
| 94 | for (int i=0;i<48*4;i++){ |
| 95 | for (int m=0;m<4;m++){ |
| 96 | if ((BcalMatrixU[i][m]) && (BcalMatrixD[i][m])){ |
| 97 | float meant = (BcalMatrixU[i][m]->t + BcalMatrixD[i][m]->t - 390./16.75)/2.; // HARD CODED VALUES |
| 98 | // this mean time has to be positive!! |
| 99 | // and the cut could be made stronger. |
| 100 | if (meant>0){ |
| 101 | MatrixUD[i][m] = 1; |
| 102 | } else { |
| 103 | BcalMatrixU[i][m] = NULL__null; |
| 104 | BcalMatrixD[i][m] = NULL__null; |
| 105 | } |
| 106 | } |
| 107 | } |
| 108 | } |
| 109 | |
| 110 | |
| 111 | // first get rid of single-ended isolated hits |
| 112 | for (int i=0;i<48*4;i++){ |
| 113 | int bef = i-1; |
| 114 | int aft = i+1; |
| 115 | if (bef<0){ |
| 116 | bef+=48*4; |
| 117 | } |
| 118 | if (aft>=48*4){ |
| 119 | aft -=48*4; |
| 120 | } |
| 121 | // layer 1; |
| 122 | if (!MatrixUD[i][0]){ |
| 123 | if (BcalMatrixU[i][0]){ |
| 124 | if ((!BcalMatrixU[bef][0]) && |
| 125 | (!BcalMatrixU[aft][0]) && |
| 126 | (!BcalMatrixU[i][1]) && |
| 127 | (!BcalMatrixU[aft][1]) && |
| 128 | (!BcalMatrixU[aft][1])){ |
| 129 | BcalMatrixU[i][0] = NULL__null; |
| 130 | } |
| 131 | } |
| 132 | if (BcalMatrixD[i][0]){ |
| 133 | if ((!BcalMatrixD[bef][0]) && |
| 134 | (!BcalMatrixD[aft][0]) && |
| 135 | (!BcalMatrixD[i][1]) && |
| 136 | (!BcalMatrixD[aft][1]) && |
| 137 | (!BcalMatrixD[aft][1])){ |
| 138 | BcalMatrixD[i][0] = NULL__null; |
| 139 | } |
| 140 | } |
| 141 | } |
| 142 | // layer 2; |
| 143 | if (!MatrixUD[i][1]){ |
| 144 | if (BcalMatrixU[i][1]){ |
| 145 | if ((!BcalMatrixU[bef][1]) && |
| 146 | (!BcalMatrixU[aft][1]) && |
| 147 | (!BcalMatrixU[i][2]) && |
| 148 | (!BcalMatrixU[aft][2]) && |
| 149 | (!BcalMatrixU[aft][2])){ |
| 150 | BcalMatrixU[i][1] = NULL__null; |
| 151 | } |
| 152 | } |
| 153 | if (BcalMatrixD[i][1]){ |
| 154 | if ((!BcalMatrixD[bef][1]) && |
| 155 | (!BcalMatrixD[aft][1]) && |
| 156 | (!BcalMatrixD[i][2]) && |
| 157 | (!BcalMatrixD[aft][2]) && |
| 158 | (!BcalMatrixD[aft][2])){ |
| 159 | BcalMatrixD[i][1] = NULL__null; |
| 160 | } |
| 161 | } |
| 162 | } |
| 163 | // layer 3; |
| 164 | if (!MatrixUD[i][2]){ |
| 165 | if (BcalMatrixU[i][2]){ |
| 166 | if ((!BcalMatrixU[bef][2]) && |
| 167 | (!BcalMatrixU[aft][2]) && |
| 168 | (!BcalMatrixU[i][3]) && |
| 169 | (!BcalMatrixU[aft][3]) && |
| 170 | (!BcalMatrixU[aft][3])){ |
| 171 | BcalMatrixU[i][2] = NULL__null; |
| 172 | } |
| 173 | } |
| 174 | if (BcalMatrixD[i][2]){ |
| 175 | if ((!BcalMatrixD[bef][2]) && |
| 176 | (!BcalMatrixD[aft][2]) && |
| 177 | (!BcalMatrixD[i][3]) && |
| 178 | (!BcalMatrixD[aft][3]) && |
| 179 | (!BcalMatrixD[aft][3])){ |
| 180 | BcalMatrixD[i][2] = NULL__null; |
| 181 | } |
| 182 | } |
| 183 | } |
| 184 | // layer 4; |
| 185 | if (!MatrixUD[i][3]){ |
| 186 | if (BcalMatrixU[i][3]){ |
| 187 | if ((!BcalMatrixU[bef][3]) && |
| 188 | (!BcalMatrixU[aft][3]) && |
| 189 | (!BcalMatrixU[i][2]) && |
| 190 | (!BcalMatrixU[aft][2]) && |
| 191 | (!BcalMatrixU[aft][2])){ |
| 192 | BcalMatrixU[i][3] = NULL__null; |
| 193 | } |
| 194 | } |
| 195 | if (BcalMatrixD[i][3]){ |
| 196 | if ((!BcalMatrixD[bef][3]) && |
| 197 | (!BcalMatrixD[aft][3]) && |
| 198 | (!BcalMatrixD[i][2]) && |
| 199 | (!BcalMatrixD[aft][2]) && |
| 200 | (!BcalMatrixD[aft][2])){ |
| 201 | BcalMatrixD[i][3] = NULL__null; |
| 202 | } |
| 203 | } |
| 204 | } |
| 205 | } |
| 206 | |
| 207 | for (int i=0;i<48*4;i++){ |
| 208 | for (int m=0;m<4;m++){ |
| 209 | |
| 210 | if (BcalMatrixU[i][m]){ |
| 211 | EU[i] += BcalMatrixU[i][m]->E; |
| 212 | EUi[i] += 1; |
| 213 | } |
| 214 | if (BcalMatrixD[i][m]){ |
| 215 | ED[i] += BcalMatrixD[i][m]->E; |
| 216 | EDi[i] += 1; |
| 217 | } |
| 218 | |
| 219 | } |
| 220 | } |
| 221 | |
| 222 | // clean up Matrix from obvious stray hits |
| 223 | for (int i=1;i<48*4-1;i++) { |
| 224 | int a = EDi[i-1]+EUi[i-1]; |
| 225 | int b = EDi[i] + EUi[i]; |
| 226 | int c = EDi[i+1]+EUi[i+1]; |
| 227 | if ((b<2)&&(a==0)&&(c==0)){ |
| 228 | EDi[i] = 0; |
| 229 | EUi[i] = 0; |
| 230 | EU[i] = 0; |
| 231 | ED[i] = 0; |
| 232 | } |
| 233 | } |
| 234 | |
| 235 | int CNT=1; |
| 236 | while (CNT) { |
| 237 | |
| 238 | //cout<<"Find Next Clump"<<endl; |
| 239 | int MaxU=999; |
| 240 | int MaxD=999; |
| 241 | |
| 242 | float mU=0.; |
| 243 | float mD=0.; |
| 244 | // find largest energy deposition in the matrix |
| 245 | for (int i=0;i<48*4;i++){ |
| 246 | if ((EU[i]>0.0)&&(ED[i]>0.0)){ |
| 247 | if (EU[i]>mU){ |
| 248 | mU = EU[i]; |
| 249 | MaxU = i; |
| 250 | } |
| 251 | if (ED[i]>mD){ |
| 252 | mD = ED[i]; |
| 253 | MaxD = i; |
| 254 | } |
| 255 | } |
| 256 | } |
| 257 | |
| 258 | if (MaxU==999 || MaxD==999){ |
| 259 | return NOERROR; |
| 260 | } |
| 261 | |
| 262 | if (mU>mD){ |
| 263 | MaxD = MaxU; |
| 264 | } else { |
| 265 | MaxU = MaxD; |
| 266 | } |
| 267 | |
| 268 | //cout<<MaxD<<" / "<<MaxU<<endl; |
| 269 | |
| 270 | int mod = MaxD; |
| 271 | float Mt = 0; |
| 272 | float cnt = 0; |
| 273 | float EmaxU = 0.; |
| 274 | float EmaxD = 0.; |
| 275 | int idxMaxU[2] = {999,999}; |
| 276 | int idxMaxD[2] = {999,999}; |
| 277 | |
| 278 | // loop over all layers in the seed sector |
| 279 | // and its neighouring sectors and use all |
| 280 | // instances with hits on both ends to determine |
| 281 | // the mean time |
| 282 | for (int k=mod-1; k<mod+2; k++) { |
| 283 | int idx = k; |
| 284 | if (k<0){ |
| 285 | idx += 48*4; |
| 286 | } else if (k>=4*48){ |
| 287 | idx -= 48*4; |
| 288 | } |
| 289 | for (int n=0;n<3;n++){ |
| 290 | const DBCALHit* hitD = BcalMatrixD[idx][n]; |
| 291 | const DBCALHit* hitU = BcalMatrixU[idx][n]; |
| 292 | if (hitD && hitU){ |
| 293 | if (hitD->E > EmaxD){ |
| 294 | EmaxD = hitD->E; |
| 295 | idxMaxD[0] = idx; // index for seed |
| 296 | idxMaxD[1] = n; |
| 297 | } |
| 298 | if (hitU->E > EmaxU){ |
| 299 | EmaxU = hitU->E; |
| 300 | idxMaxU[0] = idx; // index for seed |
| 301 | idxMaxU[1] = n; |
| 302 | } |
| 303 | float mt = hitD->t + hitU->t; |
| 304 | Mt += mt; |
| 305 | cnt +=1.; |
| 306 | } |
| 307 | } |
| 308 | } |
| 309 | if (cnt>0){ |
| 310 | Mt /= cnt; // need improve to use Energy weight |
Value stored to 'Mt' is never read | |
| 311 | } |
| 312 | |
| 313 | //cout<<Mt<<" "<<cnt<<endl; |
| 314 | |
| 315 | if (EmaxU>EmaxD){ |
| 316 | idxMaxD[0] = idxMaxU[0]; |
| 317 | idxMaxD[1] = idxMaxU[1]; |
| 318 | } else { |
| 319 | idxMaxU[0] = idxMaxD[0]; |
| 320 | idxMaxU[1] = idxMaxD[1]; |
| 321 | } |
| 322 | |
| 323 | // Note: Mt is the flight time of the patricle creating this shower! |
| 324 | // This means Mt is the same for up and down stream!. |
| 325 | // |
| 326 | // Now collect all hits related to the seed hit for up and down stream |
| 327 | // separately. Only hits that appear in time (+/- 2ns) with the seed hit |
| 328 | // are added to the lists for layers 1,2 and 3. Layer 4 is always added. |
| 329 | // This may be changed if the timing information from the FADC is good enough. |
| 330 | // All hits in the list regardless if they contribut or not are then removed |
| 331 | // from the list of hits so they can not be used twice. |
| 332 | |
| 333 | CNT = (int)cnt; |
| 334 | |
| 335 | if (cnt>0) { |
| 336 | vector <const DBCALHit*> HitListU; |
| 337 | vector <const DBCALHit*> HitListD; |
| 338 | vector <float> MeanTime; |
| 339 | vector <float> DeltaTime; |
| 340 | vector <int> sector; |
| 341 | vector <int> layer; |
| 342 | |
| 343 | // float seedTime = BcalMatrixU[idxMaxU[0]][idxMaxU[1]]->t; |
| 344 | int Idx = MaxU; |
| 345 | // collect hits for up stream shower |
| 346 | // first move upwards in sectors from the seed |
| 347 | while(EUi[Idx]) { |
| 348 | //cout<<"EUi> "<<Idx<<" "<<EUi[Idx]<<endl; |
| 349 | for (int i=0;i<4;i++){ |
| 350 | const DBCALHit* hitU = BcalMatrixU[Idx][i]; |
| 351 | const DBCALHit* hitD = BcalMatrixD[Idx][i]; |
| 352 | if ((hitU)&&(hitD)){ |
| 353 | float mt = hitD->t + hitU->t; |
| 354 | //mt -= (390./VELOCITY); |
| 355 | mt -= (390./16.75); // subract detector internal path HARD CODED VALUES!!!!!! |
| 356 | mt /= 2.; |
| 357 | MeanTime.push_back(mt); |
| 358 | |
| 359 | // Note below D-U means positive dt is upstream |
| 360 | // from bcal center |
| 361 | float dt = hitD->t - hitU->t; |
| 362 | //dt = dt/2./VELOCITY; |
| 363 | dt = dt/2.*16.75; // convert to cm HARD CODED VALUE!!!!! |
| 364 | DeltaTime.push_back(dt); |
| 365 | sector.push_back(Idx); |
| 366 | layer.push_back(i); |
| 367 | } |
| 368 | |
| 369 | if (hitU){ |
| 370 | //if ((hitU->t - seedTime)<2.){ |
| 371 | HitListU.push_back(hitU); |
| 372 | //} |
| 373 | BcalMatrixU[Idx][i] = NULL__null; |
| 374 | } |
| 375 | } |
| 376 | EUi[Idx] = 0; |
| 377 | EU[Idx] = 0.; |
| 378 | Idx++; |
| 379 | if (Idx>=4*48){ |
| 380 | Idx -= 4*48; |
| 381 | } |
| 382 | } |
| 383 | // now move downwards from the seed center |
| 384 | Idx = MaxU-1; |
| 385 | if (Idx<0){ |
| 386 | Idx += 4*48; |
| 387 | } |
| 388 | while(EUi[Idx]) { |
| 389 | //cout<<"EUi< "<<Idx<<" "<<EUi[Idx]<<endl; |
| 390 | for (int i=0;i<4;i++){ |
| 391 | const DBCALHit* hitU = BcalMatrixU[Idx][i]; |
| 392 | const DBCALHit* hitD = BcalMatrixD[Idx][i]; |
| 393 | if ((hitU)&&(hitD)){ |
| 394 | float mt = hitD->t + hitU->t; |
| 395 | //mt -= (390./VELOCITY); |
| 396 | mt -= (390./16.5); // subract detector internal drift |
| 397 | mt /= 2.; |
| 398 | MeanTime.push_back(mt); |
| 399 | float dt = hitD->t - hitU->t; |
| 400 | //dt = dt/2./VELOCITY; |
| 401 | dt = dt/2.*16.5; // convert to cm |
| 402 | DeltaTime.push_back(dt); |
| 403 | sector.push_back(Idx); |
| 404 | layer.push_back(i); |
| 405 | } |
| 406 | |
| 407 | if (hitU){ |
| 408 | //if ((hitU->t - seedTime)<2.){ |
| 409 | HitListU.push_back(hitU); |
| 410 | //} |
| 411 | BcalMatrixU[Idx][i] = NULL__null; |
| 412 | } |
| 413 | } |
| 414 | |
| 415 | EUi[Idx] = 0; |
| 416 | EU[Idx] = 0.; |
| 417 | Idx--; |
| 418 | if (Idx<0){ |
| 419 | Idx += 4*48; |
| 420 | } |
| 421 | } |
| 422 | |
| 423 | // seedTime = BcalMatrixD[idxMaxD[0]][idxMaxD[1]]->t; |
| 424 | // collect hits for down stream shower |
| 425 | // first move upwards in sectors |
| 426 | Idx = MaxD; |
| 427 | while(EDi[Idx]) { |
| 428 | //cout<<"EDi> "<<Idx<<" "<<EUi[Idx]<<endl; |
| 429 | for (int i=0;i<4;i++){ |
| 430 | const DBCALHit* hit = BcalMatrixD[Idx][i]; |
| 431 | if (hit){ |
| 432 | //if ((hit->t - seedTime)<2.){ |
| 433 | HitListD.push_back(hit); |
| 434 | //} |
| 435 | BcalMatrixD[Idx][i] = NULL__null; |
| 436 | } |
| 437 | } |
| 438 | EDi[Idx] = 0; |
| 439 | ED[Idx] = 0.; |
| 440 | Idx++; |
| 441 | if (Idx>=4*48){ |
| 442 | Idx -= 4*48; |
| 443 | } |
| 444 | } |
| 445 | // now move downwrads from center |
| 446 | Idx = MaxD-1; |
| 447 | if (Idx<0){ |
| 448 | Idx += 4*48; |
| 449 | } |
| 450 | while(EDi[Idx]) { |
| 451 | //cout<<"EDi< "<<Idx<<" "<<EUi[Idx]<<endl; |
| 452 | for (int i=0;i<4;i++){ |
| 453 | const DBCALHit* hit = BcalMatrixD[Idx][i]; |
| 454 | if (hit){ |
| 455 | //if ((hit->t - seedTime)<2.){ |
| 456 | HitListD.push_back(hit); |
| 457 | //} |
| 458 | BcalMatrixD[Idx][i] = NULL__null; |
| 459 | } |
| 460 | } |
| 461 | EDi[Idx] = 0; |
| 462 | ED[Idx] = 0.; |
| 463 | Idx--; |
| 464 | if (Idx<0){ |
| 465 | Idx += 4*48; |
| 466 | } |
| 467 | } |
| 468 | |
| 469 | if ( (HitListU.size()>0) && (HitListD.size()>0) ) { |
| 470 | //cout<<"Make New Clump:"<<cnt<<endl; |
| 471 | DBCALClump* myClump = new DBCALClump(HitListU, HitListD); |
| 472 | myClump->MeanTime = MeanTime; |
| 473 | myClump->DeltaTime = DeltaTime; |
| 474 | myClump->Sector = sector; |
| 475 | myClump->Layer = layer; |
| 476 | myClump->AnalyzeClump(); |
| 477 | |
| 478 | int oK = 1; |
| 479 | if ((HitListU.size()==1) || (HitListD.size()==1)){ |
| 480 | if (myClump->ClumpE[0]<50.) { // HARD CODED |
| 481 | oK = 0; |
| 482 | } |
| 483 | } |
| 484 | if (oK){ |
| 485 | _data.push_back(myClump); |
| 486 | } else { |
| 487 | delete myClump; |
| 488 | } |
| 489 | } else { |
| 490 | cout<<"Error no hits in this Clump!!!! Event: "<<eventnumber<<" cnt= "<<cnt <<endl; |
| 491 | } |
| 492 | } |
| 493 | |
| 494 | } |
| 495 | |
| 496 | return NOERROR; |
| 497 | } |
| 498 |