ARA ROOT v3.13 Software

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 doxygen 1.4.7