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

  // -----  Load libraries   ------------------------------------------------
//  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
 gROOT->Macro("$VMCWORKDIR/macro/mvd/Libs.C");

  gROOT->LoadMacro("../Tools.C");
  LoadPandaStyle();

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

  // -----   Get Data from framework   --------------------------------------
  // load first the file with the GeoManager
  PndMvdFileNameCreator namecreator("../data/MvdTracking.root");
  std::string inFile = namecreator.GetSimFileName(false);
  std::string recoFile = namecreator.GetRecoFileName(false);
  TString picture = namecreator.GetRecoFileName(false);
  picture.ReplaceAll(".root",".ps");

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

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

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

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

  TClonesArray* strclust_array=new TClonesArray("PndMvdCluster");
  t->SetBranchAddress("MVDStripClusterCand",&strclust_array);//Branch names

  TClonesArray* pixclust_array=new TClonesArray("PndMvdCluster");
  t->SetBranchAddress("MVDClusterCand",&pixclust_array);//Branch names

  TClonesArray* strhit_array=new TClonesArray("PndMvdHit");
  t->SetBranchAddress("MVDHitsStrip",&strhit_array);//Branch names

  TClonesArray* pixhit_array=new TClonesArray("PndMvdHit");
  t->SetBranchAddress("MVDHitsPixel",&pixhit_array);//Branch names


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

  Int_t nbins = 100; //200
  TH2D* hisxymc = new TH2D("xymc","",nbins,-15.,15.,nbins,-15.,15.);
  hisxymc->SetTitle("MVD reco Cluster, xy view;x / cm;y / cm");
  TH2D* hisrzmc = new TH2D("rzmc","",nbins,-20.,20.,nbins,-15.,25.);
  hisrzmc->SetTitle("MVD reco Cluster, rz view;z / cm;r/ cm");
  TH1D* hisde = new TH1D("de","MVD MC Points, Energyloss",nbins,0.,0.002);
  TH1D* hismom = new TH1D("mom","MVD MC Points, momentum",nbins,0.,1.5);


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


  TH2D* hisDiffXY = new TH2D("diffxy","",nbins,-0.02,0.02,nbins,-0.02,0.02);
  hisDiffXY->SetTitle("(MC - RECO) Pixel Hit coordinates xy view;#Deltax / cm;#Deltay / cm");
  TH2D* hisDiffRZ = new TH2D("diffrz","",nbins,-0.02,0.02,nbins,-0.00,0.04);
  hisDiffRZ->SetTitle("(MC - RECO) Pixel Hit coordinates rz view;#Deltaz / cm;#Deltar / cm");

  Double_t range = 0.01;
  Double_t zed = 0.01;
  Int_t bins = 100;

  TH1D* hisDiffX = new TH1D("diffx","",bins,-range,range);
  hisDiffX->SetTitle("(MC - RECO) Pixel Hit coordinate x;#Deltax / cm;");
  TH1D* hisDiffY = new TH1D("diffy","",bins,-range,range);
  hisDiffY->SetTitle("(MC - RECO) Pixel Hit coordinate y;#Deltay / cm;");
  TH1D* hisDiffZ = new TH1D("diffz","",bins,-range,range);
  hisDiffZ->SetTitle("(MC - RECO) Pixel Hit coordinate z;#Deltaz / cm;");

  double angmax = 0.02;
  TH1D* hisDiffTheta = new TH1D("DiffTheta","",100,-angmax,angmax);
  hisDiffTheta->SetTitle("#Theta resolution;#Delta#Theta/mrad;");

  TH1D* hisDiffPhi = new TH1D("DiffPhi","",100,-angmax,angmax);
  hisDiffPhi->SetTitle("#Phi resolution;#Delta#Phi/mrad;");

  TH1D* hisDiffStripX = new TH1D("diffLocalx","",bins,-range,range);
  hisDiffStripX->SetTitle("(MC - RECO) Strip Hit coordinate x;#Deltax_{L} / cm;");
  TH1D* hisDiffStripY = new TH1D("diffLocaly","",bins,-range,range);
  hisDiffStripY->SetTitle("(MC - RECO) Strip Hit coordinate y;#Deltay_{L} / cm;");
  TH1D* hisDiffStripZ = new TH1D("diffLocalz","",bins,-zed,zed);
  hisDiffStripZ->SetTitle("(MC - RECO) Strip Hit coordinate z;#Deltaz_{L} / cm;");

  TH2D* hisDiffStripXY = new TH2D("diffLocalxy","",bins,-range,range,bins,-range,range);
  hisDiffStripXY->SetTitle("(MC - RECO) Strip Hit XY;#Deltax_{L} / cm;#Deltay_{L} / cm");
  TH2D* hisDiffStripRZ = new TH2D("diffLocalrz","",bins,-range,range,bins,-0.00,2.*range);
  hisDiffStripRZ->SetTitle("(MC - RECO) Strip Hit RZ;#Deltaz_{L} / cm;#Deltar_{L} / cm");

  TH2D* hisSensorXY = new TH2D("Localxy","",bins,-5.,5.,bins,-5.,5.);
  hisSensorXY->SetTitle("Local Hit XY;x_{L} / cm;y_{L} / cm");
  TH2D* hisSensorXYMC = new TH2D("LocalxyMC","",bins,-5.,5.,bins,-5.,5.);
  hisSensorXYMC->SetTitle("Local MC Point XY;x_{L} / cm;y_{L} / cm");

  TH1I* hisClustSize = new TH1I("ClustSize","Cluster size",10,0,10);
  hisClustSize->SetTitle("Cluster size;n_{cl};");
  TH1I* hisClustSizeDisk = new TH1I("ClustSizeDisk","Cluster size",10,0,10);
  TH1I* hisClustSizeSide = new TH1I("ClustSizeSize","Cluster size",10,0,10);
  TH1I* hisClustSizeTop = new TH1I("ClustSizetop","Cluster size",10,0,10);
  TH1I* hisClustSizeBot = new TH1I("ClustSizebot","Cluster size",10,0,10);

  TVector3 vecs, vecmc, vecdiff, mommc, sensorDim;
  Double_t tmpMaster[3],tmpLocal[3];
  TVector2 locals, localmc, localdiff;
  TGeoHMatrix* currentTransMat;
  TString detname;
  Double_t difftheta, diffphi;

  for (Int_t j=0; j<nEvents && j<t->GetEntriesFast(); j++)
  {
    t->GetEntry(j);
    if(verbose) cout<<"Event No "<<j<<endl;
    else if (!(j%100)) cout <<"Event No "<<j<<endl;
    if(verbose) std::cout<<"Check sanity of the arrays:"<<std::endl;
    if(mc_array){
      if(verbose) std::cout<<"mc_array ok with "<<mc_array->GetEntriesFast()<<" entries"<<std::endl;
    }
    else std::cout<<"*** mc_array broken - check your file! ."<<std::endl;
    if(digiPixel_array){
      if(verbose) std::cout<<"digiPixel_array ok with "<<digiStrip_array->GetEntriesFast()<<" entries"<<std::endl;
    }
    else std::cout<<"*** digiPixel_array broken - check your file! ."<<std::endl;
    if(digiStrip_array){
      if(verbose) std::cout<<"digiStrip_array ok with "<<digiStrip_array->GetEntriesFast()<<" entries"<<std::endl;
    }
    else std::cout<<"*** digiStrip_array broken - check your file! ."<<std::endl;
    if(strclust_array){
      if(verbose) std::cout<<"strclust_array ok with "<<strclust_array->GetEntriesFast()<<" entries"<<std::endl;
    }
    else std::cout<<"*** strclust_array broken - check your file! ."<<std::endl;
    if(pixclust_array){
      if(verbose) std::cout<<"pixclust_array ok with "<<pixclust_array->GetEntriesFast()<<" entries"<<std::endl;
    }
    else std::cout<<"*** pixclust_array broken - check your file! ."<<std::endl;
    if(strhit_array){
      if(verbose) std::cout<<"strhit_array ok with "<<strhit_array->GetEntriesFast()<<" entries"<<std::endl;
    }
    else std::cout<<"*** strhit_array broken - check your file! ."<<std::endl;
    if(pixhit_array){
      if(verbose) std::cout<<"pixhit_array ok with "<<pixhit_array->GetEntriesFast()<<" entries"<<std::endl;
    }
    else std::cout<<"*** pixhit_array broken - check your file! ."<<std::endl;

    // ----- MC Points -----
    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());
      mommc.SetXYZ(point->GetPx(),point->GetPy(),point->GetPz());
      hisxymc->Fill(vecmc.x(),vecmc.y());
      if(vecmc.y() > 0.) hisrzmc->Fill(vecmc.z(),vecmc.Perp());
      else hisrzmc->Fill(vecmc.z(),-1.*vecmc.Perp());
      hisde->Fill(point->GetEnergyLoss());
      hismom->Fill(mommc.Mag());
      detname = fGeoH->GetPath( point->GetDetName());

      if(verbose) cout << detname<<endl;

      geoMan->cd( detname.Data() );
      currentTransMat = geoMan->GetCurrentMatrix();
      tmpMaster[0]=vecmc.x();tmpMaster[1]=vecmc.y();tmpMaster[2]=vecmc.z();
      currentTransMat->MasterToLocal(tmpMaster,tmpLocal);
      vecmc.SetXYZ(tmpLocal[0],tmpLocal[1],tmpLocal[2]);
      hisSensorXYMC->Fill(vecmc.x(),vecmc.y());

    }


    // -----  PIXEL HITS  -----
    for (Int_t ii=0; ii<pixhit_array->GetEntriesFast(); ii++)
    {

      PndMvdHit *hit=(PndMvdHit*)pixhit_array->At(ii);
	  if(verbose) cout <<ii<< ".";
	  Int_t clustid = hit->GetRefIndex();
      if(verbose)cout<<clustid<<".";

      detname = fGeoH->GetPath( hit->GetDetName());
      geoMan->cd( detname.Data() );
      currentTransMat = geoMan->GetCurrentMatrix();
      vecs.SetXYZ(hit->GetX(), hit->GetY(), hit->GetZ());
      hisxy->Fill(vecs.x(),vecs.y());
      if(vecs.y() > 0.) hisrz->Fill(vecs.z(),vecs.Perp());
      else hisrz->Fill(vecs.z(),-1.*vecs.Perp());


      PndMvdCluster *clust = (PndMvdCluster*)pixclust_array->At(clustid);
      Int_t clsize = clust->GetClusterSize();
      Int_t hclsize = hit->GetNDigiHits();
      if (hclsize < clsize)
        cout<<"Strange cluster sizes - this shall not happen!"<<endl;
      Int_t sid = clust->GetDigiIndex(0);
      if(verbose)cout<<sid<<".";
      if(sid>digiPixel_array->GetEntriesFast()) cout<<"Exceeding digi array size"<<endl;
      hisClustSize->Fill(clsize);

      PndMvdDigiPixel *apixeldigi = (PndMvdDigiPixel*)digiPixel_array->At(sid);
      if(0==apixeldigi)cout<<"no pixel digi found"<<endl;
      Int_t mcid = apixeldigi->GetIndex();
      if(verbose)cout<<mcid<<"  ";
      if(mcid<0 || mcid >= mc_array->GetEntriesFast()) continue;

      PndMvdMCPoint *point=(PndMvdMCPoint*)mc_array->At(mcid);
      Int_t mcpdg = -1;

      vecmc.SetXYZ( 0.5 * (point->GetX() + point->GetXOut()),
                    0.5 * (point->GetY() + point->GetYOut()),
                    0.5 * (point->GetZ() + point->GetZOut()));
      vecdiff = vecmc - vecs;
      difftheta = vecmc.Theta() - vecs.Theta();
      diffphi = vecmc.Phi() - vecs.Phi();
      //convert deg to mrad
      difftheta = difftheta*1000.*TMath::Pi()/180.;
      diffphi= diffphi*1000.*TMath::Pi()/180.;

      hisDiffTheta->Fill(difftheta);
      hisDiffPhi->Fill(diffphi);


      //--- Now move vectors to the local sensor system
      tmpMaster[0]=vecs.x();tmpMaster[1]=vecs.y();tmpMaster[2]=vecs.z();
      currentTransMat->MasterToLocal(tmpMaster,tmpLocal);
      vecs.SetXYZ(tmpLocal[0],tmpLocal[1],tmpLocal[2]);
      tmpMaster[0]=vecmc.x();tmpMaster[1]=vecmc.y();tmpMaster[2]=vecmc.z();
      currentTransMat->MasterToLocal(tmpMaster,tmpLocal);
      vecmc.SetXYZ(tmpLocal[0],tmpLocal[1],tmpLocal[2]);
      vecdiff = vecmc - vecs;

      hisSensorXY->Fill(vecs.x(),vecs.y());
      hisDiffX->Fill(vecdiff.X());
      hisDiffY->Fill(vecdiff.Y());
      hisDiffZ->Fill(vecdiff.Z());
      hisDiffXY->Fill(vecdiff.X(),vecdiff.Y());
      hisDiffRZ->Fill(vecdiff.Z(),vecdiff.Perp());
    }//end for ii (pixel hits in event)
	  if(verbose) cout <<endl;



    // -----  STRIP HITS  -----
    for (Int_t iii=0; iii<strhit_array->GetEntriesFast(); iii++)
    {
      PndMvdHit *hit=(PndMvdHit*)strhit_array->At(iii);
      Int_t topclustid = hit->GetRefIndex();
      Int_t botclustid = hit->GetBotIndex();
      detname = fGeoH->GetPath( hit->GetDetName());
      geoMan->cd( detname.Data() );
      currentTransMat = geoMan->GetCurrentMatrix();

      vecs.SetXYZ(hit->GetX(), hit->GetY(), hit->GetZ());
      hisxy->Fill(vecs.x(),vecs.y());
      if(vecs.y() > 0.) hisrz->Fill(vecs.z(),vecs.Perp());
      else hisrz->Fill(vecs.z(),-1.*vecs.Perp());


      PndMvdCluster *clust = (PndMvdCluster*)strclust_array->At(topclustid);
      Int_t clsize = clust->GetClusterSize();
      Int_t hclsize = hit->GetNDigiHits();
      if (hclsize < clsize)
        cout<<"Strange cluster sizes - this shall not happen!"<<endl;
      Int_t bclsize =-1;
      if(botclustid>=0){
        PndMvdCluster *botcl = (PndMvdCluster*)strclust_array->At(botclustid);
        bclsize = botcl->GetClusterSize();
      }
      Int_t sid = clust->GetDigiIndex(0);
      if(sid>digiStrip_array->GetEntriesFast()) cout<<"Exceeding digi array size"<<endl;
      PndMvdDigiStrip *astripdigi = (PndMvdDigiStrip*)digiStrip_array->At(sid);

      if(0==astripdigi)cout<<"no strip digi found"<<endl;
      Int_t mcid = astripdigi->GetIndex();
      if(mcid<0 || mcid >= mc_array->GetEntriesFast()) continue;
      PndMvdMCPoint *point=(PndMvdMCPoint*)mc_array->At(mcid);
      Int_t mcpdg = -1;

      vecmc.SetXYZ( 0.5 * (point->GetX() + point->GetXOut()),
                    0.5 * (point->GetY() + point->GetYOut()),
                    0.5 * (point->GetZ() + point->GetZOut()));
      vecdiff = vecmc - vecs;
      difftheta = vecmc.Theta() - vecs.Theta();
      diffphi = vecmc.Phi() - vecs.Phi();
      //convert deg to mrad
      difftheta = difftheta*1000.*TMath::Pi()/180.;
      diffphi= diffphi*1000.*TMath::Pi()/180.;

      hisDiffTheta->Fill(difftheta);
      hisDiffPhi->Fill(diffphi);

      hisClustSize->Fill(clsize);
      if(bclsize>0) hisClustSize->Fill(bclsize);
      hisClustSizeTop->Fill(clsize);
      if(bclsize>0) hisClustSizeBot->Fill(bclsize);

      //--- Now move vectors to the local sensor system
      tmpMaster[0]=vecs.x();tmpMaster[1]=vecs.y();tmpMaster[2]=vecs.z();
      currentTransMat->MasterToLocal(tmpMaster,tmpLocal);
      vecs.SetXYZ(tmpLocal[0],tmpLocal[1],tmpLocal[2]);
      tmpMaster[0]=vecmc.x();tmpMaster[1]=vecmc.y();tmpMaster[2]=vecmc.z();
      currentTransMat->MasterToLocal(tmpMaster,tmpLocal);
      vecmc.SetXYZ(tmpLocal[0],tmpLocal[1],tmpLocal[2]);
      vecdiff = vecmc - vecs;

      hisSensorXY->Fill(vecs.x(),vecs.y());
      hisDiffStripX->Fill(vecdiff.X());
      hisDiffStripY->Fill(vecdiff.Y());
      hisDiffStripZ->Fill(vecdiff.Z());
      hisDiffStripXY->Fill(vecdiff.X(),vecdiff.Y());
      hisDiffStripRZ->Fill(vecdiff.Z(),vecdiff.Perp());
    }//end for iii (strip hits in event)

  }// end for j (events)


//-----
Int_t pix = 300;
Int_t a = 4, b = 3;

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

can1->cd(1);
mypad=gPad;
mypad.Divide(2,2);
mypad->cd(1);hisxymc->SetStats(false);hisxymc->DrawCopy("colz");
mypad->cd(2);hisrzmc->SetStats(false);hisrzmc->DrawCopy("colz");
mypad->cd(3);hisSensorXYMC->SetStats(false);hisSensorXYMC->DrawCopy("colz");
mypad->cd(4);gPad->SetLogy();hisde->DrawCopy();
// hismom->DrawCopy();

can1->cd(2);
mypad=gPad;
mypad.Divide(2,2);
mypad->cd(1);hisxy->SetStats(false);hisxy->DrawCopy("colz");
mypad->cd(2);hisrz->SetStats(false);hisrz->DrawCopy("colz");
mypad->cd(3);hisSensorXY->SetStats(false);hisSensorXY->DrawCopy("colz");
mypad->cd(4);
  hisClustSize->Draw();
  hisClustSizeBot->SetLineColor(kRed);
  hisClustSizeTop->SetLineColor(kBlue);
  hisClustSizeBot->Draw("sames");
  hisClustSizeTop->Draw("sames");
  can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);

can1->cd(3);
mypad=gPad;
mypad.Divide(2,2);
mypad->cd(1);hisDiffXY->SetStats(false);hisDiffXY->DrawCopy("colz");
mypad->cd(2);hisDiffRZ->SetStats(false);hisDiffRZ->DrawCopy("colz");
mypad->cd(3);hisDiffStripXY->SetStats(false);hisDiffStripXY->DrawCopy("colz");
mypad->cd(4);hisDiffStripRZ->SetStats(false);hisDiffStripRZ->DrawCopy("colz");

can1->cd(4);
hisClustSizeDisk->SetLineColor(kRed);
hisClustSizeSide->SetLineColor(kBlue);
hisClustSize->DrawCopy("");
hisClustSizeDisk->DrawCopy("sames");
hisClustSizeSide->DrawCopy("sames");
can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);

can1->cd(5);
hisDiffX->GetXaxis()->SetNdivisions(-05);
hisDiffX->DrawCopy("");
can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);

can1->cd(6);
hisDiffY->GetXaxis()->SetNdivisions(-05);
hisDiffY->DrawCopy("");
can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);

can1->cd(7);
hisDiffZ->GetXaxis()->SetNdivisions(-05);
hisDiffZ->DrawCopy("");
can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);

can1->cd(8);
hisDiffTheta->GetXaxis()->SetNdivisions(-05);
hisDiffTheta->DrawCopy("");
can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);



can1->cd(9);
hisDiffStripX->GetXaxis()->SetNdivisions(-05);
hisDiffStripX->DrawCopy("");
can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);

can1->cd(10);
hisDiffStripY->GetXaxis()->SetNdivisions(-05);
hisDiffStripY->DrawCopy("");
can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);

can1->cd(11);
hisDiffStripZ->GetXaxis()->SetNdivisions(-05);
hisDiffStripZ->DrawCopy("");
can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);

can1->cd(12);
hisDiffPhi->GetXaxis()->SetNdivisions(-05);
hisDiffPhi->DrawCopy("");
can1->Update();mypad=(TPad*)gPad;BetterStatBox(mypad);

can1->Print(picture);


  // -----   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;
  // ------------------------------------------------------------------------

}