// This file deals with the outputs of both the LCMS and Cluster monitors.
// Code to re-order the Events and Hits trees for each detector, in the correct time order. (In the original file, the 5 time bins associated with each BAT trigger are not always consecutive, though each block of 5 is reliably distinct). The script automatically adds a variable "OrderTrue" to each tree, which labels each event by the time bin it falls into (0,1,2,3,4).   
// This version also adds hitId to the Event tree. For each event, the hitID tells us the number of the first corresponding entry in the Hit tree. We also have a BadStrip flag in the Hit tree to identify the noisy strips. This done simply by having an array containing the positions of all the bad strips, and using this to tag events.
// Set up for source tests (simple, only one sensor)
// This version also adds two variables related to the timing of the MIP tests. WindowDelay is the delay (relative to the input clock) of the 5ns window. SampleTime is the time each sample was taken, determined by the OrderTrue and the WindowDelay. WindowDelay is set as an argument of the function.

// THIS SCRIPT SHOULD BE COMPILED BEFORE USE

#include <iostream>
#include "TTree.h"
#include "TFile.h"

//

void sort(char rootFile[], int WindowDelay) {

  cout << "====  LHCb input file:  " << rootFile << endl;

        
        /// Code to sort our trees so that the 5 samples taken after each hit on the BAT is sorted into the right order.
        /// Want to modify so that it runs on all the trees in a file, which we give by name.

	TFile f(rootFile);

	//open new file to store the output
        
	Char_t outFile[100];
	sprintf(outFile, "Sort_%s", rootFile); // creates an output filename

	TFile f2(outFile,"recreate");

	cout << "====  Output file:  " << outFile << endl;

	// Set up the information about the sensors

	int num_sensors=1;

	// Also, the no of samples for each trigger, and the no of the initial samples used to find the pedestals (default 1100)

	int num_samples=5;
	int num_skip=1100;

	TString sensorName[num_sensors];

	// Specify the names of the sensors here

	sensorName[0] = "";

	// Also, set up the bad strip array. This is initialised to all zeroes.
	
	int sensor=0; // initialise variable to keep track of sensor we're working on
	int strip=0; // keeps track of individual noisy strip entries

	int num_noisystrips=7; // max no. of noisy strips on a sensor - choose this yourself

	int noisyStrips[num_sensors][num_noisystrips];

	// Loop to initialise entries in array to zero

        for(sensor = 0; sensor < num_sensors ; sensor++){
	  for(strip = 0; strip < num_noisystrips ; strip++){
	    noisyStrips[sensor][strip]=0;
	  }
	}

	// Then, we can specify specific noisy strips. Need to only choose the relevant sensor!
	/*
	// Gla2 - no proper noisy strips, but due to an oddity in the pattern of unbonded strips I want to ignore strip 3.
	noisyStrips[0][0]=3; 
	noisyStrips[0][1]=127; 
	*/

	// GLA1 - SINGLE MIP TESTING

	// Gla1 - avoid a single dubious strip near the start (bond pattern) and also last strip in each Alink (possibly shares signal with next strip)
	noisyStrips[0][0]=387; 
	noisyStrips[0][1]=511; 
	// Last strip in each Alink is dodgy too
	noisyStrips[0][2]=415; 
	noisyStrips[0][3]=447; 
	noisyStrips[0][4]=479; 

	// On 3D sensor, finally want to avoid a strip with no proper signal, plus the adjacent one
	noisyStrips[0][5]=450; 
	noisyStrips[0][6]=451; 

	
	/*
	// Gla0
	noisyStrips[0][0]=259; 
	noisyStrips[0][1]=383; 
	 

	// Fre2
	// The strips are around 1173 and 205, which are no use after re-ordering (fall between header strips) and also 1153, which is noisy. Note that 1152 is not used for finding hits anyway (edge of array).
	// So, block 1172-74, 1204-06, 1153
	noisyStrips[0][0]=1172;  
	noisyStrips[0][1]=1173;
	noisyStrips[0][2]=1174;
	noisyStrips[0][3]=1204;
	noisyStrips[0][4]=1205;
	noisyStrips[0][5]=1206;
 	noisyStrips[0][6]=1153;
	
	
	// On Fre0, we have a few noisy strips - 1467, 1469, 1471. Also have strips falling between header strips as above - 1416, 1448. Lastly, 1460 is unbonded.
	// So, block 1415-7, 1447-49, 1459-61, 1467-71
	noisyStrips[0][0]=1415;  
	noisyStrips[0][1]=1416; 
	noisyStrips[0][2]=1417;
	noisyStrips[0][3]=1447; 
	noisyStrips[0][4]=1448;
	noisyStrips[0][5]=1449;
	noisyStrips[0][6]=1459;
	noisyStrips[0][7]=1460;
	noisyStrips[0][8]=1461;
	noisyStrips[0][9]=1467;
	noisyStrips[0][10]=1468;
	noisyStrips[0][11]=1469;
	noisyStrips[0][12]=1470;
	noisyStrips[0][13]=1471;
	*/





	// Define variables used when adding new data to output file
	
	int OrderTrue = 0;

	int SampleTime = 0;

	int hitId=0;
	int BadStrip=0;


	int iHit=0;
	int iSens=0;
	int iSensOffset=0;

	Int_t nentriesEvent;
	Int_t nentriesHit;
	  
	Int_t *indexEvent;
	Int_t *indexHit;
	
	int BCntLoopEvent=0;
	int BCntOldEvent=0;  // Used to check for consecutive values
	int BATTrigEvent=0;
	int nHitsEvent=0;
	int BCntLoopHit=0; 
	int PosXHit=0;

        int checker=0;

	// The trees:

	TTree *treeEvent;
	TTree *treeHit;
	TTree *tsortedEvent;
	TTree *tsortedHit;

	// Strings I'll use
	TString pathEvent;
	TString pathHit;
	

	// Now, need to do the following for each sensor. Will copy and paste code. Usual syntax applies
	// During debugging, work with Fre2, then fix others after. 


   // UBER-LOOP will deal with the separate LCMS and Cluster trees. LCMS is case 0, Cluster is case 1.
   for(checker = 0; checker < 2 ; checker++)
     {

     //  First loop will loop over all the sensors, saving us from having to use new code for each sensor. sensor keeps track.
     for(sensor = 0; sensor < num_sensors ; sensor++)

       {

	//*******************************
	// Setup for each sensor
	//*******************************

	// Getting the relevant tree within the file. Here, I want to get the Event AND Hit trees. Need to use 
	// Need to use different versions for LCMS and ClusterMoni checkers
	
	if(checker==0){    
          pathEvent = "VeloLCMSMoni/LCMSEvent" + sensorName[sensor]; // Put together the path	
	  pathHit = "VeloLCMSMoni/LCMSHit" + sensorName[sensor];
	}
        else{
          pathEvent = "VeloClusterMoni/ClusterEvent" + sensorName[sensor]; // Put together the path	
	  pathHit = "VeloClusterMoni/ClusterHit" + sensorName[sensor];
	}

	treeEvent = (TTree*)f.Get(pathEvent);
	nentriesEvent = (Int_t)treeEvent->GetEntries();
	treeHit = (TTree*)f.Get(pathHit);
	nentriesHit = (Int_t)treeHit->GetEntries();


	//Drawing variable BCntLoop with no graphics option, for the event tree. 
	//Array V1 (see TTree::Draw) is then an array of these BCntLoop entries
	treeEvent->Draw("BCntLoop","","goff");

	indexEvent = new Int_t[nentriesEvent];

	//Sorts the array containing BCntLoop values (V1) and puts the sorted indices into indexEvent
	//"false" options sets increasing order sorting
	TMath::Sort(nentriesEvent,treeEvent->GetV1(),indexEvent,false);
	
	//Create an empty clone of the original tree
	//I want to be able to add extra entries to this cloned tree
	// Copy the existing tree structure, call the tree *tsorted
	
	tsortedEvent = (TTree*)treeEvent->CloneTree(0);

	// Add new branches for the true ordering of events in each batch of 5 (OrderTrue), the hitId (to make sorting the hit tree easier) and the WindowDelay and SampleTime variables
	// Note that I'm trying to add the same info to a couple of different trees here. Have to be careful.

        tsortedEvent->Branch("OrderTrue",&OrderTrue,"OrderTrue/I");
        tsortedEvent->Branch("hitId",&hitId,"hitId/I");
        tsortedEvent->Branch("WindowDelay",&WindowDelay,"WindowDelay/I");
        tsortedEvent->Branch("SampleTime",&SampleTime,"SampleTime/I");


	//****** 
	// Next, we need to do something similar for the Hit tree.

	//Drawing variable BCntLoop with no graphics option, for the hit tree. 
	//Array V1 (see TTree::Draw) is then an array of these BCntLoop entries
	treeHit->Draw("BCntLoop","","goff");
	//Then, we create an array, with the same no. of entries as the Event tree
	indexHit = new Int_t[nentriesHit];

	//Sorts the array containing BCntLoop values (V1) and puts the sorted indices into indexHit
	//"false" options sets increasing order sorting
	TMath::Sort(nentriesHit,treeHit->GetV1(),indexHit,false);
	
	//Create an empty clone of the original tree
	//I want to be able to add extra entries to this cloned tree
	// Copy the existing tree structure, call the tree *tsorted
	
	tsortedHit = (TTree*)treeHit->CloneTree(0);

	// Add new branches for the true ordering of events in each batch of 5 (OrderTrue), the raw timing info (TimeDiff), the calculated hit time in ns HitTime, and the noisy strip flag BadStrip

        tsortedHit->Branch("OrderTrue",&OrderTrue,"OrderTrue/I");
        tsortedHit->Branch("BadStrip",&BadStrip,"BadStrip/I");
        tsortedHit->Branch("WindowDelay",&WindowDelay,"WindowDelay/I");
        tsortedHit->Branch("SampleTime",&SampleTime,"SampleTime/I");

	//***



	// Have to specify the following so that we can look at specific variables after using GetEntry. Also have to be careful when dealing with different detector trees - use different variables for different detectors. I get the BCntLoop to check that there are no errors in the file.
	BCntLoopEvent=0;
	BCntOldEvent=0;  // Used to check for consecutive values
	BATTrigEvent=0;

	treeEvent->SetBranchAddress("BCntLoop",&BCntLoopEvent);
	treeEvent->SetBranchAddress("BATTrig",&BATTrigEvent);
	treeEvent->SetBranchAddress("nHits",&nHitsEvent);

	BCntLoopHit=0;
	PosXHit=0;

	treeHit->SetBranchAddress("BCntLoop",&BCntLoopHit);

        // Slight difference in data from trees - cluster tree gives positions of both start of cluster and its main strip. Use the main strip value.
	if(checker==0) treeHit->SetBranchAddress("PosX",&PosXHit);
        else treeHit->SetBranchAddress("CentrePosX",&PosXHit);

	//***************
	//
	// LOOPING OVER ALL THE ENTRIES IN THE EVENT TREE
	// Need to use iSens to loop over the sensors.
	// Within this loop, we use the array of indices to grab the events in the correct order, writing them to a new tree.


	iHit=0;  // iterator to keep track of our position within the Hit tree

	for (iSens=iSensOffset;iSens<nentriesEvent;iSens++) {
	  
	  treeEvent->GetEntry(indexEvent[iSens]); // Read in an event. Previous statements mean BCntLoop is accessible.

	  // Check to make sure that the trigger numbering is working correctly by comparing the BATTrig in the Vetra file, eventNum in the TDC, and the iTDC variable being used here
	  if (BATTrigEvent!=((iSens+num_skip)/num_samples)){
	    cout << "Event numbering problem at BATTrig=" << BATTrigEvent << " entry number" << iSens << endl;
	  }

	  
	  // Apply the numbering to the "num_samples" time bins

          OrderTrue = (iSens-iSensOffset)%num_samples;

	  // Calculate the SampleTime, based on the correct ordering and the delay of the acceptance window

          SampleTime=25*OrderTrue+(20-WindowDelay);

	  // Check the bunch counters are consecutive (if OrderTrue>0)
	  if ((OrderTrue>0)&&((BCntLoopEvent-BCntOldEvent) != 1)){
	    cout << "Non-consecutive entries at BCntLoopEvent=" << BCntLoopEvent << " BAT Trigger no=" << BATTrigEvent << " OrderTrue=" << OrderTrue << endl;
	  }

	  //		cout << "Value of BCntLoopEvent is  " << BCntLoopEvent << endl;
		  


	  // Add hitId. This is iHit, if there are any hits in this event. If not, this is set to -1.
	  if (nHitsEvent==0){hitId=-1;}
	  else {hitId=iHit;}
	  
	  BCntOldEvent=BCntLoopEvent; // Update to the old bunch count value
	  tsortedEvent->Fill(); // Fill the Event tree


	  //******************
	  // Next, deal with the Hit tree. To do this, we need to use a loop to check all the Hit events which have the same BCntLoop value. We'll have to do some checking within the loop, so we'll use a break statement to get out of it (more flexible)
	  while (iHit<nentriesHit){ // Just keep looping till we break, unless the hit tree has run out of entries
	  treeHit->GetEntry(indexHit[iHit]); // Using iHit to keep track of indices, read in an event from the sorted Hit tree
	    
	    // As a test, check to see if the BCntLoop for the current hit is lower than that for event. If this is true, something has gone wrong.
	    if (BCntLoopEvent>BCntLoopHit) cout << "Error with Hit tree ordering - missed an entry somewhere  " << endl;

	    // Check to see if the BCntLoop number in this Hit match the current event. If not, break.
	    if (BCntLoopEvent!=BCntLoopHit) break;
	    
	    // Since we have not broken, this Hit must match the Event that we still have in memory. So, we want to attach the corresponding order and time info.
	    // Note that the Hit tree also uses &OrderTrue, &TimeDiff and &HitTime as the locations where the data to be added to the tree is stored. Hopefully this means that I can just use Fill and it'll work. Note that I've read in the original Hit info.

	    // Here, we check if the strip is noisy or not, and set the BadStrip flag appropriately. Need to check the hit position, PosXHit. We loop over all the strips for the current sensor.
	    BadStrip=0;
	    for(strip = 0; strip < num_noisystrips ; strip++){
	      if (noisyStrips[sensor][strip]==PosXHit) BadStrip=1;
	    }
	    
	    tsortedHit->Fill(); // Fill the Hit tree
	    
	    ++iHit; // Move on to the next entry in the sorted hit tree. Note that if the BCntLoop comparison fails, then we'll stick with the same event till the next one. This should work fine, but be on the lookout for bugs!
	  }



	  }

	// write the histos once this is finished
	tsortedEvent->Write();
	tsortedHit->Write();

	delete [] indexEvent;
	delete [] indexHit;

	if(checker==0){    
	  cout << "LCMS trees" << endl;
 	}
        else{
	  cout << "Cluster trees" << endl;
	}
	cout << sensorName[sensor] <<" done" << endl;





     } // End loop over all sensors


     } // End loop over both checkers

     // Write some output info

     cout << "No of entries in detector file: " << nentriesEvent << endl;
	
}