![(please configure the [header_logo] section in trac.ini)](/PPS/chrome/site/your_project_logo.png){kind=logo}
Changeset 222 in Idarraga
- Timestamp:
- Jul 26, 2011, 3:34:46 PM (13 years ago)
- Location:
- mafalda
- Files:
-
- 11 edited
-
AnalysisCore/MediPixAnalysisCore.cpp (modified) (3 diffs)
-
AnalysisCore/MediPixAnalysisCore.h (modified) (4 diffs)
-
AnalysisCore/MediPixAnalysisManager.h (modified) (4 diffs)
-
BlobsFinder/BlobsFinder.cpp (modified) (3 diffs)
-
BlobsFinder/BlobsFinder.h (modified) (6 diffs)
-
FEI3/FEI3Mips/FEI3Mips.cpp (modified) (27 diffs)
-
FEI3/FEI3Mips/FEI3Mips.h (modified) (3 diffs)
-
MAFTools/MAFTools.cpp (modified) (2 diffs)
-
MAFTools/MAFTools.h (modified) (3 diffs)
-
MPXAlgo/MediPixAlgo.h (modified) (1 diff)
-
ShallowAngleFEIX/ShallowAngleFEIX.cpp (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
mafalda/AnalysisCore/MediPixAnalysisCore.cpp
r179 r222 45 45 m_outputFilename = "MPXResults.root"; // output filename by default 46 46 47 // fast matrix 48 m_frameXC_Fast = 0x0; 49 50 // TOT cut 51 m_totMaxCut = 11800; // timepix cut, settable through toplayer 47 52 } 48 53 … … 134 139 nb = fChain->GetEntry(jentry); nbytes += nb; 135 140 136 //jFrameId = jentry; 137 138 // Calling algos to execution 141 // At this point the Baskets have been opened. 142 // Reading the pixel matrix (std::map) might be slow at this point. 143 // Converting to a C-type array in order to fasten the reading 144 RemakePixelMap(); 145 146 // Calling algos for execution 139 147 for (algScheduleItr = m_algosSchedule.begin() ; algScheduleItr != m_algosSchedule.end() ; algScheduleItr++ ) { 140 148 … … 339 347 } 340 348 349 void MediPixAnalysisCore::RemakePixelMap(){ 350 351 if ( m_frameXC_Fast == 0x0 ) { // first time 352 353 Log << MSG::ALWAYS << "Creating fast matrix map" << endreq; 354 m_frameXC_Fast = new int[fWidth*fHeight]; 355 356 } else if ( fWidth != fWidth_previous || fHeight != fHeight_previous ) { // in case matrix size changes during run 357 358 delete m_frameXC_Fast; 359 m_frameXC_Fast = new int[fWidth*fHeight]; 360 361 } 362 363 // clean up always 364 for( int i = 0 ; i < fWidth*fHeight ; i++ ) m_frameXC_Fast[i] = 0; 365 366 // copy information from the map 367 std::map<int,int>::iterator itr = m_frameXC.begin(); 368 int arraylength = fWidth*fHeight; 369 for ( ; itr != m_frameXC.end() ; itr++) { 370 371 if((*itr).first > arraylength) { 372 std::cout << "[{AnalysisCore} WARNING] Attempt to fill up entries outside of the pad --> inconsistent data !" << std::endl; 373 std::cout << " Trying to recover by avoiding this entry. Position : " << (*itr).first << std::endl; 374 std::cout << " in a " << fWidth << "x" << fHeight << " matrix." << std::endl; 375 continue; 376 } 377 378 if((*itr).second < m_totMaxCut) { 379 m_frameXC_Fast[(*itr).first] = (*itr).second; 380 } 381 382 } 383 384 fWidth_previous = fWidth; 385 fHeight_previous = fHeight; 386 } 341 387 342 388 #endif -
mafalda/AnalysisCore/MediPixAnalysisCore.h
r179 r222 72 72 std::map<int,double> m_frameXC_E; 73 73 74 // Fast versions 75 int * m_frameXC_Fast; 76 74 77 std::map<int, int> m_lvl1; 75 78 //Int_t m_frameMatrix[256][256]; … … 81 84 82 85 Int_t fFormat; 86 // size of the matrix 83 87 Int_t fWidth; 84 88 Int_t fHeight; 89 // check previous size 90 Int_t fWidth_previous; 91 Int_t fHeight_previous; 92 85 93 Int_t fAcq_mode; 86 94 Double_t fAcq_time; … … 194 202 inline void ListenToTheManager(AnalysisManager * mngr){theManager = mngr;}; 195 203 204 inline void SetTOTMaxCut(int c){m_totMaxCut = c;}; 205 196 206 private: 207 208 // Copy the frameXC map to a C-array. 209 // This fasten things up 210 void RemakePixelMap(); 211 197 212 198 213 MediPixAlgoTimer * algoTimers; … … 211 226 // default branches 212 227 std::map<string, Double_t> m_timePerFrameMap; 228 229 // TOT cuts 230 int m_totMaxCut; 213 231 214 232 }; -
mafalda/AnalysisCore/MediPixAnalysisManager.h
r175 r222 59 59 MediPixStoreGate * GetStoreGate(){return storeGate;}; 60 60 61 /* Retreive data safely */ 61 /* 62 //Retreive data safely 62 63 inline Int_t GetMatrixElement(Int_t col, Int_t row) 63 64 { … … 69 70 return analysisCore->m_frameXC[xplace]; 70 71 } 71 //return analysisCore->m_frameMatrix[col][row]; 72 73 return 0; 74 }; 75 */ 76 inline Int_t GetMatrixElement(Int_t col, Int_t row) 77 { 78 int xplace = row*analysisCore->fWidth + col; 79 return analysisCore->m_frameXC_Fast[xplace]; 72 80 }; 73 81 … … 142 150 143 151 inline std::vector<Int_t> GetfCounters() 144 {return analysisCore->fCounters;};152 {return analysisCore->fCounters;}; 145 153 146 154 inline std::vector<Int_t> GetfDACs() 147 {return analysisCore->fDACs;};155 {return analysisCore->fDACs;}; 148 156 149 157 inline Int_t GetTHL(){ … … 206 214 207 215 inline TApplication * GetApplication() 208 {return analysisCore->g_theApp;};216 {return analysisCore->g_theApp;}; 209 217 210 218 private: -
mafalda/BlobsFinder/BlobsFinder.cpp
r163 r222 12 12 #include "MAFTools/MAFTools.h" 13 13 14 14 15 using namespace MSG; 15 16 … … 28 29 void AllBlobsContainer::SetBlobType(Int_t id, blobtype type){ 29 30 allBlobs[id].SetBlobType(type); 31 30 32 } 31 33 … … 38 40 } 39 41 40 void blob::SetBlobType(blobtype type) 41 { 42 void blob::SetBlobType(blobtype type) { 42 43 btype = type; 43 44 } 45 46 /* 47 template<class T> 48 void blob<T>::SetBlobType(T type) { 49 //void blob<T>::SetBlobType(blobtype type) { 50 btype = type; 51 } 52 */ 44 53 45 54 /** -
mafalda/BlobsFinder/BlobsFinder.h
r95 r222 34 34 #define __BLOB_GOOD_TO_STORE 201 35 35 36 #define MAX_BLOBS_PER_FRAME 2000 36 #define MAX_BLOBS_PER_FRAME 20000 37 37 // numbering of quadrants 38 38 #define Q_FIRST 1 … … 129 129 * One blob contents 130 130 */ 131 //typedef struct { 131 132 132 class blob { 133 133 … … 135 135 136 136 TString GetTypeAsString(); 137 //void SetBlobType(blobtype); 137 138 void SetBlobType(blobtype); 138 139 blobProperties GetBlobProperties(){return bP;}; 139 140 blobtype GetBlobType(){return btype;}; 140 141 set< pair<Int_t, Int_t> > GetContentSet(){return blobContentSet;}; 142 list< pair < pair<Int_t, Int_t>, Int_t > > GetBlobContent(){return blobContent;}; 141 143 142 144 // this will go private ... TODO !!! … … 153 155 154 156 blobProperties bP; 157 //blobtype btype; 155 158 blobtype btype; 156 159 157 160 }; 158 161 162 //template class blob<blobtype>; 159 163 160 164 /** … … 163 167 * 164 168 */ 169 165 170 /* *************************************** */ 166 171 class AllBlobsContainer : public CandidateContainer { 167 172 168 173 public: 174 169 175 AllBlobsContainer(MediPixAlgo *); 170 176 virtual ~AllBlobsContainer(){}; … … 173 179 blobtype GetBlobType(Int_t); 174 180 TString GetTypeAsString(Int_t); 181 void PushBack(blob b){ allBlobs.push_back(b); }; 175 182 176 183 public: -
mafalda/FEI3/FEI3Mips/FEI3Mips.cpp
r179 r222 12 12 #include "MPXAlgo/Highlighter.h" 13 13 14 #include "TH1F.h"15 16 14 using namespace MSG; 17 15 … … 20 18 FEI3Mips::FEI3Mips() 21 19 : m_minNPixels ( 5 ), 20 m_minMipLength ( 3.0 ), // [mm] 22 21 m_nDivisions ( 4 ), 23 22 m_pixSizeX ( 0.40 ), // FEI3 [mm] … … 25 24 m_FEI3_PIXELSIZE_YX_FRACTION ( m_pixSizeY/m_pixSizeX ), 26 25 m_totalMIPsLength ( 0.0 ), 27 m_dRayBalSearch ( 1.4 ) 26 m_dRayBalSearch ( 1.4 ), 27 m_distanceCutHot ( 3 ), 28 m_distancePerpHot ( 3 ), 29 m_guardDistanceX ( 1 ), 30 m_guardDistanceY ( 5 ) 28 31 { 29 32 … … 41 44 getMyTree()->Branch("lvl1Weights", &m_lvl1Weights); 42 45 getMyTree()->Branch("chargePerSection", &m_chargePerSection); 46 getMyTree()->Branch("deltaRayTOT", &m_deltaRayTOT); 47 getMyTree()->Branch("mipLength", &m_mipLength); 48 getMyTree()->Branch("meanChargePerSection", &m_meanChargePerSection); 49 getMyTree()->Branch("sectionCenter_x", &m_sectionCenter_x); 50 getMyTree()->Branch("sectionCenter_y", &m_sectionCenter_y); 51 getMyTree()->Branch("nDeltaRayExplicit", &m_nDeltaRayPerFrame); 52 getMyTree()->Branch("nDeltaRaySoft", &m_nDeltaRaySoftPerFrame); 43 53 44 54 getMyTree()->Branch("clusterSize", &m_clusterSize, "clusterSize/I"); … … 48 58 // A configuration value that can be tuned from the Viewer 49 59 RegisterConfigurationValue(&m_minNPixels, "minNPixels"); 60 RegisterConfigurationValue(&m_minMipLength, "minMipLength"); 50 61 RegisterConfigurationValue(&m_nDivisions, "nDivisions"); 51 62 RegisterConfigurationValue(&m_pixSizeX, "pixelSizeX"); 52 63 RegisterConfigurationValue(&m_pixSizeY, "pixelSizeY"); 53 64 RegisterConfigurationValue(&m_dRayBalSearch, "dRayBalSearch"); 54 65 RegisterConfigurationValue(&m_distanceCutHot, "distanceCutHot"); 66 RegisterConfigurationValue(&m_distancePerpHot, "distancePerpHot"); 67 RegisterConfigurationValue(&m_guardDistanceX, "guardDistanceX"); 68 RegisterConfigurationValue(&m_guardDistanceY, "guardDistanceY"); 55 69 56 70 } 57 71 58 72 void FEI3Mips::Execute(){ 73 74 // A check on parameters 75 // The number of divisions can not be smaller than the minimum amount of 76 // pixels required in a mip 77 if(m_nDivisions > m_minNPixels) { 78 Log << MSG::ERROR << "The number of divisions can not be bigger than the minimum amount of" << endreq; 79 Log << MSG::ERROR << "pixels required in a mip to process. Fix the parameters. Giving up." << endreq; 80 exit(1); 81 } 59 82 60 83 // Ask the store gate if the previous algorithm (BlobsFinder --> reponsible for clustering) … … 68 91 // AllBlobsContainer is a box full of Blobs(Clusters). Now we can iterate over all clusters inside. 69 92 vector<blob> blobsVector = m_aB->GetBlobsVector(); 93 94 // If nothing to do, return 95 if(blobsVector.empty()) return; 96 70 97 Log << MSG::INFO << "Number of blobs from clustering = " << (Int_t) blobsVector.size() << endreq; 71 98 vector<blob>::iterator blobsItr = blobsVector.begin(); //allBlobs.begin(); … … 88 115 m_clusterSize = (*blobsItr).GetBlobProperties().nPixels; 89 116 117 // Here we obtain two pairs with the pixels at the extreme ends of the Mip 118 //vector< pair<Int_t, Int_t> > limitPixels = FindLimitPixels(*blobsItr); 119 //i- 120 vector< pair<Float_t, Float_t> > limitPixels = FindLimitPixels(*blobsItr); 121 122 // Make sure the mip starts and ends inside the pad at a certain minimum 123 // distance from the borders 124 if ( MAFTools::BlobAtADistanceFromEdges(*blobsItr, limitPixels, GetWidth(), GetHeight(), 125 m_guardDistanceX, m_guardDistanceY) ) { 126 continue; 127 } 128 90 129 // create all slots to put the weighted charge and others observables 91 130 for(int ni = 0 ; ni < m_nDivisions ; ni++){ … … 93 132 m_lvl1Weights.push_back(0.0); 94 133 m_chargePerSection.push_back(vector<Int_t>()); // push back empty vectors 95 } 96 97 // Here we obtain two pairs with the pixels at the extreme ends of the Mip 98 vector< pair<Int_t, Int_t> > limitPixels = FindLimitPixels(*blobsItr); 134 m_XYPerSection.push_back(vector<pair <Int_t, Int_t> >()); 135 m_meanChargePerSection.push_back(0.0); 136 m_sectionCenter_x.push_back(0.0); 137 m_sectionCenter_y.push_back(0.0); 138 } 99 139 100 140 // If the whole mip is in the same column, don't process it. … … 102 142 // same column in FE-I3. FE-I4 too ? 103 143 if (limitPixels[__TAIL_INDX].first - limitPixels[__HEAD_INDX].first == 0) { 104 return; 144 ClearThisAlgo(); 145 continue; 105 146 } 106 147 … … 112 153 limitPixels[__TAIL_INDX].first * m_pixSizeX, 113 154 limitPixels[__TAIL_INDX].second * m_pixSizeY); 155 m_mipLength.push_back( totalLength ); // mm 114 156 m_totalMIPsLength += totalLength; 115 Log << MSG::DEBUG << "Integrated mips length = " << m_totalMIPsLength << " [mm]" << endreq; 157 Log << MSG::INFO << "Integrated mips length = " << m_totalMIPsLength << " [mm]" << endreq; 158 159 // cut by mip length 160 if ( totalLength < m_minMipLength ) { 161 ClearThisAlgo(); 162 continue; 163 } 116 164 117 165 // Fech the linear regression information for this mip … … 153 201 if(x2i < x1i) startX = x2i; 154 202 155 // build the rest of the point within the boundaries203 // Build the rest of the point within the boundaries 156 204 for (int i = 1 ; i < m_nDivisions ; i++) { 157 205 … … 159 207 Float_t newPY = newPX*slope + cut; 160 208 161 limitPixels.push_back( make_pair( TMath::FloorNint(newPX), TMath::FloorNint(newPY) ) ); 209 Log << MSG::INFO << "division point : " << newPX << " , " << newPY << endreq; 210 211 //limitPixels.push_back( make_pair( TMath::FloorNint(newPX), TMath::FloorNint(newPY) ) ); 212 //i- 213 limitPixels.push_back( make_pair( newPX, newPY ) ); 162 214 } 163 215 … … 166 218 vector< pair<Float_t, Float_t > > linePars; 167 219 220 // Flag which indicates if we are dealing with a track with slope = 0. Special case. 221 bool parallelLine = false; 168 222 for (int i = 0 ; i < m_nDivisions+1 ; i++) { // there are m_nDivisions+1 lines ! 169 linePars.push_back( GetParallelLineAt(slope, cut, limitPixels[i]) ); 223 224 // WARNING ! If slope == 0. a flag indicating it is raised 225 linePars.push_back( GetParallelLineAt(slope, cut, limitPixels[i], ¶llelLine) ); 226 170 227 } 171 228 … … 195 252 for( separatorItr = linePars.end()-1 ; separatorItr != linePars.begin()-1 ; separatorItr-- ) { 196 253 197 // calculate the perpendicular distance to each line 198 Double_t perpDistance = MAFTools::CalcPerpDistanceToLine( (*separatorItr).first, 199 (*separatorItr).second, *contItr ); 254 // deal with the slope=0 case 255 Double_t perpDistance = 0.; 256 if(parallelLine) { 257 //Log << MSG::DEBUG << x << " --> " << (*separatorItr).first << endreq; 258 perpDistance = TMath::Abs( x - (*separatorItr).first); 259 } else { 260 // calculate the perpendicular distance to each line 261 perpDistance = MAFTools::CalcPerpDistanceToLine( (*separatorItr).first, 262 (*separatorItr).second, *contItr ); 263 } 200 264 201 265 Log << MSG::DEBUG << "pixel : " << x << " , " << y << " | perpDistance = " << perpDistance << " "; … … 205 269 locationBitWord ^= 0x1; 206 270 Log << MSG::DEBUG << " 1 "; 207 } else{271 } else { 208 272 //locationBitWord ^= 0x0; 209 273 Log << MSG::DEBUG << " 0 "; … … 212 276 if(separatorItr != linePars.begin()) locationBitWord = locationBitWord << 1; 213 277 Log << MSG::DEBUG << endreq; 278 214 279 } 215 280 … … 232 297 // store charge for landau distribution per section 233 298 m_chargePerSection[__HEAD_INDX].push_back( GetMatrixElement(*contItr) ); 299 m_XYPerSection[__HEAD_INDX].push_back( make_pair(x, y) ); 300 301 // mean charge per section 302 m_meanChargePerSection[__HEAD_INDX] += GetMatrixElement(*contItr); 303 304 // mean position 305 m_sectionCenter_x[__HEAD_INDX] += x; 306 m_sectionCenter_y[__HEAD_INDX] += y; 307 234 308 235 309 }else if( firstThreeBits == __RIGHT_CORNER_MASK_1 … … 242 316 // store charge for landau distribution per section 243 317 m_chargePerSection[__TAIL_INDX].push_back( GetMatrixElement(*contItr) ); 318 m_XYPerSection[__TAIL_INDX].push_back( make_pair(x, y) ); 319 320 // mean charge per section 321 m_meanChargePerSection[__TAIL_INDX] += GetMatrixElement(*contItr); 322 323 // mean position 324 m_sectionCenter_x[__TAIL_INDX] += x; 325 m_sectionCenter_y[__TAIL_INDX] += y; 244 326 245 327 }else{ … … 259 341 nB--; 260 342 } 343 261 344 Log << MSG::DEBUG << "... column " << nB << " ... " << GetFrameId() << endreq; 262 345 m_chargeWeights[nB] += GetMatrixElement(*contItr); 263 346 m_lvl1Weights[nB] += GetLVL1(*contItr); 347 264 348 // Store charge for landau distribution per section 265 349 m_chargePerSection[nB].push_back( GetMatrixElement(*contItr) ); 350 m_XYPerSection[nB].push_back( make_pair(x, y) ); 351 352 // charge per section 353 m_meanChargePerSection[nB] += GetMatrixElement(*contItr); 354 355 // mean position 356 m_sectionCenter_x[nB] += x; 357 m_sectionCenter_y[nB] += y; 266 358 267 359 } … … 269 361 } 270 362 363 // FIXME !!! 364 // finish up the mean charge per section and position 365 for (int ix = 0 ; ix < (int)m_chargePerSection.size() ; ix++) { 366 int nEntriesSection = m_chargePerSection[ix].size(); 367 m_meanChargePerSection[ix] /= nEntriesSection; 368 m_sectionCenter_x[ix] /= nEntriesSection; 369 m_sectionCenter_y[ix] /= nEntriesSection; 370 } 271 371 272 372 // normalize by the total charge 273 373 vector<Float_t>::iterator itr = m_chargeWeights.begin(); 274 374 int cntr = 0 ; 275 for( ; itr != m_chargeWeights.end() ; itr++) {375 for( ; itr != m_chargeWeights.end() ; itr++) { 276 376 (*itr) /= (*blobsItr).GetBlobProperties().totalCharge; 277 377 (*itr) *= 100.; // percentage … … 288 388 m_nDeltaRaySoft += nDeltaSoft; 289 389 390 m_nDeltaRayPerFrame += nDeltaR; // this one is cleared on a frame-per-frame basis 391 m_nDeltaRaySoftPerFrame += nDeltaSoft; // this one is cleared on a frame-per-frame basis 392 393 // Look at the hard delta rays trying to extract an energy profile 394 if (nDeltaR > 0) { 395 m_deltaRayTOT = ExtractDeltaRayEnergy(*blobsItr); 396 } 397 290 398 // WARNING, filling once per mip !!! 291 399 m_frameId = GetFrameId(); … … 295 403 m_lvl1Weights.clear(); 296 404 m_chargePerSection.clear(); 405 m_XYPerSection.clear(); 406 m_mipLength.clear(); 407 m_meanChargePerSection.clear(); 408 m_sectionCenter_x.clear(); 409 m_sectionCenter_y.clear(); 410 297 411 m_alphaAngle = 0.; 298 412 m_clusterSize = 0; 413 m_nDeltaRayPerFrame = 0; 414 m_nDeltaRaySoftPerFrame = 0; 299 415 300 416 } else { … … 311 427 m_lvl1Weights.clear(); 312 428 m_chargePerSection.clear(); 313 314 } 315 429 m_XYPerSection.clear(); 430 m_deltaRayTOT.clear(); 431 m_mipLength.clear(); 432 433 } 434 435 void FEI3Mips::ClearThisAlgo(){ 436 437 m_chargeWeights.clear(); 438 m_lvl1Weights.clear(); 439 m_chargePerSection.clear(); 440 m_XYPerSection.clear(); 441 m_mipLength.clear(); 442 m_meanChargePerSection.clear(); 443 m_sectionCenter_x.clear(); 444 m_sectionCenter_y.clear(); 445 446 } 316 447 317 448 void FEI3Mips::Finalize() { … … 325 456 Log << MSG::INFO << "Total Number delta rays per centimeter = " << (m_nDeltaRay+m_nDeltaRaySoft)/(m_totalMIPsLength*0.1) << endreq; 326 457 458 } 459 460 vector<Float_t> FEI3Mips::ExtractDeltaRayEnergy (blob bl) { 461 462 // At this point the information in m_chargeWeights is already 463 // expressed in percentage per segment 464 vector<Float_t>::iterator itr = m_chargeWeights.begin(); 465 Int_t indxMaxSeg = 0; 466 Float_t maxSeg = 0.; 467 468 Int_t cntr = 0; 469 // Find the hotest segment 470 for ( ; itr != m_chargeWeights.end() ; itr++ ) { 471 472 if ( (*itr) > maxSeg ) { 473 maxSeg = (*itr); 474 indxMaxSeg = cntr; 475 } 476 cntr++; 477 } 478 479 // indxMaxSeg is the index of the hotest segment 480 Log << MSG::INFO << "[" << indxMaxSeg << "] max = " << maxSeg << endreq; 481 482 // Now in this segment find the hotest point ( look at indxMaxSeg ! ) 483 vector<Int_t>::iterator cItr = m_chargePerSection[indxMaxSeg].begin(); 484 vector< pair<Int_t, Int_t> >::iterator xyItr = m_XYPerSection[indxMaxSeg].begin(); 485 Int_t highC = 0; 486 Int_t highIndx = 0; 487 cntr = 0; 488 for ( ; cItr != m_chargePerSection[indxMaxSeg].end() ; ) { 489 490 if( *cItr > highC ) { 491 highC = *cItr; 492 highIndx = cntr; 493 } 494 495 cntr++; 496 cItr++; 497 xyItr++; 498 } 499 500 Log << MSG::INFO << "spot <" << m_XYPerSection[indxMaxSeg][highIndx].first << ", " 501 << m_XYPerSection[indxMaxSeg][highIndx].second << "> : " 502 << m_chargePerSection[indxMaxSeg][highIndx] << endreq; 503 504 // Once having the hotest pixel I can search for all the charge belonging to the 505 // delta ray which has been emitted at m_XYPerSection[indxMaxSeg][highIndx]. 506 // 507 // The approach is the following: 508 // 1) Pixels very close to the hotest pixel probably have a bit of a higher TOT that 509 // the mean energy in the rest of the mip track. That extra energy should be added 510 // to the delta ray spectrum. 511 // 2) If the delta ray stayed in the plane of the sensor, there will probably be a few 512 // pixels ON close to the mip track with the shape of a low energy electron (curly). 513 // This energy ought to be added too. 514 515 vector<Float_t> deltaRaysTOT; 516 517 // The mean TOT in all the track is 518 Float_t meanTOT = Float_t( bl.GetBlobProperties().totalCharge ) / 519 Float_t( bl.GetBlobProperties().nPixels ); 520 521 Log << MSG::INFO << "mean TOT in mip = " << meanTOT << endreq; 522 523 // calculate perp distance in each point and also the distance to the hot point 524 Float_t slope = bl.GetBlobProperties().fitSlope; 525 Float_t cut = bl.GetBlobProperties().fitCut; 526 // Prepare perpendicular line passing by hot pixel 527 Float_t slopePrim = -1/slope; 528 Float_t cutPrim = m_XYPerSection[indxMaxSeg][highIndx].second; 529 cutPrim -= slopePrim * m_XYPerSection[indxMaxSeg][highIndx].first; 530 531 list< pair < pair<Int_t, Int_t>, Int_t > > contents = bl.GetBlobContent(); 532 list< pair < pair<Int_t, Int_t>, Int_t > >::iterator bcItr = contents.begin(); 533 534 Double_t distP = 0., distHot = 0., distHotPerp = 0.; 535 Float_t deltaRayTotalTOT = 0.; 536 537 // draw the perp line passing by hot pixel 538 Float_t x1p = m_XYPerSection[indxMaxSeg][highIndx].first; 539 Float_t y1p = m_XYPerSection[indxMaxSeg][highIndx].second; 540 Float_t x2p = x1p*1.5; 541 Float_t y2p = slopePrim * x2p + cutPrim; 542 MAFTools::DrawLine(this, x1p, y1p, x2p, y2p, 2, 1, kBlack); 543 544 for ( ; bcItr != contents.end() ; bcItr++ ) { 545 546 // perpendicular distance to Linear Regression 547 distP = MAFTools::CalcPerpDistanceToLine(slope, cut, (*bcItr).first); 548 // distance to hot pixel (vertex) 549 distHot = MAFTools::CalcDistance(m_XYPerSection[indxMaxSeg][highIndx], (*bcItr).first); 550 // perpendicular distance to perp to Linear Regression passing by hot pixel 551 distHotPerp = MAFTools::CalcPerpDistanceToLine(slopePrim, cutPrim, (*bcItr).first); 552 553 if( (distP < 1.0 && distHot < m_distanceCutHot) // close to vertex 554 || distHotPerp < m_distancePerpHot) // parallel path 555 { 556 557 Highlighter * fp = new Highlighter((*bcItr).first.first, 558 (*bcItr).first.second, 559 "arrow", this); 560 fp->SetLineColor(kRed); 561 fp->SetLineWidth(1); 562 PullToStoreGateAccess(fp, MPXDefs::DO_NOT_SERIALIZE_ME); 563 564 Log << MSG::DEBUG << "perp dist = " << distP << ", distance to Hot = " << distHot << endreq; 565 // Finally store the excess energy if it is over the mean 566 Float_t TOTDiff = GetMatrixElement( (*bcItr).first ) - meanTOT; 567 if( TOTDiff > 0 ) { 568 // FIXME !, store the whole energy ! 569 deltaRayTotalTOT += TOTDiff; 570 } 571 } 572 } 573 574 deltaRaysTOT.push_back(deltaRayTotalTOT); 575 576 return deltaRaysTOT; 327 577 } 328 578 … … 379 629 } 380 630 381 vector<pair< Int_t, Int_t> > FEI3Mips::FindLimitPixels(blob theBlob) {631 vector<pair<Float_t, Float_t> > FEI3Mips::FindLimitPixels(blob theBlob) { 382 632 383 633 // result vector 384 vector<pair< Int_t, Int_t> > limits;634 vector<pair<Float_t, Float_t> > limits; 385 635 limits.push_back( make_pair(0, 0) ); 386 636 limits.push_back( make_pair(0, 0) ); … … 400 650 401 651 if( dist > max ) { 652 402 653 max = dist; 654 403 655 limits[0] = *contItrA; 404 656 limits[1] = *contItrB; 657 658 // The limit in the extreme right needs an adjustement of 1 pixel 659 limits[1].first += 1; // 1 pixel 660 405 661 } 406 662 … … 428 684 } 429 685 430 pair<Float_t, Float_t> FEI3Mips::GetParallelLineAt(Float_t slope, Float_t, pair<Int_t, Int_t> point) { 686 // FIXME ... shouldn't this be called "Perpendicular" ? just check the logics. 687 //pair<Float_t, Float_t> FEI3Mips::GetParallelLineAt(Float_t slope, Float_t cut, pair<Int_t, Int_t> point, bool * pL) { 688 //i- 689 pair<Float_t, Float_t> FEI3Mips::GetParallelLineAt(Float_t slope, Float_t cut, pair<Float_t, Float_t> point, bool * pL) { 690 691 if ( slope == 0. ) { 692 // In this trivial case, I basically return the pair 'point'. 693 // I will only need the x coordinate to calculate the distance later on 694 *pL = true; 695 return point; 696 697 } 431 698 432 699 Float_t cutPrim = 0.; -
mafalda/FEI3/FEI3Mips/FEI3Mips.h
r175 r222 35 35 void Finalize(); 36 36 37 vector<pair<Int_t, Int_t> > FindLimitPixels(blob); 38 pair<Float_t, Float_t> GetParallelLineAt(Float_t, Float_t, pair<Int_t, Int_t>); 37 vector<pair<Float_t, Float_t> > FindLimitPixels(blob); 38 //-i 39 //pair<Float_t, Float_t> GetParallelLineAt(Float_t, Float_t, pair<Int_t, Int_t>, bool *); 40 pair<Float_t, Float_t> GetParallelLineAt(Float_t, Float_t, pair<Float_t, Float_t>, bool *); 39 41 Int_t NumberOfBitsOnBeforeDoubleZero(UInt_t); 40 42 Int_t GuessNumberOfDeltaRaysFromUnbalance(Int_t &, Int_t &); 43 vector<Float_t> ExtractDeltaRayEnergy(blob); 44 void ClearThisAlgo(); 41 45 42 46 private: … … 44 48 AllBlobsContainer * m_aB; 45 49 Int_t m_minNPixels; 50 Float_t m_minMipLength; 46 51 Int_t m_nDivisions; 47 52 Float_t m_pixSizeX; 48 53 Float_t m_pixSizeY; 49 54 Float_t m_dRayBalSearch; 55 Int_t m_guardDistanceX; 56 Int_t m_guardDistanceY; 50 57 const Float_t m_FEI3_PIXELSIZE_YX_FRACTION; 51 58 … … 54 61 vector<Float_t> m_lvl1Weights; 55 62 63 // charge per section 56 64 vector< vector<Int_t> > m_chargePerSection; 65 vector< Float_t > m_meanChargePerSection; 66 vector< Float_t > m_sectionCenter_x; 67 vector< Float_t > m_sectionCenter_y; 57 68 69 // store the coordinates too 70 vector< vector< pair<Int_t, Int_t> > > m_XYPerSection; 71 72 vector<Float_t> m_deltaRayTOT; 73 vector<Float_t> m_mipLength; 58 74 Double_t m_totalMIPsLength; 59 Int_t m_nDeltaRay; 60 Int_t m_nDeltaRaySoft; 75 76 Int_t m_nDeltaRay; // counter over the whole run 77 Int_t m_nDeltaRaySoft; // counter over the whole run 78 Int_t m_nDeltaRayPerFrame; // to store frame per frame 79 Int_t m_nDeltaRaySoftPerFrame; // to store frame per frame 80 81 Int_t m_distanceCutHot; 82 Int_t m_distancePerpHot; 61 83 62 84 Int_t m_clusterSize; -
mafalda/MAFTools/MAFTools.cpp
r124 r222 12 12 #include "MPXAlgo/Highlighter.h" 13 13 #include "MPXAlgo/MediPixAlgo.h" 14 15 #include "BlobsFinder/BlobsFinder.h" 14 16 15 17 #include <set> … … 157 159 } 158 160 161 vector<pair<Float_t, Float_t> > FindLimitPixels(blob theBlob) { 162 163 // result vector 164 vector<pair<Float_t, Float_t> > limits; 165 limits.push_back( make_pair(0, 0) ); 166 limits.push_back( make_pair(0, 0) ); 167 168 set< pair<Int_t, Int_t> > contents = theBlob.GetContentSet(); 169 170 set< pair<Int_t, Int_t> >::iterator contItrA = contents.begin(); 171 set< pair<Int_t, Int_t> >::iterator contItrB = contents.begin(); 172 173 double max = 0., dist = 0.; 174 175 for( ; contItrA != contents.end() ; contItrA++) { 176 177 for(contItrB = contItrA ; contItrB != contents.end() ; contItrB++) { 178 179 dist = MAFTools::CalcDistance( *contItrA, *contItrB ); 180 181 if( dist > max ) { 182 max = dist; 183 limits[0] = *contItrA; 184 limits[1] = *contItrB; 185 } 186 187 } 188 189 } 190 191 return limits; 192 } 193 194 195 Bool_t BlobAtADistanceFromEdges(blob, vector<pair<Float_t, Float_t> > l, Int_t width, Int_t height, 196 Int_t cutx, Int_t cuty) { 197 198 Int_t dist = 0.; // in units of pixels 199 200 // lower border 201 for (int i = 0 ; i < (int)l.size() ; i++) { 202 dist = l[i].second; 203 if ( dist <= cuty ) return true; 204 } 205 // left border 206 for (int i = 0 ; i < (int)l.size() ; i++) { 207 dist = l[i].first; 208 if ( dist <= cutx ) return true; 209 } 210 // top border 211 for (int i = 0 ; i < (int)l.size() ; i++) { 212 dist = height - l[i].second; 213 if ( dist <= cuty ) return true; 214 } 215 // right border 216 for (int i = 0 ; i < (int)l.size() ; i++) { 217 dist = width - l[i].first; 218 if ( dist <= cutx ) return true; 219 } 220 221 return false; 222 } 223 159 224 } 160 225 -
mafalda/MAFTools/MAFTools.h
r124 r222 1 1 /** 2 * 3 * Author: John Idarraga <idarraga@cern.ch> , 2009 4 * Medipix Group, Universite de Montreal 5 * 2 * Author: John Idarraga <idarraga@cern.ch> 6 3 */ 7 4 … … 15 12 #include <TROOT.h> 16 13 14 class blob; 17 15 class MediPixAlgo; 18 16 … … 45 43 std::string DumpWordInBinary(UInt_t word); 46 44 45 // Blobs handling 46 vector<pair<Float_t, Float_t> > FindLimitPixels(blob); 47 48 /* 49 * Check if a blob appears within a certain distance from the edges 50 * Parameters: 51 * - blob 52 * - vector containing farthest pixels in the blob (see FindLimitPixels(blob)) 53 * - Width of sensor in units of pixels 54 * - Height of sensor in units of pixels 55 */ 56 Bool_t BlobAtADistanceFromEdges(blob, vector<pair<Float_t, Float_t> >, Int_t, Int_t, Int_t, Int_t); 57 47 58 } 48 59 -
mafalda/MPXAlgo/MediPixAlgo.h
r179 r222 65 65 Int_t GetMatrixXdim(); 66 66 Int_t GetMatrixWidth(){return GetMatrixXdim();}; // same 67 Int_t GetWidth(){return GetMatrixXdim();}; // same 68 67 69 Int_t GetMatrixYdim(); 68 70 Int_t GetMatrixHeight(){return GetMatrixYdim();}; // same 71 Int_t GetHeight(){return GetMatrixYdim();}; // same 69 72 70 73 Int_t GetMatrixElement(Int_t col, Int_t row); -
mafalda/ShallowAngleFEIX/ShallowAngleFEIX.cpp
r219 r222 130 130 131 131 132 Log << MSG::INFO << " angle = " << (*blobsItr).GetBlobProperties().rotAngle << endreq;132 Log << MSG::INFO << "Mip angle = " << (*blobsItr).GetBlobProperties().rotAngle << endreq; 133 133 if( TMath::Abs( (*blobsItr).GetBlobProperties().rotAngle ) > TMath::Pi()/3.) 134 134 continue; 135 135 136 // Start counting 137 Log << MSG::DEBUG << "tail first = " << limitPixels[__TAIL_INDX].first << endreq; 138 Log << MSG::DEBUG << "head first = " << limitPixels[__HEAD_INDX].first << endreq; 139 Log << MSG::DEBUG << "tail second = " << limitPixels[__TAIL_INDX].second << endreq; 140 Log << MSG::DEBUG << "head second = " << limitPixels[__HEAD_INDX].second << endreq; 141 142 int initx = limitPixels[__HEAD_INDX].first; 143 int endx = limitPixels[__TAIL_INDX].first; 144 if ( initx > endx ) { 145 initx = limitPixels[__TAIL_INDX].first; 146 endx = limitPixels[__HEAD_INDX].first; 147 } 148 int inity = limitPixels[__HEAD_INDX].second; 149 int endy = limitPixels[__TAIL_INDX].second; 150 if ( inity > endy ) { 151 inity = limitPixels[__TAIL_INDX].second; 152 endy = limitPixels[__HEAD_INDX].second; 153 } 136 // Find the limit pixels to add-up the charge per column. 137 // It turns out that the limit pixels found above are not exactly 138 // what I need here. Often a delta ray comes out of the box marked 139 // by these limits and it could remain not included. 140 // I will use the limits given by the box. 141 142 Float_t midx = ((Float_t)(*blobsItr).GetBlobProperties().width_x) /2.; 143 Float_t midy = ((Float_t)(*blobsItr).GetBlobProperties().width_y) /2.; 144 int initx = TMath::FloorNint((*blobsItr).GetBlobProperties().geoCenter_x - midx); 145 int endx = TMath::FloorNint((*blobsItr).GetBlobProperties().geoCenter_x + midx - 1); 146 int inity = TMath::FloorNint((*blobsItr).GetBlobProperties().geoCenter_y - midy); 147 int endy = TMath::FloorNint((*blobsItr).GetBlobProperties().geoCenter_y + midy - 1); 148 149 Log << MSG::INFO << "Searching in box limited by x:" << initx << " , " << endx << endreq; 150 Log << MSG::INFO << " y:" << inity << " , " << endy << endreq; 154 151 155 152 m_clusterHeight = (*blobsItr).GetBlobProperties().width_y; … … 163 160 for (int y = inity ; y <= endy ; y++) { 164 161 165 Log << MSG::DEBUG << " search : " << x << " , " << y<< endreq;162 Log << MSG::DEBUG << " --> " << x << " , " << y << " : " << GetMatrixElement(x,y) << endreq; 166 163 m_totPerSegment[cntr] += GetMatrixElement(x,y); 167 164 m_chargeWeights[cntr] += GetMatrixElement(x,y); // normalized later … … 189 186 Log << MSG::DEBUG << "[" << segCntr << "] " 190 187 << "(TOT , fraction) : " << 191 *itrT << " , " << cout.precision(2) << (*itrW)*100188 *itrT << " , " << (*itrW)*100 << "\%" 192 189 << endreq; 193 190 … … 232 229 Int_t nExplicit = 0, nSoft = 0; 233 230 231 Log << MSG::INFO << "Starting delta ray search " << endreq; 232 234 233 for ( ; itr != m_chargeWeights.end() ; itr++ ) { 235 234 … … 241 240 unbalance = (*itr)/(*itrC); 242 241 if ( unbalance > m_dRayBalSearch ) { 243 Log << MSG:: INFO<< "section " << cntr << " much hotter than " << cntrE << endreq;242 Log << MSG::DEBUG << "section " << cntr << " much hotter than " << cntrE << endreq; 244 243 nNoticeableUnbalances[cntr]++; 245 244 } … … 258 257 for( ; itrU != nNoticeableUnbalances.end() ; itrU++){ 259 258 260 if ( (*itrU) == nDiv - 1 ) { // an explicit delta ray. Un balance to all other segments259 if ( (*itrU) == nDiv - 1 ) { // an explicit delta ray. Un-balance to all other segments 261 260 nExplicit++; 262 } else if ( (*itrU) == nDiv - 2 ) { 261 } else if ( (*itrU) == nDiv - 2 ) { // a soft. Unbalance to all-1 other segments 263 262 nSoft++; 264 263 }
Note: See TracChangeset
for help on using the changeset viewer.
