// root macro to analyze the clusterization output
{
  int  nEvents = 1000;
  bool verbose = false;

  // -----  Load libraries   ------------------------------------------------
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  basiclibs();
  gSystem->Load("libGeoBase");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");
  gSystem->Load("libMCStack");
  gSystem->Load("libField");
  gSystem->Load("libGen");
  gSystem->Load("libPassive");

  gSystem->Load("libgenfit");
  gSystem->Load("libtrackrep");
  gSystem->Load("libtpc");
  gSystem->Load("libtpcreco");
  gSystem->Load("librecotasks");

  gSystem->Load("libMvd");
  gSystem->Load("libMvdReco");

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


  PndMvdFileNameCreator namecreator("/home/ralfk/MVD/mvdmacros/data/mvdevalg4.root");
  std::string inFile = namecreator.GetSimFileName(false);
  std::string digiFile = namecreator.GetDigiFileName(false);

  TFile* f = new TFile(inFile.c_str()); // the sim file you want to analyse
  TTree *t=(TTree *) f->Get("cbmsim") ;
  t->AddFriend("cbmsim",digiFile.c_str()); // the digi file

  TClonesArray* mc_array=new TClonesArray("PndMvdMCPoint");
  t->SetBranchAddress("MVDPoint",&mc_array);//Branch names

  TClonesArray* digiPixel_array=new TClonesArray("PndMvdDigiPixel");
  t->SetBranchAddress("MVDPixelHit",&digiPixel_array);//Branch names

  TClonesArray* digiStrip_array=new TClonesArray("PndMvdDigiStrip");
  t->SetBranchAddress("MVDStripHit",&digiStrip_array);//Branch names

  TGeoManager *geoMan = (TGeoManager*) gDirectory->Get("CBMGeom");

  TH2D* hisxy = new TH2D("hisxy","",400,-15.,15.,400,-15.,15.);
  hisxy->SetTitle("MVD MC Cluster, xy view;x / cm;y / cm");
  TH2D* hisrz = new TH2D("hisrz","",400,-20.,20.,400,-15.,25.);
  hisrz->SetTitle("MVD MC Cluster, rz view;z / cm;r/ cm");

  int n = 100;
  int low = 0;
  TH1I* hisPixelCol = new TH1I("hispixelcol","Pixel column channel number on FE"
    ,n,low,low+n);
  TH1I* hisPixelRow = new TH1I("hispixelrow","Pixel row channel number on FE"
    ,n,low,low+n);
  TH1I* hisPixelFE = new TH1I("hispixelfe","Pixel FE number"
    ,n,low,low+n);

  TH1I* hisCol = new TH1I("hiscol","column number",1001,0,1000);
  TH1I* hisRow = new TH1I("hisrow","row number",1001,0,1000);

  TH2I* hisColRow = new TH2I("hiscolrow","col & row number",1001,0,1000,1001,0,1000);

  TH1D* hisPixelCharge = new TH1D("hispixelcharge","Pixel Charge content",100,0.,1e5);

  TH2I* hisStripTop = new TH2I("hisstriptop","Strip Top channel&fe numbers",20,0,20,200,0,200);
  TH2I* hisStripBot = new TH2I("hisstripbot","Strip Bot channel&fe numbers",20,0,20,200,0,200);
  TH1I* hisStripStrip = new TH1I("hisstripstrip","Strip numbers",15*128,0,15*128);
  TVector3 vecmc;
  Double_t tmpx,tmpy,tmpz;
  TVector2 locals, localmc, localdiff;
  int col, row, fe;
  double x,y;

  PndMvdCalcFePixel pixelcalc(76, 84, 10);

  double pitch=0.004921;
  double orient=TMath::Pi()*(0.5);
  int nrFeChannels=128;
  int nrStrips=10*nrFeChannels;
  TVector2 Anchor(0.,0.);
  double offset=0.;
  double noise=0.;
//   PndMvdCalcStrip stripcalc(pitch,orient,nrStrips,nrFeChannels,Anchor,offset,noise)

  for (Int_t j=0; j<nEvents && j<t->GetEntriesFast(); j++)
  {
    t->GetEntry(j);
    if(verbose) cout<<"Event No "<<j<<endl;

    for (Int_t i=0; i<mc_array->GetEntriesFast(); i++)
    {
      if(verbose) cout<<"Point No "<<i<<endl;
      PndMvdMCPoint *point=(PndMvdMCPoint*)mc_array->At(i);
      vecmc.SetXYZ(point->GetX(),point->GetY(),point->GetZ());
      hisxy->Fill(vecmc.x(),vecmc.y());
      if(vecmc.y() > 0.) hisrz->Fill(vecmc.z(),vecmc.Perp());
      else hisrz->Fill(vecmc.z(),-1.*vecmc.Perp());
    }

    for (Int_t i=0; i<digiPixel_array->GetEntriesFast(); i++)
    {
      PndMvdDigiPixel *pixeldigi = digiPixel_array->At(i);
      fe = pixeldigi->GetFE();
      col = pixeldigi->GetPixelColumn();
      row = pixeldigi->GetPixelRow();
      hisPixelCol->Fill(col);
      hisPixelRow->Fill(row);
      hisPixelFE->Fill(fe);
      hisPixelCharge->Fill(pixeldigi->GetCharge());
      pixelcalc.CalcSensorColRow(col,row,fe);
      hisCol->Fill(col);
      hisRow->Fill(row);
      hisColRow->Fill(col,row);
    }

    for (Int_t i=0; i<digiStrip_array->GetEntriesFast(); i++)
    {
      PndMvdDigiStrip *stripdigi = digiStrip_array->At(i);
      fe = stripdigi->GetFE();
      col = stripdigi->GetChannel();
      int strip = fe * nrFeChannels + col;
      hisStripStrip->Fill(strip);
      if (strip <= nrStrips) hisStripTop->Fill(fe,col);
      else hisStripBot->Fill(fe,col);
    }

  }// end for j (events)

Int_t a = 4, b = 3;

TCanvas* can1 = new TCanvas("MvdTestPlot","MCHit view in MVD",0,0,a*300,b*300);
can1->Divide(a,b);

can1->cd(1);hisxy->DrawCopy("col");
can1->cd(2);hisrz->DrawCopy("col");
can1->cd(3);hisPixelCol->DrawCopy();
can1->cd(4);hisStripStrip->DrawCopy();
can1->cd(5);hisPixelRow->DrawCopy();
can1->cd(6);hisPixelFE->DrawCopy();
can1->cd(7);hisPixelCharge->DrawCopy();
can1->cd(8);hisStripTop->DrawCopy("col");
can1->cd(9);hisCol->DrawCopy();
can1->cd(10);hisRow->DrawCopy();
can1->cd(11);hisColRow->DrawCopy("col");
can1->cd(12);hisStripBot->DrawCopy("col");


// can1->Print("testoutput.ps");


  // -----   Finish   -------------------------------------------------------
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout << endl << endl;
  cout << "Macro finished succesfully." << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
  cout << endl;
  // ------------------------------------------------------------------------

}