// --------------------------------------------------------------------------
//
// CBM Computing School, Kolkata, February 2018
// Macro for Exercise 3: Analysis of reco data
// V. Friese   15/02/2018
//
// The macro compares the x position of MuchHits with the corresponding MuchPoints,
// i.e. it checks the position resolution of the MUCH detector.
//
// --------------------------------------------------------------------------

void ana_reco() {
	
  // We start by defining the time difference histogram
  TH1F* hDeltaX = new TH1F("dx", "x 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 reco data file created by reconstruction
  TFile* recoFile = new TFile("exercise.reco.root");
  
  // Get the MC tree and the raw data tree
  TTree* mcTree = (TTree*) mcFile->Get("cbmsim");
  TTree* recoTree = (TTree*) recoFile->Get("cbmsim");
  
  // For our analysis, we need the branches of MuchHit, MuchHitMatch and MuchPoint.
  
  // Get the MuchHit branch (from the reco data tree)
  TClonesArray* hits = new TClonesArray("CbmMuchPixelHit");
  recoTree->SetBranchAddress("MuchPixelHit", &hits);

  // Get the MuchHitMatch branch (from the reco data tree)
  TClonesArray* matches = new TClonesArray("CbmMatch");
  recoTree->SetBranchAddress("MuchPixelHitMatch", &matches);
  
  // Get the MuchPoint branch (from the MC tree)
  TClonesArray* points = new TClonesArray("CbmMuchPoint");
  mcTree->SetBranchAddress("MuchPoint", &points);  
  
  // Loop over events
  Int_t nEvents = recoTree->GetEntries();
  for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
  	recoTree->GetEntry(iEvent);
  	mcTree->GetEntry(iEvent);
  	
  	// Loop over MuchHits in this event
  	Int_t nHits = hits->GetEntriesFast();
  	std::cout << "Event " << iEvent << ", MUCH hits: " << nHits << std::endl;
  	for (Int_t iHit = 0; iHit < nHits; iHit++) {
  		CbmMuchPixelHit* hit = (CbmMuchPixelHit*) hits->At(iHit);
  		assert(hit);
  		
  		// Get the corresponding match object. Note that is it at the
  		// same index in its array as the hit itself.
  		CbmMatch* match = (CbmMatch*) matches->At(iHit);
  		assert(match);
  		
  		// Get the MuchPoint 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 MuchPoint object at this index
  		CbmMuchPoint* point = (CbmMuchPoint*) points->At(index);
  		assert(point);
  		
  		// Calculate the x difference of hit and point
  		Double_t deltaX = hit->GetX() - ( 0.5 * point->GetXIn() + point->GetXOut() );
  		
  		// Fill the time difference in the histogram
  		hDeltaX->Fill(deltaX);	
  		
  	} // Loop over MuchHits
  	
  } // Event loop
  	
  	// Close the input files
  	mcFile->Close();
  	recoFile->Close();
  	
  	// Plot the histogram
  	hDeltaX->Draw();
  	
  	// Fit with a Gauss function
  	hDeltaX->Fit("gaus");
  	
} // End of macro function