// --------------------------------------------------------------------------
//
// CBM Computing School, Kolkata, February 2018
// Macro for Exercise 2: Analysis of raw (digi) data
// V. Friese   15/02/2018
//
// The macro compares the time of StsDigis to te corresponding StsPoint,
// i.e. it checks the time resolution of the STS detector.
//
// --------------------------------------------------------------------------

void ana_raw() {
	
  // We start by defining the time difference histogram
  TH1F* hDeltaT = new TH1F("dt", "time difference", 100, -100., 100.);	
	
  // First we have to open the MC file created with the transport macro
  TFile* mcFile = new TFile("exercise.mc.root");
  
  // In addition, we open the raw data file created by digitization
  TFile* rawFile = new TFile("exercise.raw.root");
  
  // Get the MC tree and the raw data tree
  TTree* mcTree = (TTree*) mcFile->Get("cbmsim");
  TTree* rawTree = (TTree*) rawFile->Get("cbmsim");
  
  // For our analysis, we need the branches of StsDigi, StsDigiMatch and StsPoint:
  
  // Get the StsDigi branch (from the raw data tree)
  TClonesArray* digis = new TClonesArray("CbmStsDigi");
  rawTree->SetBranchAddress("StsDigi", &digis);

  // Get the StsDigiMatch branch (from the raw data tree)
  TClonesArray* matches = new TClonesArray("CbmMatch");
  rawTree->SetBranchAddress("StsDigiMatch", &matches);
  
  // Get the StsPoint branch (from the MC tree)
  TClonesArray* points = new TClonesArray("CbmStsPoint");
  mcTree->SetBranchAddress("StsPoint", &points);  
  
  // Loop over events
  Int_t nEvents = rawTree->GetEntries();
  for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
  	rawTree->GetEntry(iEvent);
  	mcTree->GetEntry(iEvent);
  	
  	// Loop over StsDigis in this event
  	Int_t nDigis = digis->GetEntriesFast();
  	std::cout << "Event " << iEvent << ", STS digis: " << nDigis << std::endl;
  	for (Int_t iDigi = 0; iDigi < nDigis; iDigi++) {
  		CbmStsDigi* digi = (CbmStsDigi*) digis->At(iDigi);
  		
  		// Get the corresponding match object. Note that is it at the
  		// same index in its array as the digi itself.
  		CbmMatch* match = (CbmMatch*) matches->At(iDigi);
  		
  		// Get the StsPoint index from the best link in the match object
  		// (that is the one with the largest weight)
  		Int_t index = match->GetMatchedLink().GetIndex();
  		
  		// Get the StsPoint object at this index
  		CbmStsPoint* point = (CbmStsPoint*) points->At(index);
  		
  		// Calculate the time difference of digi and point
  		Double_t deltaT = digi->GetTime() - point->GetTime();
  		
  		// Fill the time difference in the histogram
  		hDeltaT->Fill(deltaT);	
  		
  	} // Loop over StsDigis
  	
  } // Event loop
  	
  	// Close the input files
  	mcFile->Close();
  	rawFile->Close();
  	
  	// Plot the histogram
  	hDeltaT->Draw();
  	
  	// Fit with a Gauss function
  	hDeltaT->Fit("gaus");
  	
} // End of macro function