MpdFlowScripts/PicoDst/PicoDstConverter.cpp

#include <iostream>
#include <set>

#include <TROOT.h>
#include <TChain.h>
#include <TH1.h>
#include <TF1.h>
#include <TFile.h>
#include <TTree.h>
#include <TRandom1.h>
#include <TClonesArray.h>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include <TMath.h>
#include <TSystem.h>
#include <TVector3.h>
#include <TStopwatch.h>

#include "FairMCEventHeader.h"
#include "FairMCTrack.h"
#include "MpdEvent.h"
#include "MpdZdcDigi.h"
#include "MpdPid.h"
#include "MpdTrack.h"
#include "MpdKalmanTrack.h"
#include "MpdVertex.h"

#include "PicoDstMCEvent.h"
#include "PicoDstRecoEvent.h"
#include "PicoDstMCTrack.h"
#include "PicoDstRecoTrack.h"
#include "PicoDstFHCal.h"

int main(int argc, char **argv)
{

  TString iFileName, oFileName;
  if (argc < 5)
  {
    std::cerr << "./reducedTreeCreator -i inputfile -o outputfile" << std::endl;
    return 1;
  }
  for (int i = 1; i < argc; i++)
  {
    if (std::string(argv[i]) != "-i" &&
      std::string(argv[i]) != "-o")
    {
            std::cerr << "\n[ERROR]: Unknown parameter " << i << ": " << argv[i] << std::endl;
      return 1;
    }
    else
    {
      if (std::string(argv[i]) == "-i" && i != argc - 1)
      {
        iFileName = argv[++i];
        continue;
      }
      if (std::string(argv[i]) == "-i" && i == argc - 1)
      {
        std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
        return 1;
      }
      if (std::string(argv[i]) == "-o" && i != argc - 1)
      {
        oFileName = argv[++i];
        continue;
      }
      if (std::string(argv[i]) == "-o" && i == argc - 1)
      {
        std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
        return 1;
      }
    }
  }

  TStopwatch timer;
  timer.Start();

  ///////////////////////////////////////////
  //TChain *dstTree;
  TFile *fi = new TFile(iFileName.Data(),"read");
  TTree  *dstTree = (TTree*) fi->Get("mpdsim");

  TFile *outFile = new TFile(oFileName.Data(), "RECREATE");
  
  // TTree connfiguration
  TTree *outTree = new TTree("picodst","Pico Dst for flow calculations at MPD");

  const Double_t PIDsigM = 4.0;
  const Double_t PIDsigE = 4.0;
  const Double_t PIDenergy = 11.;
  const Double_t PIDkoeff = 1.;
  const TString PIDgenerator = "URQMD";
  const TString PIDtracking = "CF";
  const TString PIDparticles = "pikapr";
  const Int_t   Num_Of_Modules = 90;

  MpdPid *pid = new MpdPid(PIDsigM, PIDsigE, PIDenergy, PIDkoeff, PIDgenerator, PIDtracking, PIDparticles);

  // PicoDst-related variables
  PicoDstMCEvent *mcEvent = new PicoDstMCEvent();
  PicoDstRecoEvent *recoEvent = new PicoDstRecoEvent();
  TClonesArray *mcTracks = new TClonesArray("PicoDstMCTrack");
  TClonesArray *recoTracks = new TClonesArray("PicoDstRecoTrack");
  TClonesArray *fhcalModules = new TClonesArray("PicoDstFHCal");

  outTree->Branch("mcevent.",     &mcEvent,      32000,  99);
  outTree->Branch("recoevent.",   &recoEvent,    32000,  99);
  outTree->Branch("mctracks",     &mcTracks,     256000, 99);
  outTree->Branch("recotracks",   &recoTracks,   256000, 99);
  outTree->Branch("FHCalModules", &fhcalModules, 128000, 99);
  mcTracks->BypassStreamer();
  recoTracks->BypassStreamer();
  fhcalModules->BypassStreamer();

  FairMCEventHeader *MCHeader;
  TClonesArray *MCTracks;
  MpdEvent *MPDEvent;
  TClonesArray *FHCalHits;
  TClonesArray *MpdGlobalTracks;
  MpdZdcDigi *FHCalHit;
  //TClonesArray *mpdKalmanTracks;
  TClonesArray *vertexes;

  MCHeader = nullptr;
  MCTracks = nullptr;
  MPDEvent = nullptr;
  FHCalHits = nullptr;
  //mpdKalmanTracks = nullptr;
  vertexes = nullptr;

  dstTree->SetBranchAddress("MCEventHeader.", &MCHeader);
  dstTree->SetBranchAddress("MCTrack", &MCTracks);
  dstTree->SetBranchAddress("MPDEvent.", &MPDEvent);
  dstTree->SetBranchAddress("ZdcDigi", &FHCalHits);
  //dstTree->SetBranchAddress("TpcKalmanTrack", &mpdKalmanTracks);
  dstTree->SetBranchAddress("Vertex", &vertexes);

  // Starting event loop
  TVector3 primaryVertex;
  //std::set <Int_t> UsedMCTracks;
  std::map <Int_t,Int_t> InitMcRecoId; // map[init-mc-id] = reco-id
  std::map <Int_t,Int_t> InitMcNewMcId; // map[old-mc-id] = new-mc-id
  Float_t FHCalSumEnergy[Num_Of_Modules];
  Int_t n_entries = dstTree->GetEntriesFast();
  Bool_t isGoodPID;
  Short_t charge_mpd;
  for (int iEv = 0; iEv < n_entries; iEv++)
  {
    std::cout << "EVENT N " << iEv << std::endl;
    dstTree->GetEntry(iEv);

    mcTracks->Clear();
    recoTracks->Clear();
    InitMcRecoId.clear();
    InitMcNewMcId.clear();
    for (int i=0; i<Num_Of_Modules; i++)
    {
      FHCalSumEnergy[i] = 0.;
    }

    // Read MC Event
    mcEvent->SetB(MCHeader->GetB());
    mcEvent->SetPhiRP(MCHeader->GetRotZ());
    mcEvent->SetVertex(MCHeader->GetX(), MCHeader->GetY(), MCHeader->GetZ());

    // Reading Reco Event
    MpdVertex *vertex = (MpdVertex *)vertexes->First();
    vertex->Position(primaryVertex);
    recoEvent->SetVertex(primaryVertex.X(), primaryVertex.Y(), primaryVertex.Z());

    // Read energy reconstructed in the FHCal modules
    Int_t number_of_FHCal_hits = FHCalHits->GetEntriesFast();
    for(int ihit=0; ihit<number_of_FHCal_hits; ihit++)
    {
      FHCalHit = (MpdZdcDigi*) FHCalHits->UncheckedAt(ihit);
      Int_t DetId = FHCalHit->GetDetectorID();
      FHCalSumEnergy[(FHCalHit->GetModuleID()-1)+(Num_Of_Modules/2)*(DetId-1)] += FHCalHit->GetELoss();
    }
    // Fill fDetectorID = 1 part
    for (int imodule=0; imodule<Num_Of_Modules; imodule++)
    {
      PicoDstFHCal *fhcalModule = (PicoDstFHCal*)fhcalModules->ConstructedAt(imodule);
      fhcalModule->SetEnergy(FHCalSumEnergy[imodule]);
    }
    
    // Reading Reco Tracks
    MpdGlobalTracks = (TClonesArray*)MPDEvent->GetGlobalTracks();
    for (int itrack=0; itrack<MpdGlobalTracks->GetEntriesFast(); itrack++)
    {
      MpdTrack* mpdtrack = (MpdTrack*) MpdGlobalTracks->UncheckedAt(itrack);
      InitMcRecoId[mpdtrack->GetID()] = itrack;
      PicoDstRecoTrack *recoTrack = (PicoDstRecoTrack*)recoTracks->ConstructedAt(itrack);
      
      recoTrack->SetInitialMcId(mpdtrack->GetID());
      recoTrack->SetPx(mpdtrack->GetPx());
      recoTrack->SetPy(mpdtrack->GetPy());
      recoTrack->SetPz(mpdtrack->GetPz());
      recoTrack->SetNhits(mpdtrack->GetNofHits());
      recoTrack->SetNhitsPoss(mpdtrack->GetNofHitsPossTpc());
      recoTrack->SetTofMass2(mpdtrack->GetTofMass2());
      recoTrack->SetTofFlag(mpdtrack->GetTofFlag());
      recoTrack->SetTpcdEdx(mpdtrack->GetdEdXTPC());
      recoTrack->SetChi2(mpdtrack->GetChi2());
      recoTrack->SetDCAx(mpdtrack->GetDCAX());
      recoTrack->SetDCAy(mpdtrack->GetDCAY());
      recoTrack->SetDCAz(mpdtrack->GetDCAZ());
      charge_mpd = (mpdtrack->GetPt() < 0) ? 1 : -1;
      recoTrack->SetCharge(charge_mpd);

      // PID
      if (mpdtrack->GetTofFlag() == 2 || mpdtrack->GetTofFlag() == 6)
      {
        isGoodPID = pid->FillProbs(TMath::Abs(mpdtrack->GetPt()) * TMath::CosH(mpdtrack->GetEta()),
                                   mpdtrack->GetdEdXTPC(), mpdtrack->GetTofMass2(), charge_mpd);
      }
      else
      {
        isGoodPID = pid->FillProbs(TMath::Abs(mpdtrack->GetPt()) * TMath::CosH(mpdtrack->GetEta()),
                                   mpdtrack->GetdEdXTPC(), charge_mpd);
      }

      if (isGoodPID)
      {
        recoTrack->SetPidProbPion(pid->GetProbPi());   //mpdtrack->GetTPCPidProbPion();
        recoTrack->SetPidProbKaon(pid->GetProbKa());   //mpdtrack->GetTPCPidProbKaon();
        recoTrack->SetPidProbProton(pid->GetProbPr()); //mpdtrack->GetTPCPidProbProton();
      }
      else
      {
        recoTrack->SetPidProbPion(-999.);
        recoTrack->SetPidProbKaon(-999.);
        recoTrack->SetPidProbProton(-999.);
      }
    }

    // Reading MC Tracks
    int mcTrackCounter = 0;
    for (int imctrack=0; imctrack<MCTracks->GetEntriesFast(); imctrack++)
    {
      FairMCTrack *mctrack = (FairMCTrack*) MCTracks->UncheckedAt(imctrack);
      bool isUsed = (InitMcRecoId.count(imctrack));

      // If motherId != 1 and mc track doesn't relate to any reco track - skip
      //if (mctrack->GetMotherId() != -1) continue;
      if (mctrack->GetMotherId() != -1 && !isUsed) continue;
      if (isUsed) //std::cout << "VVVVVVV" << std::endl;
      {
        InitMcNewMcId[imctrack] = mcTrackCounter;
      }
      PicoDstMCTrack *mcPicoTrack = (PicoDstMCTrack*)mcTracks->ConstructedAt(mcTrackCounter);
      mcPicoTrack->SetInitialId(imctrack);
      mcPicoTrack->SetPx(mctrack->GetPx());
      mcPicoTrack->SetPy(mctrack->GetPy());
      mcPicoTrack->SetPz(mctrack->GetPz());
      mcPicoTrack->SetMotherId(mctrack->GetMotherId());
      mcPicoTrack->SetPdg(mctrack->GetPdgCode());
      mcPicoTrack->SetEnergy(mctrack->GetEnergy());
      mcTrackCounter++;
    }

    for (int itrack=0; itrack<recoTracks->GetEntriesFast(); itrack++)
    {
      PicoDstRecoTrack *recoTrack = (PicoDstRecoTrack*) recoTracks->UncheckedAt(itrack);
      MpdTrack* mpdtrack = (MpdTrack*) MpdGlobalTracks->UncheckedAt(itrack);
      recoTrack->SetMcId(InitMcNewMcId[mpdtrack->GetID()]);
    }

    outTree->Fill();
  } // End of the Event loop

  outFile->cd();
  outTree->Print();
  outTree->Write();
  outFile->Close();

  timer.Stop();
  timer.Print();

  return 0;
}