MCandRECTrkmatches_exe.C

Go to the documentation of this file.

#include <iostream>
#include <sstream>

#include<TApplication.h>
#include<TCanvas.h>
#include<TROOT.h>
#include<TString.h>
#include<TChain.h>
#include<TFile.h>
#include<TClonesArray.h>
#include<TSystem.h>
#include<TH1.h>
#include<TH2.h>
#include<TF1.h>
#include<TRotation.h>
#include<TVector3.h>
#include<TMath.h>
#include<TGaxis.h>
#include<TNtuple.h>
#include<TMatrixD.h>
#include<TStopwatch.h>
#include<TNtuple.h>
#include<PndMCTrack.h>
#include<PndTrack.h>
#include<PndSdsMCPoint.h>
#include "PndLmdDim.h"

//lmd track
#include<PndLinTrack.h>
#include<PndTrackCand.h>
#include<PndTrackCandHit.h>
#include<PndSdsHit.h>
#include<PndSdsClusterStrip.h>
#include<PndSdsDigiStrip.h>
#include<PndSdsDigiPixel.h>
#include<PndSdsClusterPixel.h>
#include<PndSdsMergedHit.h>

// needed for geane backtracking
#include<FairRunAna.h>
#include<FairRootManager.h>
#include<FairGeane.h>
#include<FairRtdbRun.h>
#include<FairRuntimeDb.h>
#include<FairParRootFileIo.h>
#include<FairLogger.h>
#include<FairTrackParH.h>

#include<string>
// void combitransFromLumiFrame(TVector3& hitPos, int senType=0){
//   Double_t  kTransZ, kRotUmZ, kRot, kTransX;
//   if(senType>0){ //Pixel sensors
//     //do the transformation from lab frame to LUMI frame (with z-axis perp. to lumi planes)
//     double end_seg_upstream = 360.1; // where bending starts with
//     double r_bend = 5750.; // the bending radius
//     double phi_bend = 40.068e-3; // and the angle of the circle path
//     // const Double_t kRot = phi_bend/3.141*180.;//=2.295727;//2.326; //(deg) //Rotate to dipol
//     kRot = phi_bend;//here we need abgle in rad!
//     // the point where both tangents of the straight beam pipe tubes meet is
//     kRotUmZ = end_seg_upstream + tan(phi_bend/2.)*r_bend;//476.03; //(cm) //z-point to rotate
//     kTransZ = 1130.; //(cm) //move at z-position
//     kTransX = (kTransZ - end_seg_upstream - tan(phi_bend/2.)*r_bend)*tan(phi_bend); // 25 (cm) //move at x-position
//   }
//   else{
//     /// THIS WORKS ONLY FOR DESIGN WITH STRIP SENSORS!
//     //do the transformation from lab frame to LUMI frame (with z-axis perp. to lumi planes)
//     kTransZ = 1100.; //(cm) //move at z-position
//     kRotUmZ = 476.03; //(cm) //z-point to rotate
//     kTransX = 25; //(cm) //move at x-position
//     kRot =  0.040596358401388; // 2.326 degree  = 4.05963584013881024e-02 rad
//   }
//   TVector3 LumiTrans(0,0,kTransZ-kRotUmZ);
//   hitPos +=LumiTrans;
//   hitPos.RotateY(kRot);
//   LumiTrans = TVector3(0,0,kRotUmZ);
//   hitPos +=LumiTrans;
// }
// TMatrixD rotateFromLumiFrame(TMatrixD& hitCov, int senType=0){
//   Double_t  sintheta, costheta;
//   if(senType>0){ //Pixel sensors
//     double phi_bend = 40.068e-3; // and the angle of the circle path
//     sintheta = TMath::Sin(phi_bend);
//     costheta = TMath::Cos(phi_bend);
//   }
//   else{
//     /// THIS WORKS ONLY FOR DESIGN WITH STRIP SENSORS!
//     Double_t theta=-2.326;
//     Double_t degrad = TMath::Pi()/180.;
//     sintheta = TMath::Sin(degrad*theta);
//     costheta = TMath::Cos(degrad*theta);
//   }
//   TMatrixD rot(3,3);// Rotation around Y axis
//   rot[0][0]= costheta;
//   rot[0][1]= 0;
//   rot[0][2]= sintheta;
//   rot[1][0]= 0;
//   rot[1][1]= 1;
//   rot[1][2]= 0;
//   rot[2][0]= -sintheta;
//   rot[2][1]= 0;
//   rot[2][2]= costheta;
//   TMatrixD result = rot;
//   result.T();
//   result*=hitCov;
//   hitCov = result;
//   result*=rot;
//   return result;
// }


// void rotateFromLumiFrame(TVector3& hitPos, bool err, int senType=0){
//   TMatrixD hitMtx(3,3);
//   if(err){
//     hitMtx[0][0] = hitPos[0]*hitPos[0];
//     hitMtx[1][1] = hitPos[1]*hitPos[1];
//     hitMtx[2][2] = hitPos[2]*hitPos[2];
//     TMatrixD res = rotateFromLumiFrame(hitMtx,senType);
//     hitPos = TVector3(sqrt(res(0,0)),sqrt(res(1,1)),sqrt(res(2,2)));
//   }
//   else{
//     hitMtx[0][0] = hitPos[0];
//     hitMtx[1][0] = hitPos[1];
//     hitMtx[2][0] = hitPos[2];
//     TMatrixD res = rotateFromLumiFrame(hitMtx,senType);
//     hitPos = TVector3(hitMtx(0,0),hitMtx(1,0),hitMtx(2,0));
//   }
// }

using namespace std;
int main(int __argc,char *__argv[]) {
  //gROOT->Macro("/PANDA/pandaroot/macro/lmd/Style_Imported_Style.C");
  // gROOT->SetStyle("Imported_Style");

  //TODO: read this like params!
  //  const int nEvents=500000;
  int nEvents=1000;
  int nMCtracks=1;
  int startEvent=0;
  TString storePath="/data/FAIRsorf/pandaroot/trunk/macro/lmd/tmpOutput";
  double Plab=15.;
  int verboseLevel=5;
  int mh = 0;//if>0 : use merged hit collection
  int dnu=0;//if>0 : use namespace for pixel
  std::string startStr="", momStr="", nStr="", pathStr="", verbStr="" , mcTrkStr="", useNewDStr="",usedMHStr="";
  // decode arguments
  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
                    || strcmp( __argv[1], "--help" ) == 0 ) ){

    std::cout << "This is script for comparision reconstructed and simulated tracks with parameters\n"
              <<"-s start event \n"
              <<"-n Number of events \n"
              <<"-t Number of trks per event \n"
              <<"-mom Beam Momentum \n"
              <<"-path path to the file(s) \n"
              <<"-v verbose Level (if>0, print out some information) \n"
              <<"-npx if>0: use namespace for pixel \n"
              <<"-mh if>0: use merged hits\n"
              <<"Have fun! \n"
              << std::endl;
    return 0;
  }
  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
    bool found=false;
    std::string sw = __argv[optind];
    if (sw=="-s") {
      optind++;
      startStr = __argv[optind];
      found=true;
    }
    if (sw=="-n"){
      optind++;
      nStr = __argv[optind];
      found=true;
    }
    if (sw=="-t"){
      optind++;
      mcTrkStr = __argv[optind];
      found=true;
    }
    if (sw=="-path"){
      optind++;
      pathStr = __argv[optind];
      found=true;
    }
    if (sw=="-mom"){
      optind++;
      momStr = __argv[optind];
      found=true;
    }
    if (sw=="-v"){
      optind++;
      verbStr = __argv[optind];
      found=true;
    }
    if (sw=="-npx"){
      optind++;
      useNewDStr = __argv[optind];
      found=true;
    }
    if (sw=="-mh"){
      optind++;
      usedMHStr = __argv[optind];
      found=true;
    }
    if (!found){
      std::cout<< "Unknown switch: "
               << __argv[optind] <<std::endl;
      optind++;
    }

  while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
  }

  std::stringstream startSStr(startStr), momSStr(momStr), nSStr(nStr), pathSStr(pathStr), verbSStr(verbStr), mcTrkSStr(mcTrkStr),useNewDSStr(useNewDStr), usedMHSStr(usedMHStr);

  startSStr >> startEvent;
  momSStr >> Plab;
  nSStr >> nEvents;
  pathSStr >> storePath;
  verbSStr >> verboseLevel;
  mcTrkSStr >> nMCtracks;
  useNewDSStr >> dnu;
  usedMHSStr >> mh;
  cout<<"====================================="<<endl;
  cout<<"some INFO about this macro params: "<<endl;
  cout<<"expected number of tracks per event: "<<nMCtracks<<endl;
  cout<<"expected number of events: "<<nEvents<<endl;
  cout<<"Pbeam: "<<Plab<<endl;
  cout<<"Will be used Path: "<<storePath<<endl;
  cout<<"====================================="<<endl;
  //void MCandRECTrkmatches(const int nEvents=2, const int startEvent=0, TString storePath="tmpOutput", const int verboseLevel=3, double dv=0.5, bool no4d=false)
  //{
  // ----  Load libraries   -------------------------------------------------
  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
  gSystem->Load("libLmdTrk");
  // gROOT->LoadMacro("line3Dfit.C");
  // ------------------------------------------------------------------------

  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

  // ---- Input files --------------------------------------------------------
  TString simMC=storePath+"/Lumi_MC_";
  simMC += startEvent;
  simMC += ".root";
  TChain tMC("pndsim");
  tMC.Add(simMC);

  TString DigiFile = storePath+"/Lumi_digi_";
  DigiFile += startEvent;
  DigiFile += ".root";
  TChain tdigiHits("pndsim");
  tdigiHits.Add(DigiFile);

  TString recHit=storePath+"/Lumi_reco";
  // if(dnu>0) recHit+="Merged";
  recHit +="_";
  recHit += startEvent;
  recHit += ".root";
  TChain tHits("pndsim");
  tHits.Add(recHit);

  TString recHitmerged=storePath+"/Lumi_recoMerged_";
  recHitmerged += startEvent;
  recHitmerged += ".root";
  TChain tHitsMerged("pndsim");
  if(mh)
    tHitsMerged.Add(recHitmerged);


  TString trkCand = storePath+"/Lumi_TCand_";
  trkCand += startEvent;
  trkCand += ".root";
  TChain tTrkCand("pndsim");
  tTrkCand.Add(trkCand);

  TString recTrack;
  TChain tTrkRec("pndsim");
  recTrack = storePath+"/Lumi_Track_";
  recTrack += startEvent;
  recTrack += ".root";
  tTrkRec.Add(recTrack);
  // if(kf<1){
  //   recTrack = storePath+"/Lumi_Track_";
  //   recTrack += startEvent;
  //   recTrack += ".root";
  //   tTrkRec.Add(recTrack);
  // }

  TString geaneFile = storePath+"/Lumi_Geane_";
  geaneFile += startEvent;
  geaneFile += ".root";
  TChain tgeane("pndsim");
  tgeane.Add(geaneFile);

  // ---------------------------------------------------------------------------------

  // ---- Output file ----------------------------------------------------------------
  TString out=storePath+"/Lumi_out_MC_and_REC_trks_matches_with_IDs";
  out += startEvent;
  out += ".root";
  TFile *f = new TFile(out,"RECREATE");
  // ---------------------------------------------------------------------------------

  //--- MC info -----------------------------------------------------------------
  TClonesArray* true_tracks=new TClonesArray("PndMCTrack");
  tMC.SetBranchAddress("MCTrack",&true_tracks);  //True Track to compare

  TClonesArray* true_points=new TClonesArray("PndSdsMCPoint");
  tMC.SetBranchAddress("LMDPoint",&true_points);  //True Points to compare
  //----------------------------------------------------------------------------------


  //--- Digitization info ------------------------------------------------------------
  TClonesArray* fStripClusterArray;
  if(dnu>0){
    fStripClusterArray = new TClonesArray("PndSdsClusterPixel");
    tHits.SetBranchAddress("LMDPixelClusterCand",&fStripClusterArray);
  }
  else{
    fStripClusterArray = new TClonesArray("PndSdsClusterStrip");
    tHits.SetBranchAddress("LMDStripClusterCand",&fStripClusterArray);
  }


  TClonesArray* fStripDigiArray;
  if(dnu>0){
    fStripDigiArray = new TClonesArray("PndSdsDigiPixel");
    tdigiHits.SetBranchAddress("LMDPixelDigis",&fStripDigiArray);
  }
  else{
    fStripDigiArray = new TClonesArray("PndSdsDigiStrip");
    tdigiHits.SetBranchAddress("LMDStripDigis",&fStripDigiArray);
  }
  //----------------------------------------------------------------------------------

  //--- Real Hits --------------------------------------------------------------------
  TClonesArray* rechit_array;
  if(dnu>0){
    if(mh>0){
      rechit_array = new TClonesArray("PndSdsMergedHit");
      tHitsMerged.SetBranchAddress("LMDHitsMerged",&rechit_array);  //Points for Tracks
    }
    else{
      rechit_array = new TClonesArray("PndSdsHit");
      tHits.SetBranchAddress("LMDHitsPixel",&rechit_array);  //Points for Tracks
    }
  }
  else{
    rechit_array = new TClonesArray("PndSdsHit");
    tHits.SetBranchAddress("LMDHitsStrip",&rechit_array);  //Points for Tracks
  }
  //----------------------------------------------------------------------------------


  //--- Track Candidate ---------------------------------------------------------------
  TClonesArray* trkcand_array=new TClonesArray("PndTrackCand");
  tTrkCand.SetBranchAddress("LMDTrackCand",&trkcand_array); //Points for Track Canidates
  //-----------------------------------------------------------------------------------

  //--- Real tracks -------------------------------------------------------------------
  TClonesArray* rec_trk=new TClonesArray("PndLinTrack");
  tTrkRec.SetBranchAddress("LMDTrack",&rec_trk);  //Tracks
  //  tgeane.SetBranchAddress("PndTrackLmd",&rec_trk);  //Tracks
  // if(kf>0){
  //   rec_trk = new TClonesArray("PndTrack");
  //   tgeane.SetBranchAddress("PndTrackLmd",&rec_trk);  //Tracks
  // }
  // else{
  //   rec_trk = new TClonesArray("PndLinTrack");
  //   tTrkRec.SetBranchAddress("LMDTrack",&rec_trk);  //Tracks
  // }

  cout<<"Here we go to GEANE file: "<<endl;
  //----------------------------------------------------------------------------------

  //--- Geane info ------------------------------------------------------------------
  // TClonesArray* geaneArray =new TClonesArray("FairTrackParP");
  TClonesArray* geaneArray =new TClonesArray("FairTrackParH");
  tgeane.SetBranchAddress("GeaneTrackFinal",&geaneArray);  //Tracks with parabolic parametrisation

  cout<<"And we'll make some hists"<<endl;
  //--- Output histogram -----------------------------------------------------
  TH1 *hchi2 = new TH1F("hchi2","#chi^2 for reconstructed tracks;#chi^2;",1.5e2,0,15.);
  TH2 *hnRecnMC = new TH2F("hnRecnMC","Number reconstracted tracks vs. Number simulated tracks; N_{MC}; N_{rec}",
                           100,0,100,100,0,100);
  TH1 *hDiffIDs = new TH1F("hDiffIDs","Number of track-candidates with hits from diff. MC-track;N_{IDs}",10,0,10);
  TH2 *hMCtrkvshits = new TH2F("hMCtrkvshits","Number of simulated tracks vs. Number of rec. hits",
                               1000,0,1000,100,0,100);
  TH1 *hntrkcand = new TH1F("hntrkcand","Number of track-candidates per event;N_{trk-cand}",100,0,100);
  TH1 *hntrk = new TH1F("hntrk","Number of tracks per event;N_{trk}",100,0,100);
  TH1 *hntrkmissed_I = new TH1F("hntrkmissed_I","Number of missed tracks per event (#phi, #theta trks);N_{trk}",30,0,30);
  TH1 *hntrkgood_I = new TH1F("hntrkgood_I","Number of good tracks per event (#phi, #theta trks);N_{trk}",30,0,30);
  TH1 *hntrkghost_I = new TH1F("hntrkghost_I","Number of ghost tracks per event (#phi, #theta trks);N_{trk}",30,0,30);
  TH1 *hntrkmissed_II = new TH1F("hntrkmissed_II","Number of missed tracks per event (hits matching);N_{trk}",30,0,30);
  TH1 *hntrkgood_II = new TH1F("hntrkgood_II","Number of good tracks per event (hits matching);N_{trk}",30,0,30);
  TH1 *hntrkghost_II = new TH1F("hntrkghost_II","Number of ghost tracks per event (hits matching);N_{trk}",30,0,30);
  TH2 *hntrkmissedPhiTheta = new TH2F("hntrkmissedPhiTheta",";#delta#theta/#sigma#theta;#delta#phi/#sigma#phi",1000,0,100,1000,0,100);
  TH2 *hntrkcandvsIDs = new TH2F("hntrkcandvsIDs","Number of track-candidates per event vs Number of track-candidates with hits from diff. MC-track;N_{IDs};N_{trk-cand}",10,0,10,100,0,100);
  TH2 *hntrkcandvsMC = new TH2F("hntrkcandvsMC","Number of per event vs. Number of simulated tracks; N_{MC};N_{trk-cand}",100,0,100,100,0,100);

  TH2 *hnhitsvsIDs = new TH2F("hnhitsvsIDs","Number of hits per track-candidate vs Number of track-candidates with hits from diff. MC-track;N_{IDs};Number of hits",10,0,10,10,2,12);

  TH1 * hallTheta = new TH1F("hallTheta","#theta for all MC-track;#theta, rad",50,2e-3,10e-3);
  TH1 * hgoodTheta = new TH1F("hgoodTheta","#theta for good track-candidates",50,2e-3,10e-3);
  TH1 * hallPhi = new TH1F("hallPhi","#phi for all MC-track;#phi, rad",50,-355./113,355./113);
  TH1 * hgoodPhi = new TH1F("hgoodPhi","#phi for good track-candidates",50,-355./113,355./113);
  TH2 * hgoodPhiTheta = new TH2F("hgoodPhiTheta","#phi and #thete for good track-candidates",50,2e-3,10e-3,100,-355./113,355./113);
 TH2 * hallPhiTheta = new TH2F("hallPhiTheta","#phi & #theta for all MC",50,2e-3,10e-3,100,-355./113,355./113);

  TH2 * hgoodPhichi2 = new TH2F("hgoodPhichi2","#phi for good track-candidates vs. #chi^{2} ",5e2,0,50,20,-355./113,355./113);
  TH2 * hgoodThetachi2 = new TH2F("hgoodThetachi2","#theta for good track-candidates vs. #chi^{2} ",5e2,0,50,20,3.5e-3,8.5e-3);

  //  TH1 *hMomMC = new TH1F("hMomMC","P_{MC};P,GeV/c",1e3,0,16.);
  TNtuple *nmomMC = new TNtuple("nmomMC","MCmomentum","p0:p_nearLUMI:p_afterLUMI");
  TNtuple *ntuprecTrk = new TNtuple("ntuprecTrk","Info about reconstructed trks: goodTrk [0=good,+1=ghost]","x:y:z:mom:theta:phi:goodTrk");
  TNtuple *ntupMCTrk = new TNtuple("ntupMCTrk","Info about simulated trks: goodTrk [0=good, -1=missed]","x:y:z:mom:theta:phi:goodTrk:nrechits");
  //  TH1 *hMomMCnearLUMI = new TH1F("hMomMCnear","P_{MC};P,GeV/c",1e3,0,16.);

  TH1 *hResMom = new TH1F("hResMom","P_{MC}-P_{rec};#deltaP,GeV/c",1e3,-1e-4,1e-4);
  TH1 *hErrMom = new TH1F("hErrMom","#sigma_{P};#sigmaP,GeV/c",1e3,0,1e-3);
  TH1 *hPullMom = new TH1F("hPullMom","(P_{MC}-P_{rec})/#sigma_{P};",1e3,-1e1,1e1);
  //  TH1 *hResTheta = new TH1F("hResTheta","#theta_{MC}-#theta_{rec};#delta#theta,rad",1e3,-1e-2,1e-2);
  TH1 *hResTheta = new TH1F("hResTheta","#theta_{MC}-#theta_{rec};#delta#theta,rad",1e2,-1e-3,1e-3);//TEST
  TH1 *hErrTheta = new TH1F("hErrTheta","#sigma(#theta_{rec});#sigma,rad",1e3,0,0.01);
  TH1 *hPullTheta = new TH1F("hPullTheta","(#theta_{MC}-#theta_{rec})/#sigma_{#theta};",1e2,-10,10);
  TH1 *hResPhi = new TH1F("hResPhi","#phi_{MC}-#phi_{rec};#delta#phi,rad",2e3,-1.,1.);
  TH1 *hErrPhi = new TH1F("hErrPhi","#sigma(#phi_{rec});#sigma,rad",1e3,0,0.1);
  TH1 *hPullPhi = new TH1F("hPullPhi","(#phi_{MC}-#phi_{rec})/#sigma_{#phi};",1e2,-10,10);

  TH1 *hResLumiTrkMom = new TH1F("hResLumiTrkMom","P_{MC}-P_{rec}(near Lumi);#deltaP,GeV/c",1e3,-6e-7,6e-7);
  TH1 *hResLumiTrkTheta = new TH1F("hResLumiTrkTheta","#theta_{MC}-#theta_{rec}(near Lumi);#delta#theta,rad",1e3,-6e-3,6e-3);
  TH1 *hResLumiTrkPhi = new TH1F("hResLumiTrkPhi","#phi_{MC}-#phi_{rec}(near Lumi);#delta#phi,rad",2e3,-1e-1,1e-1);

  TH1 *hangMCRec = new TH1F("hangMCRec", "Angle between Pmc and Prec;#delta#alpha, rad", 1e3,0.,0.1);
  TH1 *hResPointX = new TH1F("hResPointX","X_{MC}-X_{rec};#deltaX,cm",1e2,-2.,2.);
  TH1 *hResPointY = new TH1F("hResPointY","Y_{MC}-Y_{rec};#deltaY,cm",1e2,-2.,2.);
  TH1 *hResPointZ = new TH1F("hResPointZ","Z_{MC}-Z_{rec};#deltaZ,cm",1e3,-0.15,0.15);

  TH1 *hResLumiTrkPointX = new TH1F("hResLumiTrkPointX","X_{MC}-X_{rec}(near Lumi);#deltaX,cm",1e2,-0.02,0.02);
  TH1 *hResLumiTrkPointY = new TH1F("hResLumiTrkPointY","Y_{MC}-Y_{rec}(near Lumi);#deltaY,cm",1e2,-0.02,0.02);
  TH1 *hResLumiTrkPointZ = new TH1F("hResLumiTrkPointZ","Z_{MC}-Z_{rec}(near Lumi);#deltaZ,cm",1e2,-0.02,0.02);

  TH1 *hResLumiTrkPointPx = new TH1F("hResLumiTrkPointPx","Px_{MC}-Px_{rec}(near Lumi);#deltaPx, GeV/c",1e2,-0.001,0.001);
  TH1 *hResLumiTrkPointPy = new TH1F("hResLumiTrkPointPy","Py_{MC}-Py_{rec}(near Lumi);#deltaPy, GeV/c",1e2,-0.001,0.001);
  TH1 *hResLumiTrkPointPz = new TH1F("hResLumiTrkPointPz","Pz_{MC}-Pz_{rec}(near Lumi);#deltaPz, GeV/c",1e2,-0.0001,0.0001);

  TH1 *hResLumiTrkPointXPull = new TH1F("hResLumiTrkPointXPull","(X_{MC}-X_{rec})/#sigma;(X_{MC}-X_{rec})/#sigma",1e2,-10.,10.);
  TH1 *hResLumiTrkPointYPull = new TH1F("hResLumiTrkPointYPull","(Y_{MC}-Y_{rec})/#sigma;(Y_{MC}-Y_{rec})/#sigma",1e2,-10.,10.);
  TH1 *hResLumiTrkPointZPull = new TH1F("hResLumiTrkPointZPull","(Z_{MC}-Z_{rec})/#sigma;(Z_{MC}-Z_{rec})/#sigma",1e2,-10.,10.);

  TH1 *hResLumiTrkPointPxPull = new TH1F("hResLumiTrkPointPxPull","(Px_{MC}-Px_{rec})/#sigma;(Px_{MC}-Px_{rec})/#sigma",1e2,-10,10.);
  TH1 *hResLumiTrkPointPyPull = new TH1F("hResLumiTrkPointPyPull","(Py_{MC}-Py_{rec})/#sigma;(Py_{MC}-Py_{rec})/#sigma",1e2,-10,10);
  TH1 *hResLumiTrkPointPzPull = new TH1F("hResLumiTrkPointPzPull","(Pz_{MC}-Pz_{rec})/#sigma;(Pz_{MC}-Pz_{rec})/#sigma",1e2,-10,10);

  // TH1 *hResLumiTrkTheta = new TH1F("hResLumiTrkTheta","#theta_{MC}-#theta_{rec}(near Lumi);#delta#theta, rad",1e2,-0.001,0.001);
  //  TH1 *hResLumiTrkPhi = new TH1F("hResLumiTrkPhi","#phi_{MC}-#phi_{rec}(near Lumi);#delta#phi, rad",1e2,-0.01,0.01);
  TH1 *hResLumiTrkThetaPull = new TH1F("hResLumiTrkThetaPull","(#theta_{MC}-#theta_{rec})/#sigma (near Lumi);#delta#theta, rad",1e2,-10,10);
  TH1 *hResLumiTrkPhiPull = new TH1F("hResLumiTrkPhiPull","(#phi_{MC}-#phi_{rec})/#sigma (near Lumi);#delta#phi, rad",1e2,-10,10);


  TH1 *hErrPointX = new TH1F("hErrPointX","#sigma_{X};#sigmaX,cm",1e3,0,1.);
  TH1 *hErrPointY = new TH1F("hErrPointY","#sigma_{Y};#sigmaY,cm",1e3,0,1.);
  TH1 *hErrLinPointX = new TH1F("hErrLinPointX","#sigma_{X};#sigmaX,cm",1e3,0,0.01);
  TH1 *hErrLinPointY = new TH1F("hErrLinPointY","#sigma_{Y};#sigmaY,cm",1e3,0,0.001);
  TH1 *hErrPointZ = new TH1F("hErrPointZ","#sigma_{Z};#sigmaZ,cm",1e2,0,0.01);
  TH1 *hPullPointX = new TH1F("hPullPointX","(X_{MC}-X_{rec})/#sigma_{X};(X_{MC}-X_{rec})/#sigma_{X}",1e2,-10,10);
  TH1 *hPullPointY = new TH1F("hPullPointY","(Y_{MC}-Y_{rec})/#sigma_{Y};(Y_{MC}-Y_{rec})/#sigma_{Y}",1e2,-10,10);
  TH1 *hPullPointZ = new TH1F("hPullPointZ","(Z_{MC}-Z_{rec})/#sigma_{Z};(Z_{MC}-Z_{rec})/#sigma_{Z}",1e3,-100,100);

  TH2 *hPullPointXphi = new TH2F("hPullPointXphi","(X_{MC}-X_{rec})/#sigma_{X} vs. #phi;#phi_{MC}, rad;(X_{MC}-X_{rec})/#sigma_{X}",2e1,-355./113,355./113,1e2,-10,10);
TH2 *hPullPointYphi = new TH2F("hPullPointYphi","(Y_{MC}-Y_{rec})/#sigma_{Y} vs. #phi;#phi_{MC}, rad;(Y_{MC}-Y_{rec})/#sigma_{Y}",2e1,-355./113,355./113,1e2,-10,10);
 TH2 *hPullPointXtheta = new TH2F("hPullPointXtheta","(X_{MC}-X_{rec})/#sigma_{X} vs. #theta;#theta_{MC}, rad;(X_{MC}-X_{rec})/#sigma_{X}",2e1,0.003, 0.009,1e2,-10,10);
TH2 *hPullPointYtheta = new TH2F("hPullPointYtheta","(Y_{MC}-Y_{rec})/#sigma_{Y} vs. #theta;#theta_{MC}, rad;(Y_{MC}-Y_{rec})/#sigma_{Y}",2e1,0.003, 0.009,1e2,-10,10);
  TH1 *hResPointPx = new TH1F("hResPointPx","Px_{MC}-Px_{rec};#deltaPx, GeV/c",1e2,-0.01,0.01);
  TH1 *hErrPointPx = new TH1F("hErrPointPx","#sigma_{Px};#sigmaPx, GeV/c",1e3,0,0.01);
  TH1 *hPullPointPx = new TH1F("hPullPointPx","(Px_{MC}-Px_{rec})/#sigma_{Px};(Px_{MC}-Px_{rec})/#sigma_{Px}",1e2,-10,10);

  TH1 *hResPointPy = new TH1F("hResPointPy","Py_{MC}-Py_{rec};#deltaPy, GeV/c",1e2,-0.01,0.01);
  TH1 *hErrPointPy = new TH1F("hErrPointPy","#sigma_{Py};#sigmaPy, GeV/c",1e3,0,0.01);
  TH1 *hPullPointPy = new TH1F("hPullPointPy","(Py_{MC}-Py_{rec})/#sigma_{Py};(Py_{MC}-Py_{rec})/#sigma_{Py}",1e2,-10,10);

  TH1 *hResPointPz = new TH1F("hResPointPz","Pz_{MC}-Pz_{rec};#deltaPz, GeV/c",1e2,-1e-3,1e-3);
  TH1 *hErrPointPz = new TH1F("hErrPointPz","#sigma_{Pz};#sigmaPz, GeV/c",1e3,0,1e-1);
  TH1 *hPullPointPz = new TH1F("hPullPointPz","(Pz_{MC}-Pz_{rec})/#sigma_{Pz};(Pz_{MC}-Pz_{rec})/#sigma_{Pz}",1e2,-10,10);

  TH1 *hResDCA = new TH1F("hDCA","|DCA|;|DCA|,cm",1e3,0,10.);
  TH1 *hErrDCA = new TH1F("hErrDCA","#sigma_{DCA};#sigmaDCA,cm",1e3,0,1.);
  TH1 *hPullDCA = new TH1F("hPullDCA","|DCA|/#sigma_{DCA};",1e2,0,10);

  TH2 *hPointXY = new TH2F("hPointXY","PCA;X,cm;Y,cm",1e3,-10,10.,1e3,-10,10.);
  TH2 *hPointXYcut = new TH2F("hPointXYcut","PCA;#X,cm;#Y,cm",1e3,-1.,1.,1e3,-1.,1.);
  TH2 *hErrPointXY = new TH2F("hErrPointXY","#sigma_{PCA};#sigmaX,cm;#sigmaY,cm",1e3,0,30.,1e3,0,30.);

  TH2 *hHitIDHitX = new TH2F("hHitIDHitX",";MC_{ID} ;hit_{X}, cm;",6,-1,5,1e3,24.,45.);
  TH2 *hHitIDHitY = new TH2F("hHitIDHitY",";MC_{ID} ;hit_{Y}, cm;",6,-1,5,1e3,-10.,10.);

  TH2 *hRecGEANEX = new TH2F("hRecGEANEX",";X_{rec}, cm;X_{GEANE}, cm",1e3,-15.,15.,1e3,-100,100.);
  TH2 *hRecGEANEY = new TH2F("hRecGEANEY",";Y_{rec}, cm;Y_{GEANE}, cm",1e3,-15.,15.,1e3,-100,100.);
  TH2 *hRecGEANEZ = new TH2F("hRecGEANEZ",";Z_{rec}, cm;Z_{GEANE}, cm",1e3,-15.,15.,1e3,-1.,1.);
  TH2 *hRecGEANER = new TH2F("hRecGEANER",";R_{rec}, cm;R_{GEANE}, cm",1e3,0,10.,1e3,0,10.);
  TH2 *hRecGEANETheta = new TH2F("hRecGEANETheta",";#theta_{rec}, rad;#theta_{GEANE}, rad",1e3,0,1,1e3,0,1);
  TH2 *hRecGEANEPhi = new TH2F("hRecGEANEPhi",";#phi_{rec}, rad;#phi_{GEANE}, rad",1e3,-7,7,1e3,-7,7);
  TH2 *hRecThetaPhi = new TH2F("hRecThetaPhi",";#theta_{rec}, rad;#phi_{rec}, rad",1e3,0,1.,1e3,-7,7);

  TH2 *hSeedGEANEX = new TH2F("hSeedGEANEX",";X_{cand}, cm;X_{GEANE}, cm",1e3,-100.,100.,1e3,-100.,100.);
  TH2 *hSeedGEANEY = new TH2F("hSeedGEANEY",";Y_{cand}, cm;Y_{GEANE}, cm",1e3,-100.,100.,1e3,-100.,100.);
  TH2 *hSeedGEANEZ = new TH2F("hSeedGEANEZ",";Z_{cand}, cm;Z_{GEANE}, cm",1e3,-0.015,0.015,1e3,-0.05,0.05);
  TH2 *hSeedGEANER = new TH2F("hSeedGEANER",";R_{cand}, cm;R_{GEANE}, cm",1e3,0,100.,1e3,0,100.);
  TH2 *hSeedGEANETheta = new TH2F("hSeedGEANETheta",";#theta_{cand}, rad;#theta_{GEANE}, rad",1e3,0,0.1,1e3,0,0.1);
  TH2 *hSeedGEANEPhi = new TH2F("hSeedGEANEPhi",";#phi_{cand}, rad;#phi_{GEANE}, rad",1e3,-7,7,1e3,-7,7);
  TH2 *hSeedThetaPhi = new TH2F("hSeedThetaPhi",";#theta_{cand}, rad;#phi_{cand}, rad",1e3,0,1.,1e3,-1.,1.);

  TH2 *hSeedThetaDiffIDs = new TH2F("hSeedThetaDiffIDs",";# IDs;#theta_{cand}, rad",7,0,7,1e3,0,0.1);
  TH2 *hSeedPhiDiffIDs = new TH2F("hSeedPhiDiffIDs",";# IDs;#phi_{cand}, rad",7,0,7,1e3,-4.,4.);
  TH2 *hDiffIDsPointX = new TH2F("hDiffIDsPointX",";# IDs;X,cm",7,0,7,1e3,-10,10.);
  TH2 *hDiffIDsPointY = new TH2F("hDiffIDsPointY",";# IDs;Y,cm",7,0,7,1e3,-10,10.);
  TH2 *hDiffIDsPointZ = new TH2F("hDiffIDsPointZ",";# IDs;Z,cm",7,0,7,1e3,-1.,1.);

  TH2 *hMCThetaGEANETheta = new TH2F("hMCThetaGEANETheta",";#theta_{MC}, rad;#theta_{GEANE}, rad",1e3,0,0.012,1e3,0,0.012);
  TH2 *hMCPhiGEANEPhi = new TH2F("hMCPhiGEANEPhi",";#phi_{MC}, rad;#phi_{GEANE}, rad",1e3,-3.15,3.15,1e3,-3.15,3.15);

  TH2 *hMCThetaResTheta = new TH2F("hMCThetaResTheta",";#theta_{MC}, rad;#Delta#theta, rad",1e2,0,0.012,1e3,-1e-2,1e-2);
  TH2 *hMCPhiResPhi = new TH2F("hMCPhiResPhi",";#phi_{MC}, rad;#Delta#phi, rad",1e2,-3.15,3.15,1e2,-5e-1,5e-1);
  TH2 *hMCThetaResSeedTheta = new TH2F("hMCThetaResSeedTheta",";#theta_{MC}, rad;#Delta#theta, rad",1e2,0.,0.1,1e3,-1e-2,1e-2);
  TH2 *hMCPhiResSeedPhi = new TH2F("hMCPhiResSeedPhi",";#phi_{MC}, rad;#Delta#phi, rad",1e2,-3.15,3.15,1e2,-1e-1,1e-1);
  TH2 *hResZResPhi = new TH2F("hResZResPhi",";#Delta#phi, rad;#Delta Z, cm",1e3,-2e-1,2e-1,1e3,-0.02,0.02);
  TH2 *hntrkMCtrkID = new TH2F("hntrkMCtrkID","Number of participation this MCid in rec.tracks;ID;Number of rec.trks",10,0,10,10,0,10);
  TH2 *hchi2MCdiffID = new TH2F("hchi2MCdiffID","number diff MCid in rec.tracks vs. #chi^{2};#chi^{2};N",5e2,0,50,10,0,10);

  TH1 *hLumiTrkA = new TH1F("hLumiTrkA","param a [x=a*z+b];",1e3,-0.1,0.1);
  TH1 *hLumiTrkB = new TH1F("hLumiTrkB","param b [x=a*z+b];",1e3,-10.,10.);
  TH1 *hLumiTrkC = new TH1F("hLumiTrkC","param c [y=c*z+d];",1e3,-0.1,0.1);
  TH1 *hLumiTrkD = new TH1F("hLumiTrkD","param d [y=c*z+d];",1e3,-10.,10.);

  //-----------------------------------------------------------------------------

  TH2 *hchi2Errx = new TH2F("hchi2Errx","#sigma_{X} vs. #chi^2 for reconstructed tracks;#chi^2;#sigma_{X}, cm",
                            1e3,0,100,2e1,0,0.2);
  TH2 *hchi2Erry = new TH2F("hchi2Erry","#sigma_{Y} vs. #chi^2 for reconstructed tracks;#chi^2;#sigma_{X}, cm",
                            1e3,0,100,2e1,0,0.2);
 TH2 *hchi2nTrkCand = new TH2F("hchi2nTrkCand"," ;Number of trk-cand;#chi^2",
                            30,0,30,5e2,0,50.);
  TNtuple *nrecpointall = new TNtuple("nrecpointall","recpointAll","xrecbp:yrecbp:zrecbp:xrec:yrec:zrec:xseed:yseed:zseed:chi2");
  TNtuple *nrecdirall = new TNtuple("nrecdirall","recdirAll","dxmc:dymc:dzmc:dxrecbp:dyrecbp:dzrecbp:dirxrec:diryrec:dirzrec:dirxseed:diryseed:dirzseed");
  TNtuple *nsectors = new TNtuple("nsectors","sectors","thetares:sector");
  TH2I *hnhits = new TH2I("hnhits","# rec hits vs. # sim hits; sim; rec",100,0,100,100,0,100);
  //Load lumi geo params
  PndLmdDim *lmddim = PndLmdDim::Instance();
  // lmddim -> Read_transformation_matrices("matrices.txt", true);
  lmddim -> Read_transformation_matrices("matrices_perfect.txt", false);
  int glBADGEANE=0;
  int glBadEv = 0;
  int glNoisehit = 0;// total number of noise hits
  for (Int_t j=0; j<nEvents; j++){
    //  cout<<"Event #"<<j<<endl;
    // Read GEANE & MC info -----------------------------------------------------------------
    //    if(kf<1)
    tTrkRec.GetEntry(j);
    tgeane.GetEntry(j);
    tMC.GetEntry(j);
    tTrkCand.GetEntry(j);
    tHits.GetEntry(j);
    tdigiHits.GetEntry(j);
    if(mh>0)
      tHitsMerged.GetEntry(j);
    const int nGeaneTrks = geaneArray->GetEntriesFast();
    const int nParticles = true_tracks->GetEntriesFast();
    const int numTrk = nGeaneTrks;
    const int nRecHits = rechit_array->GetEntriesFast();
    const int nMCHits = true_points->GetEntriesFast();

    const int nTrkCandidates = trkcand_array->GetEntriesFast();
    const int nRecTrks = rec_trk->GetEntriesFast();
    if(verboseLevel>0)
      cout<<"Event #"<<j<<" has "<<nParticles<<" true particles, "<<" out of it "<<nRecHits<<" hits, "<<nTrkCandidates
          <<" trk-cands, "<<numTrk<<" tracks and "<<nGeaneTrks<<" geane Trks!"<<endl;
    hnRecnMC->Fill(nParticles,nGeaneTrks);
    if(nParticles!=nMCtracks) continue;
    hnhits->Fill(nMCHits,nRecHits);
    if(nRecHits<3*nMCtracks) glBadEv++;
    // if(nRecHits<3*nMCtracks) cout<<"Event #"<<j<<" doesn't have enough rec.hits!!!"<<endl;
    // if(nRecHits<3*nMCtracks) continue;
    double chi2Cont[5*numTrk];
    double ndiffIDCont[5*numTrk];

    hMCtrkvshits->Fill(nRecHits,nParticles);
    if(nTrkCandidates>numTrk) cout<<"Event #"<<j<<" has "<<nTrkCandidates<<" trk-cands and "<<numTrk<<" tracks!"<<endl;

    hntrkcand->Fill(nTrkCandidates);
    hntrkcandvsMC->Fill(nParticles,nTrkCandidates);

    Int_t nRecGEANEtrk = 0;
    int MCtrk[nParticles]; //Number of participation this MCid in rec.tracks
    int RECtrkMCid[nGeaneTrks];//Assignment MC id to REC trk;
    for(int nk=0;nk<nParticles;nk++)
      MCtrk[nk]=0;
    bool goodTrk[nGeaneTrks];
    bool ghostTrk[nGeaneTrks];
    for (Int_t iN=0; iN<nGeaneTrks; iN++){
      goodTrk[iN] = false;
      ghostTrk[iN] = false;
      RECtrkMCid[iN]=-1;
    }

    int goodRectrk=0;//for missed trk-search
    for (Int_t iN=0; iN<nGeaneTrks; iN++){// loop over all reconstructed trks


      vector<int> MCtrkID; //arrray of hits MCid
      Int_t diffIDs=1;
      FairTrackParH *fRes = (FairTrackParH*)geaneArray->At(iN);
      TVector3 PosRec = fRes->GetPosition();
      double pca_lim = 1.;//=10*sigma_Xpca~10*{0.093,0.11,0.12,0.22,0.55};
      if(Plab<5) pca_lim = 2.;
      if(Plab<2) pca_lim = 5.;
      // if(fabs(PosRec.X())>pca_lim && fabs(PosRec.Y())>pca_lim){
      //        cout<<"666 Event #"<<j<<" has too large X_pca ot Y_pca 666"<<endl;
      // }
      if(fabs(PosRec.X())>pca_lim && fabs(PosRec.Y())>pca_lim) continue; // PCA_x and PCA_y should be < 10sigmaX


      Double_t lyambda = fRes->GetLambda();
      if(lyambda==0){
        cout<<"GEANE didn't propagate this trk!"<<endl;
        cout<<"Event #"<<j<<" diffIDs = "<<diffIDs<<endl;
        glBADGEANE++;
      }
      if(lyambda==0) continue;
      PndLinTrack *trk;
      double linpar[6];
      Double_t errparlin[6];
      double chi2;
      PndTrack *trkpnd;
      Int_t candID;
      PndTrackCand *trkcand;
      Int_t trkSensor;
        trk = (PndLinTrack*)rec_trk->At(iN);
        trk->GetPar(linpar);
        hLumiTrkA->Fill(linpar[1]);
        hLumiTrkB->Fill(linpar[0]);
        hLumiTrkC->Fill(linpar[3]);
        hLumiTrkD->Fill(linpar[2]);
        trk->GetParErr(errparlin);
        candID = trk->GetTCandID();
        trkcand = (PndTrackCand*)trkcand_array->At(candID);
        chi2 = trk->GetChiSquare();
        chi2Cont[iN] = chi2;
        hchi2nTrkCand->Fill(nTrkCandidates,chi2);

        // PndTrackCand *trkcand = (PndTrackCand*)trkcand_array->At(iN); //TODO: find out how to check trk-cand in case then cut are applied to trk-fit results
        const int Ntrkcandhits= trkcand->GetNHits();
        PndSdsMCPoint* MCPointHit;

        //Matching between MC & Rec on hits level-----------------------------------
        if(verboseLevel>1) cout<<"MCidTOP:"<<endl;
        double momMC0,momMC1,momMC2;
        momMC0 = Plab;
        for (Int_t iHit = 0; iHit < Ntrkcandhits; iHit++){ // loop over rec.hits
          PndTrackCandHit candhit = (PndTrackCandHit)(trkcand->GetSortedHit(iHit));
        Int_t hitID = candhit.GetHitId();
        PndSdsHit* myHit = (PndSdsHit*)(rechit_array->At(hitID));

        //for pixel design
        double MCpointMom;
        if(dnu>0){
          //  if(verboseLevel>1) cout<<"Rec hit("<<myHit->GetClusterIndex()<<")";
          PndSdsClusterPixel* myCluster = (PndSdsClusterPixel*)(fStripClusterArray->At(myHit->GetClusterIndex()));
          PndSdsDigiPixel* astripdigi = (PndSdsDigiPixel*)(fStripDigiArray->At(myCluster->GetDigiIndex(0)));
          if (astripdigi->GetIndex(0) == -1){
            glNoisehit++;
            //  MCtrkID.push_back(-111);
            continue;
          }
        int sensorID = myHit->GetSensorID();
        int ihalf, iplane, imodule, iside, idie, isensor;
        // calculate the plane and sensor on this plane
        lmddim->Get_sensor_by_id(sensorID, ihalf, iplane, imodule, iside, idie, isensor);
        trkSensor = ihalf*5+imodule;
          // if (astripdigi->GetIndex(0) == -1) continue; // sort out noise
          PndSdsMCPoint* MCPoint = (PndSdsMCPoint*)(true_points->At(astripdigi->GetIndex(0)));
          MCpointMom = sqrt(MCPoint->GetPx()*MCPoint->GetPx()+MCPoint->GetPy()*MCPoint->GetPy()+MCPoint->GetPz()*MCPoint->GetPz());
          int MCidTOP = MCPoint->GetTrackID();
          if(iHit==0) MCPointHit = MCPoint;
          MCtrkID.push_back(MCidTOP);
          //      if(verboseLevel>1) cout<<"MCid("<<MCidTOP<<") ";
          if(verboseLevel>1) cout<<MCidTOP<<" ";
        }

        //for strip design
        else{
        PndSdsClusterStrip* myCluster =  (PndSdsClusterStrip*)(fStripClusterArray->At(myHit->GetClusterIndex()));
        PndSdsDigiStrip* astripdigi = (PndSdsDigiStrip*)(fStripDigiArray->At(myCluster->GetDigiIndex(0)));
        if (astripdigi->GetIndex(0) == -1) glNoisehit++;
        if (astripdigi->GetIndex(0) == -1) continue; // sort out noise
        PndSdsMCPoint* MCPoint = (PndSdsMCPoint*)(true_points->At(astripdigi->GetIndex(0)));
        MCpointMom = sqrt(MCPoint->GetPx()*MCPoint->GetPx()+MCPoint->GetPy()*MCPoint->GetPy()+MCPoint->GetPz()*MCPoint->GetPz());

        int MCidTOP = MCPoint->GetTrackID();
        if(verboseLevel>1) cout<<MCidTOP<<" ";
        if(iHit==0) MCPointHit = MCPoint;
        // double zTop = MCPoint->GetZ();
        // double xTop = MCPoint->GetX();
        // double yTop = MCPoint->GetY();
        Int_t  botIndex = myHit->GetBotIndex();
        PndSdsClusterStrip* myClusterBot =  (PndSdsClusterStrip*)(fStripClusterArray->At(botIndex));
        PndSdsDigiStrip* astripdigiBot = (PndSdsDigiStrip*)(fStripDigiArray->At(myClusterBot->GetDigiIndex(0)));
        if (astripdigiBot->GetIndex(0) == -1) glNoisehit++;
        if (astripdigiBot->GetIndex(0) == -1) continue; // sort out noise
        PndSdsMCPoint* MCPointBot = (PndSdsMCPoint*)(true_points->At(astripdigiBot->GetIndex(0)));
        // double zBot = MCPointBot->GetZ();
        // double xBot = MCPointBot->GetX();
        // double yBot = MCPointBot->GetY();
        // cout<<"zTop - zBot = "<<zTop-zBot<<endl;
        // cout<<"xTop - xBot = "<<xTop-xBot<<endl;
        // cout<<"yTop - yBot = "<<yTop-yBot<<endl;
        int MCidBOT = MCPointBot->GetTrackID();
        //      if(MCidTOP!=MCidBOT) continue;
        MCtrkID.push_back(MCidTOP);
        MCtrkID.push_back(MCidBOT);
        }
        if(iHit==0) momMC1 = MCpointMom;
        if(iHit==(Ntrkcandhits-1)) momMC2 = MCpointMom;
      }

      // TNtuple *nmomMC = new TNtuple("nmomMC","MCmomentum","p0:p_nearLUMI:p_afterLUMI");
      nmomMC->Fill(momMC0,momMC1,momMC2);
      if(verboseLevel>1) cout<<""<<endl;
      //   Sorting MC IDs ----------------------------------------
      Int_t k, x;
      bool ch=false; //Was element changed?
      Int_t nch = 0; //How many times?
      for(Int_t n=0; n<MCtrkID.size(); n++) { // n - current position
        k=n; x=MCtrkID[n];
        for(Int_t m=n+1; m<MCtrkID.size(); m++) // find the least element
          if (MCtrkID[m]<x){
            k=m; x=MCtrkID[m];          // k - index for the least element
            ch=true; nch++;
          }
        MCtrkID[k] = MCtrkID[n]; MCtrkID[n] = x; // change position between the least and current elements
      }

      Int_t prevID = MCtrkID[0];
      for(int nk=0;nk<nParticles;nk++){
        if(MCtrkID[0]==nk)
          MCtrk[nk]=MCtrk[nk]+1;
      }
      diffIDs = 1;
      for(Int_t n=1; n<MCtrkID.size(); n++){
        if(prevID<MCtrkID[n]){
          diffIDs++;
          prevID=MCtrkID[n];
          for(int nk=0;nk<nParticles;nk++){
            if(MCtrkID[n]==nk)
              MCtrk[nk]=MCtrk[nk]+1;
          }
        }
      }

      ndiffIDCont[iN]=diffIDs;
      hDiffIDs->Fill(diffIDs);
      hntrkcandvsIDs->Fill(diffIDs,nTrkCandidates);
      hnhitsvsIDs->Fill(diffIDs,Ntrkcandhits);

      if(diffIDs<2){
        RECtrkMCid[iN] = MCtrkID[0];
        goodTrk[iN] = true;
      }
      else{
        vector<int> countMC_IDs(diffIDs);
        prevID = MCtrkID[0];
        int diffCount=0;
        for(Int_t n=0; n<MCtrkID.size(); n++) {
          countMC_IDs[diffCount]++;
          if(prevID<MCtrkID[n]){
            diffCount++;
            prevID=MCtrkID[n];
          }
        }
        int maxID=countMC_IDs[0];
        int posIDmax=0;
        for(int kn=0;kn<diffIDs;kn++){
          //      cout<<"countMC_IDs["<<kn<<"]="<<countMC_IDs[kn]<<" 0.7*MCtrkID.size() = "<<0.7*MCtrkID.size()<<endl;
          if(countMC_IDs[kn]>0.65*MCtrkID.size()){ //more then 65% of hits come from the same MC id
            goodTrk[iN] = true;
            ghostTrk[iN] = false;
          }
          else{
            if(!goodTrk[iN]) ghostTrk[iN] = true;
          }

          if(countMC_IDs[kn]>maxID){
            maxID=countMC_IDs[kn];
            posIDmax = kn;
          }
        }
        prevID = MCtrkID[0];
        diffCount=0;
        for(Int_t n=0; n<MCtrkID.size(); n++) {
          if(diffCount==posIDmax) RECtrkMCid[iN] = prevID;
          if(prevID<MCtrkID[n]){
            diffCount++;
            prevID=MCtrkID[n];
          }
        }
      }

      TVector3 posSeed = trkcand->getPosSeed();
      TVector3 dirSeed = trkcand->getDirSeed();
      hSeedThetaDiffIDs->Fill(diffIDs,dirSeed.Theta());
      hSeedPhiDiffIDs->Fill(diffIDs,dirSeed.Phi());
      if(diffIDs>-1){ // All tracks
        TVector3 pos_prop_geane_trk = fRes->GetPosition();
        Double_t theta_prop_geane_trk = TMath::Pi()/2. - lyambda;
        Double_t phi_prop_geane_trk = fRes->GetPhi();
        TVector3 pos_rec_trk,dir_rec_trk;

          pos_rec_trk = trk->GetStartVec();
          dir_rec_trk = trk->GetDirectionVec();
          hRecGEANEX->Fill(pos_rec_trk.X(),pos_prop_geane_trk.X());
          hRecGEANEY->Fill(pos_rec_trk.Y(),pos_prop_geane_trk.Y());
          hRecGEANEZ->Fill(pos_rec_trk.Z(),pos_prop_geane_trk.Z());
          hRecGEANER->Fill(pos_rec_trk.Perp(),pos_prop_geane_trk.Perp());
          hRecGEANETheta->Fill(dir_rec_trk.Theta(), theta_prop_geane_trk);
          hRecGEANEPhi->Fill(dir_rec_trk.Phi(), phi_prop_geane_trk);
          hRecThetaPhi->Fill(dir_rec_trk.Theta(),dir_rec_trk.Phi());

          hSeedGEANEX->Fill(posSeed.X(),pos_prop_geane_trk.X());
          hSeedGEANEY->Fill(posSeed.Y(),pos_prop_geane_trk.Y());
          hSeedGEANEZ->Fill(posSeed.Z(),pos_prop_geane_trk.Z());
          hSeedGEANER->Fill(posSeed.Perp(),pos_prop_geane_trk.Perp());
          hSeedGEANETheta->Fill(dirSeed.Theta(), theta_prop_geane_trk);
          hSeedGEANEPhi->Fill(dirSeed.Phi(), phi_prop_geane_trk);
          hSeedThetaPhi->Fill(dirSeed.Theta(),dir_rec_trk.Phi());



        hDiffIDsPointX->Fill(diffIDs,PosRec.X());
        hDiffIDsPointY->Fill(diffIDs,PosRec.Y());
        hDiffIDsPointZ->Fill(diffIDs,PosRec.Z());
        Double_t CovGEANELAB[6][6];
        fRes->GetMARSCov(CovGEANELAB);

        Double_t errX = fRes->GetDX();
        Double_t errY = fRes->GetDY();
        Double_t errZ = fRes->GetDZ();

        Double_t errXlin = errparlin[0];
        Double_t errYlin = errparlin[1];
        hchi2Errx->Fill(chi2,errX);
        hchi2Erry->Fill(chi2,errY);
        nRecGEANEtrk++;

        int MCidforREC = RECtrkMCid[iN];
        PndMCTrack *mctrk =(PndMCTrack*) true_tracks->At(MCidforREC);
        Int_t mcID = mctrk->GetPdgCode();
        TVector3 MomMC = mctrk->GetMomentum();
        Double_t thetaMC = MomMC.Theta();
        Double_t phiMC = MomMC.Phi();

        //TODO: problem with covarance matrix in FairTrackParH???

        Double_t errPx = fRes->GetDPx();
        Double_t errPy = fRes->GetDPy();
        Double_t errPz = fRes->GetDPz();
        TVector3 errMomBP(errPx,errPy,errPz);

        Double_t thetaBP = TMath::Pi()/2. - lyambda;
        Double_t err_lyambda = fRes->GetDLambda();
        if(err_lyambda==0) err_lyambda = errMomBP.Theta();
        // Double_t err_lyambda = errMom.Theta();
        Double_t phiBP = fRes->GetPhi();
        Double_t err_phi = fRes->GetDPhi();
        if(err_phi==0) err_phi=errMomBP.Phi();
        Double_t errMomRecBP = fRes->GetDQp();


        TVector3 MomRecBP = fRes->GetMomentum();



        Double_t resTheta = thetaBP-thetaMC;

        double resMom = MomRecBP.Mag()-MomMC.Mag();
        hgoodTheta->Fill(thetaMC);
        hgoodPhi->Fill(phiMC);
        hgoodPhiTheta->Fill(thetaMC,phiMC);
        hgoodThetachi2->Fill(chi2,thetaMC);
        hgoodPhichi2->Fill(chi2,phiMC);

        hMCThetaGEANETheta->Fill(thetaMC,thetaBP);
        hMCPhiGEANEPhi->Fill(phiMC,phiBP);
        hMCThetaResTheta->Fill(thetaMC,thetaMC-thetaBP);
        hMCPhiResPhi->Fill(phiMC,phiMC-phiBP);
        hResPointX->Fill(PosRec.X());
        hResPointY->Fill(PosRec.Y());
        hResDCA->Fill(PosRec.Mag());
        Double_t errPCA = TMath::Sqrt(errX*errX+errY*errY+errZ*errZ);
        hErrDCA->Fill(errPCA);
        hPullDCA->Fill(PosRec.Mag()/errPCA);
        hPointXY->Fill(-PosRec.X(), -PosRec.Y());
        hPointXYcut->Fill(-PosRec.X(), -PosRec.Y());
        hResPointZ->Fill(-PosRec.Z());
        hErrPointX->Fill(errX);
        hErrPointY->Fill(errY);
        hErrLinPointX->Fill(errXlin);
        hErrLinPointY->Fill(errYlin);
        hErrPointXY->Fill(errX,errY);
        hErrPointZ->Fill(errZ);
        hPullPointX->Fill(-PosRec.X()/errX);
        hPullPointY->Fill(-PosRec.Y()/errY);
        hPullPointZ->Fill(-PosRec.Z()/errZ);
        hPullPointXphi->Fill(phiMC,(-PosRec.X()/errX));
        hPullPointYphi->Fill(phiMC,(-PosRec.Y()/errY));
        hPullPointXtheta->Fill(thetaMC,(-PosRec.X()/errX));
        hPullPointYtheta->Fill(thetaMC,(-PosRec.Y()/errY));

        hResPointPx->Fill(MomMC.X()-MomRecBP.X());
        hResPointPy->Fill(MomMC.Y()-MomRecBP.Y());
        hResPointPz->Fill(MomMC.Z()-MomRecBP.Z());
        hErrPointPx->Fill(errPx);
        hPullPointPx->Fill((MomMC.X()-MomRecBP.X())/errPx);
        hErrPointPy->Fill(errPy);
        hPullPointPy->Fill((MomMC.Y()-MomRecBP.Y())/errPy);
        hErrPointPz->Fill(errPz);
        hPullPointPz->Fill((MomMC.Z()-MomRecBP.Z())/errPz);

        Double_t resPhi = phiBP-phiMC;
        hResMom->Fill(resMom);
        hErrMom->Fill(errMomRecBP);
        hPullMom->Fill(resMom/errMomRecBP);
        hResTheta->Fill(resTheta);
        nsectors->Fill(resTheta,trkSensor);
        hErrTheta->Fill(err_lyambda);
        hPullTheta->Fill(resTheta/err_lyambda);
        hResPhi->Fill(resPhi);
        hErrPhi->Fill(err_phi);
        hPullPhi->Fill(resPhi/err_phi);
        hResZResPhi->Fill(resPhi,PosRec.Z());

        Double_t angMCRec = MomMC.Angle(MomRecBP); // Angle between two vectors
        hangMCRec->Fill(angMCRec);
        hchi2->Fill(chi2);


        for(Int_t jk=0; jk< nParticles;jk++){
          PndMCTrack *mctrk =(PndMCTrack*) true_tracks->At(jk);
          Int_t mcID = mctrk->GetPdgCode();
          if(mcID==-2212){
            TVector3 MomMC = mctrk->GetMomentum();
            double diffTheta=fabs(MomMC.Theta()-thetaBP)/err_lyambda;
            double diffPhi=fabs(MomMC.Phi()-phiBP)/err_phi;

            if(fabs(MomMC.Phi())<0.14 || fabs(MomMC.Phi())>3){ //for #phi near edge reconstructed angle can differ
              double phiMC1 = -MomMC.Phi();
              double diffPhi2=fabs(phiMC1-phiBP)/err_phi;
              if(diffPhi2<diffPhi) diffPhi = diffPhi2;
            }

            hntrkmissedPhiTheta->Fill(diffTheta,diffPhi);
            if(diffTheta<4. && diffPhi<4.) goodRectrk++;
          }
        }
        nrecpointall->Fill(pos_prop_geane_trk.X(),pos_prop_geane_trk.Y(),pos_prop_geane_trk.Z(),
                           pos_rec_trk.X(),pos_rec_trk.Y(),pos_rec_trk.Z(),
                           posSeed.X(),posSeed.Y(),posSeed.Z(),chi2);
        nrecdirall->Fill(MomMC.X()/MomMC.Mag(),MomMC.Y()/MomMC.Mag(),MomMC.Z()/MomMC.Mag(),MomRecBP.X()/MomRecBP.Mag(),MomRecBP.Y()/MomRecBP.Mag(),MomRecBP.Z()/MomRecBP.Mag(),
                         dir_rec_trk.X(),dir_rec_trk.Y(),dir_rec_trk.Z(),
                         dirSeed.X(),dirSeed.Y(),dirSeed.Z());


        // cout<<" "<<endl;
        // cout<<"=============================="<<endl;
        TVector3 vtxLumi = trk->GetStartVec();
        TVector3 vtxLumiErr = trk->GetStartErrVec();
        TVector3 dirLumi = trk->GetDirectionVec();
        TVector3 dirLumiErr = trk->GetDirectionErrVec();
        // // //do the transformation from LUMI frame (with z-axis perp. to lumi planes) to lab frame
        // if(dnu>0){
        //   vtxLumi = lmddim->Transform_lmd_local_to_global(vtxLumi, false, false);
        //   dirLumi = lmddim->Transform_lmd_local_to_global(dirLumi, true, false);
        //   vtxLumiErr = lmddim->Transform_lmd_local_to_global(vtxLumiErr, true, false);
        //   dirLumiErr = lmddim->Transform_lmd_local_to_global(dirLumiErr, true, false);
        // }
        // else{
        //   combitransFromLumiFrame(vtxLumi,dnu);
        //   rotateFromLumiFrame(dirLumi, false, dnu);
        //   rotateFromLumiFrame(dirLumiErr, true, dnu);
        //   rotateFromLumiFrame(vtxLumiErr, true, dnu);
        // }
        double xTrue = MCPointHit->GetX();
        double yTrue = MCPointHit->GetY();
        double zTrue = MCPointHit->GetZ();
        TVector3 vtxLumiMC = TVector3(xTrue,yTrue,zTrue);


        double pxTrue =  MCPointHit->GetPx();
        double pyTrue =  MCPointHit->GetPy();
        double pzTrue =  MCPointHit->GetPz();
        TVector3 dirLumiMC = TVector3(pxTrue,pyTrue,pzTrue);
        //      hMomMC->Fill(dirLumiMC.Mag());
        dirLumiMC *=1./dirLumiMC.Mag();

        double dz = -zTrue+vtxLumi.Z();//Correct definition Z coord for comparision
        //      double dz = 0;//TEST
        double dx = dirLumiMC.X()*dz;
        double dy = dirLumiMC.Y()*dz;
        vtxLumiMC += TVector3(dx,dy,dz);

        hResLumiTrkMom->Fill(dirLumiMC.Mag()-dirLumi.Mag());
        hResLumiTrkTheta->Fill(dirLumiMC.Theta()-dirLumi.Theta());
        hResLumiTrkPhi->Fill(dirLumiMC.Phi()-dirLumi.Phi());
        hResLumiTrkThetaPull->Fill((dirLumiMC.Theta()-dirLumi.Theta())/dirLumiErr.Theta());
        //      hResLumiTrkPhiPull->Fill(dirLumiMC.Phi()-dirLumi.Phi());
        hResLumiTrkPhiPull->Fill((dirLumiMC.Phi()-dirLumi.Phi())/dirLumiErr.Phi());

        hResLumiTrkPointX->Fill(vtxLumiMC.X()-vtxLumi.X());
        hResLumiTrkPointY->Fill(vtxLumiMC.Y()-vtxLumi.Y());
        hResLumiTrkPointZ->Fill(vtxLumiMC.Z()-vtxLumi.Z());

        hResLumiTrkPointPx->Fill(dirLumiMC.X()-dirLumi.X());
        hResLumiTrkPointPy->Fill(dirLumiMC.Y()-dirLumi.Y());
        hResLumiTrkPointPz->Fill(dirLumiMC.Z()-dirLumi.Z());

        hResLumiTrkPointXPull->Fill((vtxLumiMC.X()-vtxLumi.X())/vtxLumiErr.X());
        hResLumiTrkPointYPull->Fill((vtxLumiMC.Y()-vtxLumi.Y())/vtxLumiErr.Y());
        hResLumiTrkPointZPull->Fill((vtxLumiMC.Z()-vtxLumi.Z())/vtxLumiErr.Z());

        hResLumiTrkPointPxPull->Fill((dirLumiMC.X()-dirLumi.X())/dirLumiErr.X());
        hResLumiTrkPointPyPull->Fill((dirLumiMC.Y()-dirLumi.Y())/dirLumiErr.Y());
        hResLumiTrkPointPzPull->Fill((dirLumiMC.Z()-dirLumi.Z())/dirLumiErr.Z());

        hMCThetaResSeedTheta->Fill(dirLumiMC.Theta(),(dirLumiMC.Theta()-dirLumi.Theta()));
        hMCPhiResSeedPhi->Fill(dirLumiMC.Phi(),(dirLumiMC.Phi()-dirLumi.Phi()));
      }
      MCtrkID.clear();
    }
    for(int nk=0;nk<numTrk;nk++){
      hntrkMCtrkID->Fill(nk,MCtrk[nk]);
      hchi2MCdiffID->Fill(chi2Cont[nk],ndiffIDCont[nk]);
    }
    hntrkgood_I->Fill(goodRectrk);
    if((nMCtracks-goodRectrk)>0) hntrkmissed_I->Fill(nMCtracks-goodRectrk);
    if((nGeaneTrks-goodRectrk)>0) hntrkghost_I->Fill(nGeaneTrks-goodRectrk);
    if(verboseLevel>0){
    cout<<nMCtracks<<" trks per event were simulated and "<<nGeaneTrks<<" were reconstructed"<<endl;
    cout<<"-- (I) dPhi\dTheta -------------------------"<<endl;
    cout<<"Good trks: "<<goodRectrk<<" missed: "<<nMCtracks-goodRectrk<<" ghost: "<<nGeaneTrks-goodRectrk<<endl;
    cout<<"--------------------------------"<<endl;
    }
    int goodRecII=0, ghostRecII=0;
    for (Int_t iN=0; iN<nGeaneTrks; iN++){
      int trkType = -1;
      if(goodTrk[iN]){
        goodRecII++;
        trkType = 0;
      }
      if(ghostTrk[iN]){
        ghostRecII++;
        trkType = +1;
      }

      FairTrackParH *fRes = (FairTrackParH*)geaneArray->At(iN);
      Double_t lyambda = fRes->GetLambda();
      Double_t thetaBP = TMath::Pi()/2. - lyambda;
      Double_t phiBP = fRes->GetPhi();
      TVector3 MomRecBP = fRes->GetMomentum();
      TVector3 PosBP = fRes->GetPosition();
      ntuprecTrk->Fill(PosBP.X(),PosBP.Y(),PosBP.Z(),MomRecBP.Mag(),thetaBP,phiBP,trkType);
    }
    hntrkgood_II->Fill(goodRecII);
    if(ghostRecII>0) hntrkghost_II->Fill(ghostRecII);
    if((nMCtracks-goodRecII)>0){
      hntrkmissed_II->Fill(nMCtracks-goodRecII);
    }

    for(int imc=0;imc<nParticles;imc++){//MC trks
      bool missTrk=true;
      for(int irec=0;irec<nGeaneTrks;irec++){//RECids assigment
        int mc_comp = RECtrkMCid[irec];
        if(mc_comp==imc) missTrk=false;
      }

      PndMCTrack *mctrk =(PndMCTrack*) true_tracks->At(imc);
      TVector3 MomMC = mctrk->GetMomentum();
      TVector3 PosMC = mctrk->GetStartVertex();
      int trkQ;
      if(missTrk) trkQ=-1;
      else trkQ=0;
      ntupMCTrk->Fill(PosMC.X(),PosMC.Y(),PosMC.Z(),MomMC.Mag(),MomMC.Theta(),MomMC.Phi(),trkQ,nRecHits);
    }

    if(verboseLevel>0){
      cout<<"-- (II) Hits matching -------------------------"<<endl;
      cout<<"Good trks: "<<goodRecII<<" missed: "<<nMCtracks-goodRecII<<" ghost: "<<ghostRecII<<endl;
      cout<<"--------------------------------"<<endl;
      cout<<" "<<endl;
    }
    //  cout<<"number of good: "<<goodRecII<<" number of missed: "<<nMCtracks-goodRecII<<" number of ghost: "<<ghostRecII<<endl;

    //    cout<<"nGeaneTrks = "<<nGeaneTrks<<" nRecGEANEtrk = "<<nRecGEANEtrk<<endl;
    hntrk->Fill(nRecGEANEtrk);
    //  cout<<"========="<<endl;
    for(Int_t jk=0; jk< nParticles;jk++){
      // if(nParticles!=4) continue;
      PndMCTrack *mctrk =(PndMCTrack*) true_tracks->At(jk);
      Int_t mcID = mctrk->GetPdgCode();
      if(mcID==-2212){
        //  if(nParticles!=4) cout<<"For event #"<<j<<" mcID="<<mcID<<endl;
        TVector3 MomMC = mctrk->GetMomentum();
        hallTheta->Fill(MomMC.Theta());
        hallPhi->Fill(MomMC.Phi());
        hallPhiTheta->Fill(MomMC.Theta(),MomMC.Phi());
      }
    }
  }

  // TCanvas *c1 = new TCanvas("IDsINFO");
  // c1->Divide(2,2);
  // c1->cd(1);
  // hnRecnMC->Draw();
  // c1->cd(2);
  // hntrkcand->Draw();
  // hntrk->SetLineColor(2);
  // hntrk->Draw("same");
  // c1->cd(3);
  // hDiffIDs->Draw();
  // c1->cd(4);
  // // hntrkcandvsIDs->Draw();
  // hchi2->Draw();
  // c1->Write();
  // c1->Close();
  TCanvas *c1 = new TCanvas("IDsINFO");
  c1->Divide(2,2);
  c1->cd(1);
  hnRecnMC->Draw("colz");
  c1->cd(2);
  hntrkcand->Draw();
  hntrk->SetLineColor(2);
  hntrk->Draw("same");
  c1->cd(3);
  hDiffIDs->Draw();
  c1->cd(4);
  //  hntrkcandvsIDs->Draw("colz");
  // hntrkcandvsMC->Draw("colz");
  hnhitsvsIDs->Draw("colz");
  //hchi2->Draw();
  c1->Write();
  c1->Close();

  TH1F *heffTheta = (TH1F*)hgoodTheta->Clone("heffTheta");
  heffTheta->SetTitle("Efficiency vs. #theta");
  heffTheta->Divide(hallTheta);
  TH1F *heffPhi = (TH1F*)hgoodPhi->Clone("heffPhi");
  heffPhi->SetTitle("Efficiency vs. #phi");
  heffPhi->Divide(hallPhi);

  TH2F *heffPhiTheta = (TH2F*)hgoodPhiTheta->Clone("heffPhiTheta");
  heffPhiTheta->SetTitle("Efficiency #phi & #theta");
  heffPhiTheta->Divide(hallPhiTheta);

  TCanvas *c2 = new TCanvas("EffINFO");
  c2->Divide(3,2);
  c2->cd(1);
  hallTheta->SetLineWidth(2);
  hallTheta->SetFillStyle(3354);
  hallTheta->SetFillColor(kBlack);
  hallTheta->Draw();
  hgoodTheta->SetLineWidth(2);
  hgoodTheta->SetLineColor(2);
  hgoodTheta->SetFillStyle(3690);
  hgoodTheta->SetFillColor(kRed);
  hgoodTheta->Draw("same");
  c2->cd(4);
  heffTheta->SetMinimum(0.9);
  heffTheta->Draw();
  c2->cd(3);
  hallPhiTheta->Draw("colz");
  c2->cd(2);
  hallPhi->SetLineWidth(2);
  hallPhi->SetFillStyle(3354);
  hallPhi->SetFillColor(kBlack);
  hallPhi->Draw();
  hgoodPhi->SetLineWidth(2);
  hgoodPhi->SetLineColor(2);
  hgoodPhi->SetFillColor(kRed);
  hgoodPhi->SetFillStyle(3690);
  hgoodPhi->Draw("same");
  c2->cd(5);
  heffPhi->SetMinimum(0.9);
  heffPhi->Draw();
  c2->cd(6);
  heffPhiTheta->Draw("colz");

  c2->Write();
  c2->Close();

  TF1 *funrp = new TF1("fitrp","gaus",-1e-5,1e-5);
  funrp->SetParameters(100,0,3e-3);
  funrp->SetParNames("Constant","Mean","Sigma");

  TF1 *funrth = new TF1("fitrth","gaus",-0.01,0.01);
  funrth->SetParameters(100,0,3e-3);
  funrth->SetParNames("Constant","Mean","Sigma");

  TF1 *funrphi = new TF1("fitrphi","gaus",-1,1);
  funrphi->SetParameters(100,0,3e-3);
  funrphi->SetParNames("Constant","Mean","Sigma");

  TF1 *funp = new TF1("fitp","gaus",-20,20);
  funp->SetParameters(100,0,3e-3);
  funp->SetParNames("Constant","Mean","Sigma");

  TCanvas *c3 = new TCanvas("ResMomINFO");
  c3->Divide(3,3);
  c3->cd(1);
  hResMom->Fit(funrp,"r");
  hResMom->Draw();
  c3->cd(2);
  hResTheta->Fit(funrth,"r");
  hResTheta->Draw();
  c3->cd(3);
  hResPhi->Fit(funrphi,"r");
  hResPhi->Draw();
  c3->cd(4);
  //  hPullMom->Fit(funp,"r");
  hPullMom->Draw();
  c3->cd(5);
  hPullTheta->Fit(funp,"r");
  hPullTheta->Draw();
  c3->cd(6);
  hPullPhi->Fit(funp,"r");
  hPullPhi->Draw();
  c3->cd(7);
  hErrMom->Draw();
  c3->cd(8);
  hErrTheta->Draw();
  c3->cd(9);
  hErrPhi->Draw();

  c3->Write();
  c3->Close();


  TF1 *funcoord = new TF1("fitcoord","gaus",-1,1.);
  funcoord->SetParameters(1e4,0,1e-1);
  funcoord->SetParNames("Constant","Mean","Sigma");

  // TF1 *funp = new TF1("fitp","gaus",-10,10);
  funp->SetParameters(100,0,1);
  funp->SetParNames("Constant","Mean","Sigma");

  TCanvas *c4 = new TCanvas("ResPCAINFO");
  c4->Divide(3,3);
  c4->cd(1);
  hResPointX->Fit(funcoord,"r");
  hResPointX->Draw();
  c4->cd(2);
  hResPointY->Fit(funcoord,"r");
  hResPointY->Draw();
  c4->cd(3);
  hResPointZ->Fit(funcoord,"r");
  hResPointZ->Draw();
  c4->cd(4);
  hErrPointX->Draw();
  c4->cd(5);
  hErrPointY->Draw();
  c4->cd(6);
  hErrPointZ->Draw();
  c4->cd(7);
  hPullPointX->Fit(funp,"r");
  hPullPointX->Draw();
  c4->cd(8);
  hPullPointY->Fit(funp,"r");
  hPullPointY->Draw();
  c4->cd(9);
  hPullPointZ->Fit(funp,"r");
  hPullPointZ->Draw();
  c4->Write();
  c4->Close();

  TCanvas *c42 = new TCanvas("ResMomPointINFO");
  c42->Divide(3,3);
  c42->cd(1);
  hResPointPx->Fit(funcoord,"r");
  hResPointPx->Draw();
  c42->cd(2);
  hResPointPy->Fit(funcoord,"r");
  hResPointPy->Draw();
  c42->cd(3);
  hResPointPz->Fit(funcoord,"r");
  hResPointPz->Draw();
  c42->cd(4);
  hErrPointPx->Draw();
  c42->cd(5);
  hErrPointPy->Draw();
  c42->cd(6);
  hErrPointPz->Draw();
  c42->cd(7);
  hPullPointPx->Fit(funp,"r");
  hPullPointPx->Draw();
  c42->cd(8);
  hPullPointPy->Fit(funp,"r");
  hPullPointPy->Draw();
  c42->cd(9);
  hPullPointPz->Fit(funp,"r");
  hPullPointPz->Draw();
  c42->Write();
  c42->Close();



  TCanvas *c41 = new TCanvas("ResLinINFO");
  c41->Divide(2,2);
  c41->cd(1);
  hErrLinPointX->Draw();
  c41->cd(2);
  hErrLinPointY->Draw();
  c41->Write();
  c41->Close();

  // TCanvas *c5 = new TCanvas("HitIdINFO");
  // c5->Divide(1,2);
  // c5->cd(1);
  // hHitIDHitX->SetMarkerStyle(24);
  // hHitIDHitX->Draw();
  // c5->cd(2);
  // hHitIDHitY->SetMarkerStyle(24);
  // hHitIDHitY->Draw();
  // c5->Write();
  // c5->Close();

  TCanvas *c6 = new TCanvas("VertexINFO");
  c6->Divide(2,2);
  c6->cd(1);
  hPointXY->Draw();
  c6->cd(2);
  hPointXYcut->Draw();
  c6->cd(3);
  hErrPointXY->Draw();
  c6->Write();
  c6->Close();

  TCanvas *c7 = new TCanvas("GEANEvsRecINFO");
  c7->Divide(3,2);
  c7->cd(1);
  hRecGEANEX->Draw();
  c7->cd(2);
  hRecGEANEY->Draw();
  c7->cd(3);
  hRecGEANEZ->Draw();
  c7->cd(4);
  hRecGEANETheta->Draw();
  c7->cd(5);
  hRecGEANEPhi->Draw();
  c7->cd(6);
  hRecThetaPhi->Draw();
  c7->Write();
  c7->Close();

  TCanvas *c8 = new TCanvas("GEANEvsCandINFO");
  c8->Divide(3,2);
  c8->cd(1);
  hSeedGEANEX->Draw();
  c8->cd(2);
  hSeedGEANEY->Draw();
  c8->cd(3);
  hSeedGEANEZ->Draw();
  c8->cd(4);
  hSeedGEANETheta->Draw();
  TLine *uplim = new TLine(0.05, 0, 0.05, 0.1);
  TLine *downlim = new TLine(0.03, 0, 0.03, 0.1);
  uplim->SetLineColor(kRed);
  downlim->SetLineColor(kRed);
  uplim->SetLineWidth(2);
  downlim->SetLineWidth(2);
  uplim->Draw();
  downlim->Draw();
  // ---- Count number events in range [3-8] mrad ------------------------
  Int_t binx1, binx2;
  Int_t binmax = hSeedGEANETheta->GetNbinsX();
  for(int bin=0;bin<binmax;bin++){
    Double_t binCenter = hSeedGEANETheta->GetBinCenter(bin);
    if(binCenter>0.05){
    //    if(binCenter>5){
      binx2 = bin;
      break;
    }
  }
  for(int bin=0;bin<binmax;bin++){
    Double_t binCenter = hSeedGEANETheta->GetBinCenter(bin);
    if(binCenter>0.03){
      binx1 = bin;
      break;
    }
  }
  //  cout<<"binmax = "<<binmax<<" binx1 = "<<binx1<<" binx2 = "<<binx2<<endl;
  Double_t nall = hSeedGEANETheta->Integral(0, binmax);
  Double_t neff = hSeedGEANETheta->Integral(binx1, binx2);
  // cout<<"Number of track-candidate with #theta from [30-50] mrad = "<<neff<<endl;
  cout<<"Total number of track-candidate = "<<nall<<endl;

  //----------------------------------------------------------------------
  c8->cd(5);
  hSeedGEANEPhi->Draw();
  c8->cd(6);
  hSeedThetaPhi->Draw();
  c8->Write();
  c8->Close();

  TCanvas *c9 = new TCanvas("DCA");
  c9->Divide(3,2);
  c9->cd(1);
  hResDCA->Draw();
  c9->cd(2);
  hErrDCA->Draw();
  c9->cd(3);
  hPullDCA->Draw();
  c9->Write();
  c9->Close();

  hangMCRec->Write();
  hchi2->Write();
  heffTheta->Write();
  heffPhi->Write();

  hnRecnMC->Write();
  hDiffIDs->Write();
  hntrkcand->Write();
  hntrk->Write();
  hntrkcandvsIDs->Write();

  hResMom->Write();
  hErrMom->Write();
  hPullMom->Write();
  hResPhi->Write();
  hErrPhi->Write();
  hPullPhi->Write();
  hResTheta->Write();
  hErrTheta->Write();
  hPullTheta->Write();
  hangMCRec->Write();
  hPointXY->Write();
  hPointXYcut->Write();
  hSeedGEANEX->Write();
  hSeedGEANEY->Write();
  hSeedGEANEZ->Write();
  hSeedGEANER->Write();
  hSeedGEANETheta->Write();
  hSeedGEANEPhi->Write();
  hSeedThetaPhi->Write();
  hRecGEANEX->Write();
  hRecGEANEY->Write();
  hRecGEANEZ->Write();
  hchi2Errx->Write();
  hchi2Erry->Write();

  TCanvas *c10 = new TCanvas("ThetaPhiDistr");
  c10->Divide(4,2);
  // c10->cd(1);
  //  hSeedThetaDiffIDs->Draw("colz");
  // c10->cd(2);
  //  hSeedPhiDiffIDs->Draw("colz");
  // c10->cd(1);
  // hMCThetaGEANETheta->Draw("colz");
  // c10->cd(2);
  // hMCPhiGEANEPhi->Draw("colz");
  c10->cd(1);
  hMCThetaResTheta->Draw("colz");
  c10->cd(5);
  hMCPhiResPhi->Draw("colz");

  TH1F *fdiffThetaMean;
  TH1F *fdiffThetaSigma;
  TH1F *fdiffThetaChi2;
  hMCThetaResTheta->FitSlicesY();
  fdiffThetaMean = (TH1F*)gDirectory->Get("hMCThetaResTheta_1");
  fdiffThetaSigma = (TH1F*)gDirectory->Get("hMCThetaResTheta_2");
  fdiffThetaChi2= (TH1F*)gDirectory->Get("hMCThetaResTheta_chi2");
  TH1F *fdiffPhiMean;
  TH1F *fdiffPhiSigma;
  TH1F *fdiffPhiChi2;
  hMCPhiResPhi->FitSlicesY();
  fdiffPhiMean = (TH1F*)gDirectory->Get("hMCPhiResPhi_1");
  fdiffPhiSigma = (TH1F*)gDirectory->Get("hMCPhiResPhi_2");
  fdiffPhiChi2 = (TH1F*)gDirectory->Get("hMCPhiResPhi_chi2");
  c10->cd(2);
  fdiffThetaMean->Draw();
  c10->cd(6);
  fdiffPhiMean->Draw();
  c10->cd(3);
  fdiffThetaSigma->Draw();
  c10->cd(7);
  fdiffPhiSigma->Draw();
  c10->cd(4);
  fdiffThetaChi2->Draw();
  c10->cd(8);
  fdiffPhiChi2->Draw();
  c10->Write();
  c10->Close();

   TCanvas *c11 = new TCanvas("HitsUsing");
   c11->Divide(2,2);
   c11->cd(1);
   hntrkMCtrkID->Draw("colz");
   c11->cd(2);
   // hntrkMCtrkID->ProjectionY();
   // hntrkMCtrkID_py->Draw();
   c11->cd(3);
   // hntrkMCtrkID->ProjectionX();
   // hntrkMCtrkID_px->Draw();
   c11->Write();
   c11->Close();

  TCanvas *c12 = new TCanvas("TrkLin_and_MC");
  c12->Divide(3,3);
  c12->cd(1);
  hResLumiTrkMom->Fit(funrp,"r");
  hResLumiTrkMom->Draw();
  c12->cd(2);
  hResLumiTrkTheta->Fit(funrth,"r");
  hResLumiTrkTheta->Draw();
  c12->cd(3);
  hResLumiTrkPhi->Fit(funrphi,"r");
  hResLumiTrkPhi->Draw();
  c12->cd(4);
  hResLumiTrkPointX->Fit(funcoord,"r");
  hResLumiTrkPointX->Draw();
  c12->cd(5);
  hResLumiTrkPointY->Fit(funcoord,"r");
  hResLumiTrkPointY->Draw();
  c12->cd(6);
  hResLumiTrkPointZ->Fit(funcoord,"r");
  hResLumiTrkPointZ->Draw();
  c12->cd(7);
  hResLumiTrkPointPx->Fit(funp,"r");
  hResLumiTrkPointPx->Draw();
  c12->cd(8);
  hResLumiTrkPointPy->Fit(funp,"r");
  hResLumiTrkPointPy->Draw();
  c12->cd(9);
  hResLumiTrkPointPz->Fit(funp,"r");
  hResLumiTrkPointPz->Draw();
  c12->Write();
  c12->Close();

 TCanvas *c13 = new TCanvas("TrkLin_and_MC_pulls");
 c13->Divide(3,2);
 c13->cd(1);
 hResLumiTrkPointXPull->Fit(funp,"r");
 hResLumiTrkPointXPull->Draw();
 c13->cd(2);
 hResLumiTrkPointYPull->Fit(funp,"r");
 hResLumiTrkPointYPull->Draw();
 c13->cd(3);
 hResLumiTrkPointZPull->Fit(funp,"r");
 hResLumiTrkPointZPull->Draw();
 c13->cd(4);
 hResLumiTrkPointPxPull->Fit(funp,"r");
 hResLumiTrkPointPxPull->Draw();
 c13->cd(5);
 hResLumiTrkPointPyPull->Fit(funp,"r");
 hResLumiTrkPointPyPull->Draw();
 c13->cd(6);
 hResLumiTrkPointPzPull->Fit(funp,"r");
 hResLumiTrkPointPzPull->Draw();
 c13->Write();
 c13->Close();

 TCanvas *c14 = new TCanvas("TrkIP_andLin_and_MC_pulls");
 c14->Divide(3,2);
 c14->cd(1);
 //Grey = Trk near IP
 //Blue = Trk near Lumi
 hPullPointX->SetLineWidth(2);
 hPullPointX->SetLineColor(12);
 hPullPointX->SetFillColor(12);
 hPullPointX->SetFillStyle(3003);
 hPullPointX->Draw();
 hResLumiTrkPointXPull->SetLineColor(4);
 hResLumiTrkPointXPull->SetLineWidth(2);
 hResLumiTrkPointXPull->Draw("same");
 c14->cd(2);
 hPullPointY->SetLineWidth(2);
 hPullPointY->SetLineColor(12);
 hPullPointY->SetFillColor(12);
 hPullPointY->SetFillStyle(3003);
 hPullPointY->Draw();
 hResLumiTrkPointYPull->SetLineColor(4);
 hResLumiTrkPointYPull->SetLineWidth(2);
 hResLumiTrkPointYPull->Draw("same");
 c14->cd(3);
 hPullPointZ->SetLineWidth(2);
 hPullPointZ->SetLineColor(12);
 hPullPointZ->SetFillColor(12);
 hPullPointZ->SetFillStyle(3003);
 hPullPointZ->Draw();
 hResLumiTrkPointZPull->SetLineColor(4);
 hResLumiTrkPointZPull->SetLineWidth(2);
 hResLumiTrkPointZPull->Draw("same");
 c14->cd(4);
 hPullPointPx->SetLineWidth(2);
 hPullPointPx->SetLineColor(12);
 hPullPointPx->SetFillColor(12);
 hPullPointPx->SetFillStyle(3003);
 hPullPointPx->Draw();
 hResLumiTrkPointPxPull->SetLineColor(4);
 hResLumiTrkPointPxPull->SetLineWidth(2);
 hResLumiTrkPointPxPull->Draw("same");
 c14->cd(5);
 hPullPointPy->SetLineWidth(2);
 hPullPointPy->SetLineColor(12);
 hPullPointPy->SetFillColor(12);
 hPullPointPy->SetFillStyle(3003);
 hPullPointPy->Draw();
 hResLumiTrkPointPyPull->SetLineColor(4);
 hResLumiTrkPointPyPull->SetLineWidth(2);
 hResLumiTrkPointPyPull->Draw("same");
 c14->cd(6);
 hPullPointPz->SetLineWidth(2);
 hPullPointPz->SetLineColor(12);
 hPullPointPz->SetFillColor(12);
 hPullPointPz->SetFillStyle(3003);
 hPullPointPz->Draw();
 hResLumiTrkPointPzPull->SetLineColor(4);
 hResLumiTrkPointPzPull->SetLineWidth(2);
 hResLumiTrkPointPzPull->Draw("same");
 c14->Write();
 c14->Close();


  hSeedThetaDiffIDs->Write();
  hSeedPhiDiffIDs->Write();
  hMCThetaGEANETheta->Write();
  hMCPhiGEANEPhi->Write();
  hMCThetaResTheta->Write();
  hMCPhiResPhi->Write();

  hDiffIDsPointX->Write();
  hDiffIDsPointY->Write();
  hDiffIDsPointZ->Write();

  hResPointPx->Write();
  hErrPointPx->Write();
  hPullPointPx->Write();
  hResPointPy->Write();
  hErrPointPy->Write();
  hPullPointPy->Write();
  hResPointPz->Write();
  hErrPointPz->Write();
  hPullPointPz->Write();

  hResPointX->Write();
  hErrPointX->Write();
  hPullPointX->Write();
  hResPointY->Write();
  hErrPointY->Write();
  hPullPointY->Write();
  hResPointZ->Write();
  hErrPointZ->Write();
  hPullPointZ->Write();

  hMCtrkvshits->Write();
  hntrkMCtrkID->Write();
  hchi2nTrkCand->Write();
  hgoodPhichi2->Write();
  hgoodThetachi2->Write();
  hchi2MCdiffID->Write();
  hnhitsvsIDs->Write();
  hResLumiTrkMom->Write();
  hResLumiTrkTheta->Write();
  hResLumiTrkPhi->Write();
  hResLumiTrkThetaPull->Write();
  hResLumiTrkPhiPull->Write();
  hResLumiTrkPointX->Write();
  hResLumiTrkPointY->Write();
  hResLumiTrkPointZ->Write();
  hResLumiTrkPointPx->Write();
  hResLumiTrkPointPy->Write();
  hResLumiTrkPointPz->Write();
  hResLumiTrkPointXPull->Write();
  hResLumiTrkPointYPull->Write();
  hResLumiTrkPointZPull->Write();
  hResLumiTrkPointPxPull->Write();
  hResLumiTrkPointPyPull->Write();
  hResLumiTrkPointPzPull->Write();
  //  hMomMC->Write();
  hResZResPhi->Write();
  nrecpointall->Write();
  nrecdirall->Write();

  TCanvas *c134 = new TCanvas("TrkLinparams");
 c134->Divide(2,2);
 c134->cd(1);
 hLumiTrkA->Draw();
 c134->cd(2);
 hLumiTrkB->Draw();
 c134->cd(3);
 hLumiTrkC->Draw();
 c134->cd(4);
 hLumiTrkD->Draw();
 c134->Write();
 c134->Close();
 hLumiTrkA->Write();
 hLumiTrkB->Write();
 hLumiTrkC->Write();
 hLumiTrkD->Write();
 hMCThetaResSeedTheta->Write();
 hMCPhiResSeedPhi->Write();
 hntrkmissedPhiTheta->Write();
 hntrkmissed_I->Write();
 hntrkghost_I->Write();
 hntrkmissed_II->Write();
 hntrkghost_II->Write();
 hntrkgood_I->Write();
 hntrkgood_II->Write();
 heffPhiTheta->Write();
 nmomMC->Write();
 ntuprecTrk->Write();
 ntupMCTrk->Write();
 hnhits->Write();
 nsectors->Write();
 hPullPointXphi->Write();
 hPullPointYphi->Write();
 hPullPointXtheta->Write();
 hPullPointYtheta->Write();
 f->Close();
 cout<<"Number of events with low number of hits (less then 3 per trk): "<<glBadEv<<endl;
 cout<<"Number of trks where GEANE failed: "<<glBADGEANE<<endl;
 cout<<"Total number of noise hits = "<<glNoisehit<<endl;
}