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

Generated on Mon Mar 18 16:00:44 2013 for ARA ROOT v3.4 Software by doxygen 1.4.7