{
  gROOT->Reset();
  gSystem.Load("libBlast.so");
  TBLRaw raw("electronics.map");
  TBLGeometry geom;
  geom.GetData("blast.geom");
  TBLScRecon scRecon;

  scRecon.Init(geom, "/home/blast/blast/commis/phase1/blast.sc_cal");//for large TOF

  //  TString script("show_p.C");
  TString datadir(gSystem.Getenv("DATADIR"));
  TString analdir(gSystem.Getenv("ANALDIR"));
  TString data(datadir+"commis-");
  int outFlag;
  char outFile[100], title[125];

  raw.ReadEpicsMap("/home/blast/blast/commis/epics/Epics.map");
  int  beam_index = raw.GetEpicsID("ldcct:ma");

  // read in the data file numbers, put them in list[]
  int arc=gApplication->Argc();
  char** arv=gApplication->Argv();
  int arg=0;

  int n1, n2, list[100], nl, il;
  int ntot, n3;
  int i, k, j, kk;

  outFlag = 0;
  for (i=0;i<arc;i++){
    k = sscanf(arv[i],"%d-%d",&n1,&n2);
    if (k == 1 && n1 > 0) list[nl++] = n1;
    if (k == 2 && n1 > 0 && n2 > 0 && n2>n1) 
      for (j=n1;j<n2+1;j++) list[nl++] = j;
    if (strcmp(arv[i],"-o")==0 && i+1 < arc) {
      strcpy(outFile,arv[i+1]); 
      outFlag = 1;   }
  }

  int iPadLeft, iPadRight;
  float timeLeft, timeRight, rtOR;
  float adcLeft,  adcRight, posLeft, posRight;

  TH1F *tdiff_L[6][6];
  TH1F *tdiff_R[6][6];

//  TH1F *copla_L[6][6];
//  TH1F *copla_R[6][6];

  TH2F *copla_2d_L[6][6];
  TH2F *copla_2d_R[6][6];

  for (k=10;k<16;k++) {
      kk = k - 10;
    for (j=0;j<6;j++) {
      sprintf(title,"tdiff_L[%d]_R[%d]",k,j);
      TH1F *tdiff_L[kk][j]  = new TH1F(title,title,400,-800,800);
      sprintf(title,"tdiff_R[%d]_L[%d]",k,j);
      TH1F *tdiff_R[kk][j]  = new TH1F(title,title,400,-800,800);

//      sprintf(title,"copla_L[%d]_R[%d]",k,j);
//      TH2F *copla_L[kk][j]  = new TH2F(title,title,150,-300,300);
//      sprintf(title,"copla_R[%d]_L[%d]",k,j);
//      TH2F *copla_R[kk][j]  = new TH2F(title,title,150,-300,300);

      sprintf(title,"copla_L[%d]_R[%d]",k,j);
      TH2F *copla_2d_L[kk][j]  = new TH2F(title,title,50,-400,400,50,-400,400);
      sprintf(title,"copla_R[%d]_L[%d]",k,j);
      TH2F *copla_2d_R[kk][j]  = new TH2F(title,title,50,-400,400,50,-400,400);

    }
  }

  char tmp[50];

  
  for (il=0;il<nl;il++){
    sprintf(tmp,"%d",list[il]);
    TString run(data+tmp+".dat");
    cout<<"Processing "<<run<<" ..."<<endl;
    raw.Open(run.Data()); 
    raw.ReadNext();
    raw.goTime = 0; 


    // global statistics
    int phys, clean, good_Lpad, good_Rpad, passed_Lcut, passed_Rcut, total_passed;
    phys = 0; clean = 0; 
    passed_Lcut = passed_Rcut = total_passed = 0 ;
    
    int nCC, nNC, good_LCC, good_RCC, beam;


    for (int ctr=0; !raw.ReadNext() && ctr < 1e9 ; ctr++) {
      if (ctr%50000==0) printf("\n%9d ", ctr);
      else if (ctr%10000==0) cout<<":"<<flush;
      else if (ctr%1000==0) cout<<"."<<flush;


      if (raw.GetEventType() == 7) beam = raw.epics[beam_index].value;


      if (raw.GetEventType() == 1) {

      // Physics Events - now go select good coincs (event by event)
      good_Lpad = 0; good_Rpad = 0;
      good_LCC = 0; good_RCC = 0;

	phys++;
	scRecon << raw;
	scRecon.PaddleFit();
	iPadLeft =  iPadRight = -1;


 	for (j = 0;j<15;j++)
	  if (scRecon.Data[0][j].quality > 0 ) {
	      good_Lpad++;
	      iPadLeft = j; }
	//  printf("\n found good LEFT paddle %9d %9d",j,scRecon.Data[0][j]);

	for (j = 0;j<15;j++)
	  if (scRecon.Data[1][j].quality > 0 ) {
	      good_Rpad++;
	      iPadRight = j; }
	// printf("\n found good RIGHT paddle %9d %9d",j,scRecon.Data[1][j]);


	 if ( good_Rpad == 1 && good_Lpad == 1 ) {
	  clean++;

	  timeLeft  = (scRecon.Data[0][iPadLeft].tdc1  
		       + scRecon.Data[0][iPadLeft].tdc2)/2.0;
	  timeRight = (scRecon.Data[1][iPadRight].tdc1 
		       + scRecon.Data[1][iPadRight].tdc2)/2.0;
	  posLeft  = (scRecon.Data[0][iPadLeft].tdc1  
		       - scRecon.Data[0][iPadLeft].tdc2)/2.0;
	  posRight = (scRecon.Data[1][iPadRight].tdc1 
		       - scRecon.Data[1][iPadRight].tdc2)/2.0;
	  adcLeft  = (scRecon.Data[0][iPadLeft].adc1  
		      + scRecon.Data[0][iPadLeft].adc2);
	  adcRight = (scRecon.Data[1][iPadRight].adc1 
		      + scRecon.Data[1][iPadRight].adc2);
	  
	  int sect, pad, lr ;
	  nNC = raw.GetNCnhits();
	     for (i=0;i<nNC;i++){
              sect = raw.GetNCsector(i);
	      pad  = raw.GetNCpaddle(i);
	      lr   = raw.GetNClr(i);
	      if (sect == 1 && pad == 0 && lr == 0) rtOR = raw.GetNCtdc(i);
	     }


	  nCC = raw.GetCCnhits();
	  //	  printf("\n %9d ----------------------------------------------",nCC);
	     for (i=0;i<nCC;i++){
	  //	  printf("\n %3d %3d %6d",raw.GetCCsector(i),raw.GetCCbox(i),raw.GetCCtdc(i) );
	       if (raw.GetCCsector(i) == 0){
		 if ( iPadLeft < 4 && raw.GetCCbox(i) == 0 && raw.GetCCtdc(i) > -1) good_LCC++ ;
		 else if ( iPadLeft < 6 && raw.GetCCbox(i) == 1 && raw.GetCCtdc(i) > -1) good_LCC++ ;

	       }
	       if (raw.GetCCsector(i) == 1){
		 if ( iPadRight < 4 && raw.GetCCbox(i) == 0 && raw.GetCCtdc(i) > -1) good_RCC++ ;
		 else if ( iPadRight < 6 && raw.GetCCbox(i) == 1 && raw.GetCCtdc(i) > -1) good_RCC++ ;
	       }
             }




	  if (iPadLeft > 9 && iPadRight < 6 
	                                      && beam < 20         // cut on beam current
	      //                              && good_RCC == 1     // cut R-forward CC hit
	      //	      	              && adcRight < 2000   // Electron ADC  (R)     
	      //	      	              && adcLeft > 3000    // Hadron ADC (L)  
	      //                              && rtOR < 1670 && rtOR > 1660 // retiming-OR self timing!
         ) {
	    k = iPadLeft - 10;
	    j = iPadRight;
            passed_Lcut++;
	    tdiff_L[k][j]->Fill(timeLeft-timeRight); 
	    //      copla_L[k][j]->Fill(posLeft + posRight);
	    	    copla_2d_L[k][j]->Fill(posLeft , posRight); 
           }


	  else if (iPadLeft < 6 && iPadRight > 9 
	                                      && beam < 20         // cut on beam current
	      //                              && good_LCC == 1     // cut L-forward CC hit 
	      //		   	      && adcLeft < 2000    // Electron ADC (L)      
	      //		   	      && adcRight > 3000   // Hadron ADC (R)
	      //                              && rtOR < 1670 && rtOR > 1660 // retiming-OR self timing!
         ) {
	    k = iPadRight - 10;
	    j = iPadLeft;
            passed_Rcut++;
	    tdiff_R[k][j]->Fill(timeRight-timeLeft);
	    //	    copla_R[k][j]->Fill(posLeft + posRight);
	    	    copla_2d_R[k][j]->Fill(posLeft , posRight); 
	   }


    	  }  // end of coinc was good (only 1 L and 1 R valid hit)

        } // end of physics event

      }	// loop on events

    /*
    total_passed = passed_Lcut+passed_Rcut;
    if (phys >0 ) {
	printf("\n Total physics events         %9d",phys);
	printf("\n Clean events (1Lhit + 1Rhit) %9d  (%4d \%)",clean,100.*clean/phys);      
	printf("\n Clean that passed proton cut %9d  (%4d \%)",total_passed,100.*total_passed/phys);     
    if (clean >0 ) {
	printf("\n       ... of which in left  sector %9d  (%4d \%)",passed_Lcut,100.*passed_Lcut/total_passed);      
	printf("\n       ... of which in right sector %9d  (%4d \%)",passed_Rcut,100.*passed_Rcut/total_passed);      
	printf("\n");
	printf("\n");
    } }
    */
          raw.CodaRead.closedatafile();

  }	// end of loop on runs - reset phys and clean counters



    TCanvas *c1 = new TCanvas("c1","L-R coincs", 10, 10, 1200, 900);
    c1->Divide(6,5); // columns, rows
    TCanvas *c2 = new TCanvas("c2","R-L coincs", 10, 10, 1200, 900);
    c2->Divide(6,5); 

//  TCanvas *c1 = new TCanvas("c1","coincs", 10, 10, 1200, 900);
//  c1->Divide(3,3); // columns, rows
//  TCanvas *c2 = new TCanvas("c2","coincs", 10, 10, 1200, 900);
//  c2->Divide(3,3); // columns, rows

//    TCanvas *c3 = new TCanvas("c3","L-R coplanarity", 10, 10, 1200, 900);
//    c3->Divide(6,5); 
//    TCanvas *c4 = new TCanvas("c4","R-L coplanarity", 10, 10, 1200, 900);
//    c4->Divide(6,5); 

    for (k=0; k<5; k++) { // the backward paddle (rows)
      for (j=0; j<6; j++) { // the forward paddle (columns)

//   for (k=1; k<4; k++) { // the backward paddle (rows)
//      for (j=3; j<6; j++) { // the forward paddle (columns)

  c1->cd(6*k+j+1);
  tdiff_L[k][j]->Draw();
  c2->cd(6*k+j+1);
  tdiff_R[k][j]->Draw();

	//c3->cd(6*k+j+1);
        //copla_2d_L[k][j]->Draw();
        //c4->cd(6*k+j+1);
        //copla_2d_R[k][j]->Draw(); 

//  c1->cd(j-3+1+(k-1)*3);
//        copla_2d_L[k][j]->Draw();
//  c2->cd(j-3+1+(k-1)*3);
//        copla_2d_R[k][j]->Draw();

      } 
    }


}