// --------------------------------------------------------------------------
//
// CBM Computing School, Kolkata, February 2018
// Macro for Exercise 1: Analysis of MC data
// V. Friese   15/02/2018
//
// --------------------------------------------------------------------------

void ana_mc() {
	
  // We start by defining the momentum histogram
  TH1F* hMom = new TH1F("momentum", "muon momentum", 100, 0., 11.);	
	
  // First we have to open the MC file created with the transport macro
  TFile* inFile = new TFile("exercise.mc.root");
  
  // Get the tree
  TTree* tree = (TTree*) inFile->Get("cbmsim");
  
  // Get the MCTrack branch
  TClonesArray* tracks = new TClonesArray("CbmMCTrack");
  tree->SetBranchAddress("MCTrack", &tracks);

  // Get the TrdPoint branch
  TClonesArray* trdpoints = new TClonesArray("CbmTrdPoint");
  tree->SetBranchAddress("TrdPoint", &trdpoints);
  
  // Loop over events
  Int_t nEvents = tree->GetEntries();
  for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
  	tree->GetEntry(iEvent);
  	
  	// Loop over MCTracks in this event
  	Int_t nTracks = tracks->GetEntriesFast();
  	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
  		CbmMCTrack* track = (CbmMCTrack*) tracks->At(iTrack);
  		
  		// Check if it is a mu-. Else next track.
  		if ( track->GetPdgCode() != 13 ) continue;
  		
  		// Check whether the track has points in the TRD
  		if ( track->GetNPoints(kTrd) == 0 ) continue;
  		
  		// Reaching here, the track has reached the TRD, i.e.
  		// it has passed through all muon absorbers. We can enter
  		// its momentum into the histogram.
  		hMom->Fill(track->GetP());	
  		
  	} // Loop over MCTracks	
  	
  } // Event loop
  	
  	// Close the input file
  	inFile->Close();
  	
  	// Plot the histogram
  	hMom->Draw();
  	
} // End of macro function