///// Simple code which parses Hrec, extracts the FrSimEvent and
// also recomputes their approximate SNR
{
  char buffer[256];
  int numInjections = 2;

  FrFile * iFile = FrFileINew("/virgoData/ffl/WSR7/hrec-v1.ffl");
  //FrFile * iFile = FrFileINew("/virgoData/ffl/hrecOnline.ffl");

  double GPSlist[2] = {852683756, 852849240};
  double duration[2] = {10000, 10000};
  double calGainBurst[2] = {1.0, 1.0};
  double calGainCB[2] = {1.0, 1.0};

  double originalDistanceCB  = 1.0;
  double originalDistanceDFM = 0.01;

  double originalSnrCB[2]  = {26.2448, 26.2448};
  double originalHrssCB[2] = {8.4933e-21, 8.4933e-21};

  double originalSnrDFM[2] = {0.529793,0.529793};
  double originalHrssDFM[2] = {3.79892e-22,3.79892e-22};

  double originalSnrGaussian[2] = {2.28999e+21,2.28999e+21};

  double originalSnrSineGaussian1000q5[2] = {6.5425e+21, 6.5425e+21};
  double originalSnrSineGaussian1000q15[2] = {7.01737e+21,7.01737e+21};
  double originalSnrSineGaussian1300q5[2] = {5.91894e+21,5.91894e+21};
  double originalSnrSineGaussian1300q15[2] = {5.90462e+21,5.90462e+21};
  double originalSnrSineGaussian1600q5[2] = {5.03796e+21,5.03796e+21};
  double originalSnrSineGaussian1600q15[2] = {4.73465e+21,4.73465e+21};

  double scaleCB[2] = {0.457234,0.457234};
  double scaleDFM[2] = {28.3129,28.3129};
  double scaleGaussian[2] = {6.55024e-21,6.55024e-21};
  double scaleSineGaussian1000q5[2] = {2.2927e-21, 2.2927e-21};
  double scaleSineGaussian1000q15[2] = {2.13755e-21, 2.13755e-21};
  double scaleSineGaussian1300q5[2] = {2.53454e-21, 2.53454e-21};
  double scaleSineGaussian1300q15[2] = {2.54038e-21, 2.54038e-21};
  double scaleSineGaussian1600q5[2] = {2.9774e-21, 2.9774e-21};
  double scaleSineGaussian1600q15[2] = {3.16813e-21, 3.16813e-21};

  double amplitudeCB = 2.4e-20;
  double amplitudeDFM = 4.85e-21;
  double amplitudeGaussian = 23.7527;
  double amplitudeSineGaussian1000q5 = 37.5;
  double amplitudeSineGaussian1000q15 = 21.5;
  double amplitudeSineGaussian1300q5 = 42.5;
  double amplitudeSineGaussian1300q15 = 24.5;
  double amplitudeSineGaussian1600q5 = 47.0;
  double amplitudeSineGaussian1600q15 = 28.7;

  double gpsInspiralScales[2] = {0.66666666666666, 1.0};
  double inspiralScale = 0.;

  double gpsBurstStart[2] = {852687665,852853596};
  double gpsBurstScales[3] = {0.5, 1.0, 1.666666666};
  double burstScale = 0.;

  char * day[2] = {"13-Jan-2007","14-Jan-2007"};

  FrSimEvent * simEvent;
  FrSimEvent * copyEvent;

  for(int block = 0; block < numInjections; block++) {
    sprintf(buffer,"V-WSR7-Injections-%s",day[block]);
    FrFile * oFile = FrFileONewM(buffer,6,"hrec",duration[block]);
    FrameH * frameH = FrameHNew("InjectionsFrame");
    frameH->dt = duration[block];
    frameH->GTimeS = GPSlist[block];
    frameH->GTimeN = 0;
    frameH->ULeapS = 32;
    simEvent = FrSimEventReadTF(iFile,
				"HI*",
				GPSlist[block],
				duration[block],
				0, 0);
    while(simEvent != NULL) {
      double gtime = ((double) simEvent->GTimeS + 
		      1.0e-9 * ((double) simEvent->GTimeN));
      /* print about the event found */
      FrSimEventDump(simEvent, stderr, 2);
      /* Set a special scale in case it is an inspiral */
      if(strcmp(simEvent->name,"HI_CB") == 0) {
	if(gtime < GPSlist[block] + 1800.) {
	  inspiralScale = gpsInspiralScales[0];
	} else {
	  inspiralScale = gpsInspiralScales[1];
	}
      }
      /* Set a special scale in case it is a burst */
      if(strcmp(simEvent->name,"HI_Burst") == 0) {
	if(gtime < gpsBurstStart[block] + 900.) {
	  burstScale = gpsBurstScales[0];
	} else if(gtime < gpsBurstStart[block] + 1800.) {
	  burstScale = gpsBurstScales[1];
	} else {
	  burstScale = gpsBurstScales[2];
	}
      }
      /* Act differently depending on the event */
      if(strcmp(simEvent->name,"HI_Line") == 0 ||
	 ((strcmp(simEvent->name,"HI_CB") == 0) &&
	  (simEvent->timeBefore > 90.0))) {
	copyEvent = FrSimEventNew(frameH,
				  "HI_Line",
				  "Lines injected",
				  "Sc_NI_LoopIn",
				  gtime + (simEvent->timeAfter - simEvent->timeBefore)/2.,
				  (simEvent->timeBefore + simEvent->timeAfter)/2.,
				  (simEvent->timeBefore + simEvent->timeAfter)/2.,
				  simEvent->amplitude,
				  0);
	  FrSimEventAddParam(copyEvent,"freqLine1",102.);
	  FrSimEventAddParam(copyEvent,"voltLine1",4.e-6);
	  FrSimEventAddParam(copyEvent,"freqLine2",352.);
	  FrSimEventAddParam(copyEvent,"voltLine2",8.e-6);
      } else if(strcmp(simEvent->name,"HI_CB") == 0) {
	copyEvent = FrSimEventNew(frameH,
				  "HI_CB",
				  "NS-NS injected event",
				  "Sc_NI_LoopIn",
				  gtime,
				  (float) simEvent->timeBefore,
				  (float) simEvent->timeAfter,
				  inspiralScale * 
				  amplitudeCB * 
				  scaleCB[block] * 
				  calGainCB[block],
				  0);
	FrSimEventAddParam(copyEvent,"m1",1.39);
	FrSimEventAddParam(copyEvent,"m2",1.47);
	FrSimEventAddParam(copyEvent,"nuLow",50.0);
	FrSimEventAddParam(copyEvent,"distance",
			   originalDistanceCB/
			   (inspiralScale * scaleCB[block] * 
			    calGainCB[block]));
	FrSimEventAddParam(copyEvent,"fAverage",136.395);
	FrSimEventAddParam(copyEvent,"deltaf",163.042);
	FrSimEventAddParam(copyEvent,"deltat",4.03731);
	FrSimEventAddParam(copyEvent,"hrss",
			   inspiralScale *
			   originalHrssCB[block] * 
			   scaleCB[block] * 
			   calGainCB[block]);
	FrSimEventAddParam(copyEvent,"SNR",
			   inspiralScale *
			   originalSnrCB[block] * 
			   scaleCB[block] * 
			   calGainCB[block]);
      } else if(strcmp(simEvent->name,"HI_Burst") == 0 && 
		((simEvent->timeBefore + simEvent->timeAfter) < 0.006)) {
	copyEvent = 
	  FrSimEventNew(frameH,
			"HI_Burst_SineGaussian_f1600_Q5",
			"SineGaussian f1600 Q5 injected event",
			"Sc_NI_LoopIn",
			gtime,
			simEvent->timeBefore,
			simEvent->timeAfter,
			burstScale * 
			amplitudeSineGaussian1600q5 * 
			scaleSineGaussian1600q5[block] *
			calGainBurst[block],
			0);
	FrSimEventAddParam(copyEvent,"Q",5);
	FrSimEventAddParam(copyEvent,"f0",1600);
	FrSimEventAddParam(copyEvent,"hrss",
			   burstScale * 
			   scaleSineGaussian1600q5[block] *
			   calGainBurst[block]);
	FrSimEventAddParam(copyEvent,"sigma",0.0004974);
	FrSimEventAddParam(copyEvent,"width",0.0009947);
	FrSimEventAddParam(copyEvent,"deltat",0.000351686);
	FrSimEventAddParam(copyEvent,"fAverage",1600);
	FrSimEventAddParam(copyEvent,"deltaf",226.274);
	FrSimEventAddParam(copyEvent,"SNR",
			   burstScale *
			   originalSnrSineGaussian1600q5[block] * 
			   scaleSineGaussian1600q5[block] *
			   calGainBurst[block]);
      } else if(strcmp(simEvent->name,"HI_Burst") == 0 && 
		((simEvent->timeBefore + simEvent->timeAfter) < 0.007)) {
	copyEvent = 
	  FrSimEventNew(frameH,
			"HI_Burst_SineGaussian_f1300_Q5",
			"SineGaussian f1300 Q5 injected event",
			"Sc_NI_LoopIn",
			gtime,
			simEvent->timeBefore,
			simEvent->timeAfter,
			burstScale * 
			amplitudeSineGaussian1300q5 * 
			scaleSineGaussian1300q5[block] *
			calGainBurst[block],
			0);
	FrSimEventAddParam(copyEvent,"Q",5);
	FrSimEventAddParam(copyEvent,"f0",1300);
	FrSimEventAddParam(copyEvent,"hrss",
			   burstScale * 
			   scaleSineGaussian1300q5[block] *
			   calGainBurst[block]);
	FrSimEventAddParam(copyEvent,"sigma",0.0006121);
	FrSimEventAddParam(copyEvent,"width",0.001224);
	FrSimEventAddParam(copyEvent,"deltat",0.000432844);
	FrSimEventAddParam(copyEvent,"fAverage",1300);
	FrSimEventAddParam(copyEvent,"deltaf",183.848);
	FrSimEventAddParam(copyEvent,"SNR",
			   burstScale *
			   originalSnrSineGaussian1300q5[block] * 
			   scaleSineGaussian1300q5[block] *
			   calGainBurst[block]);
      } else if(strcmp(simEvent->name,"HI_Burst") == 0 && 
		((simEvent->timeBefore + simEvent->timeAfter) < 0.009)) {
	copyEvent = 
	  FrSimEventNew(frameH,
			"HI_Burst_SineGaussian_f1000_Q5",
			"SineGaussian f1000 Q5 injected event",
			"Sc_NI_LoopIn",
			gtime,
			simEvent->timeBefore,
			simEvent->timeAfter,
			burstScale * 
			amplitudeSineGaussian1000q5 * 
			scaleSineGaussian1000q5[block] *
			calGainBurst[block],
			0);
	FrSimEventAddParam(copyEvent,"Q",5);
	FrSimEventAddParam(copyEvent,"f0",1000);
	FrSimEventAddParam(copyEvent,"hrss",
			   burstScale * 
			   scaleSineGaussian1000q5[block] *
			   calGainBurst[block]);
	FrSimEventAddParam(copyEvent,"sigma",0.0007958);
	FrSimEventAddParam(copyEvent,"width",0.001592);
	FrSimEventAddParam(copyEvent,"deltat",0.000562698);
	FrSimEventAddParam(copyEvent,"fAverage",1000);
	FrSimEventAddParam(copyEvent,"deltaf",141.421);
	FrSimEventAddParam(copyEvent,"SNR",
			   burstScale *
			   originalSnrSineGaussian1000q5[block] * 
			   scaleSineGaussian1000q5[block] *
			   calGainBurst[block]);
      } else if(strcmp(simEvent->name,"HI_Burst") == 0 && 
		(simEvent->timeBefore + simEvent->timeAfter) < 0.011) {
	copyEvent = 
	  FrSimEventNew(frameH,
			"HI_Burst_Gaussian",
			"Gaussian injected event",
			"Sc_NI_LoopIn",
			gtime,
			simEvent->timeBefore,
			simEvent->timeAfter,
			burstScale * 
			amplitudeGaussian * 
			scaleGaussian[block] * 
			calGainBurst[block],
			0);
	FrSimEventAddParam(copyEvent,"sigma",0.001);
	FrSimEventAddParam(copyEvent,"width",0.002);
	FrSimEventAddParam(copyEvent,"hrss",
			   burstScale * 
			   scaleGaussian[block] * 
			   calGainBurst[block]);
	FrSimEventAddParam(copyEvent,"deltat",0.000707107);
	FrSimEventAddParam(copyEvent,"fAverage",89.79);
	FrSimEventAddParam(copyEvent,"deltaf",67.91);
	FrSimEventAddParam(copyEvent,"SNR",
			   burstScale * 
			   originalSnrGaussian[block] * 
			   scaleGaussian[block] * 
			   calGainBurst[block]);
      } else if(strcmp(simEvent->name,"HI_Burst") == 0 && 
		((simEvent->timeBefore + simEvent->timeAfter) < 0.016)) {
	copyEvent = 
	  FrSimEventNew(frameH,
			"HI_Burst_SineGaussian_f1600_Q15",
			"SineGaussian f1600 Q15 injected event",
			"Sc_NI_LoopIn",
			gtime,
			simEvent->timeBefore,
			simEvent->timeAfter,
			burstScale * 
			amplitudeSineGaussian1600q15 * 
			scaleSineGaussian1600q15[block] *
			calGainBurst[block],
			0);
	FrSimEventAddParam(copyEvent,"Q",15);
	FrSimEventAddParam(copyEvent,"f0",1600);
	FrSimEventAddParam(copyEvent,"hrss",
			   burstScale * 
			   scaleSineGaussian1600q15[block] *
			   calGainBurst[block]);
	FrSimEventAddParam(copyEvent,"sigma",0.001492);
	FrSimEventAddParam(copyEvent,"width",0.002984);
	FrSimEventAddParam(copyEvent,"deltat",0.00105506);
	FrSimEventAddParam(copyEvent,"fAverage",1600);
	FrSimEventAddParam(copyEvent,"deltaf",75.4247);
	FrSimEventAddParam(copyEvent,"SNR",
			   burstScale *
			   originalSnrSineGaussian1600q15[block] * 
			   scaleSineGaussian1600q15[block] *
			   calGainBurst[block]);
      } else if(strcmp(simEvent->name,"HI_Burst") == 0 && 
		((simEvent->timeBefore + simEvent->timeAfter) < 0.019)) {
	copyEvent = 
	  FrSimEventNew(frameH,
			"HI_Burst_SineGaussian_f1300_Q15",
			"SineGaussian f1600 Q15 injected event",
			"Sc_NI_LoopIn",
			gtime,
			simEvent->timeBefore,
			simEvent->timeAfter,
			burstScale * 
			amplitudeSineGaussian1300q15 * 
			scaleSineGaussian1300q15[block] *
			calGainBurst[block],
			0);
	FrSimEventAddParam(copyEvent,"Q",15);
	FrSimEventAddParam(copyEvent,"f0",1300);
	FrSimEventAddParam(copyEvent,"hrss",
			   burstScale * 
			   scaleSineGaussian1300q15[block] *
			   calGainBurst[block]);
	FrSimEventAddParam(copyEvent,"sigma",0.001836);
	FrSimEventAddParam(copyEvent,"width",0.003673);
	FrSimEventAddParam(copyEvent,"deltat",0.00129854);
	FrSimEventAddParam(copyEvent,"fAverage",1300);
	FrSimEventAddParam(copyEvent,"deltaf",61.2826);
	FrSimEventAddParam(copyEvent,"SNR",
			   burstScale *
			   originalSnrSineGaussian1300q15[block] * 
			   scaleSineGaussian1300q15[block] *
			   calGainBurst[block]);
      } else if(strcmp(simEvent->name,"HI_Burst") == 0 && 
		((simEvent->timeBefore + simEvent->timeAfter) < 0.025)) {
	copyEvent = 
	  FrSimEventNew(frameH,
			"HI_Burst_SineGaussian_f1000_Q15",
			"SineGaussian f1000 Q15 injected event",
			"Sc_NI_LoopIn",
			gtime,
			simEvent->timeBefore,
			simEvent->timeAfter,
			burstScale * 
			amplitudeSineGaussian1000q15 * 
			scaleSineGaussian1000q15[block] *
			calGainBurst[block],
			0);
	FrSimEventAddParam(copyEvent,"Q",15);
	FrSimEventAddParam(copyEvent,"f0",1000);
	FrSimEventAddParam(copyEvent,"hrss",
			   burstScale * 
			   scaleSineGaussian1000q15[block] *
			   calGainBurst[block]);
	FrSimEventAddParam(copyEvent,"sigma",0.002387);
	FrSimEventAddParam(copyEvent,"width",0.004775);
	FrSimEventAddParam(copyEvent,"deltat",0.00168809);
	FrSimEventAddParam(copyEvent,"fAverage",1000);
	FrSimEventAddParam(copyEvent,"deltaf",47.1405);
	FrSimEventAddParam(copyEvent,"SNR",
			   burstScale *
			   originalSnrSineGaussian1000q15[block] * 
			   scaleSineGaussian1000q15[block] *
			   calGainBurst[block]);
      } else if(strcmp(simEvent->name,"HI_Burst") == 0 && 
		((simEvent->timeBefore + simEvent->timeAfter) < 0.2)) {
	copyEvent = 
	  FrSimEventNew(frameH,
			"HI_Burst_DFM_a:2_b:4_g:1",
			"DFM a:2 b:4 g:1 injected event",
			"Sc_NI_LoopIn",
			gtime,
			simEvent->timeBefore,
			simEvent->timeAfter,
			burstScale * 
			amplitudeDFM * 
			scaleDFM[block] *
			calGainBurst[block],
			0);
	FrSimEventAddParam(copyEvent,"a",2);
	FrSimEventAddParam(copyEvent,"b",4);
	FrSimEventAddParam(copyEvent,"g",1);
	FrSimEventAddParam(copyEvent,"hrss",
			   burstScale * 
			   originalHrssDFM[block] * 
			   scaleDFM[block] *
			   calGainBurst[block]);
	FrSimEventAddParam(copyEvent,"deltat",0.0255508);
	FrSimEventAddParam(copyEvent,"fAverage",55.8127);
	FrSimEventAddParam(copyEvent,"deltaf",65.5145);
	FrSimEventAddParam(copyEvent,"SNR",
			   burstScale * 
			   originalSnrDFM[block] * 
			   scaleDFM[block] *
			   calGainBurst[block]);
	FrSimEventAddParam(copyEvent,"distance",
			   originalDistanceDFM/
			   (burstScale * 
			    scaleDFM[block] *
			    calGainCB[block]));
      } else {
	copyEvent = 
	  FrSimEventNew(frameH,
			simEvent->name,
			simEvent->comment,
			simEvent->inputs,
			gtime,
			simEvent->timeBefore,
			simEvent->timeAfter,
			simEvent->amplitude,
			0);
      }
      simEvent = simEvent->next;
    }
    FrameWrite(frameH,oFile);
    FrameFree(frameH);
    FrFileOEnd(oFile);
  }
}