AraEvent/AraEventCalibrator.cxx
00001 00002 00003 00004 00005 00006 00007 00008 00009 #include "AraEventCalibrator.h" 00010 #include "UsefulIcrrStationEvent.h" 00011 #include "UsefulAtriStationEvent.h" 00012 #include "AraGeomTool.h" 00013 #include "TMath.h" 00014 #include "TGraph.h" 00015 #include "FFTtools.h" 00016 #include <iostream> 00017 #include <fstream> 00018 #include <cstring> 00019 #include <zlib.h> 00020 #include <cstdlib> 00021 00022 00023 Bool_t AraCalType::hasCableDelays(AraCalType::AraCalType_t calType) 00024 { 00025 if(calType==kFirstCalibPlusCables || calType==kSecondCalibPlusCables || calType==kSecondCalibPlusCablesUnDiplexed) 00026 return kTRUE; 00027 return kFALSE; 00028 } 00029 00030 Bool_t AraCalType::hasInterleaveCalib(AraCalType::AraCalType_t calType) 00031 { 00032 if(calType==kFirstCalibPlusCables || calType==kSecondCalibPlusCables || 00033 calType==kFirstCalib || calType==kSecondCalib || calType==kSecondCalibPlusCablesUnDiplexed) 00034 return kTRUE; 00035 return kFALSE; 00036 } 00037 00038 Bool_t AraCalType::hasBinWidthCalib(AraCalType::AraCalType_t calType) 00039 { 00040 if(calType>=kFirstCalib) 00041 return kTRUE; 00042 return kFALSE; 00043 } 00044 00045 00046 Bool_t AraCalType::hasClockAlignment(AraCalType::AraCalType_t calType) 00047 { 00048 if(calType==kSecondCalibPlusCables || 00049 calType==kSecondCalib || calType==kSecondCalibPlusCablesUnDiplexed) 00050 return kTRUE; 00051 return kFALSE; 00052 } 00053 00054 00055 00056 Bool_t AraCalType::hasPedestalSubtraction(AraCalType::AraCalType_t calType) 00057 { 00058 if(calType==kNoCalib) return kFALSE; 00059 return kTRUE; 00060 } 00061 00062 Bool_t AraCalType::hasCommonMode(AraCalType::AraCalType_t calType) 00063 { 00064 return kFALSE; 00065 if(calType==kNoCalib) return kFALSE; 00066 return kTRUE; 00067 } 00068 00069 Bool_t AraCalType::hasUnDiplexing(AraCalType::AraCalType_t calType) 00070 { 00071 if(calType==kSecondCalibPlusCablesUnDiplexed) return kTRUE; 00072 return kFALSE; 00073 } 00074 00075 00076 00077 ClassImp(AraEventCalibrator); 00078 00079 AraEventCalibrator * AraEventCalibrator::fgInstance=0; 00080 00081 00082 AraEventCalibrator::AraEventCalibrator() 00083 { 00084 00085 //Loop through and set the got ped / calib flags to zero. This assumes stationId's first n elements 00086 //are for ICRR stations and the remaining N-n are ATRI 00087 for(AraStationId_t station=0;station<ICRR_NO_STATIONS+ATRI_NO_STATIONS;station++){ 00088 // fprintf(stderr, "AraEventCalibrator::AraEventCalibrator() -- station %i\n", station); 00089 if(AraGeomTool::isIcrrStation(station)){ 00090 gotIcrrPedFile[AraGeomTool::getStationCalibIndex(station)]=0; 00091 gotIcrrCalibFile[AraGeomTool::getStationCalibIndex(station)]=0; 00092 } 00093 if(AraGeomTool::isAtriStation(station)){ 00094 fGotAtriPedFile[AraGeomTool::getStationCalibIndex(station)]=0; 00095 fGotAtriCalibFile[AraGeomTool::getStationCalibIndex(station)]=0; 00096 } 00097 00098 } 00099 00100 } 00101 00102 AraEventCalibrator::~AraEventCalibrator() { 00103 //Default Destructor 00104 } 00105 00106 //______________________________________________________________________________ 00107 AraEventCalibrator* AraEventCalibrator::Instance() 00108 { 00109 00110 // printf("AraEventCalibrator::Instance() creating an instance of the calibrator\n"); 00111 //static function 00112 if(fgInstance) 00113 return fgInstance; 00114 00115 fgInstance = new AraEventCalibrator(); 00116 00117 return fgInstance; 00118 } 00119 00120 00121 void AraEventCalibrator::setPedFile(char fileName[], AraStationId_t stationId) 00122 { 00123 00124 strncpy(IcrrPedFile[stationId],fileName,FILENAME_MAX); 00125 if(stationId==ARA_TESTBED){ 00126 // fprintf(stdout, "AraEventCalibrator::setPedFile() setting TestBed IcrrPedFile to %s\n", IcrrPedFile[stationId]); 00127 gotIcrrPedFile[0]=1; //Protects from loading the default pedestal File 00128 } 00129 00130 if(stationId==ARA_STATION1){ 00131 // fprintf(stdout, "AraEventCalibrator::setPedFile() setting Station1 IcrrPedFile to %s\n", IcrrPedFile[stationId]); 00132 gotIcrrPedFile[1]=1; //Protects from loading the default pedestal File 00133 } 00134 loadIcrrPedestals(stationId); 00135 } 00136 00137 void AraEventCalibrator::loadIcrrPedestals(AraStationId_t stationId) 00138 { 00139 00140 if(stationId==ARA_TESTBED&&gotIcrrPedFile[0]==0){ 00141 //populate the IcrrPedFile[0] variable 00142 char *pedFileEnv = getenv( "ARA_PEDESTAL_FILE" ); 00143 if ( pedFileEnv == NULL ) { 00144 char calibDir[FILENAME_MAX]; 00145 char *calibEnv=getenv("ARA_CALIB_DIR"); 00146 if(!calibEnv) { 00147 char *utilEnv=getenv("ARA_UTIL_INSTALL_DIR"); 00148 if(!utilEnv) { 00149 sprintf(calibDir,"calib"); 00150 } else { 00151 sprintf(calibDir,"%s/share/araCalib",utilEnv); 00152 } 00153 } 00154 else { 00155 strncpy(calibDir,calibEnv,FILENAME_MAX); 00156 } 00157 sprintf(IcrrPedFile[0],"%s/ICRR/TestBed/peds_1294924296.869787.run001202.dat",calibDir); 00158 } // end of IF-block for pedestal file not specified by environment variable 00159 else { 00160 strncpy(IcrrPedFile[0],pedFileEnv,FILENAME_MAX); 00161 } // end of IF-block for pedestal file specified by environment variable 00162 00163 00164 } 00165 if(stationId==ARA_STATION1&&gotIcrrPedFile[1]==0){ 00166 //populate the IcrrPedFile[1] variable 00167 char *pedFileEnv = getenv( "ARA_PEDESTAL_FILE" ); 00168 if ( pedFileEnv == NULL ) { 00169 char calibDir[FILENAME_MAX]; 00170 char *calibEnv=getenv("ARA_CALIB_DIR"); 00171 if(!calibEnv) { 00172 char *utilEnv=getenv("ARA_UTIL_INSTALL_DIR"); 00173 if(!utilEnv) { 00174 sprintf(calibDir,"calib"); 00175 } else { 00176 sprintf(calibDir,"%s/share/araCalib",utilEnv); 00177 } 00178 } 00179 else { 00180 strncpy(calibDir,calibEnv,FILENAME_MAX); 00181 } 00182 sprintf(IcrrPedFile[1],"%s/ICRR/Station1/peds_1326108401.602169.run003747.dat",calibDir); 00183 } // end of IF-block for pedestal file not specified by environment variable 00184 else { 00185 strncpy(IcrrPedFile[1],pedFileEnv,FILENAME_MAX); 00186 } // end of IF-block for pedestal file specified by environment variable 00187 00188 } 00189 //Finished populating the IcrrPedFile variables 00190 00191 FullLabChipPedStruct_t peds; 00192 00193 //Now we load in the pedestals 00194 if(stationId==ARA_TESTBED){ 00195 gzFile inPedTB = gzopen(IcrrPedFile[0],"r"); 00196 fprintf(stdout, "%s : loading IcrrPedFile TestBed %s\n", __FUNCTION__, IcrrPedFile[0]);//DEBUG 00197 00198 if( !inPedTB ){ 00199 fprintf(stderr,"ERROR - Failed to open pedestal file for TestBed %s.\n",IcrrPedFile[0]); 00200 return; 00201 } 00202 int nRead = gzread(inPedTB,&peds,sizeof(FullLabChipPedStruct_t)); 00203 if( nRead != sizeof(FullLabChipPedStruct_t)){ 00204 int numErr; 00205 fprintf(stderr,"ERROR - Error reading pedestal file %s; %s\n",IcrrPedFile[0],gzerror(inPedTB,&numErr)); 00206 gzclose(inPedTB); 00207 return; 00208 } 00209 00210 int chip,chan,samp; 00211 for(chip=0;chip<LAB3_PER_ICRR;++chip) { 00212 for(chan=0;chan<CHANNELS_PER_LAB3;++chan) { 00213 int chanIndex = chip*CHANNELS_PER_LAB3+chan; 00214 for(samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;++samp) { 00215 pedestalData[0][chip][chan][samp]=peds.chan[chanIndex].pedMean[samp]; 00216 } 00217 } 00218 } 00219 gzclose(inPedTB); 00220 00221 gotIcrrPedFile[0]=1; 00222 gotIcrrPedFile[1]=0; 00223 //last thing is to set the gotIcrrPedFile flags appropriately 00224 } 00225 if(stationId==ARA_STATION1){ 00226 gzFile inPedAra1 = gzopen(IcrrPedFile[1],"r"); 00227 fprintf(stdout, "%s : loading IcrrPedFile Station1 %s\n", __FUNCTION__, IcrrPedFile[1]);//DEBUG 00228 if( !inPedAra1 ){ 00229 fprintf(stderr,"ERROR - Failed to open pedestal file for Station1 %s\n",IcrrPedFile[1]); 00230 return; 00231 } 00232 00233 int nRead = gzread(inPedAra1,&peds,sizeof(FullLabChipPedStruct_t)); 00234 if( nRead != sizeof(FullLabChipPedStruct_t)){ 00235 int numErr; 00236 fprintf(stderr,"ERROR - Error reading pedestal file %s; %s\n",IcrrPedFile[1],gzerror(inPedAra1,&numErr)); 00237 gzclose(inPedAra1); 00238 return; 00239 } 00240 00241 int chip, chan, samp; 00242 for(chip=0;chip<LAB3_PER_ICRR;++chip) { 00243 for(chan=0;chan<CHANNELS_PER_LAB3;++chan) { 00244 int chanIndex = chip*CHANNELS_PER_LAB3+chan; 00245 for(samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;++samp) { 00246 pedestalData[1][chip][chan][samp]=peds.chan[chanIndex].pedMean[samp]; 00247 } 00248 } 00249 } 00250 gzclose(inPedAra1); 00251 gotIcrrPedFile[0]=0; 00252 gotIcrrPedFile[1]=1; 00253 //last thing is to set the gotIcrrPedFile flags appropriately 00254 } 00255 //Finished loading pedestals 00256 00257 00258 00259 } 00260 00261 int AraEventCalibrator::doBinCalibration(UsefulIcrrStationEvent *theEvent, int chanIndex,int overrideRCO) 00262 { 00263 AraStationId_t stationId=theEvent->stationId; 00264 00265 int nChip=theEvent->chan[chanIndex].chanId/CHANNELS_PER_LAB3; 00266 int nChan=theEvent->chan[chanIndex].chanId%CHANNELS_PER_LAB3; 00267 int rco=overrideRCO; 00268 if(overrideRCO!=0 && overrideRCO!=1) 00269 rco=theEvent->getRCO(chanIndex); 00270 int hbwrap = theEvent->chan[chanIndex].chipIdFlag&0x08; 00271 char hbextra = (theEvent->chan[chanIndex].chipIdFlag&0xf0)>>4; 00272 short hbstart=theEvent->chan[chanIndex].firstHitbus; 00273 short hbend=theEvent->chan[chanIndex].lastHitbus+hbextra; 00274 00275 double calTime=0; 00276 for(int samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;++samp){ 00277 rawadc[samp]=theEvent->chan[chanIndex].data[samp]; 00278 if(theEvent->chan[chanIndex].data[samp]==0){ 00279 calwv[samp]=0; 00280 pedsubadc[samp]=0; 00281 }else{ 00282 pedsubadc[samp]=rawadc[samp]-pedestalData[stationId][nChip][nChan][samp]; 00283 calwv[samp]=pedsubadc[samp]*ADCMV; 00284 if(calwv[samp]>SATURATION) calwv[samp]=SATURATION; 00285 if(calwv[samp]<-1*SATURATION) calwv[samp]=-1*SATURATION; 00286 } 00287 } 00288 00289 00290 if(nChan==0 || (overrideRCO==0 || overrideRCO==1)) { 00291 //Do the calibration for each chip 00292 calTime=0; 00293 int ir=0; 00294 if(hbwrap){ // Wrapped hitbus 00295 for(int samp=hbstart+1;samp<hbend;++samp) { 00296 tempTimeNums[ir++]=calTime; 00297 calTime+=binWidths[stationId][nChip][rco][samp]; 00298 } 00299 } 00300 else{ 00301 for(int samp=hbend+1;samp<MAX_NUMBER_SAMPLES_LAB3;++samp) { 00302 tempTimeNums[ir++]=calTime; 00303 calTime+=binWidths[stationId][nChip][1-rco][samp]; 00304 } 00305 //Now add epsilon 00306 calTime+=epsilonVals[stationId][nChip][rco]; //Need to check if this is rco or 1-rco 00307 for(int samp=0;samp<hbstart;++samp) { 00308 tempTimeNums[ir++]=calTime; 00309 calTime+=binWidths[stationId][nChip][rco][samp]; 00310 } 00311 } 00312 for(int samp=ir;samp<MAX_NUMBER_SAMPLES_LAB3;++samp) { 00313 tempTimeNums[samp]=calTime; 00314 calTime+=1; 00315 } 00316 TMath::Sort(MAX_NUMBER_SAMPLES_LAB3,tempTimeNums,indexNums,kFALSE); 00317 } 00318 // ... and rotate 00319 int ir=0; 00320 int numValid=0; 00321 if(hbwrap){ // Wrapped hitbus 00322 for(int samp=hbstart+1;samp<hbend;++samp) 00323 v[ir++]=calwv[samp]; 00324 }else{ 00325 for(int samp=hbend+1;samp<MAX_NUMBER_SAMPLES_LAB3;++samp) 00326 v[ir++]=calwv[samp]; 00327 for(int samp=0;samp<hbstart;++samp) 00328 v[ir++]=calwv[samp]; 00329 } 00330 numValid=ir; 00331 // Fill in remaining bins with zeros 00332 for(int samp=ir;samp<MAX_NUMBER_SAMPLES_LAB3;++samp) 00333 v[samp]=0; 00334 00335 //Now we just have to make sure the times are monotonically increasing 00336 for(int i=0;i<MAX_NUMBER_SAMPLES_LAB3;i++) { 00337 calTimeNums[i]=tempTimeNums[indexNums[i]]; 00338 calVoltNums[i]=v[indexNums[i]]; 00339 // std::cout << i << "\t" << indexNums[i] << "\t" << calTimeNums[i] << "\t" << calVoltNums[i] << "\n"; 00340 } 00341 return numValid; 00342 } 00343 00344 00345 void AraEventCalibrator::calibrateEvent(UsefulIcrrStationEvent *theEvent, AraCalType::AraCalType_t calType) 00346 { 00347 00348 AraStationId_t stationId=theEvent->stationId; 00349 if(gotIcrrPedFile[stationId]==0) 00350 { 00351 loadIcrrPedestals(stationId); //This will only load the pedestals once 00352 } 00353 00354 if(gotIcrrCalibFile[stationId]==0){ 00355 loadIcrrCalib(stationId); //This will only load the calib values once 00356 } 00357 00358 if(AraCalType::hasBinWidthCalib(calType)) 00359 theEvent->guessRCO(0); //Forces the calculation of the RCO phase from the clock 00360 //FIXME -- should this be a zero or do we want the actual channel number? 00361 00362 00363 //set the number of rfchans in the event 00364 00365 if(stationId==ARA_TESTBED) theEvent->numRFChans=RFCHANS_TESTBED; 00366 else theEvent->numRFChans=RFCHANS_STATION1; 00367 00368 // Calibrate waveforms 00369 for(int samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;++samp){ 00370 sampNums[samp]=samp; 00371 timeNums[samp]=samp*NSPERSAMP; 00372 } 00373 00374 for(int chanIndex = 0; chanIndex < NUM_DIGITIZED_ICRR_CHANNELS; ++chanIndex ){ 00375 int numValid=doBinCalibration(theEvent,chanIndex); 00376 00377 //Now we stuff it back into the UsefulIcrrStationEvent object 00378 00379 if(calType==AraCalType::kNoCalib) { 00380 theEvent->fNumPoints[chanIndex]=MAX_NUMBER_SAMPLES_LAB3; 00381 for(int samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;samp++) { 00382 theEvent->fVolts[chanIndex][samp]=rawadc[samp]; 00383 theEvent->fTimes[chanIndex][samp]=sampNums[samp]; 00384 } 00385 } 00386 if(calType==AraCalType::kJustUnwrap || calType==AraCalType::kADC) { 00387 theEvent->fNumPoints[chanIndex]=numValid; 00388 for(int samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;samp++) { 00389 theEvent->fTimes[chanIndex][samp]=sampNums[samp]; 00390 if(samp<numValid) { 00391 theEvent->fVolts[chanIndex][samp]=v[samp]; 00392 } 00393 else { 00394 theEvent->fVolts[chanIndex][samp]=0; 00395 } 00396 } 00397 } 00398 if(calType==AraCalType::kVoltageTime) { 00399 theEvent->fNumPoints[chanIndex]=numValid; 00400 for(int samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;samp++) { 00401 theEvent->fTimes[chanIndex][samp]=timeNums[samp]; 00402 if(samp<numValid) { 00403 theEvent->fVolts[chanIndex][samp]=v[samp]; 00404 } 00405 else { 00406 theEvent->fVolts[chanIndex][samp]=0; 00407 } 00408 } 00409 } 00410 if(calType==AraCalType::kJustPed) { 00411 theEvent->fNumPoints[chanIndex]=MAX_NUMBER_SAMPLES_LAB3; 00412 for(int samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;samp++) { 00413 theEvent->fVolts[chanIndex][samp]=pedsubadc[samp]; 00414 theEvent->fTimes[chanIndex][samp]=sampNums[samp]; 00415 } 00416 } 00417 if(AraCalType::hasBinWidthCalib(calType)) { 00418 //Almost always want this 00419 theEvent->fNumPoints[chanIndex]=numValid; 00420 for(int samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;samp++) { 00421 theEvent->fTimes[chanIndex][samp]=calTimeNums[samp]; 00422 if(samp<numValid) { 00423 theEvent->fVolts[chanIndex][samp]=calVoltNums[samp]; 00424 } 00425 else { 00426 theEvent->fVolts[chanIndex][samp]=0; 00427 } 00428 } 00429 } 00430 } 00431 00432 //Now we have done the initial bin-by-bin and epsilon calibrations 00433 //Next up is to do the clock alignment 00434 if(AraCalType::hasClockAlignment(calType)) { 00435 //All of the higher calibrations do some form of clock alignment 00436 calcClockAlignVals(theEvent,calType); 00437 for(int chanIndex = 0; chanIndex < NUM_DIGITIZED_ICRR_CHANNELS; ++chanIndex ){ 00438 int nChip=theEvent->chan[chanIndex].chanId/CHANNELS_PER_LAB3; 00439 for(int samp=0;samp<MAX_NUMBER_SAMPLES_LAB3;samp++) { 00440 theEvent->fTimes[chanIndex][samp]+=clockAlignVals[stationId][nChip]; 00441 } 00442 } 00443 } 00444 00445 //For now we just have the one calibration type for interleaving 00446 AraGeomTool *tempGeom = AraGeomTool::Instance(); 00447 00448 for(int rfchan = 0; rfchan < (theEvent->numRFChans); ++rfchan ){ 00449 // printf("AraEventCalibrator::calibrateEvent() rfchan %i stationId %i\n", rfchan, stationId); 00450 00451 memset(theEvent->fVoltsRF[rfchan],0,sizeof(Double_t)*2*MAX_NUMBER_SAMPLES_LAB3); 00452 memset(theEvent->fTimesRF[rfchan],0,sizeof(Double_t)*2*MAX_NUMBER_SAMPLES_LAB3); 00453 00454 00455 if(tempGeom->getNumLabChansForChan(rfchan, stationId)==2) { 00456 // printf("interleaved channel\n"); 00457 // std::cout << rfchan << "\t" 00458 // << tempGeom->getFirstLabChanIndexForChan(rfchan, stationId) << "\t" 00459 // << tempGeom->getSecondLabChanIndexForChan(rfchan, stationId) << "\n"; 00460 00461 //printf("rfchan %i stationId %i labChip %i FirstLabChan %i SecondLabChan %i numLabChans %i labChans[0] %i labChans[1] %i\n", rfchan, stationId, tempGeom->getLabChipForChan(rfchan, stationId), tempGeom->getFirstLabChanForChan(rfchan, stationId), tempGeom->getSecondLabChanForChan(rfchan, stationId),tempGeom->getNumLabChansForChan(rfchan, stationId), tempGeom->getFirstLabChanIndexForChan(rfchan, stationId), tempGeom->getSecondLabChanIndexForChan(rfchan, stationId)); 00462 00463 int ci1=tempGeom->getFirstLabChanIndexForChan(rfchan, stationId); 00464 int ci2=tempGeom->getSecondLabChanIndexForChan(rfchan, stationId); 00465 theEvent->fNumPointsRF[rfchan]=theEvent->fNumPoints[ci1]+theEvent->fNumPoints[ci2]; 00466 //Need to zero mean, maybe should do this in each half seperately? 00467 //RJN hack 13/01/11 00468 double mean=0; 00469 for(int i=0;i<theEvent->fNumPoints[ci1];i++) { 00470 mean+=theEvent->fVolts[ci1][i]; 00471 } 00472 for(int i=0;i<theEvent->fNumPoints[ci2];i++) { 00473 mean+=theEvent->fVolts[ci2][i]; 00474 } 00475 mean/=theEvent->fNumPointsRF[rfchan]; 00476 00477 //Only the interleaved types have any special calib here 00478 if(AraCalType::hasInterleaveCalib(calType)) { 00479 int i1=0; 00480 int i2=0; 00481 for(int i=0;i<theEvent->fNumPointsRF[rfchan];i++) { 00482 if(i1<theEvent->fNumPoints[ci1] && i2<theEvent->fNumPoints[ci2]) { 00483 //Both in play 00484 if(theEvent->fTimes[ci1][i1]<(theEvent->fTimes[ci2][i2]+interleaveVals[stationId][rfchan])) { 00485 theEvent->fTimesRF[rfchan][i]=theEvent->fTimes[ci1][i1]; 00486 theEvent->fVoltsRF[rfchan][i]=theEvent->fVolts[ci1][i1]-mean; 00487 // std::cout << "A: " << i << "\t" << theEvent->fTimesRF[rfchan][i] << "\n"; 00488 00489 i1++; 00490 continue; 00491 } 00492 else { 00493 theEvent->fTimesRF[rfchan][i]=theEvent->fTimes[ci2][i2]+interleaveVals[stationId][rfchan]; 00494 theEvent->fVoltsRF[rfchan][i]=theEvent->fVolts[ci2][i2]-mean; 00495 // std::cout << "B: " << i << "\t" << theEvent->fTimesRF[rfchan][i] << "\n"; 00496 i2++; 00497 continue; 00498 } 00499 } 00500 else if(i1<theEvent->fNumPoints[ci1]) { 00501 theEvent->fTimesRF[rfchan][i]=theEvent->fTimes[ci1][i1]; 00502 theEvent->fVoltsRF[rfchan][i]=theEvent->fVolts[ci1][i1]-mean; 00503 i1++; 00504 // std::cout << "C: " << i << "\t" << theEvent->fTimesRF[rfchan][i] << "\n"; 00505 continue; 00506 } 00507 else if(i2<theEvent->fNumPoints[ci2]) { 00508 00509 theEvent->fTimesRF[rfchan][i]=theEvent->fTimes[ci2][i2]+interleaveVals[stationId][rfchan]; 00510 theEvent->fVoltsRF[rfchan][i]=theEvent->fVolts[ci2][i2]-mean; 00511 // std::cout << "D: " << i << "\t" << theEvent->fTimesRF[rfchan][i] << "\n"; 00512 i2++; 00513 continue; 00514 } 00515 } 00516 }//interleaved 00517 00518 else { 00519 for(int i=0;i<theEvent->fNumPointsRF[rfchan];i++) { 00520 if(i%2==0) { 00521 //ci2 00522 theEvent->fVoltsRF[rfchan][i]=theEvent->fVolts[ci2][i/2]; 00523 theEvent->fTimesRF[rfchan][i]=theEvent->fTimes[ci2][i/2]; 00524 } 00525 else { 00526 //ci1 00527 theEvent->fVoltsRF[rfchan][i]=theEvent->fVolts[ci1][i/2]; 00528 theEvent->fTimesRF[rfchan][i]=theEvent->fTimes[ci1][i/2]+0.5*NSPERSAMP; 00529 } 00530 00531 } 00532 } 00533 00534 }//numLabChans==2 00535 00536 //Perform undiplexing of channels (ARA1 has surface antennae diplexed in to HPOL channels) 00537 else if(tempGeom->isDiplexed(rfchan, stationId)==0||!(AraCalType::hasUnDiplexing(calType))){ 00538 // printf("non-interleaved non-diplexed channel\n"); 00539 // std::cout << rfchan << "\t" 00540 // << tempGeom->getFirstLabChanIndexForChan(rfchan) << "\n"; 00541 int ci=tempGeom->getFirstLabChanIndexForChan(rfchan, stationId); 00542 theEvent->fNumPointsRF[rfchan]=theEvent->fNumPoints[ci]; 00543 00544 //Need to zero mean the data 00545 double mean=0; 00546 for(int i=0;i<theEvent->fNumPoints[ci];i++) { 00547 mean+=theEvent->fVolts[ci][i]; 00548 } 00549 mean/=theEvent->fNumPoints[ci]; 00550 00551 for(int i=0;i<theEvent->fNumPointsRF[rfchan];i++) { 00552 theEvent->fVoltsRF[rfchan][i]=theEvent->fVolts[ci][i]-mean; 00553 theEvent->fTimesRF[rfchan][i]=theEvent->fTimes[ci][i]; 00554 00555 } 00556 }//numLabChannels!=2&&isDiplexed==0 00557 00558 else if(tempGeom->isDiplexed(rfchan, stationId)==1&&(AraCalType::hasUnDiplexing(calType))){ 00559 // printf("non-interleaved diplexed channel\n"); 00560 //now need to work out how to do a diplexed channel 00561 //find lab channel from which to derived rfchan 00562 int ci=tempGeom->getFirstLabChanIndexForChan(rfchan,stationId); 00563 int noPoints=theEvent->fNumPoints[ci]; 00564 Double_t lowPassFreq=tempGeom->getLowPassFilter(rfchan,stationId); 00565 Double_t highPassFreq=tempGeom->getHighPassFilter(rfchan,stationId); 00566 00567 TGraph *gUnfiltered = new TGraph(noPoints, theEvent->fTimes[ci], theEvent->fVolts[ci]); 00568 00569 TGraph *gFiltered = FFTtools::simplePassBandFilter(gUnfiltered, highPassFreq, lowPassFreq); 00570 Double_t *filteredT = gFiltered->GetX(); 00571 Double_t *filteredV = gFiltered->GetY(); 00572 00573 00574 for(int i=0;i<noPoints;i++){ 00575 theEvent->fTimesRF[rfchan][i]=filteredT[i]; 00576 theEvent->fVoltsRF[rfchan][i]=filteredV[i]; 00577 } 00578 theEvent->fNumPointsRF[rfchan]=noPoints; 00579 00580 delete gFiltered; 00581 delete gUnfiltered; 00582 00583 }//isDiplexed==1 00584 00585 if(AraCalType::hasCableDelays(calType)) { 00586 Double_t delay=tempGeom->fStationInfo[stationId].fAntInfo[rfchan].cableDelay; 00587 // delay-=180; //Just an arbitrary offset 00588 for(int i=0;i<theEvent->fNumPointsRF[rfchan];i++) { 00589 theEvent->fTimesRF[rfchan][i]-=delay; 00590 } 00591 } 00592 00593 00594 } 00595 00596 } 00597 00598 void AraEventCalibrator::loadIcrrCalib(AraStationId_t stationId) 00599 { 00600 char calibDir[FILENAME_MAX]; 00601 char binWidthFileName[FILENAME_MAX]; 00602 char epsilonFileName[FILENAME_MAX]; 00603 char interleaveFileName[FILENAME_MAX]; 00604 char *calibEnv=getenv("ARA_CALIB_DIR"); 00605 if(!calibEnv) { 00606 char *utilEnv=getenv("ARA_UTIL_INSTALL_DIR"); 00607 if(!utilEnv) 00608 sprintf(calibDir,"calib"); 00609 else 00610 sprintf(calibDir,"%s/share/araCalib",utilEnv); 00611 } 00612 else { 00613 strncpy(calibDir,calibEnv,FILENAME_MAX); 00614 } 00615 //Populate the binWidthFileName, epsilonFileName and interleaveFileName variables 00616 if(stationId==ARA_TESTBED&&gotIcrrCalibFile[0]==0){ 00617 sprintf(binWidthFileName,"%s/ICRR/TestBed/binWidths.txt",calibDir); 00618 sprintf(epsilonFileName,"%s/ICRR/TestBed/epsilonFile.txt",calibDir); 00619 sprintf(interleaveFileName,"%s/ICRR/TestBed/interleaveFile.txt",calibDir); 00620 } 00621 if(stationId==ARA_STATION1&&gotIcrrCalibFile[1]==0){ 00622 sprintf(binWidthFileName,"%s/ICRR/Station1/binWidths.txt",calibDir); 00623 sprintf(epsilonFileName,"%s/ICRR/Station1/epsilonFile.txt",calibDir); 00624 sprintf(interleaveFileName,"%s/ICRR/Station1/interleaveFile.txt",calibDir); 00625 } 00626 //load the calibration files 00627 00628 std::ifstream binWidthFile(binWidthFileName); 00629 std::ifstream epsilonFile(epsilonFileName); 00630 std::ifstream interleaveFile(interleaveFileName); 00631 int chip, rco, chan; 00632 double width, epsilon, interleave; 00633 while(binWidthFile >> chip >> rco){ 00634 for(int samp=0; samp<MAX_NUMBER_SAMPLES_LAB3;++samp){ 00635 binWidthFile >> width; 00636 binWidths[stationId][chip][rco][samp]=width; 00637 } 00638 } 00639 binWidthFile.close(); 00640 while(epsilonFile >> chip >> rco >> epsilon){ 00641 epsilonVals[stationId][chip][rco]=epsilon; 00642 } 00643 epsilonFile.close(); 00644 while(interleaveFile >> chip >> chan >> interleave){ 00645 interleaveVals[stationId][chan+4*chip]=interleave; 00646 } 00647 interleaveFile.close(); 00648 00649 //Set the gotIcrrCalibFile flags 00650 if(stationId==ARA_TESTBED){ 00651 gotIcrrCalibFile[0]=1; 00652 gotIcrrCalibFile[1]=0; 00653 } 00654 if(stationId==ARA_STATION1){ 00655 gotIcrrCalibFile[0]=0; 00656 gotIcrrCalibFile[1]=1; 00657 } 00658 00659 00660 } 00661 00662 00663 void AraEventCalibrator::calcClockAlignVals(UsefulIcrrStationEvent *theEvent, AraCalType::AraCalType_t calType) 00664 { 00665 AraStationId_t stationId=theEvent->stationId; 00666 00667 if(!AraCalType::hasClockAlignment(calType)) return; 00668 TGraph *grClock[LAB3_PER_ICRR]={0}; 00669 Double_t lag[LAB3_PER_ICRR]={0}; 00670 for(int chip=0;chip<LAB3_PER_ICRR;chip++) { 00671 clockAlignVals[stationId][chip]=0; 00672 int chanIndex=ICRR1_CLOCK_CHANNEL+CHANNELS_PER_LAB3*chip; 00673 grClock[chip]=theEvent->getGraphFromElecChan(chanIndex); 00674 lag[chip]=estimateClockLag(grClock[chip]); 00675 delete grClock[chip]; 00676 00677 if(chip>0) { 00678 //Then can actually do some alignment 00679 clockAlignVals[stationId][chip]=lag[0]-lag[chip]; 00680 //The below fudge factors were "tuned" using pulser data 00681 // to try and remove period ambiguities resulting from wrong cycle lag 00682 if(lag[chip]<8 && lag[0]>9) 00683 clockAlignVals[stationId][chip]-=25; 00684 if(lag[chip]>9 && lag[0]<7) 00685 clockAlignVals[stationId][chip]+=25; 00686 // std::cout << "clockAlignVals[ " << chip << "] = " << clockAlignVals[chip] << "\n"; 00687 } 00688 00689 } 00690 } 00691 00692 00693 Double_t AraEventCalibrator::estimateClockLag(TGraph *grClock) 00694 { 00695 // This funciton estimates the clock lag (i.e. phase but expressed in terms of a deltaT between 0 and 1 period) by just using all the negative-positive zero crossing 00696 00697 Double_t period=ICRR1_CLOCK_PERIOD; 00698 Double_t mean=grClock->GetMean(2); 00699 Int_t numPoints=grClock->GetN(); 00700 if(numPoints<3) return 0; 00701 Double_t *tVals=grClock->GetX(); 00702 Double_t *vVals=grClock->GetY(); 00703 00704 Double_t zc[1000]={0}; 00705 Double_t rawZc[1000]={0}; 00706 int countZC=0; 00707 for(int i=2;i<numPoints;i++) { 00708 if(vVals[i-1]<0 && vVals[i]>0) { 00709 Double_t x1=tVals[i-1]; 00710 Double_t x2=tVals[i]; 00711 Double_t y1=vVals[i-1]-mean; 00712 Double_t y2=vVals[i]-mean; 00713 // std::cout << i << "\t" << y2 << "\t" << y1 << "\t" << (y2-y1) << "\n"; 00714 zc[countZC]=(((0-y1)/(y2-y1))*(x2-x1))+x1; 00715 rawZc[countZC]=zc[countZC]; 00716 countZC++; 00717 // if(countZC==1) 00718 // break; 00719 } 00720 00721 } 00722 00723 Double_t firstZC=zc[0]; 00724 while(firstZC>period) firstZC-=period; 00725 Double_t meanZC=0; 00726 Double_t meanZC2=0; 00727 for(int i=0;i<countZC;i++) { 00728 while((zc[i])>period) zc[i]-=period; 00729 if(TMath::Abs((zc[i]-period)-firstZC)<TMath::Abs(zc[i]-firstZC)) 00730 zc[i]-=period; 00731 if(TMath::Abs((zc[i]+period)-firstZC)<TMath::Abs(zc[i]-firstZC)) 00732 zc[i]+=period; 00733 meanZC+=zc[i]; 00734 meanZC2+=zc[i]*zc[i]; 00735 // std::cout << i << "\t" << zc[i] << "\t" << rawZc[i] << "\n"; 00736 } 00737 meanZC/=countZC; 00738 meanZC2/=countZC; 00739 // Double_t rms=TMath::Sqrt(meanZC2-meanZC*meanZC); 00740 // std::cout << meanZC << "\t" << rms << "\n"; 00741 return meanZC; 00742 00743 } 00744 00745 void AraEventCalibrator::fillRCOGuessArray(UsefulIcrrStationEvent *theEvent, int rcoGuess[LAB3_PER_ICRR]) 00746 { 00747 for(int chip=0;chip<LAB3_PER_ICRR;chip++) { 00748 int chanIndex=ICRR1_CLOCK_CHANNEL+CHANNELS_PER_LAB3*chip; 00749 rcoGuess[chip]=theEvent->getRawRCO(chanIndex); 00750 00751 Double_t period[2]={0}; 00752 Double_t rms[2]={0}; 00753 for(int rcoTest=0;rcoTest<2;rcoTest++) { 00754 int numValid=doBinCalibration(theEvent,chanIndex,rcoTest); 00755 period[rcoTest]=estimateClockPeriod(numValid,rms[rcoTest]); 00756 } 00757 00758 Double_t periodTest=(TMath::Abs(period[0]-ICRR1_CLOCK_PERIOD)-TMath::Abs(period[1]-ICRR1_CLOCK_PERIOD)); 00759 Double_t rmsTest=(rms[0]-rms[1]); 00760 if(periodTest<0) { 00761 rcoGuess[chip]=0; 00762 if(TMath::Abs(periodTest)<0.5 && rmsTest>0) 00763 rcoGuess[chip]=1; 00764 } 00765 else { 00766 rcoGuess[chip]=1; 00767 if(TMath::Abs(periodTest)<0.5 && rmsTest<0) { 00768 rcoGuess[chip]=0; 00769 } 00770 } 00771 if(rms[rcoGuess[chip]]>4) { 00772 rcoGuess[chip]=theEvent->getRawRCO(chanIndex); 00773 } 00774 00775 // std::cout << "AraEventCalibrator:\t" << period[0] << "\t" << period[1] << "\n"; 00776 // std::cout << "AraEventCalibrator:\t" << periodTest << "\t" << rmsTest << "\t" << rcoGuess[chip] << "\n"; 00777 } 00778 } 00779 00780 00781 00782 Double_t AraEventCalibrator::estimateClockPeriod(Int_t numPoints, Double_t &rms) 00783 { 00784 // This funciton estimates the period by just using all the negative-positive zero crossing 00785 if(numPoints<3) return 0; 00786 Double_t mean=0; 00787 Double_t vVals[MAX_NUMBER_SAMPLES_LAB3]={0}; 00788 for(int i=0;i<numPoints;i++) { 00789 mean+=calVoltNums[i]; 00790 } 00791 mean/=numPoints; 00792 for(int i=0;i<numPoints;i++) { 00793 vVals[i]=calVoltNums[i]-mean; 00794 } 00795 00796 00797 Double_t zc[1000]={0}; 00798 Double_t periods[1000]={0}; 00799 int countZC=0; 00800 for(int i=2;i<numPoints;i++) { 00801 if(vVals[i-1]<0 && vVals[i]>0) { 00802 Double_t x1=calTimeNums[i-1]; 00803 Double_t x2=calTimeNums[i]; 00804 Double_t y1=vVals[i-1]; 00805 Double_t y2=vVals[i]; 00806 // std::cout << i << "\t" << y2 << "\t" << y1 << "\t" << (y2-y1) << "\n"; 00807 Double_t zcTime=(((0-y1)/(y2-y1))*(x2-x1))+x1; 00808 if(countZC>0) { 00809 if((zcTime-zc[countZC-1])<10) 00810 continue; 00811 } 00812 zc[countZC]=zcTime; 00813 countZC++; 00814 // if(countZC==1) 00815 // break; 00816 } 00817 } 00818 00819 if(countZC<2) return 0; 00820 00821 00822 int countPeriods=0; 00823 Double_t meanPeriod=0; 00824 Double_t meanPeriodSq=0; 00825 for(int i=1;i<countZC;i++) { 00826 periods[countPeriods]=zc[i]-zc[i-1]; 00827 meanPeriod+=periods[countPeriods]; 00828 meanPeriodSq+=periods[countPeriods]*periods[countPeriods]; 00829 countPeriods++; 00830 } 00831 meanPeriod/=countPeriods; 00832 meanPeriodSq/=countPeriods; 00833 rms=TMath::Sqrt(meanPeriodSq-meanPeriod*meanPeriod); 00834 00835 return meanPeriod; 00836 } 00837 00838 void AraEventCalibrator::calibrateEvent(UsefulAtriStationEvent *theEvent, AraCalType::AraCalType_t calType) 00839 { 00840 AraStationId_t thisStationId = theEvent->stationId; 00841 Int_t calibIndex = AraGeomTool::getStationCalibIndex(thisStationId); 00842 //Only load the pedestals / calib if they are not already loaded 00843 if(fGotAtriCalibFile[calibIndex]==0){ 00844 loadAtriCalib(thisStationId); 00845 } 00846 if(fGotAtriPedFile[calibIndex]==0){ 00847 loadAtriPedestals(thisStationId); 00848 } 00849 00850 //Step one is loop over the blocks 00851 00852 std::vector<RawAtriStationBlock>::iterator blockIt; 00853 std::vector< std::vector<UShort_t> >::iterator vecVecIt; 00854 std::map< Int_t, std::vector <Double_t> >::iterator timeMapIt; 00855 std::map< Int_t, std::vector <Double_t> >::iterator voltMapIt; 00856 std::map< Int_t, std::vector <Double_t> >::iterator voltMapIt2; 00857 std::vector< UShort_t >::iterator shortIt; 00858 for(blockIt = theEvent->blockVec.begin(); 00859 blockIt!=theEvent->blockVec.end(); 00860 blockIt++) { 00861 //Step two is determine the channel Ids 00862 Int_t irsChan[8]; 00863 Int_t numChans=0; 00864 for(Int_t bit=0;bit<8;bit++) { 00865 Int_t mask=(1<<bit); 00866 if((blockIt->channelMask)&mask) { 00867 irsChan[numChans]=bit; 00868 numChans++; 00869 } 00870 } 00871 // std::cout << "Got numChans " << numChans << "\n"; 00872 00873 //Step three is loop over the channels within a block 00874 Int_t uptoChan=0; 00875 for(vecVecIt=blockIt->data.begin(); 00876 vecVecIt!=blockIt->data.end(); 00877 vecVecIt++) { 00878 Int_t chanId=irsChan[uptoChan] | ((blockIt->channelMask&0x300)>>5); 00879 Int_t chan=irsChan[uptoChan]; 00880 Int_t dda=blockIt->getDda(); 00881 Int_t block=blockIt->getBlock(); 00882 Int_t capArray=blockIt->getCapArray(); 00883 // std::cout << "Got chanId " << chanId << "\t" << irsChan[uptoChan] << "\t" 00884 // << dda << "\t" << block << "\t" << RawAtriStationEvent::getPedIndex(dda,block,chan,0) << "\n"; 00885 uptoChan++; 00886 00887 //Step four is to check if we have already got this chanId 00888 00889 timeMapIt=theEvent->fTimes.find(chanId); 00890 Double_t time=0; 00891 Int_t firstTime=1; 00892 if(timeMapIt==theEvent->fTimes.end()) { 00893 //First time round for this channel 00894 std::vector <Double_t> tempTimes; 00895 std::vector <Double_t> tempVolts; 00896 //Now need to insert empty vector into map 00897 theEvent->fTimes.insert( std::pair< Int_t, std::vector <Double_t> >(chanId,tempTimes)); 00898 theEvent->fVolts.insert( std::pair< Int_t, std::vector <Double_t> >(chanId,tempVolts)); 00899 theEvent->fNumChannels++; 00900 } 00901 else { 00902 //Just get the last time 00903 time=timeMapIt->second.back(); 00904 // if(dda==1 &&chan==1) 00905 // std::cout << "Last time " << time << "\t" << timeMapIt->second.size() << "\n"; 00906 firstTime=0; 00907 } 00908 00909 //Set the iterators to point to the correct channel in the map 00910 timeMapIt=theEvent->fTimes.find(chanId); 00911 voltMapIt=theEvent->fVolts.find(chanId); 00912 00913 int samp=0; 00914 int index=0; 00915 //Now loop over the 64 samples 00916 Double_t tempTimes[SAMPLES_PER_BLOCK]; 00917 Double_t tempVolts[SAMPLES_PER_BLOCK]; 00918 Int_t voltIndex[SAMPLES_PER_BLOCK]; 00919 00920 //Here is the Epsilon calibration 00921 if(AraCalType::hasBinWidthCalib(calType)) { 00922 if(!firstTime) { 00923 //Add on the time between blocks 00924 time+=fAtriEpsilonTimes[dda][chan][capArray]; 00925 // if(dda==1 && chan==1) 00926 // std::cout << "Block " << time << "\t" << fAtriEpsilonTimes[dda][chan][capArray] << "\n"; 00927 00928 } 00929 } 00930 00931 for(shortIt=vecVecIt->begin(); 00932 shortIt!=vecVecIt->end(); 00933 shortIt++) { 00934 00935 if(AraCalType::hasBinWidthCalib(calType)) { 00936 00937 //Now need to work out where to place the sample 00938 voltIndex[samp]=fAtriSampleIndex[dda][chan][capArray][samp]; 00939 00940 //Now get the time 00941 tempTimes[samp]=time+fAtriSampleTimes[dda][chan][capArray][samp]; 00942 00943 // if(dda==1 && chan==1) { 00944 // // std::cout << dda << "\t" << chan << "\t" << capArray << "\t" << samp << "\t" << fAtriSampleTimes[dda][chan][capArray][samp] << "\n"; 00945 // std::cout << index << "\t" << tempTimes[samp] << "\t" << time << "\n"; 00946 // } 00947 00948 } 00949 else { 00950 time+=1; 00951 tempTimes[samp]=time; 00952 voltIndex[samp]=samp; 00953 } 00954 00955 //Now the voltage part 00956 if(!AraCalType::hasPedestalSubtraction(calType)) 00957 tempVolts[samp]=(*shortIt); 00958 else { 00959 tempVolts[samp]=(*shortIt)-(Int_t)fAtriPeds[RawAtriStationEvent::getPedIndex(dda,block,chan,samp)]; 00960 // if(dda==3 && samp<2) { 00961 // std::cout << (*shortIt) << "\t" << (Int_t)fAtriPeds[RawAtriStationEvent::getPedIndex(dda,block,chan,samp)] << "\n"; 00962 // } 00963 } 00964 samp++; 00965 } 00966 00967 for(samp=0;samp<SAMPLES_PER_BLOCK;samp++) { 00968 timeMapIt->second.push_back(tempTimes[samp]); 00969 voltMapIt->second.push_back(tempVolts[voltIndex[samp]]); //Filling with volts 00970 00971 } 00972 00973 00974 } 00975 } 00976 00977 if(hasCommonMode(calType)) { 00978 //Then we need to do a common mode correction 00979 //loop over dda 00980 //loop over chan 00981 //loop over times and subtract one 00982 00983 for(int dda=0;dda<DDA_PER_ATRI;dda++) { 00984 for(Int_t chan=0;chan<5;chan++) { 00985 Int_t chanId=chan+RFCHAN_PER_DDA*dda; 00986 voltMapIt=theEvent->fVolts.find(chanId); 00987 Int_t chanId2=5+RFCHAN_PER_DDA*dda; 00988 voltMapIt2=theEvent->fVolts.find(chanId2); 00989 00990 Int_t numPoints=(voltMapIt->second).size(); 00991 for(int samp=0;samp<numPoints;samp++) { 00992 voltMapIt->second[samp]-=voltMapIt2->second[samp]; 00993 } 00994 } 00995 } 00996 00997 } 00998 00999 // fprintf(stderr, "AraEventCalibrator::CalibrateEvent() -- finished calibrating event\n");//DEBUG 01000 01001 } 01002 01003 01004 01005 void AraEventCalibrator::setAtriPedFile(char *filename, AraStationId_t stationId) 01006 { 01007 Int_t calibIndex = AraGeomTool::getStationCalibIndex(stationId); 01008 strncpy(fAtriPedFile[calibIndex],filename,FILENAME_MAX); 01009 fGotAtriPedFile[calibIndex]=1; //Protects us from loading the default pedfile 01010 loadAtriPedestals(stationId); 01011 01012 } 01013 01014 void AraEventCalibrator::loadAtriPedestals(AraStationId_t stationId) 01015 { 01016 Int_t calibIndex = AraGeomTool::getStationCalibIndex(stationId); 01017 if(calibIndex==-1){ 01018 fprintf(stderr, "AraEventCalibrator::loadAtriPedestals -- ERROR Unknown stationId %i\n", stationId); 01019 exit(0); 01020 } 01021 01022 if(fGotAtriPedFile[calibIndex]==0){ 01023 char *pedFileEnv = getenv( "ARA_ATRI_PEDESTAL_FILE" ); //FIXME 01024 if ( pedFileEnv == NULL ) { 01025 char calibDir[FILENAME_MAX]; 01026 char *calibEnv=getenv("ARA_CALIB_DIR"); 01027 if(!calibEnv) { 01028 char *utilEnv=getenv("ARA_UTIL_INSTALL_DIR"); 01029 if(!utilEnv) { 01030 sprintf(calibDir,"calib"); 01031 // fprintf(stdout,"AraEventCalibrator::loadAtriPedestals(): INFO - Pedestal file [from ./calib]"); 01032 } else { 01033 sprintf(calibDir,"%s/share/araCalib",utilEnv); 01034 // fprintf(stdout,"AraEventCalibrator::loadAtriPedestals(): INFO - Pedestal file [from ARA_UTIL_INSTALL_DIR/share/calib]"); 01035 } 01036 } 01037 else { 01038 strncpy(calibDir,calibEnv,FILENAME_MAX); 01039 // fprintf(stdout,"AraEventCalibrator::loadAtriPedestals(): INFO - Pedestal file [from ARA_CALIB_DIR]"); 01040 } 01041 sprintf(fAtriPedFile[calibIndex],"%s/ATRI/araAtriStation%iPedestals.txt",calibDir, stationId); 01042 fprintf(stdout," = %s\n",fAtriPedFile[calibIndex]); 01043 } // end of IF-block for pedestal file not specified by environment variable 01044 else { 01045 strncpy(fAtriPedFile[calibIndex],pedFileEnv,FILENAME_MAX); 01046 // fprintf(stdout,"AraEventCalibrator::loadAtriPedestals(): INFO - Pedestal file [from ARA_ONE_PEDESTAL_FILE] = %s\n",fAtriPedFile[calibIndex]); 01047 } // end of IF-block for pedestal file specified by environment variable 01048 } 01049 01050 //Pedestal file 01051 01052 fprintf(stdout, "%s : Loading fAtriPedFile - %s\n", __FUNCTION__, fAtriPedFile[calibIndex]); 01053 std::ifstream PedFile(fAtriPedFile[calibIndex]); 01054 Int_t dda,block,chan; 01055 UShort_t pedVal; 01056 01057 //If there are pedestals in memory already then we need to delete them and replace with this station's 01058 if(fAtriPeds) delete fAtriPeds; 01059 01060 fAtriPeds = new UShort_t [DDA_PER_ATRI*BLOCKS_PER_DDA*RFCHAN_PER_DDA*SAMPLES_PER_BLOCK]; 01061 if(!fAtriPeds) { 01062 std::cerr << "Can not allocate memory for pedestal file\n"; 01063 exit(0); 01064 } 01065 01066 while(PedFile >> dda >> block >> chan) { 01067 for(int samp=0;samp<SAMPLES_PER_BLOCK;samp++) { 01068 PedFile >> pedVal; 01069 fAtriPeds[RawAtriStationEvent::getPedIndex(dda,block,chan,samp)]=pedVal; 01070 } 01071 } 01072 01073 01074 //Now we set the gotPedFile flags to indicate which station we have in memory 01075 01076 for(AraStationId_t station=0;station<ICRR_NO_STATIONS+ATRI_NO_STATIONS;station++){ 01077 if(AraGeomTool::isAtriStation(station)){ 01078 calibIndex=AraGeomTool::getStationCalibIndex(station); 01079 if(station==stationId) fGotAtriPedFile[calibIndex]=1; 01080 else fGotAtriPedFile[calibIndex]=0; 01081 } 01082 } 01083 01084 01085 } 01086 01087 01088 void AraEventCalibrator::loadAtriCalib(AraStationId_t stationId) 01089 { 01090 Int_t calibIndex = AraGeomTool::getStationCalibIndex(stationId); 01091 if(calibIndex==-1){ 01092 fprintf(stderr, "AraEventCalibrator::loadAtriCalib -- ERROR Unknown stationId %i\n", stationId); 01093 return; 01094 } 01095 01096 char calibFile[FILENAME_MAX]; 01097 char calibDir[FILENAME_MAX]; 01098 char *calibEnv=getenv("ARA_CALIB_DIR"); 01099 if(!calibEnv) { 01100 char *utilEnv=getenv("ARA_UTIL_INSTALL_DIR"); 01101 if(!utilEnv) 01102 sprintf(calibDir,"calib"); 01103 else 01104 sprintf(calibDir,"%s/share/araCalib",utilEnv); 01105 } 01106 else { 01107 strncpy(calibDir,calibEnv,FILENAME_MAX); 01108 } 01109 01110 int dda,chan,sample,capArray; 01111 sprintf(calibFile,"%s/ATRI/araAtriStation%iSampleTiming.txt",calibDir, stationId); 01112 // fprintf(stdout, "AraEventCalibrator::loadAtriCalib(): INFO - Calibration file = %s\n", calibFile);//DEBUG 01113 01114 std::ifstream SampleFile(calibFile); 01115 for(dda=0;dda<DDA_PER_ATRI;dda++) { 01116 for(chan=0;chan<RFCHAN_PER_DDA;chan++) { 01117 for(capArray=0;capArray<2;capArray++) { 01118 for(sample=0;sample<SAMPLES_PER_BLOCK;sample++) { 01119 fAtriSampleTimes[dda][chan][capArray][sample]=sample/3.2; 01120 fAtriSampleIndex[dda][chan][capArray][sample]=sample; 01121 } 01122 fAtriEpsilonTimes[dda][chan][capArray]=1/3.2; 01123 } 01124 } 01125 } 01126 01127 Double_t value; 01128 Int_t index; 01129 while(SampleFile >> dda >> chan >> capArray) { 01130 // std::cout << dda << "\t" << chan << "\t" << capArray << "\t"; 01131 for(sample=0;sample<SAMPLES_PER_BLOCK;sample++) { 01132 SampleFile >> index; 01133 fAtriSampleIndex[dda][chan][capArray][sample]=index; 01134 // std::cout << fAtriSampleIndex[dda][chan][ << " "; 01135 } 01136 01137 SampleFile >> dda >> chan >> capArray; 01138 for(sample=0;sample<SAMPLES_PER_BLOCK;sample++) { 01139 SampleFile >> value; 01140 fAtriSampleTimes[dda][chan][capArray][sample]=value; 01141 // std::cout << fAtriSampleTimes[dda][chan][capArray][sample] << " "; 01142 } 01143 // std::cout << "\n"; 01144 } 01145 SampleFile.close(); 01146 01147 char epsilonFileName[100]; 01148 sprintf(epsilonFileName,"%s/ATRI/araAtriStation%iEpsilon.txt",calibDir, stationId); 01149 // fprintf(stdout, "AraEventCalibrator::loadAtriCalib(): INFO - Epsilon file = %s\n", epsilonFileName);//DEBUG 01150 std::ifstream epsilonFile(epsilonFileName); 01151 while(epsilonFile >> dda >> chan >> capArray){ 01152 epsilonFile >> value; 01153 fAtriEpsilonTimes[dda][chan][capArray]=value; 01154 // printf("%s : dda %i channel %i capArray %f\n", __FUNCTION__, dda, chan, value); 01155 } 01156 epsilonFile.close(); 01157 01158 01159 //Now we set the gotCalibFile flags to indicate which station we have in memory 01160 01161 for(AraStationId_t station=0;station<ICRR_NO_STATIONS+ATRI_NO_STATIONS;station++){ 01162 if(AraGeomTool::isAtriStation(station)){ 01163 calibIndex=AraGeomTool::getStationCalibIndex(station); 01164 if(station==stationId) fGotAtriCalibFile[calibIndex]=1; 01165 else fGotAtriCalibFile[calibIndex]=0; 01166 } 01167 } 01168 01169 } 01170
Generated on Mon Mar 18 16:00:44 2013 for ARA ROOT v3.4 Software by
