/*
This macro creates for each pattern four 2D histograms: distance between MChit and COG vs impact angle. This information are taken
from the output of compClusHits.C. Information about the pattern are taken from the database generated by BuildClTop.C.
The macro fits these histograms with gaussians along y axis, in order to study the dependence of dZ and dX wrt
the imact angles. Two 1D histograms with dZ and dX are generated as a projection on to the y axis and are used to calculate mean value-sigma
for each pattern.
The output is a PDF with 4 2D histograms (mean-value of dZ(dX) vs alpha, sigma of dZ(dX) vs alpha,mean-value of dZ(dX) vs beta,
sigma of dZ(dX) vs beta) and 2 histograms, dZ and dX, for each pattern, a drawing of the pattern and the information from the fit and the DB and
information from DB and fits.
You can set the range of patterns you want to study by setting LowerLimit/UpperLimit. If the range coincide with the the whole ragne of patterns,
the DB is updated with the information from these fits.
It is possible to loop over few clusters instead of all the cluster setting ntotCl: usefull to check the macro. (-1 means all the clusters).
It is possilble toprint a PDF with the 2D histograms, pattern drawing and info: setting number2d histos you decide the number of pattern for which
you want 2D histograms.
*/


#if !defined(__CINT__) || defined(__MAKECINT__)
#include "TObjArray.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TArrayI.h"
#include "TArrayF.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TPad.h"
#include "TLatex.h"
#include "TBits.h"
#include "TStopwatch.h"
#include "TLine.h"
#include "TMarker.h"
#include "../ITS/UPGRADE/AliITSUClusterPix.h"
#include "../ITS/UPGRADE/AliITSURecoLayer.h"
#include "../ITS/UPGRADE/AliITSURecoDet.h"
#include "../ITS/UPGRADE/AliITSUHit.h"
#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
#include "AliITSsegmentation.h"
#include "AliGeomManager.h"
#include "TROOT.h"
#include "TStyle.h"

#endif

TObjArray* pattDB=0; //it is an array with all the patterns in TBits format (forom the most to the least frequent)
TVectorF* pattFR=0; //frequency of the pattern in the database
TVectorF* xCentrePix=0; //coordinate of the centre of the pixel containing the COG for the down-left corner in fracion of the pitch
TVectorF* zCentrePix=0;
TVectorF* xCentreShift=0; //Difference between x coordinate fo COG and the centre of the pixel containing the COG, in fraction of the pitch
TVectorF* zCentreShift=0;
TVectorF* NPix=0; //Number of fired pixels
TVectorF* NRow=0; //Number of rows of the pattern
TVectorF* NCol=0; //Number of columns of the pattern

TArrayF DeltaZmean; //mean value of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
TArrayF DeltaZmeanErr;
TArrayF DeltaXmean; //mean value of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
TArrayF DeltaXmeanErr;
TArrayF DeltaZsigma; //sigma of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
TArrayF DeltaZsigmaErr;
TArrayF DeltaXsigma; //sigma of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
TArrayF DeltaXsigmaErr;
TArrayF Chi2x; //Chi2 of the gaussian fit over x
TArrayF Chi2z; //Chi2 of the aussian fut over z
TArrayI NDFx; 
TArrayI NDFz;

TObjArray histoArr;

TStopwatch timer;
TStopwatch timer1;

Float_t pitchX; //pitch of the pixel in x direction
Float_t pitchZ; //pitch of the pixel in z direction

enum{kDXAlpha=0, kDZAlpha=1, kDXBeta=2, kDZBeta=3};//to select the kind of TH2 we want

void LoadDB(const char* fname);
void CreateHistos(TObjArray* harr);
TH2* GetHisto(Int_t pattID, Int_t kind, TObjArray* harr);
void FitHistos(Int_t pattID, Float_t Freq, Int_t nPix, Int_t nRow, Int_t nCol,
 TObjArray* harr, TCanvas* canvas, TCanvas* canvasbis);
void UpDateDB(const char* outDBFile);

typedef struct {
  Int_t evID;
  Int_t volID;
  Int_t lrID;
  Int_t clID;
  Int_t nPix;
  Int_t nX;
  Int_t nZ;
  Int_t q;
  Float_t pt;
  Float_t eta;
  Float_t phi;
  Float_t xyz[3];
  Float_t dX;
  Float_t dY;
  Float_t dZ;  
  Bool_t split;  
  Bool_t prim;
  Int_t  pdg;
  Int_t  ntr;
  Float_t alpha; // alpha is the angle in y-radius plane in local frame
  Float_t beta;  // beta is the angle in xz plane, taken from z axis, growing counterclockwise
  Int_t nRowPatt;
  Int_t nColPatt;
  Int_t pattID;
  Float_t freq;
  Float_t xCen;
  Float_t zCen;
  Float_t zShift;
  Float_t xShift;
} clSumm;

static clSumm Summ;

Int_t nPatterns=0;
Int_t LowerLimit=0;//set the first cluster of the range of interest, INCLUDED (to update DB LowerLimit=0)
Int_t UpperLimit=4000; // EXCLUDED; -1 means number of Patterns. In this case the database is updated with the information about DeltaxMean, DeltaXerr, etc.
Int_t nStudied=0;
Int_t ntotCl=-1;//Set -1 to loop over all the clusters, otherwise set the number of clusters
Int_t number2dhisto=10;//Set the number pattern from 0 to number2dhisto-1 for which you want to print the TH2 histos (dZ or dX vs alpha or beta)

void complete1(){

  //importing data

  LoadDB("clusterTopology.root");

  AliGeomManager::LoadGeometry("geometry.root");
  AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
  AliITSUClusterPix::SetGeom(gm);
  const AliITSsegmentation* segm = gm->GetSegmentation(0);
  pitchX = segm->Dpx(0)*10000; //pitch in X in micron
  printf("pitchX: %f\n",pitchX);
  pitchZ = segm->Dpz(0)*10000; //pitch in Z in micron
  printf("pitchX: %f\n",pitchX);
  delete gm;

  nPatterns = (pattDB->GetEntriesFast());

  if(UpperLimit==-1){
    UpperLimit=nPatterns;
  }

  nStudied=UpperLimit-LowerLimit+1;

  printf("\n\n nPatterns %d\n\n",nPatterns);

  DeltaZmean.Set(nPatterns+100);
  DeltaZmeanErr.Set(nPatterns+100);
  DeltaXmean.Set(nPatterns+100);
  DeltaXmeanErr.Set(nPatterns+100);
  DeltaZsigma.Set(nPatterns+100);
  DeltaZsigmaErr.Set(nPatterns+100);
  DeltaXsigma.Set(nPatterns+100);
  DeltaXsigmaErr.Set(nPatterns+100);
  Chi2x.Set(nPatterns+100);
  Chi2z.Set(nPatterns+100);
  NDFx.Set(nPatterns+100);
  NDFz.Set(nPatterns+100);

  printf("Loading input tree... ");

  TFile *input = new TFile("clInfo.root","read");
  TTree *clitsu = (TTree*) input->Get("clitsu");
  
  if(ntotCl==-1){
    ntotCl = (Int_t) clitsu->GetEntries();
  }
  TBranch *brdZ = clitsu->GetBranch("dZ");
  brdZ->SetAddress(&Summ.dZ);
  TBranch *brdX = clitsu->GetBranch("dX");
  brdX->SetAddress(&Summ.dX);
  TBranch *brpattID = clitsu->GetBranch("pattID");
  brpattID->SetAddress(&Summ.pattID);
  TBranch *brfreq = clitsu->GetBranch("freq");
  brfreq->SetAddress(&Summ.freq);
  TBranch *bralpha = clitsu->GetBranch("alpha");
  bralpha->SetAddress(&Summ.alpha);
  TBranch *brbeta = clitsu->GetBranch("beta");
  brbeta->SetAddress(&Summ.beta);
  printf("done!!\n\n");

  CreateHistos(&histoArr);
  printf("Starting to process clusters\n\n");
  timer.Start();
  for(Int_t iCl=0; iCl<ntotCl; iCl++){
    clitsu->GetEntry(iCl);
    Int_t ID = Summ.pattID;
    printf("Processing cluster %d... ",iCl);
    if(ID<LowerLimit || ID>=UpperLimit) {
      printf("skipped!\n\n");
      continue;
    }
    GetHisto(ID,kDXAlpha,&histoArr)->Fill(Summ.alpha,Summ.dX);
    GetHisto(ID,kDZAlpha,&histoArr)->Fill(Summ.alpha,Summ.dZ);
    GetHisto(ID,kDXBeta,&histoArr)->Fill(Summ.beta,Summ.dX);
    GetHisto(ID,kDZBeta,&histoArr)->Fill(Summ.beta,Summ.dZ);
    printf("done!\n\n");
  }
  printf("All clusters processed!\n");
  timer.Stop();
  timer.Print();

  TCanvas* cnv1 = new TCanvas("cnv1","cnv1",900,700);
  TCanvas* cnv2 = new TCanvas("cnv2","cnv2",900,700);

  printf("Printing PDF...\n\n");
  timer1.Start();
  
  cnv1->Print("angles.pdf[");
  cnv2->Print("angles2D.pdf[");

  for(Int_t i=LowerLimit; i<UpperLimit; i++){
    
    FitHistos(i,(*pattFR)[i],(*NPix)[i],(*NRow)[i],(*NCol)[i],&histoArr, cnv1, cnv2);
  }

  cnv1->Print("angles.pdf]");
  cnv2->Print("angles2D.pdf]");
  timer1.Stop();
  printf("\n\ndone!!\n\n");
  timer1.Print();

  delete cnv1;
  delete cnv2;

  if(LowerLimit==0 && UpperLimit==nPatterns){
    UpDateDB("clusterTopology2.root");
  }

}

void LoadDB(const char* fname){ //load the data base built by BuildClTopDB.C

  // load database
  TFile* fl = TFile::Open(fname);
  if(!fl){printf("Could not find %s",fname); exit(1);}
  pattDB = (TObjArray*)fl->Get("TopDB");
  pattFR = (TVectorF*)fl->Get("TopFreq");
  xCentrePix =(TVectorF*)fl->Get("xCOG");
  zCentrePix =(TVectorF*)fl->Get("zCOG");
  xCentreShift =(TVectorF*)fl->Get("xShift");
  zCentreShift =(TVectorF*)fl->Get("zShift");
  NPix =(TVectorF*)fl->Get("NPix");
  NCol =(TVectorF*)fl->Get("NCol");
  NRow =(TVectorF*)fl->Get("NRow");
}

void CreateHistos(TObjArray* harr){ //creates four 2D histograms (dX or dZ vs alpha or beta) for each pattern included in the study range

  printf("Creating histograms... ");

  //for each pattern we crate 4 2D histograms corresponding the cathegories within the enumeration
  if(!harr) harr=&histoArr;
    
  for(Int_t i=0; i<nStudied; i++){

    TH2* h0 = new TH2F(Form("hXalpha%d", i),
      Form("#DeltaX vs #alpha for pattID %d", i),10,-0.1, 1.1* TMath::Pi()/2,100,-30,30);
    h0->SetDirectory(0);
    harr->AddAtAndExpand(h0, i*4);//This histogram is the first of the 4-histos series
    
    TH2* h1 = new TH2F(Form("hZalpha%d", i),
      Form("#DeltaZ vs #alpha for pattID %d", i),10, -0.1, 1.1* TMath::Pi()/2,100,-30,30);
    h1->SetDirectory(0);
    harr->AddAtAndExpand(h1, i*4+1);//This histogram is the second of the 4-histos series and so on...

    TH2* h2 = new TH2F(Form("hXbeta%d", i),
      Form("#DeltaX vs #beta for pattID %d", i),10,-0.1, 1.1* TMath::Pi()/2,100,-30,30);
    h2->SetDirectory(0);
    harr->AddAtAndExpand(h2, i*4+2);

    TH2* h3 = new TH2F(Form("hZbeta%d", i),
      Form("#DeltaZ vs #beta for pattID %d", i),10,-0.1, 1.1* TMath::Pi()/2,100,-30,30);
    h3->SetDirectory(0);
    harr->AddAtAndExpand(h3, i*4+3);

  }

  printf("done!!\n\n");
}

TH2* GetHisto(Int_t pattID, Int_t kind, TObjArray* harr){

  if(!harr) harr=&histoArr ;
  TH2*  h=0;
  h=(TH2*)harr->At((pattID-LowerLimit)*4+kind);
  if (!h) {
    printf("Unknown histo id=%d\n",pattID); 
    exit(1);
  }
  return h;

}

void FitHistos(Int_t pattID, Float_t Freq, Int_t nPix, Int_t nRow, 
 Int_t nCol, TObjArray* harr, TCanvas* canvas, TCanvas* canvasbis){ //Fit the 2D histograms with gaussian along the Y and creates two 1D histograms,
                                                                    // with dX and dZ, bay a projection to the y axis and the fit with gaussian  

  static TF1* gs = new TF1("gs","gaus",-50,50);
  
  if(!harr) harr=&histoArr;

  //Creating the pads for the output
  
  canvas->Clear();
  canvas->cd();

  Float_t width = (1 - canvas->GetRightMargin() - canvas->GetLeftMargin() )/3;

  TPad* titlepad = new TPad("titlepad","",canvas->GetLeftMargin(),0.825,1-canvas->GetRightMargin(),0.95);
  titlepad->SetFillColor(kYellow);
  titlepad->Draw();
  TPad* pad1 = new TPad("pad1","",canvas->GetLeftMargin(),(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin(),canvas->GetLeftMargin()+width,0.8);
  pad1->Draw();
  TPad* pad2 = new TPad("pad2","",canvas->GetLeftMargin()+width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin(),canvas->GetLeftMargin()+2*width,0.8);
  pad2->Draw();
  TPad* pad3 = new TPad("pad3","",canvas->GetLeftMargin()+2*width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin(),1-canvas->GetRightMargin(),0.8);
  pad3->Draw();
  TPad* pad4 = new TPad("pad3","",canvas->GetLeftMargin(),canvas->GetBottomMargin(),canvas->GetLeftMargin()+width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin());
  pad4->Draw();
  TPad* pad5 = new TPad("pad5","",canvas->GetLeftMargin()+width,canvas->GetBottomMargin(),canvas->GetLeftMargin()+2*width,(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin());
  pad5->Draw();
  TPad* pad6 = new TPad("pad6","",canvas->GetLeftMargin()+2*width,canvas->GetBottomMargin(),1-canvas->GetRightMargin(),(0.8-canvas->GetBottomMargin())/2+canvas->GetBottomMargin());
  pad6->Draw();
  TPad* padPattern = new TPad("padPattern","",0.75,0.825,0.9,0.95);
  padPattern->Draw();
  
  //...................................................................................................

  // dX vs alpha
  
  TObjArray harrFXA;
  TH2* hXA = (TH2*)harr->At((pattID-LowerLimit)*4);
  hXA->GetYaxis()->SetTitle("#DeltaX, #mum");
  hXA->GetXaxis()->SetTitle("#alpha");
  hXA->SetStats(0);
  hXA->FitSlicesY(gs,0,-1,0,"qnr",&harrFXA);
  TH1* hXA_1 = (TH1*)harrFXA.At(1);
  hXA_1->SetStats(0);
  hXA_1->SetDirectory(0);
  hXA_1->GetYaxis()->SetTitle("#DeltaX_{#mu}(#DeltaZ_{#mu}), #mum");
  TH1* hXA_2 = (TH1*)harrFXA.At(2);
  hXA_2->SetStats(0);
  hXA_2->SetDirectory(0);
  hXA_2->GetYaxis()->SetTitle("#DeltaX_{#sigma}(#DeltaZ_{#sigma}), #mum");
  
  //............................................................................................................

  // dX

  Int_t fitStatusX=0;

  TH1* hdx = hXA->ProjectionY("hdx");
  hdx->GetXaxis()->SetTitle("#DeltaX, #mum");
  hdx->SetTitle(Form("#DeltaX for pattID %d",pattID));
  hdx->SetStats(0);
  pad3->cd();
  gs->SetParameters(hdx->GetMaximum(),hdx->GetMean(),hdx->GetRMS());
  if((hdx->GetEntries())<100) fitStatusX = hdx->Fit("gs","ql");
  else fitStatusX = hdx->Fit("gs","q");

  if(fitStatusX==0){
    DeltaXmean[pattID]= gs->GetParameter(1);
    DeltaXmeanErr[pattID]= gs->GetParError(1);
    DeltaXsigma[pattID]= gs->GetParameter(2);
    DeltaXsigmaErr[pattID]= gs->GetParError(2);
    Chi2x[pattID] = gs->GetChisquare();
    Int_t varNDFx = gs->GetNDF();
    if(varNDFx>=0)
    NDFx[pattID] = varNDFx;
    else{
      NDFx[pattID]=-1;
    }
  }
  else{
    DeltaXmean[pattID]=0;
    DeltaXmeanErr[pattID]=0;
    DeltaXsigma[pattID]=0;
    DeltaXsigmaErr[pattID]=0;
    Chi2x[pattID] = -1;   
  }

  pad3->cd();
  hdx->Draw();


  //.............................................................................................................

  
  // dZ vs alpha

  TObjArray harrFZA;
  TH2* hZA = (TH2*)harr->At((pattID-LowerLimit)*4+1);
  hZA->GetYaxis()->SetTitle("#DeltaZ, #mum");
  hZA->GetXaxis()->SetTitle("#alpha");
  hZA->SetStats(0);
  hZA->FitSlicesY(gs,0,-1,0,"qnr",&harrFZA);
  TH1* hZA_1 = (TH1*)harrFZA.At(1);
  hZA_1->SetMarkerColor(8);
  hZA_1->SetLineColor(8);
  hZA_1->SetStats(0);
  TH1* hZA_2 = (TH1*)harrFZA.At(2);
  hZA_2->SetMarkerColor(8);
  hZA_2->SetLineColor(8);
  hZA_2->SetStats(0);
  
  Double_t maxA_1 = TMath::Max(hXA_1->GetMaximum(),hZA_1->GetMaximum());
  Double_t minA_1 = TMath::Min(hXA_1->GetMinimum(),hZA_1->GetMinimum());
  Double_t addtorangeA_1 = (maxA_1-minA_1)*0.1;
  Double_t rangeA_1 = addtorangeA_1*10;

  Double_t maxA_2 = TMath::Max(hXA_2->GetMaximum(),hZA_2->GetMaximum());
  Double_t minA_2 = TMath::Min(hXA_2->GetMinimum(),hZA_2->GetMinimum());
  Double_t addtorangeA_2 = (maxA_2-minA_2)*0.1;
  Double_t rangeA_2 = addtorangeA_2*10;

  hXA_1->SetMaximum(maxA_1+addtorangeA_1);
  hXA_1->SetMinimum(minA_1-addtorangeA_1);

  hXA_2->SetMaximum(maxA_2+addtorangeA_2);
  hXA_2->SetMinimum(minA_2-addtorangeA_2);

  pad1->cd();
  hXA_1->Draw();
  hZA_1->Draw("same");
  TLegend* legA_1 = new TLegend(0.8*1.1* TMath::Pi()/2, maxA_1-0.2*rangeA_1, TMath::Pi()/2, maxA_1,"","");
  legA_1->AddEntry(hXA_1, "   #DeltaX_{#mu}", "l");
  legA_1->AddEntry(hZA_1, "   #DeltaZ_{#mu}", "l");
  legA_1->SetFillColor(0);
  legA_1->SetBorderSize(0);
  legA_1->Draw();
  pad2->cd();
  hXA_2->Draw();
  hZA_2->Draw("same");
  TLegend* legA_2 = new TLegend(0.8*1.1* TMath::Pi()/2, maxA_2-0.2*rangeA_2, TMath::Pi()/2, maxA_2,"","");
  legA_2->AddEntry(hXA_2, "   #DeltaX_{#sigma}", "l");
  legA_2->AddEntry((TObject*)0, "", "");
  legA_2->AddEntry(hZA_2, "   #DeltaZ_{#sigma}", "l");
  legA_2->SetFillColor(0);
  legA_2->SetBorderSize(0);
  legA_2->Draw();
 
  //...........................................................................................


  // dZ

  Int_t fitStatusZ=0;

  TH1* hdz = hZA->ProjectionY("hdz");
  hdz->GetXaxis()->SetTitle("#DeltaZ, #mum");
  hdz->SetTitle(Form("#DeltaZ for pattID %d",pattID));
  hdz->SetStats(0);
  pad6->cd();
  gs->SetParameters(hdz->GetMaximum(),hdz->GetMean(),hdz->GetRMS());
  if((hdz->GetEntries())<100) fitStatusZ = hdz->Fit("gs","ql");
  else fitStatusZ = hdz->Fit("gs","q");

  if(fitStatusZ==0){
    DeltaZmean[pattID]= gs->GetParameter(1);
    DeltaZmeanErr[pattID]= gs->GetParError(1);
    DeltaZsigma[pattID]= gs->GetParameter(2);
    DeltaZsigmaErr[pattID]= gs->GetParError(2);
    Chi2z[pattID] = gs->GetChisquare();
    Int_t varNDFz = gs->GetNDF();
    if(varNDFz>=0)
    NDFz[pattID] = varNDFz;
    else{
    NDFz[pattID]=-1;
    }
  }
  else{
    DeltaZmean[pattID]=0;
    DeltaZmeanErr[pattID]=0;
    DeltaZsigma[pattID]=0;
    DeltaZsigmaErr[pattID]=0;
    Chi2z[pattID] = -1;   
  }

  pad6->cd();
  hdz->Draw();

  //.............................................................................................................

  // dX vs beta

  TObjArray harrFXB;
  TH2* hXB = (TH2*)harr->At((pattID-LowerLimit)*4+2);
  hXB->GetYaxis()->SetTitle("#DeltaX, #mum");
  hXB->GetXaxis()->SetTitle("#beta");
  hXB->SetStats(0);
  hXB->FitSlicesY(gs,0,-1,0,"qnr",&harrFXB);
  TH1* hXB_1 = (TH1*)harrFXB.At(1);
  hXB_1->SetDirectory(0);
  hXB_1->SetStats(0);
  hXB_1->GetYaxis()->SetTitle("#DeltaX_{#mu}(#DeltaZ_{#mu}), #mum");
  TH1* hXB_2 = (TH1*)harrFXB.At(2);
  hXB_2->SetDirectory(0);
  hXB_2->SetStats(0);
  hXB_2->GetYaxis()->SetTitle("#DeltaX_{#sigma}(#DeltaZ_{#sigma}), #mum");
  
  //............................................................................................

  //dZ vs beta

  TObjArray harrFZB;
  TH2* hZB = (TH2*)harr->At((pattID-LowerLimit)*4+3);
  hZB->GetYaxis()->SetTitle("#DeltaZ, #mum");
  hZB->GetXaxis()->SetTitle("#alpha");
  hZB->SetStats(0);
  hZB->FitSlicesY(gs,0,-1,0,"qnr",&harrFZB);
  TH1* hZB_1 = (TH1*)harrFZB.At(1);
  hZB_1->SetMarkerColor(8);
  hZB_1->SetLineColor(8);
  hZB_1->SetStats(0);
  TH1* hZB_2 = (TH1*)harrFZB.At(2);
  hZB_2->SetMarkerColor(8);
  hZB_2->SetLineColor(8);
  hZB_2->SetStats(0);

  Double_t maxB_1 = TMath::Max(hXB_1->GetMaximum(),hZB_1->GetMaximum());
  Double_t minB_1 = TMath::Min(hXB_1->GetMinimum(),hZB_1->GetMinimum());
  Double_t addtorangeB_1 = (maxB_1-minB_1)*0.1;
  Double_t rangeB_1 = 10*addtorangeB_1;

  Double_t maxB_2 = TMath::Max(hXB_2->GetMaximum(),hZB_2->GetMaximum());
  Double_t minB_2 = TMath::Min(hXB_2->GetMinimum(),hZB_2->GetMinimum());
  Double_t addtorangeB_2 = (maxB_2-minB_2)*0.1;
  Double_t rangeB_2 = 10*addtorangeB_2;


  hXB_1->SetMaximum(maxB_1+addtorangeB_1);
  hXB_1->SetMinimum(minB_1-addtorangeB_1);

  hXB_2->SetMaximum(maxB_2+addtorangeB_2);
  hXB_2->SetMinimum(minB_2-addtorangeB_2);

  pad4->cd();
  hXB_1->Draw();
  hZB_1->Draw("same");
  TLegend* legB_1 = new TLegend(0.8*1.1* TMath::Pi()/2, maxB_1-0.2*rangeB_1, TMath::Pi()/2, maxB_1,"","");
  legB_1->AddEntry(hXA_1, "   #DeltaX_{#mu}", "l");
  legB_1->AddEntry((TObject*)0, "", "");
  legB_1->AddEntry(hZA_1, "   #DeltaZ_{#mu}", "l");
  legB_1->SetFillColor(0);
  legB_1->SetBorderSize(0);
  legB_1->Draw();
  pad5->cd();
  hXB_2->Draw();
  hZB_2->Draw("same");
  TLegend* legB_2 = new TLegend(0.8*1.1* TMath::Pi()/2, maxB_2-0.2*rangeB_2, TMath::Pi()/2, maxB_2,"","");
  legB_2->AddEntry(hXB_2, "   #DeltaX_{#sigma}", "l");
  legB_2->AddEntry((TObject*)0, "", "");
  legB_2->AddEntry(hZB_2, "   #DeltaZ_{#sigma}", "l");
  legB_2->SetFillColor(0);
  legB_2->SetBorderSize(0);
  legB_2->Draw();
  
  //cnv->Print("patt0.gif");

  //..................................................................................................................

  // Drawing the title

  const char* text = Form("pattID: %d   freq: %f   nPix: %d   nRow: %d   nCol: %d  xCOG (fraction): %0.2f + (%0.2f)   zCOG (fraction): %0.2f + (%0.2f)",
    pattID,Freq,nPix,nRow,nCol, (*xCentrePix)[pattID],(*xCentreShift)[pattID],(*zCentrePix)[pattID],
    (*zCentreShift)[pattID]);
  const char* text1 = Form("#DeltaX_{#mu}: (%.2f #pm %.2f) #mum     #DeltaX_{#sigma}: (%.2f #pm %.2f) #mum     #DeltaX_{MC-centre}: %.2f #mum     #chi^{2}/NDF: %.2f / %d",
    DeltaXmean[pattID],DeltaXmeanErr[pattID],DeltaXsigma[pattID],DeltaXsigmaErr[pattID], DeltaXmean[pattID]+(*xCentreShift)[pattID]*pitchX, Chi2x[pattID], NDFx[pattID]);
  const char* text2 = Form("#DeltaZ_{#mu}: (%.2f #pm %.2f) #mum     #DeltaZ_{#sigma}: (%.2f #pm %.2f) #mum     #DeltaZ_{MC-centre}: %.2f #mum     #chi^{2}/NDF: %.2f / %d",
    DeltaZmean[pattID],DeltaZmeanErr[pattID],DeltaZsigma[pattID],DeltaZsigmaErr[pattID], DeltaZmean[pattID]+(*zCentreShift)[pattID]*pitchZ, Chi2z[pattID], NDFz[pattID]);
  TLatex title;
  title.SetTextSize(0.1);
  title.SetTextAlign(12);
  titlepad->cd();
  title.DrawLatex(0.085,0.75,text);
  title.DrawLatex(0.085,0.5,text1);
  title.DrawLatex(0.085,0.25,text2);


  //Drawing the pattern

  TH2* patternIm = new TH2F("patternIm","", nCol,-0.5,nCol-0.5,nRow,-0.5,nRow-0.5);
  for (Int_t ir=0; ir<nRow; ir++) {
    for (Int_t ic=0; ic<nCol; ic++) {
      if(((TBits*)(*pattDB)[pattID])->TestBitNumber(ir*nCol+ic)){
        patternIm->Fill(ic,ir);
      }
    }
  }

  patternIm->SetStats(0);
  padPattern->cd();
  patternIm->Draw("Acol");

  //Drawing vertical lines on the pattern

  TObjArray vertlinesarray;
  for(int i=2; i<=nCol; i++){
    TLine* vline = new TLine(patternIm->GetXaxis()->GetBinLowEdge(i),patternIm->GetYaxis()->GetBinLowEdge(1),
      patternIm->GetXaxis()->GetBinLowEdge(i),patternIm->GetYaxis()->GetBinLowEdge(nRow+1));
    vertlinesarray.AddAtAndExpand(vline,i);
  }

  for(int i=2; i<=nCol; i++){
    TLine* v1 = (TLine*)vertlinesarray.At(i);
    padPattern->cd();
    v1->Draw();
  }

  //Drawing horizontal lines on the pattern
  
  TObjArray horlinesarray;
  for(int i=2; i<=nRow; i++){
    TLine* hline = new TLine(patternIm->GetXaxis()->GetBinLowEdge(1),patternIm->GetYaxis()->GetBinLowEdge(i),
      patternIm->GetXaxis()->GetBinLowEdge(nCol+1),patternIm->GetYaxis()->GetBinLowEdge(i));
    horlinesarray.AddAtAndExpand(hline,i);
  }

  for(int i=2; i<=nRow; i++){
    TLine* hor1 = (TLine*)horlinesarray.At(i);
    padPattern->cd();
    hor1->Draw();
  }

  //Drawing MChit and COG

  TMarker* COG = new TMarker((*zCentrePix)[pattID]+(*zCentreShift)[pattID]-0.5,
    (*xCentrePix)[pattID]+(*xCentreShift)[pattID]-0.5,20);
  COG->Draw();

  TMarker* MChit = new TMarker((*zCentrePix)[pattID]+(*zCentreShift)[pattID]+(DeltaZmean[pattID])/pitchZ-0.5, //20 is the pixel pitch
    (*xCentrePix)[pattID]+(*xCentreShift)[pattID]+(DeltaXmean[pattID])/pitchX-0.5,29);
  MChit->SetMarkerColor(kWhite);
  MChit->Draw();

  canvas->Print("angles.pdf");
  canvas->Clear();

  if(pattID<number2dhisto){

    canvasbis->cd();
    TPad* titlepadbis = new TPad("titlepadbis","",canvas->GetLeftMargin(),0.825,1-canvas->GetRightMargin(),0.95);
    titlepadbis->SetFillColor(kYellow);
    titlepadbis->Draw();
    TPad* pad1bis = new TPad("pad1bis","",canvasbis->GetLeftMargin(),(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin(),0.5,0.8);
    pad1bis->Draw();
    TPad* pad2bis = new TPad("pad2bis","",0.5,(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin(),1-canvasbis->GetRightMargin(),0.8);
    pad2bis->Draw();
    TPad* pad3bis = new TPad("pad3bis","",canvasbis->GetLeftMargin(),canvasbis->GetBottomMargin(),0.5,(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin());
    pad3bis->Draw();
    TPad* pad4bis = new TPad("pad4bis","",0.5,canvasbis->GetBottomMargin(),1-canvasbis->GetRightMargin(),(0.8-canvasbis->GetBottomMargin())/2+canvasbis->GetBottomMargin());
    pad4bis->Draw();
    TPad* padPatternbis = new TPad("padPatternbis","",0.75,0.825,0.9,0.95);
    padPatternbis->Draw();

    titlepadbis->cd();
    title.DrawLatex(0.085,0.75,text);
    title.DrawLatex(0.085,0.5,text1);
    title.DrawLatex(0.085,0.25,text2);
    pad1bis->cd();
    hXA->Draw("colz");
    pad2bis->cd();
    hXB->Draw("colz");
    pad3bis->cd();
    hZA->Draw("colz");
    pad4bis->cd();
    hZB->Draw("colz");
    padPatternbis->cd();
    patternIm->Draw("Acol");
    for(int i=2; i<=nCol; i++){
      TLine* v1 = (TLine*)vertlinesarray.At(i);
      padPatternbis->cd();
      v1->Draw();
    }
    for(int i=2; i<=nRow; i++){
      TLine* hor1 = (TLine*)horlinesarray.At(i);
      padPatternbis->cd();
      hor1->Draw();
    }
    COG->Draw();
    MChit->Draw();
    canvasbis->Print("angles2D.pdf");
    canvasbis->Clear();
  }

  delete patternIm;
  delete COG;
  delete MChit;

  horlinesarray.Delete();
  vertlinesarray.Delete();

}

void UpDateDB(const char* outDBFile){

  printf("\n\n\nStoring data in the DataBase\n\n\n");

  TString outDBFileName = outDBFile;
  if (outDBFileName.IsNull()) outDBFileName = "clusterTopology2.root";

  TVectorF arrDeltaZmean(nPatterns);
  TVectorF arrDeltaZmeanErr(nPatterns);
  TVectorF arrDeltaXmean(nPatterns);
  TVectorF arrDeltaXmeanErr(nPatterns);
  TVectorF arrDeltaZsigma(nPatterns);
  TVectorF arrDeltaZsigmaErr(nPatterns);
  TVectorF arrDeltaXsigma(nPatterns);
  TVectorF arrDeltaXsigmaErr(nPatterns);
  TVectorF arrChi2x(nPatterns);
  TVectorF arrNDFx(nPatterns);
  TVectorF arrChi2z(nPatterns);
  TVectorF arrNDFz(nPatterns);

  for(Int_t ID=0; ID<nPatterns; ID++){

    printf("Processing pattern %d... ", ID);

    arrDeltaZmean[ID]=DeltaZmean[ID];
    arrDeltaZmeanErr[ID]=DeltaZmeanErr[ID];
    arrDeltaXmean[ID]=DeltaXmean[ID];
    arrDeltaXmeanErr[ID]=DeltaXmeanErr[ID];
    arrDeltaZsigma[ID]=DeltaZsigma[ID];
    arrDeltaZsigmaErr[ID]=DeltaZsigmaErr[ID];
    arrDeltaXsigma[ID]=DeltaXsigma[ID];
    arrDeltaXsigmaErr[ID]=DeltaXsigmaErr[ID];
    arrChi2x[ID]=Chi2x[ID];
    arrNDFx[ID]=NDFx[ID];
    arrChi2z[ID]=Chi2z[ID];
    arrNDFz[ID]=NDFz[ID];

    printf("done!!\n\n");
  }

  TFile* flDB = TFile::Open(outDBFileName.Data(), "recreate");
  flDB->WriteObject(pattDB,"TopDB","kSingleKey");
  flDB->WriteObject(NPix,"NPix","kSingleKey");
  flDB->WriteObject(NRow,"NRow","kSingleKey");
  flDB->WriteObject(NCol,"NCol","kSingleKey");
  flDB->WriteObject(pattFR, "TopFreq","kSingleKey");
  flDB->WriteObject(xCentrePix,"xCOG","kSingleKey");
  flDB->WriteObject(xCentreShift,"xShift","kSingleKey");
  flDB->WriteObject(zCentrePix,"zCOG","kSingleKey");
  flDB->WriteObject(zCentreShift,"zShift","kSingleKey");
  flDB->WriteObject(&arrDeltaZmean,"DeltaZmean","kSingleKey");
  flDB->WriteObject(&arrDeltaZmeanErr,"DeltaZmeanErr","kSingleKey");
  flDB->WriteObject(&arrDeltaXmean,"DeltaXmean","kSingleKey");
  flDB->WriteObject(&arrDeltaXmeanErr,"DeltaXmeanErr","kSingleKey");
  flDB->WriteObject(&arrDeltaZsigma,"DeltaZsigma","kSingleKey");
  flDB->WriteObject(&arrDeltaZsigmaErr,"DeltaZsigmaErr","kSingleKey");
  flDB->WriteObject(&arrDeltaXsigma,"DeltaXsigma","kSingleKey");
  flDB->WriteObject(&arrDeltaXsigmaErr,"DeltaXsigmaErr","kSingleKey");
  flDB->WriteObject(&arrChi2x,"Chi2x","kSingleKey");
  flDB->WriteObject(&arrChi2z,"Chi2z","kSingleKey");
  flDB->WriteObject(&arrNDFx,"NDFx","kSingleKey");
  flDB->WriteObject(&arrNDFz,"NDFz","kSingleKey");
  //
  flDB->Close();
  delete flDB;

  printf("\n\nDB Complete!!\n\n");