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