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

  // -----  Load libraries   ------------------------------------------------
  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.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/MvdResDpm6GeV10k.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);


  TH1D* hisDiff = new TH1D("diff","",100,-0.008,0.008);
  hisDiff->SetTitle("MVD Hit Resolution (6 GeV/c DPM);#Deltax / 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;");

  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;


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


      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;

      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;

      hisDiff->Fill(vecdiff.X());
    }//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());

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

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

      hisDiff->Fill(vecdiff.X());
    }//end for iii (strip hits in event)

  }// end for j (events)


////-----
Int_t pix = 400;
Int_t a = 2, b = 2;

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

TCanvas* can2 = new TCanvas("MvdResPlot","MVD point resolution");

can2->cd();

  Double_t par[6];
  g1    = new TF1("g1","gaus",-0.008,0.008);
  g2    = new TF1("g2","gaus",-0.002,0.002);
  total = new TF1("total","gaus(0)+gaus(3)",-0.008,0.008);
  total->SetParameters(0., 0., 0.0005, 0., 0., 0.002);
  total->SetLineColor(4);
  total->SetLineWidth(2);
  total->SetLineStyle(7);
  hisDiff->Fit(g1,"R");
  hisDiff->Fit(g2,"R+");
  g1->GetParameters(&par[0]);
  g2->GetParameters(&par[3]);
  total->SetParameters(par);
  hisDiff->Fit(total,"R");

//  TF1* resTotalFit    = new TF1("resTotalFit", FitFunction, -0.01, 0.01, 6);
//  resTotalFit->SetParameters(0., 0., 0.004, 0., 0., 0.001);
//  hisDiff->Fit(resTotalFit,"R+");
  hisDiff->SetMarkerStyle(3);
  hisDiff->SetMarkerSize(0.4);
  hisDiff->DrawCopy("pe");

  gStyle->SetOptTitle(kFALSE);
  gStyle->SetOptStat(kFALSE);

  TString str = "#sigma_{1} = ";
  str += (10000*par[5]); for (int i=0;i<6;i++) str.Chop();
  str += " #mum";
  DrawText( 0.2, 0.8, str.Data(),0.05,1);
  str = "#sigma_{2} = ";
  str += (10000*par[2]); for (int i=0;i<6;i++) str.Chop();
  str += " #mum";
  DrawText( 0.2, 0.7, str.Data(),0.05,1);
  LoadPandaStyle();


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



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

can1->cd(4);
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;
  // ------------------------------------------------------------------------

}







Double_t Background(Double_t *x, Double_t *par)
{
  Bool_t reject= false;
  if ( reject && x[0] > 1.99 && x[0] < 2.03) {
    TF1::RejectPoint();
    return 0;
  }
  return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
}

Double_t Signal(Double_t *x, Double_t *par) {
  return par[0]*exp( -0.5*(pow((x[0] - par[1])/par[2],2)) );
}


Double_t FitFunction(Double_t *x, Double_t *par) {
//  return Background(x,par) + Signal(x,&par[3]) + Signal(x,&par[6]);
	return Signal(x,par) + Signal(x,&par[3]);
}