macros developed by Luca Barioglio for patterns analysis
authorshahoian <ruben.shahoyan@cern.ch>
Fri, 12 Sep 2014 09:18:10 +0000 (11:18 +0200)
committershahoian <ruben.shahoyan@cern.ch>
Fri, 12 Sep 2014 09:18:24 +0000 (11:18 +0200)
ITS/UPGRADE/compression/PatternAnalysis/BuildClTopDB.C [new file with mode: 0644]
ITS/UPGRADE/compression/PatternAnalysis/README [new file with mode: 0644]
ITS/UPGRADE/compression/PatternAnalysis/SecondcompClusHits.C [new file with mode: 0644]
ITS/UPGRADE/compression/PatternAnalysis/compClusHits.C [new file with mode: 0644]
ITS/UPGRADE/compression/PatternAnalysis/complete1.C [new file with mode: 0644]
ITS/UPGRADE/compression/PatternAnalysis/errClus.C [new file with mode: 0644]
ITS/UPGRADE/compression/PatternAnalysis/groupdef.C [new file with mode: 0644]

diff --git a/ITS/UPGRADE/compression/PatternAnalysis/BuildClTopDB.C b/ITS/UPGRADE/compression/PatternAnalysis/BuildClTopDB.C
new file mode 100644 (file)
index 0000000..c8f9102
--- /dev/null
@@ -0,0 +1,338 @@
+/*
+This macro loops over all the cluster (output of compClusHits.C) and creates a database containing the topologies of the cluster and infromation about:
+pattern in TBits format, coordinate of the centre of the pixel contining the COG wrt down-left corner of the osculating rectangle (fraction of the pixel pithc),
+distance between COG and centre of the pixel (fraction), number of fired pixels, number of colums, number of rows.
+*/
+
+#if !defined(__CInt_t__) || defined(__MAKECInt_t__)
+#include "AliESDEvent.h"
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TBits.h"
+#include "TArrayI.h"
+#include "TArrayF.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"
+
+#endif
+
+TTree* treeInp = 0;
+TTree* cluTree = 0;
+AliRunLoader* runLoader = 0;
+AliESDEvent* esd=0;
+
+AliITSUGeomTGeo* gm=0;
+Int_t nLayers = 0;
+AliITSURecoDet *its = 0;
+TObjArray patterns;
+TArrayF pattCount;
+Int_t nPatterns=0,nClusters=0;
+TString outDBFileName = "clusterTopology.root";
+
+void ProcessEvent(Int_t iev);
+void ProcChunk(const char* path);
+void BuildClTopDB(const char* inpData="AliESDs.root", const char* outDBFile="clusterTopology.root");
+void AccountTopology(TBits& patt);
+void CreateDB();
+
+void BuildClTopDB(const char* inpData, const char* outDBFile)
+{
+  //
+  // parse input list of single root file and process it
+  outDBFileName = outDBFile;
+  if (outDBFileName.IsNull()) outDBFileName = "clusterTopology.root";
+  //
+  AliGeomManager::LoadGeometry("geometry.root");
+  gm = new AliITSUGeomTGeo(kTRUE);
+  AliITSUClusterPix::SetGeom(gm);
+  nLayers = gm->GetNLayers();
+  its = new AliITSURecoDet(gm, "ITSInt_terface");
+  its->CreateClusterArrays();
+  esd = new AliESDEvent();
+  //
+  //
+  TString inpDtStr = inpData;
+  if (inpDtStr.EndsWith(".root") || inpDtStr.EndsWith(".zip#")) {
+    ProcChunk(inpDtStr.Data());
+  }
+  else {
+    ifstream inpf(inpData);
+    if (!inpf.good()) {
+      printf("Failed on input filename %s\n",inpData);
+      return;
+    }
+    //
+    inpDtStr.ReadLine(inpf);
+    while ( !inpDtStr.IsNull() ) {
+      inpDtStr = inpDtStr.Strip(TString::kBoth,' ');
+      if (inpDtStr.BeginsWith("//") || inpDtStr.BeginsWith("#")) {inpDtStr.ReadLine(inpf); continue;}
+      inpDtStr = inpDtStr.Strip(TString::kBoth,',');
+      inpDtStr = inpDtStr.Strip(TString::kBoth,'"');
+      ProcChunk(inpDtStr.Data());
+      inpDtStr.ReadLine(inpf);
+    }
+  }
+  //
+  CreateDB();
+  //
+}
+
+void CreateDB()
+{
+  // sort patterns according to frequencies and store them
+  printf("Storing %d patterns in %s\n",nPatterns, outDBFileName.Data());
+  if (nPatterns<1) {
+    printf("No patterns found\n");
+    return;
+  }
+  TArrayI sortIndex;
+  sortIndex.Set(nPatterns);
+  TMath::Sort(nPatterns,pattCount.GetArray(),sortIndex.GetArray());
+  //
+  //create a file.txt containing topologies and corresponding frequecies
+  ofstream top("topologies.txt");
+
+  TObjArray arrStore(nPatterns);
+  TVectorF arrFreq(nPatterns);
+  TVectorF arrzCOGPix(nPatterns); //It's the cordinate of the pixel containing along the rows direction
+  TVectorF arrxCOGPix(nPatterns); //It's the COG cordinate along the columns direction
+  TVectorF arrzCOGshift(nPatterns); //It's the distance between COG and the centre of the pixel containing COG
+  TVectorF arrxCOGshift(nPatterns);
+  TVectorF arrNpix(nPatterns);
+  TVectorF arrNcol(nPatterns);
+  TVectorF arrNrow(nPatterns);
+
+  TCanvas* cnv = new TCanvas("cnv","Freq",900,600);
+  TH1F* hfreq = new TH1F("hfreq","Pattern frequency",300/2,-0.5, 300/2-0.5);
+  hfreq->SetStats(0);
+  hfreq->GetXaxis()->SetTitle("pattern ID");
+  hfreq->SetFillColor(kRed);
+  hfreq->SetFillStyle(3005);
+  TLegend* leg = new TLegend(70, 300000, 140,450000 ,"","");
+  leg->AddEntry((TObject*)0, Form("Number of patterns: %d",nPatterns), "");
+  leg->AddEntry((TObject*)0, "Number of events: 50", "");
+  leg->SetFillColor(0);
+  leg->SetBorderSize(0);
+  
+  for (Int_t i=0;i<nPatterns;i++) {
+    Int_t ind = sortIndex[i];
+    Float_t fr = arrFreq[i] = Float_t(pattCount[ind])/nClusters;
+    arrStore.AddAt(patterns[ind],i);
+    UInt_t idd = patterns[ind]->GetUniqueID();
+
+    hfreq->Fill(i,pattCount[ind]);
+
+    Int_t firedPixels=0;
+    Int_t tempxCOG=0;
+    Int_t tempzCOG=0;
+
+    Int_t rc = idd>>16;//It's the rows span
+    Int_t cl = idd&0xffff;//It's the columns span
+    top << Form("\nPatt %4d R:%d C: %d Freq:%f \n\n",i,rc,cl,fr);
+    for (Int_t ir=rc;ir--;) {
+      top << "|"; 
+      for (Int_t ic=0; ic<cl; ic++) {
+        top << Form("%c",((TBits*)patterns[ind])->TestBitNumber(ir*cl+ic) ? '+':' ');
+        if(((TBits*)patterns[ind])->TestBitNumber(ir*cl+ic)){
+          firedPixels++;
+          tempxCOG+=ir;
+          tempzCOG+=ic;
+        }
+      }
+      top << ("|\n");
+    }
+    Float_t xsh=Float_t((tempxCOG%firedPixels))/firedPixels; //distance between COG end centre of the pixel containing COG
+    Float_t zsh=Float_t((tempzCOG%firedPixels))/firedPixels;
+    tempxCOG/=firedPixels;
+    tempzCOG/=firedPixels;
+    if(xsh>0.5){
+      tempxCOG+=1;
+      xsh-=1;
+    }
+    if(zsh>0.5){
+      tempzCOG+=1;
+      zsh-=1;
+    }
+    Float_t xc=(Float_t) tempxCOG+0.5;
+    Float_t zc=(Float_t) tempzCOG+0.5;
+    arrxCOGPix[i]=xc;
+    arrxCOGshift[i]=xsh;
+    arrzCOGPix[i]=zc;
+    arrzCOGshift[i]=zsh;
+    arrNpix[i]=firedPixels;
+    arrNrow[i]=rc;
+    arrNcol[i]=cl;
+    
+    top << "\nxCentre: " <<  xc << " + " << xsh
+    <<" zCentre: " << zc << " + " << zsh <<"\n\n...............................................\n";
+  }
+
+  cnv->cd();
+  hfreq->Draw();
+  leg->Draw();
+  cnv->Print("Pattern_freq.jpg");
+
+  top.close();
+
+  printf("\nWe have %d clusters corresponding to %d patterns!!!\nEnjoy!!\n\n", nClusters, nPatterns);
+
+  /*
+  //Controlling if frequencies have been stored correctly
+  ofstream a("frequency.txt");
+  a << setw(15) << "i" << setw(20) << "Frequency" << endl;
+  for(Int_t i=0; i<nPatterns; i++) a << setw(15) << i << setw(20) << arrFreq[i] << endl;
+  a.close();
+  */
+  //
+  TFile* flDB = TFile::Open(outDBFileName.Data(), "recreate");
+  flDB->WriteObject(&arrStore,"TopDB","kSingleKey");
+  flDB->WriteObject(&arrFreq, "TopFreq","kSingleKey");
+  flDB->WriteObject(&arrxCOGPix,"xCOG","kSingleKey");
+  flDB->WriteObject(&arrxCOGshift,"xShift","kSingleKey");
+  flDB->WriteObject(&arrzCOGPix,"zCOG","kSingleKey");
+  flDB->WriteObject(&arrzCOGshift,"zShift","kSingleKey");
+  flDB->WriteObject(&arrNpix,"NPix","kSingleKey");
+  flDB->WriteObject(&arrNrow,"NRow","kSingleKey");
+  flDB->WriteObject(&arrNcol,"NCol","kSingleKey");
+  //
+  flDB->Close();
+  delete flDB;
+}
+
+void ProcChunk(const char* path)
+{
+  //
+  TString dir=path;
+  //
+  printf("Processing %s\n",dir.Data());
+  //
+  if (dir.EndsWith("AliESDs.root")) {
+    dir.ReplaceAll("AliESDs.root","");
+  }
+  //
+  //
+  runLoader = AliRunLoader::Open(Form("%sgalice.root",dir.Data()));
+  if (!runLoader) {
+    printf("galice not found\n");
+    return;
+  }
+  runLoader->LoadgAlice();
+  gAlice = runLoader->GetAliRun();
+  //  runLoader->LoadHeader();
+  //  runLoader->LoadKinematics();
+  runLoader->LoadRecPoints("ITS");
+  // ESD
+  TFile* esdFile = TFile::Open(Form("%sAliESDs.root",dir.Data()));
+  if (!esdFile || !esdFile->IsOpen()) {
+    printf("Error in opening ESD file\n");
+    //    runLoader->UnloadKinematics();
+    //    runLoader->UnloadHeader();
+    runLoader->UnloadgAlice();
+    delete runLoader;
+    return;
+  }
+
+  treeInp = (TTree*) esdFile->Get("esdTree");
+  if (!treeInp) {
+    printf("Error: no ESD tree found\n");
+    //    runLoader->UnloadKinematics();
+    //    runLoader->UnloadHeader();
+    runLoader->UnloadgAlice();
+    delete runLoader;
+    return;
+  }
+  esd->ReadFromTree(treeInp);
+  //
+  for (Int_t iEv=0;iEv<runLoader->GetNumberOfEvents(); iEv++) {
+    ProcessEvent(iEv);
+  }
+  runLoader->UnloadRecPoints("all");
+  //  runLoader->UnloadKinematics();
+  //  runLoader->UnloadHeader();
+  runLoader->UnloadgAlice();
+  delete runLoader; 
+  runLoader = 0;
+  esdFile->Close();
+  delete esdFile;
+}
+
+
+void ProcessEvent(Int_t iEv) 
+{
+  // process single event clusters
+  runLoader->GetEvent(iEv);
+  treeInp->GetEvent(iEv);
+  cluTree = runLoader->GetTreeR("ITS",kFALSE);
+  if (!cluTree) {
+    printf("Did not get cluster tree for event %d\n",iEv);
+    return;
+  }
+  //
+  for (Int_t ilr=0;ilr<nLayers;ilr++) {
+    TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",ilr));
+    if (!br) {printf("Did not find cluster branch for lr %d\n",ilr); exit(1);}
+    br->SetAddress(its->GetLayerActive(ilr)->GetClustersAddress());
+  }
+  cluTree->GetEntry(0); 
+  its->ProcessClusters();
+  //
+  TBits currTop;
+  //
+  for (Int_t ilr=0;ilr<nLayers;ilr++) {
+    AliITSURecoLayer* lr = its->GetLayerActive(ilr);
+    TClonesArray* clr = lr->GetClusters();
+    Int_t nClu = clr->GetEntries();
+    printf("Ev%d Lr:%d Ncl:%d\n",iEv,ilr,nClu);
+    //    printf("Layer %d : %d clusters\n",ilr,nClu);
+    //
+    for (Int_t icl=0;icl<nClu;icl++) {
+      AliITSUClusterPix *cl = (AliITSUClusterPix*)clr->At(icl);
+      currTop.Clear();
+      UShort_t rs = (UShort_t)cl->GetPatternRowSpan();
+      UShort_t cs = (UShort_t)cl->GetPatternColSpan();
+      for (Int_t ir=0;ir<rs;ir++)
+            for (Int_t ic=0;ic<cs;ic++) 
+              if (cl->TestPixel(ir,ic)) currTop.SetBitNumber(ir*cs+ic);
+      //
+      currTop.SetUniqueID( (rs<<16) + (cs));
+      //
+      AccountTopology(currTop);
+    }
+  }
+  //
+}
+
+//___________________________________________
+void AccountTopology(TBits& patt)
+{
+  // account topology in the database
+  Bool_t newPatt = kTRUE;
+  nClusters++;
+  for (Int_t ip=0;ip<nPatterns;ip++) {
+    TBits* pattOld = (TBits*)patterns.At(ip);
+    if (*pattOld == patt && pattOld->GetUniqueID()==patt.GetUniqueID()) {
+      // require that the matrix size and the bit patterns are equal
+      newPatt = kFALSE;
+      pattCount[ip]++;
+      break;
+    }
+  }
+  //
+  if (newPatt) {
+    TBits* pt = new TBits(patt);
+    patterns.AddLast(pt);
+    if (pattCount.GetSize()<nPatterns+100) pattCount.Set(nPatterns+100);
+    pattCount[nPatterns++] = 1;
+  }
+  //
+}
\ No newline at end of file
diff --git a/ITS/UPGRADE/compression/PatternAnalysis/README b/ITS/UPGRADE/compression/PatternAnalysis/README
new file mode 100644 (file)
index 0000000..5396b4e
--- /dev/null
@@ -0,0 +1,9 @@
+Macros from Luca Barioglio for claster patterns analysis
+Presented in 
+https://indico.cern.ch/event/339513/contribution/3/material/slides/0.pdf
+- BuildClTopDB.C
+- compClusHits.C (modified with impact angles)
+- SecondcompClusHits.C (with DB Loading)
+- complete1.C (all the histograms)
+- errClus.C (just 1D histograms)
+- groupdef.C (grouping)
diff --git a/ITS/UPGRADE/compression/PatternAnalysis/SecondcompClusHits.C b/ITS/UPGRADE/compression/PatternAnalysis/SecondcompClusHits.C
new file mode 100644 (file)
index 0000000..4b8db84
--- /dev/null
@@ -0,0 +1,652 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TBits.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 "AliStack.h"
+#include "AliLoader.h"
+#include "AliCDBManager.h"
+
+#include "TROOT.h"
+#include "TStyle.h"
+#include "TGeoMatrix.h"
+#include "TParticle.h"
+#include "TCanvas.h"
+#include "TPaveStats.h"
+#include "TClonesArray.h"
+
+#endif
+
+TObjArray *pattDB=0;
+TVectorF*  pattFR=0;
+TVectorF*  xCentrePix=0;
+TVectorF*  zCentrePix=0;
+TVectorF*  xCentreShift=0;
+TVectorF*  zCentreShift=0;
+TObjArray histoArr;
+enum {kNPixAll=0,kNPixSPL=1,kDR=0,kDTXodd,kDTXeven,kDTZ, kDTXoddSPL,kDTXevenSPL,kDTZSPL};
+
+TPaveStats* GetStPad(TH1* hst);
+TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl=-1,Int_t col=-1);
+TCanvas* DrawNP(int np, TObjArray* harr=0, TCanvas* cnv=0);
+TH1* GetHistoClSize(int npix,int id,TObjArray* harr=0);
+void DrawReport(const char* psname, TObjArray* harr=0);
+void LoadDB(const char* fname);
+
+
+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;
+
+TBits clTop;
+
+void SecondcompClusHits(int nev=-1)
+{
+  const int kSplit=0x1<<22;
+  const int kSplCheck=0x1<<23;
+  //
+  gSystem->Load("libITSUpgradeBase");
+  gSystem->Load("libITSUpgradeSim");
+  gSystem->Load("libITSUpgradeRec");
+  gROOT->SetStyle("Plain");
+
+  AliCDBManager* man = AliCDBManager::Instance();
+  man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+  man->SetSpecificStorage("GRP/GRP/Data",
+                         Form("local://%s",gSystem->pwd()));
+  man->SetSpecificStorage("ITS/Align/Data",
+                         Form("local://%s",gSystem->pwd()));
+  man->SetSpecificStorage("ITS/Calib/RecoParam",
+                         Form("local://%s",gSystem->pwd()));
+  man->SetRun(0);
+
+
+  gAlice=NULL;
+  AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
+  runLoader->LoadgAlice();
+
+  gAlice = runLoader->GetAliRun();
+
+  runLoader->LoadHeader();
+  runLoader->LoadKinematics();
+  runLoader->LoadRecPoints();
+  runLoader->LoadSDigits();
+  runLoader->LoadHits();
+
+  AliLoader *dl = runLoader->GetDetectorLoader("ITS");
+
+  AliGeomManager::LoadGeometry("geometry.root");
+  TObjArray algITS;
+  AliGeomManager::LoadAlignObjsFromCDBSingleDet("ITS",algITS);
+  AliGeomManager::ApplyAlignObjsToGeom(algITS);
+  //
+  AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
+  AliITSUClusterPix::SetGeom(gm);
+  //
+  AliITSURecoDet *its = new AliITSURecoDet(gm, "ITSinterface");
+  its->CreateClusterArrays();
+  //
+  Double_t xg1,yg1,zg1=0.,xg0,yg0,zg0=0.,tg0;
+  Double_t xExit,yExit,zExit,xEnt,yEnt,zEnt,tof1;
+  //
+  TTree * cluTree = 0x0;
+  TTree *hitTree = 0x0;
+  TClonesArray *hitList=new TClonesArray("AliITSUHit");
+  //
+  TObjArray arrMCTracks; // array of hit arrays for each particle
+  //
+  Float_t xyzClGloF[3];
+  Double_t xyzClGlo[3],xyzClTr[3];
+  Int_t labels[3];
+  int nLab = 0;
+  int nlr=its->GetNLayersActive();
+  int ntotev = (Int_t)runLoader->GetNumberOfEvents();
+  printf("N Events : %i \n",ntotev);
+  if (nev>0) ntotev = TMath::Min(nev,ntotev);
+  //
+  // load topology database
+  LoadDB("clusterTopology.root");
+
+  Int_t nPatterns = pattDB->GetEntriesFast();
+
+  // output tree
+  TFile* flOut = TFile::Open("clInfo.root","recreate");
+  TTree* trOut = new TTree("clitsu","clitsu");
+  clSumm cSum;
+  trOut->Branch("evID", &cSum.evID ,"evID/I");
+  trOut->Branch("volID",&cSum.volID,"volID/I");
+  trOut->Branch("lrID", &cSum.lrID ,"lrID/I");  
+  trOut->Branch("clID", &cSum.clID ,"clID/I");  
+  trOut->Branch("nPix", &cSum.nPix ,"nPix/I");
+  trOut->Branch("nX"  , &cSum.nX   ,"nX/I");
+  trOut->Branch("nZ"  , &cSum.nZ   ,"nZ/I");
+  trOut->Branch("q"   , &cSum.q    ,"q/I");
+  trOut->Branch("pt"  , &cSum.pt   ,"pt/F");  
+  trOut->Branch("eta"  ,&cSum.eta  ,"eta/F");  
+  trOut->Branch("phi"  , &cSum.phi  ,"phi/F");  
+  trOut->Branch("xyz",   cSum.xyz,  "xyz[3]/F");  
+  trOut->Branch("dX"  , &cSum.dX   ,"dX/F");
+  trOut->Branch("dY"  , &cSum.dY   ,"dY/F");
+  trOut->Branch("dZ"  , &cSum.dZ   ,"dZ/F");  
+  trOut->Branch("split",&cSum.split,"split/O");
+  trOut->Branch("prim", &cSum.prim, "prim/O");
+  trOut->Branch("pdg",  &cSum.pdg,  "pdg/I");
+  trOut->Branch("ntr",  &cSum.ntr,  "ntr/I");
+  trOut->Branch("alpha", &cSum.alpha, "alpha/F");
+  trOut->Branch("beta", &cSum.beta, "beta/F");
+  trOut->Branch("nRowPatt", &cSum.nRowPatt, "nRowPatt/I");
+  trOut->Branch("nColPatt", &cSum.nColPatt, "nColPatt/I");
+  trOut->Branch("pattID", &cSum.pattID,"pattID/I");
+  trOut->Branch("freq", &cSum.freq, "freq/F");
+  trOut->Branch("xCen", &cSum.xCen, "xCen/F");
+  trOut->Branch("zCen", &cSum.zCen, "zCen/F");
+  trOut->Branch("xShift", &cSum.xShift, "xShift/F");
+  trOut->Branch("zShift", &cSum.zShift, "zShift/F");
+  //
+  for (Int_t iEvent = 0; iEvent < ntotev; iEvent++) {
+    printf("\n Event %i \n",iEvent);
+    runLoader->GetEvent(iEvent);
+    AliStack *stack = runLoader->Stack();
+    cluTree=dl->TreeR();
+    hitTree=dl->TreeH();
+    hitTree->SetBranchAddress("ITS",&hitList);
+    // 
+    // read clusters
+    for (int ilr=nlr;ilr--;) {
+      TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",ilr));
+      if (!br) {printf("Did not find cluster branch for lr %d\n",ilr); exit(1);}
+      br->SetAddress(its->GetLayerActive(ilr)->GetClustersAddress());
+    }
+    cluTree->GetEntry(0); 
+    its->ProcessClusters();
+    //
+    // read hits
+    for(Int_t iEnt=0;iEnt<hitTree->GetEntries();iEnt++){//entries loop degli hits
+      hitTree->GetEntry(iEnt);
+      int nh = hitList->GetEntries();      
+      for(Int_t iHit=0; iHit<nh;iHit++){
+            AliITSUHit *pHit = (AliITSUHit*)hitList->At(iHit);
+            int mcID = pHit->GetTrack();
+            TClonesArray* harr = arrMCTracks.GetEntriesFast()>mcID ? (TClonesArray*)arrMCTracks.At(mcID) : 0;
+            if (!harr) {
+              harr = new TClonesArray("AliITSUHit"); // 1st encounter of the MC track
+              arrMCTracks.AddAtAndExpand(harr,mcID);
+            }
+              //
+            new ( (*harr)[harr->GetEntriesFast()] ) AliITSUHit(*pHit);
+      }
+    }
+    //
+    // compare clusters and hits
+    //
+    printf(" tree entries: %lld\n",cluTree->GetEntries());
+    //
+    for (int ilr=0;ilr<nlr;ilr++) {
+      AliITSURecoLayer* lr = its->GetLayerActive(ilr);
+      TClonesArray* clr = lr->GetClusters();
+      int nClu = clr->GetEntries();
+      printf("Layer %d : %d clusters\n",ilr,nClu);
+      //
+      for (int icl=0;icl<nClu;icl++) {
+             AliITSUClusterPix *cl = (AliITSUClusterPix*)clr->At(icl);
+             int modID = cl->GetVolumeId();
+
+        //Compare clusters with DB
+
+        clTop.Clear();
+        Int_t tempPattID=-1;
+        Float_t tempFreq=0;
+        Float_t tempxCen=0;
+        Float_t tempzCen=0;
+        Float_t tempxShift=0;
+        Float_t tempzShift=0;
+        UShort_t rs = cl->GetPatternRowSpan();
+        UShort_t cs = cl->GetPatternColSpan();
+
+        for(int ir=0; ir<rs; ir++)
+          for(int ic=0; ic<cs; ic++)
+            if(cl->TestPixel(ir,ic)) clTop.SetBitNumber(ir*cs+ic);//Set cluster pattern
+        clTop.SetUniqueID((rs<<16)+cs);//Set cluster ID
+
+        for(int iDB=0; iDB<nPatterns; iDB++){
+          TBits* refPatt = (TBits*) pattDB->At(iDB);
+          if (*refPatt == clTop && refPatt->GetUniqueID()== clTop.GetUniqueID()){
+            tempPattID=iDB;
+            tempFreq=(*pattFR)[iDB];
+            tempxCen=(*xCentrePix)[iDB];
+            tempzCen=(*zCentrePix)[iDB];
+            tempxShift=(*xCentreShift)[iDB];
+            tempzShift=(*zCentreShift)[iDB];
+            break;
+          }
+        }   
+       
+            //------------ check if this is a split cluster
+            if (!cl->TestBit(kSplCheck)) {
+              cl->SetBit(kSplCheck);
+              // check if there is no other cluster with same label on this module
+              AliITSURecoSens* sens = lr->GetSensorFromID(modID);
+              int nclSn = sens->GetNClusters();
+              int offs = sens->GetFirstClusterId();
+              //       printf("To check for %d (mod:%d) N=%d from %d\n",icl,modID,nclSn,offs);
+              for (int ics=0;ics<nclSn;ics++) {
+                AliITSUClusterPix* clusT = (AliITSUClusterPix*)lr->GetCluster(offs+ics); // access to clusters
+                if (clusT==cl) continue;
+                for (int ilb0=0;ilb0<3;ilb0++) {
+                  int lb0 = cl->GetLabel(ilb0); if (lb0<=-1) break;
+                  for (int ilb1=0;ilb1<3;ilb1++) {
+                           int lb1 = clusT->GetLabel(ilb1); if (lb1<=-1) break;
+                           if (lb1==lb0) {
+                             cl->SetBit(kSplit);
+                             clusT->SetBit(kSplit);
+                             /*
+                             printf("Discard clusters of module %d:\n",modID);
+                             cl->Print();
+                             clusT->Print();
+                             */
+                             break;
+                           }
+                   }
+                 }
+               }
+             }
+            //------------
+            const AliITSsegmentation* segm = gm->GetSegmentation(ilr);
+            //
+            cl->GetGlobalXYZ(xyzClGloF);
+            int clsize = cl->GetNPix();
+            for (int i=3;i--;) xyzClGlo[i] = xyzClGloF[i];
+            const TGeoHMatrix* mat = gm->GetMatrixSens(modID);
+            if (!mat) {printf("failed to get matrix for module %d\n",cl->GetVolumeId());}
+            mat->MasterToLocal(xyzClGlo,xyzClTr);
+            //
+            int col,row;
+            segm->LocalToDet(xyzClTr[0],xyzClTr[2],row,col); // effective col/row
+            nLab = 0;
+            for (int il=0;il<3;il++) {
+              if (cl->GetLabel(il)>=0) labels[nLab++] = cl->GetLabel(il);
+              else break;
+             }
+             // find hit info
+        for (int il=0;il<nLab;il++) {
+          TClonesArray* htArr = (TClonesArray*)arrMCTracks.At(labels[il]);
+          if (!htArr) {printf("did not find MChits for label %d ",labels[il]); cl->Print(); continue;}
+          //
+          int nh = htArr->GetEntriesFast();
+          AliITSUHit *pHit=0;
+          double dst2Max = 1e33;
+          for (int ih=nh;ih--;) {
+            AliITSUHit* tHit = (AliITSUHit*)htArr->At(ih);
+            if (tHit->GetChip()!=modID) continue;
+            tHit->GetPositionG(xg1,yg1,zg1);
+            tHit->GetPositionG0(xg0,yg0,zg0,tg0);
+            double gxyzHDif[3] = { (xg1+xg0)/2 - xyzClGlo[0], (yg1+yg0)/2 - xyzClGlo[1], (zg1+zg0)/2 - xyzClGlo[2] };
+            double dst2 = gxyzHDif[0]*gxyzHDif[0] + gxyzHDif[1]*gxyzHDif[1] + gxyzHDif[2]*gxyzHDif[2];
+            if (dst2<dst2Max) {
+              pHit = tHit;
+              dst2Max = dst2;
+             }
+           }
+               if (!pHit) {
+                printf("did not find MChit for label %d on module %d ",il,modID); 
+                cl->Print(); 
+                htArr->Print();
+                continue;
+               }
+               //
+               pHit->GetPositionG(xg1,yg1,zg1);
+               pHit->GetPositionG0(xg0,yg0,zg0,tg0);
+               //
+               double txyzH[3],gxyzH[3] = { (xg1+xg0)/2, (yg1+yg0)/2, (zg1+zg0)/2 };
+               mat->MasterToLocal(gxyzH,txyzH);
+               double rcl = TMath::Sqrt(xyzClTr[0]*xyzClTr[0]+xyzClTr[1]*xyzClTr[1]);
+               double rht = TMath::Sqrt(txyzH[0]*txyzH[0]+txyzH[1]*txyzH[1]);
+               //
+          //
+          //Angles determination
+
+          pHit->GetPositionL(xExit,yExit,zExit);
+          pHit->GetPositionL0(xEnt,yEnt,zEnt,tof1);
+
+          Double_t dirHit[3]={(xExit-xEnt),(yExit-yEnt),(zExit-zEnt)};
+
+          Double_t alpha1 = TMath::ACos(TMath::Abs(dirHit[1])/TMath::Sqrt(dirHit[0]*dirHit[0]+dirHit[1]*dirHit[1]+dirHit[2]*dirHit[2])); //Polar Angle
+          Float_t alpha2 = (Float_t) alpha1; //convert to float
+          cSum.alpha = alpha2;
+
+          Double_t beta1;
+          beta1 = TMath::ATan2(TMath::Abs(dirHit[0]),TMath::Abs(dirHit[2])); //Azimuthal angle, values from -Pi to Pi
+          Float_t beta2 = (Float_t) beta1;
+          cSum.beta = beta2;
+
+               GetHistoClSize(clsize,kDR,&histoArr)->Fill((rht-rcl)*1e4);
+               if (cl->TestBit(kSplit)) {
+                 if (col%2) GetHistoClSize(clsize,kDTXoddSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
+                 else       GetHistoClSize(clsize,kDTXevenSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
+                 GetHistoClSize(clsize,kDTZSPL,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
+                 GetHistoClSize(0,kNPixSPL,&histoArr)->Fill(clsize);
+               }
+               if (col%2) GetHistoClSize(clsize,kDTXodd,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
+               else       GetHistoClSize(clsize,kDTXeven,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
+               GetHistoClSize(clsize,kDTZ,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
+               GetHistoClSize(0,kNPixAll,&histoArr)->Fill(clsize);
+               //
+               cSum.evID = iEvent;
+               cSum.volID = cl->GetVolumeId();
+         cSum.lrID = ilr;
+         cSum.clID = icl;
+         cSum.nPix = cl->GetNPix();
+         cSum.nX   = cl->GetNx();
+         cSum.nZ   = cl->GetNz();
+         cSum.q    = cl->GetQ();
+         cSum.split = cl->TestBit(kSplit);
+         cSum.dX = (txyzH[0]-xyzClTr[0])*1e4;
+         cSum.dY = (txyzH[1]-xyzClTr[1])*1e4;
+         cSum.dZ = (txyzH[2]-xyzClTr[2])*1e4;
+          cSum.nRowPatt = cl-> GetPatternRowSpan();
+          cSum.nColPatt = cl-> GetPatternColSpan();
+          cSum.pattID = tempPattID;
+          cSum.freq = tempFreq;
+          cSum.xCen = tempxCen;
+          cSum.zCen = tempzCen;
+          cSum.xShift = tempxShift;
+          cSum.zShift = tempzShift;
+
+         int label = cl->GetLabel(0);
+         TParticle* part = 0;
+         if (label>=0 && (part=stack->Particle(label)) ) {
+           cSum.pdg = part->GetPdgCode();
+           cSum.eta = part->Eta();
+           cSum.pt  = part->Pt();
+           cSum.phi = part->Phi();
+                 cSum.prim = stack->IsPhysicalPrimary(label);
+               }
+               cSum.ntr = 0;
+               for (int ilb=0;ilb<3;ilb++) if (cl->GetLabel(ilb)>=0) cSum.ntr++;
+               for (int i=0;i<3;i++) cSum.xyz[i] = xyzClGloF[i];
+               //
+               trOut->Fill();
+         /*
+         if (clsize==5) {
+           printf("\nL%d(%c) Mod%d, Cl:%d | %+5.1f %+5.1f (%d/%d)|H:%e %e %e | C:%e %e %e\n",ilr,cl->TestBit(kSplit) ? 'S':'N',
+                  modID,icl,(txyzH[0]-xyzClTr[0])*1e4,(txyzH[2]-xyzClTr[2])*1e4, row,col,
+                  gxyzH[0],gxyzH[1],gxyzH[2],xyzClGlo[0],xyzClGlo[1],xyzClGlo[2]);
+           cl->Print();
+           pHit->Print();
+           //
+           double a0,b0,c0,a1,b1,c1,e0;
+           pHit->GetPositionL0(a0,b0,c0,e0);
+           pHit->GetPositionL(a1,b1,c1);
+           float cloc[3];
+           cl->GetLocalXYZ(cloc);
+           printf("LocH: %e %e %e | %e %e %e\n",a0,b0,c0,a1,b1,c1);
+           printf("LocC: %e %e %e | %e %e %e\n",cloc[0],cloc[1],cloc[2],xyzClTr[0],xyzClTr[1],xyzClTr[2]);
+         }
+         */
+         //
+       }
+      }
+    }
+    
+    //    layerClus.Clear();
+    //
+    arrMCTracks.Delete();
+  }//event loop
+  //
+  flOut->cd();
+  trOut->Write();
+  delete trOut;
+  flOut->Close();
+  flOut->Delete();
+  DrawReport("clinfo.ps",&histoArr);
+  //
+}
+
+void DrawReport(const char* psname, TObjArray* harr) 
+{
+  gStyle->SetOptFit(1);
+  if (!harr) harr = &histoArr;
+  TCanvas* cnv = new TCanvas("cl","cl",900,600);
+  //
+  TString psnm1 = psname;
+  if (psnm1.IsNull()) psnm1 = "clusters.ps";
+  TString psnm0 = psnm1.Data(); 
+  psnm0 += "[";
+  TString psnm2 = psnm1.Data(); 
+  psnm2 += "]";
+  cnv->Print(psnm0.Data());
+  //
+  TH1* clall = GetHistoClSize(0,kNPixAll,harr);
+  clall->SetLineColor(kRed);
+  clall->Draw();
+  TH1* clSpl = GetHistoClSize(0,kNPixSPL,harr);
+  clSpl->SetLineColor(kBlue);
+  clSpl->Draw("sames");
+  gPad->Modified();
+  gPad->Update();
+  SetStPadPos(clall,0.75,0.97,0.8,1.,-1,clall->GetLineColor());
+  SetStPadPos(clSpl,0.75,0.97,0.6,0.8,-1,clSpl->GetLineColor());
+  gPad->Modified();
+  gPad->Update();
+  gPad->SetLogy(1);
+  //
+  cnv->cd();
+  cnv->Print(psnm1.Data());
+  //
+  // plot cluster sized from 1 to 10
+  for (int i=1;i<=10;i++) {
+    if (clall->GetBinContent(clall->FindBin(i))<100) continue;
+    DrawNP(i,harr,cnv);
+    cnv->Print(psnm1.Data());
+  }
+  cnv->Print(psnm2.Data());
+}
+
+//------------------------------------------
+TH1* GetHistoClSize(int npix,int id,TObjArray* harr)
+{
+  // book histos
+  TH1* h = 0;
+  if (!harr) harr = &histoArr;
+  //
+  if (npix<1) {
+    if (harr->GetEntriesFast()>=id && (h=(TH1*)harr->At(id))) return h;
+    h = new TH1F("npixAll","npixAll",150,0.5,54.5); 
+    h->SetDirectory(0);
+    h->SetLineColor(kRed);
+    harr->AddAtAndExpand(h, kNPixAll);
+    //
+    h = new TH1F("npixSpl","npixSpl",150,0.5,54.5);
+    h->SetLineColor(kBlue);
+    h->SetDirectory(0);
+    harr->AddAtAndExpand(h, kNPixSPL);
+    //
+    h = (TH1*)harr->At(id);
+    if (!h) {printf("Unknown histo id=%d\n",id); exit(1);}
+    return h;
+  }
+  //
+  int idh = npix*10+id;
+  if (harr->GetEntriesFast()>=idh && (h=(TH1*)harr->At(idh))) return h;
+  //
+  const int nbin=100;
+  const double kdiff=80;
+  // need to create set of histos
+  //
+  h = new TH1F(Form("dxy_npix%d",npix),Form("dr_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetDirectory(0);
+  harr->AddAtAndExpand(h, npix*10 + kDR);
+  //
+  h  = new TH1F(Form("dtxODD_npix%d",npix),Form("dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetDirectory(0);
+  h->SetLineColor(kRed);
+  harr->AddAtAndExpand(h, npix*10 + kDTXodd);
+  h  = new TH1F(Form("dtxEVN_npix%d",npix),Form("dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetDirectory(0);
+  h->SetLineColor(kBlue);
+  harr->AddAtAndExpand(h, npix*10 + kDTXeven);
+  //
+  h  = new TH1F(Form("dtz_npix%d",npix),Form("dtz_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetLineColor(kGreen);
+  h->SetDirectory(0);
+  harr->AddAtAndExpand(h, npix*10 + kDTZ);
+  //
+  //
+  h  = new TH1F(Form("SPL_dtxODD_npix%d",npix),Form("SPL_dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetLineColor(kMagenta);
+  h->SetFillColor(kMagenta);
+  h->SetFillStyle(3001);  
+  h->SetLineStyle(2);
+  h->SetDirectory(0);
+
+  harr->AddAtAndExpand(h, npix*10 + kDTXoddSPL);
+  h  = new TH1F(Form("SPL_dtxEVN_npix%d",npix),Form("SPL_dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetLineColor(kCyan);
+  h->SetFillColor(kCyan);
+  h->SetFillStyle(3006);  
+  h->SetLineStyle(2);
+  h->SetDirectory(0);
+  harr->AddAtAndExpand(h, npix*10 + kDTXevenSPL);
+  //
+  h  = new TH1F(Form("SPL_dtz_npix%d",npix),Form("SPLdtz_npix%d",npix),nbin,-kdiff,kdiff);
+  harr->AddAtAndExpand(h, npix*10 + kDTZSPL);
+  h->SetDirectory(0);
+  //
+  h->SetLineColor(kGreen+2);
+  h->SetFillColor(kGreen+2);
+  h->SetLineStyle(2);
+  h->SetFillStyle(3001);
+  h = (TH1*)harr->At(idh);
+  if (!h) {printf("Unknown histo id=%d\n",idh); exit(1);}
+  return h;
+}
+
+
+
+TCanvas* DrawNP(int np, TObjArray* harr, TCanvas* cnv)
+{
+  if (!harr) harr = &histoArr;
+  if (!cnv) cnv = new TCanvas(Form("cnv%d",np),Form("cnv%d",np),900,700);
+  cnv->Clear();
+  cnv->Divide(2,1);
+  cnv->cd(1);
+  //
+  TH1* dxodd = (TH1*)harr->At(np*10+kDTXodd);
+  TH1* dxevn = (TH1*)harr->At(np*10+kDTXeven);
+  TH1* dxoddS =(TH1*)harr->At(np*10+kDTXoddSPL);
+  TH1* dxevnS =(TH1*)harr->At(np*10+kDTXevenSPL);
+  double max = TMath::Max(dxodd->GetMaximum(),dxevn->GetMaximum());
+  dxodd->SetMaximum(1.1*max);
+  dxodd->GetXaxis()->SetTitle("#DeltaX, #mum");
+  dxodd->SetTitle(Form("#DeltaX for clSize=%d",np));
+  dxodd->Fit("gaus","","");
+  dxevn->Fit("gaus","","sames");
+  //
+  dxoddS->Draw("sames");
+  dxevnS->Draw("sames");
+  //
+  gPad->Modified();
+  gPad->Update();
+  SetStPadPos(dxodd,0.75,0.97,0.8,1., -1,dxodd->GetLineColor());
+  SetStPadPos(dxevn,0.75,0.97,0.6,0.8, -1,dxevn->GetLineColor());
+  SetStPadPos(dxoddS,0.75,0.97,0.4,0.6, -1,dxoddS->GetLineColor());
+  SetStPadPos(dxevnS,0.75,0.97,0.2,0.4, -1,dxevnS->GetLineColor());
+  //
+  cnv->cd(2);
+  TH1* dz  = (TH1*)harr->At(np*10+kDTZ);
+  dz->SetTitle(Form("#DeltaZ for clSize=%d",np));
+  dz->GetXaxis()->SetTitle("#DeltaZ, #mum");
+  dz->Fit("gaus");
+  TH1* dzS = (TH1*)harr->At(np*10+kDTZSPL);
+  dz->Draw("sames");
+  gPad->Modified();
+  gPad->Update();
+  SetStPadPos(dz,0.75,0.97,0.8,1., -1, dz->GetLineColor());
+  SetStPadPos(dzS,0.75,0.97,0.5,0.7, -1, dzS->GetLineColor());
+  gPad->Modified();
+  gPad->Update();
+  //
+  cnv->cd();
+  return cnv;
+}
+
+
+TPaveStats* GetStPad(TH1* hst)
+{
+  TList *lst = hst->GetListOfFunctions();
+  if (!lst) return 0;
+  int nf = lst->GetSize();
+  for (int i=0;i<nf;i++) {
+    TPaveStats *fnc = (TPaveStats*) lst->At(i);
+    if (fnc->InheritsFrom("TPaveStats")) return fnc;
+  }
+  return 0;
+  //
+}
+
+
+TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl, Int_t col)
+{
+  TPaveStats* pad = GetStPad(hst);
+  if (!pad) return 0;
+  pad->SetX1NDC( x1 );
+  pad->SetX2NDC( x2 );
+  pad->SetY1NDC( y1 );
+  pad->SetY2NDC( y2 );
+  if (stl>=0) pad->SetFillStyle(stl);
+  if (col>=0) pad->SetTextColor(col);
+  pad->SetFillColor(0);
+  //
+  gPad->Modified();
+  return pad;
+}
+
+void LoadDB(const char* fname)
+{
+  // load database
+  TFile* fl = TFile::Open(fname);
+  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");
+
+}
+
diff --git a/ITS/UPGRADE/compression/PatternAnalysis/compClusHits.C b/ITS/UPGRADE/compression/PatternAnalysis/compClusHits.C
new file mode 100644 (file)
index 0000000..25daa40
--- /dev/null
@@ -0,0 +1,581 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.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 "AliStack.h"
+#include "AliLoader.h"
+#include "AliCDBManager.h"
+#include "TROOT.h"
+#include "TStyle.h"
+#include "TGeoMatrix.h"
+#include "TParticle.h"
+#include "TCanvas.h"
+#include "TPaveStats.h"
+#include "TClonesArray.h"
+
+#endif
+
+TObjArray histoArr;
+enum {kNPixAll=0,kNPixSPL=1,kDR=0,kDTXodd,kDTXeven,kDTZ, kDTXoddSPL,kDTXevenSPL,kDTZSPL};
+
+TPaveStats* GetStPad(TH1* hst);
+TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl=-1,Int_t col=-1);
+TCanvas* DrawNP(int np, TObjArray* harr=0, TCanvas* cnv=0);
+TH1* GetHistoClSize(int npix,int id,TObjArray* harr=0);
+void DrawReport(const char* psname, TObjArray* harr=0);
+
+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;
+} clSumm;
+
+void compClusHits(int nev=-1)
+{
+
+  const int kSplit=0x1<<22;
+  const int kSplCheck=0x1<<23;
+  //
+  gSystem->Load("libITSUpgradeBase");
+  gSystem->Load("libITSUpgradeSim");
+  gSystem->Load("libITSUpgradeRec");
+  gROOT->SetStyle("Plain");
+
+  AliCDBManager* man = AliCDBManager::Instance();
+  man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+  man->SetSpecificStorage("GRP/GRP/Data",
+        Form("local://%s",gSystem->pwd()));
+  man->SetSpecificStorage("ITS/Align/Data",
+        Form("local://%s",gSystem->pwd()));
+  man->SetSpecificStorage("ITS/Calib/RecoParam",
+        Form("local://%s",gSystem->pwd()));
+  man->SetRun(0);
+
+  gAlice=NULL;
+  AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
+  runLoader->LoadgAlice();
+
+  gAlice = runLoader->GetAliRun();
+
+  runLoader->LoadHeader();
+  runLoader->LoadKinematics();
+  runLoader->LoadRecPoints();
+  runLoader->LoadSDigits();
+  runLoader->LoadHits();
+
+  AliLoader *dl = runLoader->GetDetectorLoader("ITS");
+
+  AliGeomManager::LoadGeometry("geometry.root");
+  TObjArray algITS;
+  AliGeomManager::LoadAlignObjsFromCDBSingleDet("ITS",algITS);
+  AliGeomManager::ApplyAlignObjsToGeom(algITS);
+  //
+  AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
+  AliITSUClusterPix::SetGeom(gm);
+  //
+  AliITSURecoDet *its = new AliITSURecoDet(gm, "ITSinterface");
+  its->CreateClusterArrays();
+  //
+  Double_t xg1,yg1,zg1=0.,xg0,yg0,zg0=0.,tg0;
+  Double_t xExit,yExit,zExit,xEnt,yEnt,zEnt,tof1;
+
+  //
+  TTree * cluTree = 0x0;
+  TTree *hitTree = 0x0;
+  TClonesArray *hitList=new TClonesArray("AliITSUHit");
+  //
+  TObjArray arrMCTracks; // array of hit arrays for each particle
+  //
+  Float_t xyzClGloF[3];
+  Double_t xyzClGlo[3],xyzClTr[3];
+  Int_t labels[3];
+  int nLab = 0;
+  int nlr=its->GetNLayersActive();
+  int ntotev = (Int_t)runLoader->GetNumberOfEvents();
+  printf("N Events : %i \n",ntotev);
+  if (nev>0) ntotev = TMath::Min(nev,ntotev);
+  //
+  
+  // output tree
+  TFile* flOut = TFile::Open("clInfo.root","recreate");
+  TTree* trOut = new TTree("clitsu","clitsu");
+  clSumm cSum;
+  trOut->Branch("evID", &cSum.evID ,"evID/I");
+  trOut->Branch("volID",&cSum.volID,"volID/I");
+  trOut->Branch("lrID", &cSum.lrID ,"lrID/I");  
+  trOut->Branch("clID", &cSum.clID ,"clID/I");  
+  trOut->Branch("nPix", &cSum.nPix ,"nPix/I");
+  trOut->Branch("nX"  , &cSum.nX   ,"nX/I");
+  trOut->Branch("nZ"  , &cSum.nZ   ,"nZ/I");
+  trOut->Branch("q"   , &cSum.q    ,"q/I");
+  trOut->Branch("pt"  , &cSum.pt   ,"pt/F");  
+  trOut->Branch("eta"  ,&cSum.eta  ,"eta/F");  
+  trOut->Branch("phi"  , &cSum.phi  ,"phi/F");  
+  trOut->Branch("xyz",   cSum.xyz,  "xyz[3]/F");  
+  trOut->Branch("dX"  , &cSum.dX   ,"dX/F");
+  trOut->Branch("dY"  , &cSum.dY   ,"dY/F");
+  trOut->Branch("dZ"  , &cSum.dZ   ,"dZ/F");  
+  trOut->Branch("split",&cSum.split,"split/O");
+  trOut->Branch("prim", &cSum.prim, "prim/O");
+  trOut->Branch("pdg",  &cSum.pdg,  "pdg/I");
+  trOut->Branch("ntr",  &cSum.ntr,  "ntr/I");
+  trOut->Branch("alpha", &cSum.alpha, "alpha/F");
+  trOut->Branch("beta", &cSum.beta, "beta/F");
+  trOut->Branch("nRowPatt", &cSum.nRowPatt, "nRowPatt/I");
+  trOut->Branch("nColPatt", &cSum.nColPatt, "nColPatt/I");
+  
+  for (Int_t iEvent = 0; iEvent < ntotev; iEvent++) {
+    printf("\n Event %i \n",iEvent);
+    runLoader->GetEvent(iEvent);
+    AliStack *stack = runLoader->Stack();
+    cluTree=dl->TreeR();
+    hitTree=dl->TreeH();
+    hitTree->SetBranchAddress("ITS",&hitList);
+    // 
+    // read clusters
+    for (int ilr=nlr;ilr--;) {
+      TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",ilr));
+      if (!br) {printf("Did not find cluster branch for lr %d\n",ilr); exit(1);}
+      br->SetAddress(its->GetLayerActive(ilr)->GetClustersAddress());
+    }
+    cluTree->GetEntry(0); 
+    its->ProcessClusters();
+    //
+    // read hits
+    for(Int_t iEnt=0;iEnt<hitTree->GetEntries();iEnt++){//entries loop degli hits
+      hitTree->GetEntry(iEnt);
+      int nh = hitList->GetEntries();      
+      for(Int_t iHit=0; iHit<nh;iHit++){
+        AliITSUHit *pHit = (AliITSUHit*)hitList->At(iHit);
+        int mcID = pHit->GetTrack();
+        TClonesArray* harr = arrMCTracks.GetEntriesFast()>mcID ? (TClonesArray*)arrMCTracks.At(mcID) : 0;
+        if (!harr) {
+          harr = new TClonesArray("AliITSUHit"); // 1st encounter of the MC track
+          arrMCTracks.AddAtAndExpand(harr,mcID);
+        }
+        //
+        new ( (*harr)[harr->GetEntriesFast()] ) AliITSUHit(*pHit);
+      }
+    }
+
+    //
+    // compare clusters and hits
+    //
+    printf(" tree entries: %lld\n",cluTree->GetEntries());
+    //
+    for (int ilr=0;ilr<nlr;ilr++) {
+      AliITSURecoLayer* lr = its->GetLayerActive(ilr);
+      TClonesArray* clr = lr->GetClusters();
+      int nClu = clr->GetEntries();
+      printf("Layer %d : %d clusters\n",ilr,nClu);
+      //
+      for (int icl=0;icl<nClu;icl++) {
+        AliITSUClusterPix *cl = (AliITSUClusterPix*)clr->At(icl);
+        int modID = cl->GetVolumeId();
+
+        //------------ check if this is a split cluster
+        int sInL = modID - gm->GetFirstChipIndex(ilr);
+        if (!cl->TestBit(kSplCheck)) {
+          cl->SetBit(kSplCheck);
+          // check if there is no other cluster with same label on this module
+          AliITSURecoSens* sens = lr->GetSensor(sInL);
+          int nclSn = sens->GetNClusters();
+          int offs = sens->GetFirstClusterId();
+          //  printf("To check for %d (mod:%d) N=%d from %d\n",icl,modID,nclSn,offs);
+          for (int ics=0;ics<nclSn;ics++) {
+            AliITSUClusterPix* clusT = (AliITSUClusterPix*)lr->GetCluster(offs+ics); // access to clusters
+            if (clusT==cl) continue;
+            for (int ilb0=0;ilb0<3;ilb0++) {
+              int lb0 = cl->GetLabel(ilb0); if (lb0<=-1) break;
+              for (int ilb1=0;ilb1<3;ilb1++) {
+                int lb1 = clusT->GetLabel(ilb1); if (lb1<=-1) break;
+                if (lb1==lb0) {
+                  cl->SetBit(kSplit);
+                  clusT->SetBit(kSplit);
+                  /*
+                  printf("Discard clusters of module %d:\n",modID);
+                  cl->Print();
+                  clusT->Print();
+                  */
+                  break;
+                }
+              }
+            }
+          }
+        }
+        //------------
+        const AliITSsegmentation* segm = gm->GetSegmentation(ilr);
+        //
+        cl->GetGlobalXYZ(xyzClGloF);
+        int clsize = cl->GetNPix();
+        for (int i=3;i--;) xyzClGlo[i] = xyzClGloF[i];
+        const TGeoHMatrix* mat = gm->GetMatrixSens(modID);
+        if (!mat) {printf("failed to get matrix for module %d\n",cl->GetVolumeId());}
+        mat->MasterToLocal(xyzClGlo,xyzClTr);
+        //
+        int col,row;
+        segm->LocalToDet(xyzClTr[0],xyzClTr[2],row,col); // effective col/row
+        nLab = 0;
+        for (int il=0;il<3;il++) {
+          if (cl->GetLabel(il)>=0) labels[nLab++] = cl->GetLabel(il);
+          else break;
+        }
+        // find hit info
+        for (int il=0;il<nLab;il++) {
+          TClonesArray* htArr = (TClonesArray*)arrMCTracks.At(labels[il]);
+          if (!htArr) {printf("did not find MChits for label %d ",labels[il]); cl->Print(); continue;}
+          //
+          int nh = htArr->GetEntriesFast();
+          AliITSUHit *pHit=0;
+          for (int ih=nh;ih--;) {
+            AliITSUHit* tHit = (AliITSUHit*)htArr->At(ih);
+            if (tHit->GetChip()!=modID) continue;
+            pHit = tHit;
+            break;
+          }
+          if (!pHit) {
+            printf("did not find MChit for label %d on module %d ",il,modID); 
+            cl->Print(); 
+            htArr->Print();
+            continue;
+          }
+          //
+          pHit->GetPositionG(xg1,yg1,zg1);
+          pHit->GetPositionG0(xg0,yg0,zg0,tg0);
+          //
+          double txyzH[3],gxyzH[3] = { (xg1+xg0)/2, (yg1+yg0)/2, (zg1+zg0)/2 };
+          mat->MasterToLocal(gxyzH,txyzH);
+
+          double rcl = TMath::Sqrt(xyzClTr[0]*xyzClTr[0]+xyzClTr[1]*xyzClTr[1]);
+          double rht = TMath::Sqrt(txyzH[0]*txyzH[0]+txyzH[1]*txyzH[1]);
+          //
+          //Angles determination
+
+          pHit->GetPositionL(xExit,yExit,zExit);
+          pHit->GetPositionL0(xEnt,yEnt,zEnt,tof1);
+
+          Double_t dirHit[3]={(xExit-xEnt),(yExit-yEnt),(zExit-zEnt)};
+
+          /*double PG[3] = {(double)pHit->GetPXG(), (double)pHit->GetPYG(), (double)pHit->GetPZG()}; //Momentum at hit-point in Global Frame
+          double PL[3];
+          if (TMath::Abs(PG[0])<10e-7 && TMath::Abs(PG[1])<10e-7) {
+            pHit->Dump();
+            int lb = pHit->GetTrack();
+            stack->Particle(lb)->Print();
+            continue;
+          }
+          mat->MasterToLocalVect(PG,PL); //Momentum in local Frame
+          //printf(">> %e %e   %e %e   %e %e\n",PG[0],PL[0],PG[1],PL[1],PG[2],PL[2]);*/
+
+          Double_t alpha1 = TMath::ACos(TMath::Abs(dirHit[1])/TMath::Sqrt(dirHit[0]*dirHit[0]+dirHit[1]*dirHit[1]+dirHit[2]*dirHit[2])); //Polar Angle
+          Float_t alpha2 = (Float_t) alpha1; //convert to float
+          cSum.alpha = alpha2;
+
+          Double_t beta1;
+          beta1 = TMath::ATan2(dirHit[0],dirHit[2]); //Azimuthal angle, values from -Pi to Pi
+          Float_t beta2 = (Float_t) beta1;
+          cSum.beta = beta2;
+          
+          GetHistoClSize(clsize,kDR,&histoArr)->Fill((rht-rcl)*1e4);
+          if (cl->TestBit(kSplit)) {
+            if (col%2) GetHistoClSize(clsize,kDTXoddSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
+            else       GetHistoClSize(clsize,kDTXevenSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
+            GetHistoClSize(clsize,kDTZSPL,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
+            GetHistoClSize(0,kNPixSPL,&histoArr)->Fill(clsize);
+          }
+          if (col%2) GetHistoClSize(clsize,kDTXodd,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
+          else       GetHistoClSize(clsize,kDTXeven,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
+          GetHistoClSize(clsize,kDTZ,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
+          GetHistoClSize(0,kNPixAll,&histoArr)->Fill(clsize);
+          //
+          cSum.evID = iEvent;
+          cSum.volID = cl->GetVolumeId();
+          cSum.lrID = ilr;
+          cSum.clID = icl;
+          cSum.nPix = cl->GetNPix();
+          cSum.nX   = cl->GetNx();
+          cSum.nZ   = cl->GetNz();
+          cSum.q    = cl->GetQ();
+          cSum.split = cl->TestBit(kSplit);
+          cSum.dX = (txyzH[0]-xyzClTr[0])*1e4;
+          cSum.dY = (txyzH[1]-xyzClTr[1])*1e4;
+          cSum.dZ = (txyzH[2]-xyzClTr[2])*1e4;
+          cSum.nRowPatt = cl-> GetPatternRowSpan();
+          cSum.nColPatt = cl-> GetPatternColSpan();
+                    
+          int label = cl->GetLabel(0);
+          TParticle* part = 0;
+          if (label>=0 && (part=stack->Particle(label)) ) {
+            cSum.pdg = part->GetPdgCode();
+            cSum.eta = part->Eta();
+            cSum.pt  = part->Pt();
+            cSum.phi = part->Phi();
+            cSum.prim = stack->IsPhysicalPrimary(label);
+          } 
+          cSum.ntr = 0;
+          for (int ilb=0;ilb<3;ilb++) if (cl->GetLabel(ilb)>=0) cSum.ntr++;
+          for (int i=0;i<3;i++) cSum.xyz[i] = xyzClGloF[i];
+          //
+          trOut->Fill();
+          /*
+          if (clsize==5) {
+            printf("\nL%d(%c) Mod%d, Cl:%d | %+5.1f %+5.1f (%d/%d)|H:%e %e %e | C:%e %e %e\n",ilr,cl->TestBit(kSplit) ? 'S':'N',
+             modID,icl,(txyzH[0]-xyzClTr[0])*1e4,(txyzH[2]-xyzClTr[2])*1e4, row,col,
+             gxyzH[0],gxyzH[1],gxyzH[2],xyzClGlo[0],xyzClGlo[1],xyzClGlo[2]);
+            cl->Print();
+            pHit->Print();
+            //
+            double a0,b0,c0,a1,b1,c1,e0;
+            pHit->GetPositionL0(a0,b0,c0,e0);
+            pHit->GetPositionL(a1,b1,c1);
+            float cloc[3];
+            cl->GetLocalXYZ(cloc);
+            printf("LocH: %e %e %e | %e %e %e\n",a0,b0,c0,a1,b1,c1);
+            printf("LocC: %e %e %e | %e %e %e\n",cloc[0],cloc[1],cloc[2],xyzClTr[0],xyzClTr[1],xyzClTr[2]);
+          }
+          */
+          //
+        }
+      }
+    }
+    
+    //    layerClus.Clear();
+    //
+    arrMCTracks.Delete();
+  }//event loop
+  //
+  flOut->cd();
+  trOut->Write();
+  delete trOut;
+  flOut->Close();
+  flOut->Delete();
+  DrawReport("clinfo.ps",&histoArr);
+  //
+}
+
+void DrawReport(const char* psname, TObjArray* harr) 
+{
+  gStyle->SetOptFit(1);
+  if (!harr) harr = &histoArr;
+  TCanvas* cnv = new TCanvas("cl","cl",900,600);
+  //
+  TString psnm1 = psname;
+  if (psnm1.IsNull()) psnm1 = "clusters.ps";
+  TString psnm0 = psnm1.Data(); 
+  psnm0 += "[";
+  TString psnm2 = psnm1.Data(); 
+  psnm2 += "]";
+  cnv->Print(psnm0.Data());
+  //
+  TH1* clall = GetHistoClSize(0,kNPixAll,harr);
+  clall->SetLineColor(kRed);
+  clall->Draw();
+  TH1* clSpl = GetHistoClSize(0,kNPixSPL,harr);
+  clSpl->SetLineColor(kBlue);
+  clSpl->Draw("sames");
+  gPad->Modified();
+  gPad->Update();
+  SetStPadPos(clall,0.75,0.97,0.8,1.,-1,clall->GetLineColor());
+  SetStPadPos(clSpl,0.75,0.97,0.6,0.8,-1,clSpl->GetLineColor());
+  gPad->Modified();
+  gPad->Update();
+  gPad->SetLogy(1);
+  //
+  cnv->cd();
+  cnv->Print(psnm1.Data());
+  //
+  // plot cluster sized from 1 to 10
+  for (int i=1;i<=10;i++) {
+    if (clall->GetBinContent(clall->FindBin(i))<100) continue;
+    DrawNP(i,harr,cnv);
+    cnv->Print(psnm1.Data());
+  }
+  cnv->Print(psnm2.Data());
+}
+
+//------------------------------------------
+TH1* GetHistoClSize(int npix,int id,TObjArray* harr)
+{
+  // book histos
+  TH1* h = 0;
+  if (!harr) harr = &histoArr;
+  //
+  if (npix<1) {
+    if (harr->GetEntriesFast()>=id && (h=(TH1*)harr->At(id))) return h;
+    h = new TH1F("npixAll","npixAll",150,0.5,54.5); 
+    h->SetDirectory(0);
+    h->SetLineColor(kRed);
+    harr->AddAtAndExpand(h, kNPixAll);
+    //
+    h = new TH1F("npixSpl","npixSpl",150,0.5,54.5);
+    h->SetLineColor(kBlue);
+    h->SetDirectory(0);
+    harr->AddAtAndExpand(h, kNPixSPL);
+    //
+    h = (TH1*)harr->At(id);
+    if (!h) {printf("Unknown histo id=%d\n",id); exit(1);}
+    return h;
+  }
+  //
+  int idh = npix*10+id;
+  if (harr->GetEntriesFast()>=idh && (h=(TH1*)harr->At(idh))) return h;
+  //
+  const int nbin=100;
+  const double kdiff=80;
+  // need to create set of histos
+  //
+  h = new TH1F(Form("dxy_npix%d",npix),Form("dr_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetDirectory(0);
+  harr->AddAtAndExpand(h, npix*10 + kDR);
+  //
+  h  = new TH1F(Form("dtxODD_npix%d",npix),Form("dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetDirectory(0);
+  h->SetLineColor(kRed);
+  harr->AddAtAndExpand(h, npix*10 + kDTXodd);
+  h  = new TH1F(Form("dtxEVN_npix%d",npix),Form("dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetDirectory(0);
+  h->SetLineColor(kBlue);
+  harr->AddAtAndExpand(h, npix*10 + kDTXeven);
+  //
+  h  = new TH1F(Form("dtz_npix%d",npix),Form("dtz_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetLineColor(kGreen);
+  h->SetDirectory(0);
+  harr->AddAtAndExpand(h, npix*10 + kDTZ);
+  //
+  //
+  h  = new TH1F(Form("SPL_dtxODD_npix%d",npix),Form("SPL_dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetLineColor(kMagenta);
+  h->SetFillColor(kMagenta);
+  h->SetFillStyle(3001);  
+  h->SetLineStyle(2);
+  h->SetDirectory(0);
+
+  harr->AddAtAndExpand(h, npix*10 + kDTXoddSPL);
+  h  = new TH1F(Form("SPL_dtxEVN_npix%d",npix),Form("SPL_dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
+  h->SetLineColor(kCyan);
+  h->SetFillColor(kCyan);
+  h->SetFillStyle(3006);  
+  h->SetLineStyle(2);
+  h->SetDirectory(0);
+  harr->AddAtAndExpand(h, npix*10 + kDTXevenSPL);
+  //
+  h  = new TH1F(Form("SPL_dtz_npix%d",npix),Form("SPLdtz_npix%d",npix),nbin,-kdiff,kdiff);
+  harr->AddAtAndExpand(h, npix*10 + kDTZSPL);
+  h->SetDirectory(0);
+  //
+  h->SetLineColor(kGreen+2);
+  h->SetFillColor(kGreen+2);
+  h->SetLineStyle(2);
+  h->SetFillStyle(3001);
+  h = (TH1*)harr->At(idh);
+  if (!h) {printf("Unknown histo id=%d\n",idh); exit(1);}
+  return h;
+}
+
+
+
+TCanvas* DrawNP(int np, TObjArray* harr, TCanvas* cnv)
+{
+  if (!harr) harr = &histoArr;
+  if (!cnv) cnv = new TCanvas(Form("cnv%d",np),Form("cnv%d",np),900,700);
+  cnv->Clear();
+  cnv->Divide(2,1);
+  cnv->cd(1);
+  //
+  TH1* dxodd = (TH1*)harr->At(np*10+kDTXodd);
+  TH1* dxevn = (TH1*)harr->At(np*10+kDTXeven);
+  TH1* dxoddS =(TH1*)harr->At(np*10+kDTXoddSPL);
+  TH1* dxevnS =(TH1*)harr->At(np*10+kDTXevenSPL);
+  double max = TMath::Max(dxodd->GetMaximum(),dxevn->GetMaximum());
+  dxodd->SetMaximum(1.1*max);
+  dxodd->GetXaxis()->SetTitle("#DeltaX, #mum");
+  dxodd->SetTitle(Form("#DeltaX for clSize=%d",np));
+  dxodd->Fit("gaus","","");
+  dxevn->Fit("gaus","","sames");
+  //
+  dxoddS->Draw("sames");
+  dxevnS->Draw("sames");
+  //
+  gPad->Modified();
+  gPad->Update();
+  SetStPadPos(dxodd,0.75,0.97,0.8,1., -1,dxodd->GetLineColor());
+  SetStPadPos(dxevn,0.75,0.97,0.6,0.8, -1,dxevn->GetLineColor());
+  SetStPadPos(dxoddS,0.75,0.97,0.4,0.6, -1,dxoddS->GetLineColor());
+  SetStPadPos(dxevnS,0.75,0.97,0.2,0.4, -1,dxevnS->GetLineColor());
+  //
+  cnv->cd(2);
+  TH1* dz  = (TH1*)harr->At(np*10+kDTZ);
+  dz->SetTitle(Form("#DeltaZ for clSize=%d",np));
+  dz->GetXaxis()->SetTitle("#DeltaZ, #mum");
+  dz->Fit("gaus");
+  TH1* dzS = (TH1*)harr->At(np*10+kDTZSPL);
+  dz->Draw("sames");
+  gPad->Modified();
+  gPad->Update();
+  SetStPadPos(dz,0.75,0.97,0.8,1., -1, dz->GetLineColor());
+  SetStPadPos(dzS,0.75,0.97,0.5,0.7, -1, dzS->GetLineColor());
+  gPad->Modified();
+  gPad->Update();
+  //
+  cnv->cd();
+  return cnv;
+}
+
+
+TPaveStats* GetStPad(TH1* hst)
+{
+  TList *lst = hst->GetListOfFunctions();
+  if (!lst) return 0;
+  int nf = lst->GetSize();
+  for (int i=0;i<nf;i++) {
+    TPaveStats *fnc = (TPaveStats*) lst->At(i);
+    if (fnc->InheritsFrom("TPaveStats")) return fnc;
+  }
+  return 0;
+  //
+}
+
+
+TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl, Int_t col)
+{
+  TPaveStats* pad = GetStPad(hst);
+  if (!pad) return 0;
+  pad->SetX1NDC( x1 );
+  pad->SetX2NDC( x2 );
+  pad->SetY1NDC( y1 );
+  pad->SetY2NDC( y2 );
+  if (stl>=0) pad->SetFillStyle(stl);
+  if (col>=0) pad->SetTextColor(col);
+  pad->SetFillColor(0);
+  //
+  gPad->Modified();
+  return pad;
+}
\ No newline at end of file
diff --git a/ITS/UPGRADE/compression/PatternAnalysis/complete1.C b/ITS/UPGRADE/compression/PatternAnalysis/complete1.C
new file mode 100644 (file)
index 0000000..90fce92
--- /dev/null
@@ -0,0 +1,770 @@
+/*
+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");
+}
\ No newline at end of file
diff --git a/ITS/UPGRADE/compression/PatternAnalysis/errClus.C b/ITS/UPGRADE/compression/PatternAnalysis/errClus.C
new file mode 100644 (file)
index 0000000..9a46d56
--- /dev/null
@@ -0,0 +1,394 @@
+/*
+This macro creates for each pattern in the database (generated by BuildClTop.C) two histograms, for MChit-COG difference in x and z direction respectively
+This histograms are filled with the data from compClusHits.C. With a gaussian fit we get mean value, sigma, chi2, NDF and update the database of patterns
+with this information.
+*/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TArrayI.h"
+#include "TArrayF.h"
+#include "TCanvas.h"
+#include "TBits.h"
+#include "TF1.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"
+#include "TCanvas.h"
+#include "TPaveStats.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; //frequecy 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;
+
+/*Int_t err1counter=0;
+Int_t err2counter=0;
+Int_t err3counter=0;*/
+
+enum{kXhisto=0,kZhisto=1};
+
+void LoadDB(const char* fname);
+void CreateHistos(TObjArray* harr);
+TH1* GetHistoX(Int_t pattID, TObjArray* harr);
+TH1* GetHistoZ(Int_t pattID, TObjArray* harr);
+void DrawReport(TObjArray* harr);
+void DrawNP(int pattID, TObjArray* harr);
+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;
+
+void errClus(){
+
+       //importing data
+
+  LoadDB("clusterTopology.root");
+
+  nPatterns = (pattDB->GetEntriesFast());
+
+  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);
+
+  TFile *input = new TFile("clInfo.root","read");
+       TTree *clitsu = (TTree*) input->Get("clitsu");
+
+       Int_t 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);
+
+  CreateHistos(&histoArr);
+
+  for(Int_t iCl=0; iCl<ntotCl; iCl++){
+    clitsu->GetEntry(iCl);
+    GetHistoX(Summ.pattID, &histoArr)->Fill(Summ.dX);
+    GetHistoZ(Summ.pattID, &histoArr)->Fill(Summ.dZ);  
+  }
+  
+  DrawReport(&histoArr);
+
+  //printf("\n\nNDF is 0 : %d times\n\nNDF and fitstati not corresponding: %d times\n\nNDF not corresponding : %d times\n\n", err3counter,err1counter, err2counter);
+
+  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 two histograms for each pattern in the DB
+
+  printf("\n\n Creating histograms...   ");
+
+  Int_t histonumber=0;
+
+  if(!harr) harr = &histoArr;
+
+  for(Int_t i=0; i<nPatterns; i++){
+
+    TH1* h = new TH1F(Form("hX%d", i),Form("#DeltaX for pattID %d", i),100,-30,30);
+    h->SetDirectory(0);
+    harr->AddAtAndExpand(h, i);
+    histonumber++;
+  }
+
+  for(Int_t i=0; i<nPatterns; i++){
+
+    TH1* h = new TH1F(Form("hZ%d", i),Form("#DeltaZ for pattID %d", i),100,-30,30);
+    h->SetDirectory(0);
+    harr->AddAtAndExpand(h, i+nPatterns);
+    histonumber++;
+  }
+
+  printf(" %d histograms created, corresponding to %d patterns\n\n", histonumber, histonumber/2);
+}
+
+TH1* GetHistoX(Int_t pattID, TObjArray* harr){
+
+  TH1* h=0;
+  if(!harr) harr = &histoArr;
+
+  h=(TH1*)harr->At(pattID);
+  if (!h) {printf("Unknown histo id=%d\n",pattID); exit(1);}
+  return h;
+}
+
+TH1* GetHistoZ(Int_t pattID, TObjArray* harr){
+
+  Int_t zID=nPatterns+pattID;
+
+  TH1* h=0;
+  if(!harr) harr = &histoArr;
+  
+  h=(TH1*)harr->At(zID);
+  if (!h) {printf("Unknown histo id=%d\n",zID); exit(1);}
+  return h;
+}
+
+void DrawReport(TObjArray* harr) 
+{
+  // plot all the nPatterns cluster
+  for (int i=0;i<nPatterns;i++) {
+    DrawNP(i,harr);
+  }
+}
+
+void DrawNP(int pattID, TObjArray* harr)
+{
+  static TF1* gs = new TF1("gs","gaus",-500,500);
+  if (!harr) harr = &histoArr;
+
+  printf("\n\nProcessing #DeltaX of pattern %d...",pattID);
+
+  Int_t fitStatusX=0;
+
+  TH1* hdx = (TH1*)harr->At(pattID);
+  hdx->GetXaxis()->SetTitle("#DeltaX, #mum");
+  hdx->SetTitle(Form("#DeltaX for pattID %d",pattID));
+  gs->SetParameters(hdx->GetMaximum(),hdx->GetMean(),hdx->GetRMS());
+  if((hdx->GetEntries())<100) fitStatusX = hdx->Fit("gs","Ql0");
+  else fitStatusX = hdx->Fit("gs","Q0");
+
+  if(fitStatusX==0){
+    Double_t px1 = gs->GetParameter(1);
+    DeltaXmean[pattID]=px1;
+    Double_t ex1 = gs->GetParError(1);
+    DeltaXmeanErr[pattID]=ex1;
+    Double_t px2 = gs->GetParameter(2);
+    DeltaXsigma[pattID]=px2;
+    Double_t ex2 = gs->GetParError(2);
+    DeltaXsigmaErr[pattID]=ex2;
+    Double_t ChiSqx = gs->GetChisquare();
+    Chi2x[pattID] = ChiSqx;
+    Int_t varNDFx = gs->GetNDF();
+    if(varNDFx>=0)
+    NDFx[pattID] = varNDFx;
+    else{
+      /*if(varNDFx==0){
+        printf("\n\nNDF = 0, Chi2 = %e for pattern %d\n\n",ChiSqx,pattID);
+        err3counter++;
+      }*/
+    NDFx[pattID]=-1;
+    }
+  }
+  else{
+    DeltaXmean[pattID]=0;
+    DeltaXmeanErr[pattID]=0;
+    DeltaXsigma[pattID]=0;
+    DeltaXsigmaErr[pattID]=0;
+    Chi2x[pattID] = -1;   
+  }
+
+  printf("done!!\n\n");
+
+  printf("Processing #DeltaZ of pattern %d... ",pattID);
+
+  Int_t fitStatusZ=0;
+
+  TH1* hdz = (TH1*)harr->At(pattID+nPatterns);
+  hdz->SetTitle(Form("#DeltaZ for pattID %d",pattID));
+  hdz->GetXaxis()->SetTitle("#DeltaZ, #mum");
+  gs->SetParameters(hdz->GetMaximum(),hdz->GetMean(),hdz->GetRMS());
+  if((hdz->GetEntries())<100) fitStatusZ = hdz->Fit("gs","Ql0");
+  else fitStatusZ = hdz->Fit("gs","Q0");
+
+  if(fitStatusZ==0){
+    Double_t pz1=gs->GetParameter(1);
+    DeltaZmean[pattID]=pz1;
+    Double_t ez1 = gs->GetParError(1);
+    DeltaZmeanErr[pattID]=ez1;
+    Double_t pz2 = gs->GetParameter(2);
+    DeltaZsigma[pattID]=pz2;
+    Double_t ez2 = gs->GetParError(2);
+    DeltaZsigmaErr[pattID]=ez2;
+    Double_t ChiSqz = gs->GetChisquare();
+    Chi2z[pattID]=ChiSqz;
+    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;
+  }
+
+  /*if((NDFz[pattID]!=NDFx[pattID])&&fitStatusZ!=fitStatusX){
+    printf("\n\nERR: neither NDF nor FitStatus correspond!\n\n");
+    err1counter++;
+  }
+  else if(NDFz[pattID]!=NDFx[pattID]){
+    printf("\n\nERR: NDFx and NDFz donnot correspond!\n\n");
+    err2counter++;
+  }*/
+
+  printf("done!!\n\n");
+
+}
+
+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");
+}
diff --git a/ITS/UPGRADE/compression/PatternAnalysis/groupdef.C b/ITS/UPGRADE/compression/PatternAnalysis/groupdef.C
new file mode 100644 (file)
index 0000000..218601e
--- /dev/null
@@ -0,0 +1,477 @@
+/*
+In this macro we create groups of patterns with similar characteristics. First we load the pattern database with information about MC-hit.
+We do not group patterns with a frequency below a treshold (frequencyTreshold).
+We define 7 different ID. 2 for the sigma of the MChit-COG, x and z direction (sigmaXID and sigma ZID);
+2 for the distance between COG and centre of the pixel, x and z (shiftXID and shiftZID);
+2 for the distance between MChit and centre of the pixel, x and z direction (biasXID ad biasZID).
+1 ID for the number of fired pixels.
+We define the binning (number/width of bins) and assign the previous IDs.
+We decide which the grouping method we want to use: sigma, shift and number of pixels (kShift) or sigma, bias and number of pixels (kBias)
+Finally we assign group ID. Frequent patterns (above the treshold), form a one-pattern group. The tohers are in the same group if they have
+the same IDs.
+A fil.txt is print qith he inormation for each apttern, including the IDs and the pattID of the patterns in the same group.
+
+*/
+
+#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 "TPad.h"
+#include "TPavesText.h"
+#include "TLatex.h"
+#include "TBits.h"
+#include "TGraph.h"
+#include "TStopwatch.h"
+#include "TMath.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"
+
+#endif
+
+TObjArray histoArr;
+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
+TVectorF* DeltaZmean=0; //mean value of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
+TVectorF* DeltaXmean=0; //mean value of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
+TVectorF* DeltaZsigma=0; //sigma of the difference between MChit and COG z coordinate in micron derived from a gaussian fit
+TVectorF* DeltaXsigma=0; //sigma of the difference between MChit and COG x coordinate in micron derived from a gaussian fit
+TVectorF* DeltaZmeanErr=0;
+TVectorF* DeltaXmeanErr=0;
+TVectorF* DeltaZsigmaErr=0;
+TVectorF* DeltaXsigmaErr=0;
+
+//Defining the ID to create the groups.
+TArrayI sigmaXID;
+TArrayI sigmaZID;
+TArrayI shiftXID;
+TArrayI shiftZID;
+TArrayI biasXID;
+TArrayI biasZID;
+TArrayI NPixID;
+TArrayI groupID;
+Float_t totalGroupFreq=0;
+
+void LoadDB(const char* namefile);
+
+Int_t nPatterns=0;
+
+//Defining the frequency treshold under which ti group patterns
+Float_t frequencyTreshold = 0.001;
+//Defining the bins
+Int_t NumberofSigmaXbins=1000;
+Int_t NumberofSigmaZbins=1000;
+Int_t NumberofShiftXbins=20;
+Int_t NumberofShiftZbins=20;
+Int_t NumberofBiasXbins=2000;
+Int_t NumberofBiasZbins=2000;
+//Defining the boundaries of the bins concerning the number of pixel
+const Int_t    limitsnumber = 4;
+//Defining the width of the bins to store patterns with similar sigma/shift/bias (in micron)
+Float_t sigmaXwidth=3;
+Float_t sigmaZwidth=3;
+Float_t biasXwidth=5;
+Float_t biasZwidth=5;
+Float_t shiftXwidth=1./NumberofShiftXbins; //(fraction of the pitch)
+Float_t shiftZwidth=1./NumberofShiftZbins; //(fraction of the pitch)
+
+Int_t nPixlimits[limitsnumber] = {4,10,25,50};
+
+enum{kShift=0, kBias=1};
+//Select kShift to group wrt sigma & COG-centreOfThePixel distance
+//Select kBias to group wrt sigma & MChit-centreOfThePixel distance
+
+Int_t groupingMethod = kShift;
+
+Int_t tempgroupID=0;
+
+void groupdef(){
+
+       //Importing Data
+       LoadDB("clusterTopology2.root");
+
+       //Define moudule segmentation
+
+       AliGeomManager::LoadGeometry("geometry.root");
+       AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
+       AliITSUClusterPix::SetGeom(gm);
+       const AliITSsegmentation* segm = gm->GetSegmentation(0);
+       Float_t pitchX = segm->Dpx(0)*10000; // for pitch in X in micron
+       printf("pitchX: %f\n",pitchX);
+       Float_t pitchZ = segm->Dpz(0)*10000; // for pitch in Z in micron
+       printf("pitchX: %f\n",pitchX);
+       delete gm;
+
+       //Setting the number of patterns
+       nPatterns = (pattDB->GetEntriesFast());
+
+       sigmaXID.Set(nPatterns);
+       sigmaZID.Set(nPatterns);
+
+       //Assign -2 to frequent patterns, not to group, and -1 to rare clusters, -3 to patterns with sigma set to zero
+       for(Int_t ID=0; ID<nPatterns; ID++){
+               printf("Assign temporary sigma ID to pattern %d... ",ID);
+               if((*pattFR)[ID]>frequencyTreshold){
+                       sigmaXID[ID]=-2;//In order not to consider it within the groups
+                       sigmaZID[ID]=-2;//In order not to consider it within the groups
+               }
+               else{
+                       sigmaXID[ID]=-1;
+                       sigmaZID[ID]=-1;
+                       totalGroupFreq+=(*pattFR)[ID];
+               }
+               printf("done\n\n");
+       }
+
+       //Assign to similar patterns the same sigmaXID
+       for(Int_t i=0; i<nPatterns; i++){
+               printf("Assign sigmaXID to pattern %d... ",i);
+               if(sigmaXID[i]==-1){
+                       if((*DeltaXsigma)[i]==0){
+                               sigmaXID[i]=-3; // In order not to cosider it within the groups
+                       }
+                       else{
+                               for(Int_t j=0; j<NumberofSigmaXbins; j++){
+                                       if(j*sigmaXwidth < (*DeltaXsigma)[i] && (*DeltaXsigma)[i]<= (j+1)*sigmaXwidth){
+                                               sigmaXID[i]=j+1;
+                                               break;
+                                       }
+                               }
+                       }
+               }
+               printf("done!!\n\n");
+       }
+
+       //Assign to similar patterns the same sigmaZID
+       for(Int_t i=0; i<nPatterns; i++){
+               printf("Assign sigmaZID to pattern %d... ",i);
+               if(sigmaZID[i]==-1){
+                       if((*DeltaZsigma)[i]==0){
+                               sigmaZID[i]=-3; // In order not to cosider it within the groups
+                       }
+                       else{
+                               for(int j=0; j<NumberofSigmaZbins ; j++){
+                                       if(j*sigmaZwidth < (*DeltaZsigma)[i] && (*DeltaZsigma)[i]<= (j+1)*sigmaZwidth){
+                                               sigmaZID[i]=j+1;
+                                               break;
+                                       }
+                               }
+                       }
+               }
+               printf("done!!\n\n");
+       }
+
+       //assigning shiftID
+
+       shiftXID.Set(nPatterns);
+       shiftZID.Set(nPatterns);
+
+       for(Int_t i=0; i<nPatterns; i++){
+               printf("Assign shiftXID to pattern %d... ",i);
+               
+               for(int j=-NumberofShiftXbins/2; j<NumberofShiftXbins/2; j++){
+                       
+                       if(j*shiftXwidth < (*xCentreShift)[i] && (*xCentreShift)[i]<= (j+1)*shiftXwidth){
+                               shiftXID[i]=j+1;
+                               printf("done!!\n\n");
+                               break;
+                       }       
+               }
+       }
+
+       for(Int_t i=0; i<nPatterns; i++){
+               printf("Assign shiftZID to pattern %d... ",i);
+               
+               for(int j=-NumberofShiftZbins/2; j<NumberofShiftZbins/2; j++){
+                       if(j*shiftZwidth < (*zCentreShift)[i] && (*zCentreShift)[i]<= (j+1)*shiftZwidth){
+                               shiftZID[i]=j+1;
+                               printf("done!!\n\n");
+                               break;
+                       }
+               }       
+       }
+       
+       //assigning BiasID
+
+       biasXID.Set(nPatterns);
+       biasZID.Set(nPatterns);
+
+       //Setting all the bias ID to zero
+
+       for(Int_t i=0; i<nPatterns; i++){
+               biasXID[i]=0;
+               biasZID[i]=0;
+       }
+
+       for(Int_t i=0; i<nPatterns; i++){
+               printf("Assign biasXID to pattern %d... ",i);
+               for(Int_t j=-NumberofBiasXbins/2; j<NumberofBiasXbins/2; j++){
+                       if(j*biasXwidth < ((*DeltaXmean)[i]+((*xCentreShift)[i]*pitchX))
+                               && ((*DeltaXmean)[i]+((*xCentreShift)[i]*pitchX))<= (j+1)*biasXwidth){
+                       biasXID[i]=j+1;
+                       break;
+                       }
+               }
+               printf("done!!\n\n");
+       }
+
+       for(Int_t i=0; i<nPatterns; i++){
+               biasXID[i]=0;
+               biasZID[i]=0;
+       }
+
+       for(Int_t i=0; i<nPatterns; i++){
+               printf("Assign biasZID to pattern %d... ",i);
+               for(Int_t j=-NumberofBiasZbins/2; j<NumberofBiasZbins/2; j++){
+                       if(j*biasZwidth < ((*DeltaZmean)[i]+((*zCentreShift)[i]*pitchZ)) 
+                               && ((*DeltaZmean)[i]+((*zCentreShift)[i]*pitchZ))<= (j+1)*biasZwidth){
+                       biasZID[i]=j+1;
+                       break;
+                       }
+               }
+               printf("done!!\n\n");
+       }
+
+       
+
+       //Assigning NPixID
+
+       NPixID.Set(nPatterns);
+
+       for(Int_t i=0; i<nPatterns; i++){
+               printf("Assigning NPixID to pattern %d...", i);
+               if((*NPix)[i]<=nPixlimits[0]){
+                       NPixID[i]=0;
+                       printf("done!!\n\n");
+               }
+               else if(nPixlimits[0]<(*NPix)[i] && (*NPix)[i]<=nPixlimits[limitsnumber-1]){
+                       for(Int_t j=0; j<nPatterns; j++ ){
+                               if(nPixlimits[j]<(*NPix)[i] && (*NPix)[i]<=nPixlimits[j+1]){
+                                       NPixID[i]=j+1;
+                                       printf("done!!\n\n");
+                                       break;
+                               }
+                       }
+               }
+               else if((*NPix)[i]>nPixlimits[limitsnumber-1]){
+                       NPixID[i]=limitsnumber+1;
+                       printf("done!!\n\n");
+               }
+       }
+
+       //Assigning groupID
+
+       groupID.Set(nPatterns);
+
+       //Assign -2 to frequent patterns, not to group, and -1 to rare clusters
+       for(Int_t ID=0; ID<nPatterns; ID++){
+               printf("Assign temporary group ID to pattern %d... ",ID);
+               groupID[ID]=-1;
+               printf("done\n\n");
+       }
+       Int_t k=0;
+       while((*pattFR)[k]>frequencyTreshold){
+               groupID[k]=tempgroupID;
+               tempgroupID++;
+               k++;
+       }
+
+       if(groupingMethod==kShift){
+               for(Int_t i=0; i<nPatterns; i++){
+                       if(groupID[i]!=-1) continue;    
+                       groupID[i]=tempgroupID;
+                       printf("Assigning group ID %d... ",tempgroupID);
+                       for(Int_t j=i+1; j<nPatterns; j++){
+                               if(sigmaXID[j]==-3){
+                                       groupID[j]=-1;
+                                       continue;
+                               }
+                               else if(sigmaXID[j]==sigmaXID[i] && sigmaZID[j]==sigmaZID[i] && 
+                                       shiftXID[j]==shiftXID[i] && shiftZID[j]==shiftZID[i] &&
+                                       NPixID[i]==NPixID[j]) groupID[j]=tempgroupID;
+                       }
+                       printf("done!!\n\n");
+                       tempgroupID++;
+               }
+       }
+       else if(groupingMethod==kBias){
+               for(Int_t i=0; i<nPatterns; i++){
+                       if(groupID[i]!=-1) continue;
+                       groupID[i]=tempgroupID;
+                       printf("Assigning group ID %d... ",tempgroupID);
+                       for(Int_t j=i+1; j<nPatterns; j++){
+                               if(sigmaZID[j]==-3){
+                                       groupID[j]=-1;
+                                       continue;
+                               }
+                               else if(sigmaXID[j]==sigmaXID[i] && sigmaZID[j]==sigmaZID[i]
+                                       && biasXID[j]==biasXID[i] && biasZID[j]==biasZID[i]
+                                       && NPixID[i]==NPixID[j]) groupID[j]=tempgroupID;
+                       }
+                       printf("done!!\n\n");
+                       tempgroupID++;
+               }
+       }
+
+       ofstream a("groupdef.txt");
+               
+       //setw(55) << "patterns in the group\n" << endl;
+
+       a << Form("A NEGATIVE ID means that the pattern is not in a group.") << endl << endl <<
+               Form("EXCEPTION: biasID CAN be negative") << endl <<
+               "\n\n......................................................................................." << 
+               "................................................................................................\n\n";
+
+       for(int i=0; i<nPatterns; i++){
+
+               printf("Writing info about pattern %d ...", i);
+
+               a <<  setw(30) << Form("pattID: %d",i) <<  setw(30) << Form("freq: %f", (*pattFR)[i]) << setw(30) << Form("NPix: %d",Int_t((*NPix)[i]))<< setw(45) << Form("NRow: %d ",Int_t((*NRow)[i])) << setw(45) <<
+               Form("NCol: %d", Int_t((*NCol)[i])) << endl << endl
+               << setw(45) << Form("DeltaXmean: %f (%f)",(*DeltaXmean)[i],(*DeltaXmeanErr)[i]) << setw(45) << 
+               Form("DeltaZmean: %f (%f)",(*DeltaZmean)[i],(*DeltaZmeanErr)[i])<<
+               setw(45) << Form("DeltaXsigma: %f (%f)",(*DeltaXsigma)[i],(*DeltaXsigmaErr)[i]) << setw(45) <<
+               Form("DeltaZsigma: %f (%f)",(*DeltaZsigma)[i],(*DeltaZsigmaErr)[i]) << endl << endl << setw(30) << 
+               Form("xShift: %f",(*xCentreShift)[i]) <<  setw(45) << Form("zShift: %f",(*zCentreShift)[i]) << setw(45) << Form("BiasX(MChit-centrePix): %f",(*DeltaXmean)[i]+((*xCentreShift)[i]*pitchX)) << setw(45) <<
+               Form("BiasZ(MChit-centrePix): %f",(*DeltaZmean)[i]+((*zCentreShift)[i]*pitchZ)) << endl << endl <<
+               setw(30) << Form("sigmaXID: %d",sigmaXID[i]) << setw(15) << Form("sigmaZID: %d",sigmaZID[i]) << setw(15) <<
+               Form("shiftXID: %d",shiftXID[i]) << setw(15) << Form("shiftZID: %d",shiftZID[i]) << setw(15)<< Form("biasXID: %d",biasXID[i]) <<
+               setw(15) << Form("biasZID: %d",biasZID[i]) << setw(15) << Form("NPixID: %d", NPixID[i]) << endl << endl << setw(30) << Form("groupID: %d", groupID[i])
+               << endl << endl << setw(30) << "patterns in this group:\n\n";
+
+               Int_t matchcounter=0;
+
+               for(int j=0; j<nPatterns;  j++){
+                       if(matchcounter!=0 && (matchcounter%20)==0){
+                               a << endl;
+                               matchcounter=0;
+                       }
+                       if(groupID[j]==groupID[i]) {
+                               a << j << ", ";
+                               matchcounter++;
+                       }
+               }
+               
+               a << "\n\n......................................................................................." << 
+               "................................................................................................\n\n";
+
+
+               printf("done!\n\n");
+       }
+
+       a.close();
+
+       printf("%d groups found!!!!\n\n",tempgroupID);
+
+       printf("\n\nThe total frequency of the patterns in group is %f\n\n",totalGroupFreq);
+
+       /*
+
+       TVectorF groupDeltaX(tempGID+100);
+       TVectorF groupShiftX(tempGID+100);
+       TVectorF groupDeltaZ(tempGID+100);
+       TVectorF groupShiftZ(tempGID+100);
+       TVectorF patternNum(tempGID+100);
+       groupDeltaX.Zero();
+       groupShiftX.Zero();
+       groupDeltaZ.Zero();
+       groupShiftZ.Zero();
+       patternNum.Zero();
+
+       ofstream b("groupDB.txt");
+
+       b<<setw(15)<<"groupID"<<setw(25)<<"patterns in the group"<<setw(15)<<
+       "DeltaX"<<setw(15)<<"ShiftX"<<setw(15)<<"DeltaZ"<<setw(15)<<"ShiftZ\n"<<endl;
+
+       TCanvas* c = new TCanvas("c","Patterns groups",900,600);
+       TH1F* h = new TH1F("h","Patterns groups",tempGID,-0.5, tempGID-0.5);
+       h->GetXaxis()->SetTitle("group ID");
+       h->SetStats(0);
+
+       for(Int_t gid=0; gid<tempGID; gid++){
+
+               Int_t PatternNumberInGroup=0;
+               Float_t freqSum=0.;//It is the sum of the frequencies of the patterns in the same group
+
+               for(Int_t pattID=0; pattID<nPatterns; pattID++){
+                       if (groupID[pattID]==gid)
+                       {
+                               groupDeltaX[gid]+=(*DeltaXsigma)[pattID]*(*pattFR)[pattID];
+                               groupShiftX[gid]+=(*xCentreShift)[pattID]*(*pattFR)[pattID];
+                               groupDeltaZ[gid]+=(*DeltaZsigma)[pattID]*(*pattFR)[pattID];
+                               groupShiftZ[gid]+=(*zCentreShift)[pattID]*(*pattFR)[pattID];
+                               freqSum+=(*pattFR)[pattID];
+                               PatternNumberInGroup++;
+
+                               h->Fill(gid);
+                       }
+               }
+
+               groupDeltaX[gid]=groupDeltaX[gid]/freqSum;
+               groupShiftX[gid]=groupShiftX[gid]/freqSum;
+               groupDeltaZ[gid]=groupDeltaZ[gid]/freqSum;
+               groupShiftZ[gid]=groupShiftZ[gid]/freqSum;
+               patternNum[gid]=PatternNumberInGroup;
+
+               b<<setw(15)<<gid<<setw(25)<<patternNum[gid]<<setw(15)<<
+               groupDeltaX[gid]<<setw(15)<<groupShiftX[gid]<<setw(15)<<
+               groupDeltaZ[gid]<<setw(15)<<groupShiftZ[gid]<<"\n"<<endl;
+       }
+
+       c->cd();
+       h->Draw();
+
+       TPavesText* info = new TPavesText(0.5,1,0.5,1,1,"nb");
+       info->AddText(Form("Number of groups: %d",tempGID));
+       info->AddText("#delta_{X}<5 & #delta_{Z}<5");
+       info->AddText("#DeltaX<0.05 & #DeltaZ<0.05");
+       c->cd();
+       info->Draw();
+
+       b.close();
+       */
+}
+
+void LoadDB(const char* fname){
+
+  printf("\n\nLoading DB... ");
+  // 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");
+  DeltaXmean =(TVectorF*)fl->Get("DeltaXmean");
+  DeltaZmean =(TVectorF*)fl->Get("DeltaZmean");
+  DeltaXsigma =(TVectorF*)fl->Get("DeltaXsigma");
+  DeltaZsigma =(TVectorF*)fl->Get("DeltaZsigma");
+  DeltaXmeanErr =(TVectorF*)fl->Get("DeltaXmeanErr");
+  DeltaZmeanErr =(TVectorF*)fl->Get("DeltaZmeanErr");
+  DeltaZsigmaErr =(TVectorF*)fl->Get("DeltaXsigmaErr");
+  DeltaXsigmaErr =(TVectorF*)fl->Get("DeltaZsigmaErr");
+  printf("done!!\n\n");
+}
\ No newline at end of file