/* ======================================================================== *\
!   Take different MC input files and combine them into single one
!   Author(s): Oscar Blanch & Javier Rico
\* ======================================================================== */

void mergeClean()
{
  const TString fileout = "/local_disk/jrico/mc/prueba.root";

  const Int_t nfiles = 11;
  const Int_t nlevels= 5;

  const Char_t* levelname[nlevels]={"noclean","3015","4025","5035","6045"};

  // input files -- no noise
  const Char_t* NonoiseOutFilename[nfiles] = 
    {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_nonoise.root",
     "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_nonoise.root"
    };

  // input files -- different cleaning levels
  const Char_t* NoiseOutFilename[nlevels][nfiles]   = 
    {{"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_noclean.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_noclean.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_noclean.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_noclean.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_noclean.root",
      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_noclean.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_noclean.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_noclean.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_noclean.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_noclean.root",
      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_noclean.root"},
     {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_3015.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_3015.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_3015.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_3015.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_3015.root",
      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_3015.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_3015.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_3015.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_3015.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_3015.root",
      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_3015.root"},
     {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_4025.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_4025.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_4025.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_4025.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_4025.root",
      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_4025.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_4025.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_4025.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_4025.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_4025.root",
      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_4025.root"},   
     {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_5035.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_5035.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_5035.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_5035.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_5035.root",
      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_5035.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_5035.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_5035.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_5035.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_5035.root",
      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_5035.root"},
     {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_6045.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_6045.root",   
      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_6045.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_6045.root",
      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_6045.root",
      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_6045.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_6045.root",
      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_6045.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_6045.root",
      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_6045.root",
      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_6045.root"}
    };   


  // Getting input braches 
  TChain* tno  = new TChain("Events");
  TChain* tyes[nlevels];
  for(Int_t i=0;i<nlevels;i++)
    tyes[i] = new TChain("Events");

  for(Int_t i=0;i<nfiles;i++)
    {
      tno->Add(NonoiseOutFilename[i]);
      for(Int_t j=0;j<nlevels;j++)
	tyes[j]->Add(NoiseOutFilename[j][i]);
    }

  // input braches for no noise file
  MCerPhotEvt*   nphotNo;
  MHillas*       hillasNo;
  MHillasSrc*    hillassrcNo;
  MMcEvt*        mcevtNo;
  MRawEvtHeader* evtheaderNo;
  MRawRunHeader* runheaderNo;

  tno->SetBranchAddress("MCerPhotEvt.",&nphotNo);
  tno->SetBranchAddress("MHillas.",&hillasNo);
  tno->SetBranchAddress("MHillasSrc.",&hillassrcNo);
  tno->SetBranchAddress("MMcEvt.",&mcevtNo);
  tno->SetBranchAddress("MRawEvtHeader.",&evtheaderNo);
  tno->SetBranchAddress("MRawRunHeader.",&runheaderNo);

  
  // input branches for noise files
  MCerPhotEvt*   nphotYes[nlevels];
  MHillas*       hillasYes[nlevels];
  MHillasSrc*    hillassrcYes[nlevels];
  MMcEvt*        mcevtYes[nlevels];
  MRawEvtHeader* evtheaderYes[nlevels];
  MRawRunHeader* runheaderYes[nlevels];

  for(Int_t i=0;i<nlevels;i++)
    {      
      tyes[i]->SetBranchAddress("MCerPhotEvt.",&nphotYes[i]);
      tyes[i]->SetBranchAddress("MHillas.",&hillasYes[i]);
      tyes[i]->SetBranchAddress("MHillasSrc.",&hillassrcYes[i]);
      tyes[i]->SetBranchAddress("MMcEvt.",&mcevtYes[i]);
      tyes[i]->SetBranchAddress("MRawEvtHeader.",&evtheaderYes[i]);
      tyes[i]->SetBranchAddress("MRawRunHeader.",&runheaderYes[i]);
    }

  // Setting ouput branches
  TFile* fout = new TFile(fileout,"RECREATE");
  TTree* to   = new TTree("Events","Events");   

  // common containers
  TBranch* bmcevt        = to->Branch("MMcEvt.",       "MMcEvt",       &mcevtNo);
  TBranch* brawevtheader = to->Branch("MRawEvtHeader.","MRawEvtHeader",&evtheaderNo);
  TBranch* brawrunheader = to->Branch("MRawRunHeader.","MRawRunHeader",&runheaderNo);

  // noise/cleaning dependent containers
  to->Branch("MCerPhotEvt0.",  "MCerPhotEvt",  &nphotNo);
  to->Branch("MHillas0.",      "MHillas",      &hillasNo);
  to->Branch("MHillasSrc0.",   "MHillasSrc",   &hillassrcNo);

  //  cout << "llego1" << endl;
  for(Int_t i=0;i<nlevels;i++)
    {
      //      cout << "i: " << i << endl;
      
      char mcerphot[100];
      char mhillas[100];
      char mhillassrc[100];


      sprintf(mcerphot,  "MCerPhotEvt%s.",levelname[i]);
      sprintf(mhillas,   "MHillas%s."    ,levelname[i]);
      sprintf(mhillassrc,"MHillasSrc%s." ,levelname[i]);

      to->Branch(mcerphot,  "MCerPhotEvt",  &nphotYes[i]);
      to->Branch(mhillas,   "MHillas",      &hillasYes[i]);
      to->Branch(mhillassrc,"MHillasSrc",   &hillassrcYes[i]);      
    }
  //  cout << "llego2" << endl;
  
      
  Int_t nentries1 = (Int_t)tno->GetEntries();
  cout << "Entries: No noise file: " << nentries1 << endl;
  Int_t nentries2[nlevels];
  Int_t latestk[nlevels];
  Int_t latestevt[nlevels];
  Int_t nfound[nlevels];
  for(Int_t i=0;i<nlevels;i++)
    {
      nfound[i]=0;
      latestk[i]=0;
      latestevt[i]=99999999;
      nentries2[i] = (Int_t)tyes[i]->GetEntries();
      cout << "Entries: cleaning level [" << levelname[i] << "] :" << nentries2[i] << endl;
    }

  // loop over all events in no-noise file and look for them in the other files
  for(Int_t i=0;i<nentries1;i++)
    //  for(Int_t i=0;i<100;i++)
    {      
      tno->GetEntry(i);
      Int_t runno = runheaderNo->GetRunNumber();
      Int_t evtno = evtheaderNo->GetDAQEvtNumber();

      cout << "Checking run " << runno<< ", event " << evtno << endl;

      //check if the no-noise event is present in all of the other files
      for(Int_t j=0;j<nlevels;j++)
	{
	  Int_t runyes;
	  Int_t evtyes;
	  cout << "Level " << levelname[j] << endl;
	  Int_t k=latestk[j];	  
	  while(k<nentries2[j])
	    {
	      //	      cout << k << endl;
	      tyes[j]->GetEntry(k++);
	      runyes = runheaderYes[j]->GetRunNumber();
	      evtyes = evtheaderYes[j]->GetDAQEvtNumber();
	      if(runyes==runno && evtyes==evtno)
		break;
	      // next condition speeds up selection when events are grouped by runs and
	      // in event increasing order
	      if(runyes==runno && latestevt[j]<evtno && evtno<evtyes)
		break;
	      // next condition speeds up selection when events are grouped by runs
	      // and in run increasing order
	      //	      if(runyes>runno)
	      //		break;
	    }

	  if(k>=nentries2[j])
	    {
	      k=0;
	      while(k<latestk[j])
		{
		  //		  cout << k << endl;
		  tyes[j]->GetEntry(k++);
		  runyes = runheaderYes[j]->GetRunNumber();
		  evtyes = evtheaderYes[j]->GetDAQEvtNumber();
		  if(runyes==runno && evtyes==evtno)
		    break;
		  // next condition speeds up selection when events are grouped by runs and
		  // in event increasing order
		  if(runyes==runno && latestevt[j]<evtno && evtno<evtyes)
		    break;
		  // next condition speeds up selection when events are grouped by runs
		  // and in run increasing order
		  //		  if(runyes>runno)
		  //		    break;
		}
	    }
	  
	  // the event has not been found, assign dummy containers
	  if(runyes=!runno || evtyes!=evtno)
	    {
	      cout << "Not Found!" << endl;	      
	      for(Int_t l=j;l<nlevels;l++)
		{
		  nphotYes[l]->Reset();
		  hillasYes[l]->Reset();
		  hillassrcYes[l]->Reset();		  
		}
	      break;
	    }
	  else
	    {
	      cout << "Found!" << endl;
	      latestevt[j] = evtyes;
	      latestk[j]=k;	  
	      nfound[j]++;
	    }
	}
      to->Fill();
    }

  // loop over all events in first noise file and look for them in the no-noise file in order to save also purely noisy events  
  cout << "***************************************" << endl;
  cout << "*** Looking for purely noisy events ***" << endl;
  cout << "***************************************" << endl;
  Int_t nfound2[nlevels];
  for(Int_t i=0;i<nlevels;i++)
    {
      nfound2[i]=0;
      latestk[i]=0;
      latestevt[i]=99999999;
    }
  
  // readdress the common container addresses
  bmcevt->SetAddress(&mcevtYes[0]);
  brawevtheader->SetAddress(&evtheaderYes[0]);
  brawrunheader->SetAddress(&runheaderYes[0]);

  for(Int_t i=0;i<nentries2[0];i++)
  //  for(Int_t i=0;i<100;i++)
    {      
      tyes[0]->GetEntry(i);
      Int_t run = runheaderYes[0]->GetRunNumber();
      Int_t evt = evtheaderYes[0]->GetDAQEvtNumber();
      Int_t runno;
      Int_t evtno;

      cout << "Checking run " << run << ", event " << evt << endl;

      //check if the event is present in the no-noise file
      Int_t kk=latestk[0];	  
      while(kk<nentries1)
	{
	  //	  cout << kk << endl;
	  tno->GetEntry(kk++);
	  runno = runheaderNo->GetRunNumber();
	  evtno = evtheaderNo->GetDAQEvtNumber();
	  if(run==runno && evt==evtno)
	    break;
	  if(run==runno && latestevt[0]<evt && evt<evtno)
	    break;
	  //	  if(runno>run)
	  //	    break;
	}

      if(kk>=nentries1)
	{
	  kk=0;
	  while(kk<latestk[0])
	    {
	      //		  cout << kk << endl;
	      tno->GetEntry(kk++);
	      runno = runheaderNo->GetRunNumber();
	      evtno = evtheaderNo->GetDAQEvtNumber();
	      if(run==runno && evt==evtno)
		break;
	      if(run==runno && latestevt[0]<evt && evt<evtno)
		break;
	      //	      if(runno>run)
	      //		break;
	    }
	}

      // the event is already included because it is also a no-noise event, dismiss it
      if(run==runno && evt==evtno)
	{
	  cout << "Found already in the no-noise sample, continuing... " << endl;
	  latestevt[0] = evtno;
	  latestk[0]=kk;	  
	  continue;
	}

      nfound2[0]++;
      nphotNo->Reset();
      hillasNo->Reset();
      hillassrcNo->Reset();

      // look for the event in the rest of the noisy samples
      for(Int_t j=1;j<nlevels;j++)
	{
	  Int_t runyes;
	  Int_t evtyes;
	  cout << "Level " << levelname[j] << endl;
	  Int_t k=latestk[j];	  
	  while(k<nentries2[j])
	    {
	      //	      cout << k << endl;
	      tyes[j]->GetEntry(k++);
	      runyes = runheaderYes[j]->GetRunNumber();
	      evtyes = evtheaderYes[j]->GetDAQEvtNumber();
	      if(runyes==run && evtyes==evt)
		break;
	      if(runyes==run && latestevt[j]<evt && evt<evtyes)
		break;
	      //	      if(runyes>run)
	      //		break;
	    }

	  if(k>=nentries2[j])
	    {
	      k=0;
	      while(k<latestk[j])
		{
		  //		  cout << k << endl;
		  tyes[j]->GetEntry(k++);
		  runyes = runheaderYes[j]->GetRunNumber();
		  evtyes = evtheaderYes[j]->GetDAQEvtNumber();
		  if(runyes==run && evtyes==evt)
		    break;
		  if(runyes==run && latestevt[j]<evt && evt<evtyes)
		    break;
		  //		  if(runyes>run)
		  //		    break;
		}
	    }
	  
	  // the event has not been found, assign dummy containers
	  if(runyes=!run || evtyes!=evt)
	    {
	      cout << "Not Found!" << endl;	      
	      for(Int_t l=j;l<nlevels;l++)
		{
		  nphotYes[l]->Reset();
		  hillasYes[l]->Reset();
		  hillassrcYes[l]->Reset();		  
		}
	      break;
	    }
	  else
	    {
	      cout << "Found!" << endl;
	      latestevt[j] = evtyes;
	      latestk[j]=k;	  
	      nfound2[j]++;
	    }
	}
      to->Fill();    
    }

  fout->Write();
  fout->Close();

  for(Int_t i=0;i<nlevels;i++)
    cout << "Found " << nfound[i] << "+" << nfound2[i] << " for cleaning level " << levelname[i] << endl;

}