ROOT logo
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "TObjArray.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.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

TFile* fl=0;
TTree* tr=0;
TObjArray patterns;
TArrayF pattCount;
TArrayI sortIndex,sortColRow;
TClonesArray *modClStore;
Int_t nPatterns;

TH1F* hPattFreq = 0;
TH1F* hVolDiffFreq = 0;
TH2F* hColDiffFreq = 0;
TH2F* hColFreq = 0;
TH2F* hRowDiffFreq = 0;
TH1F* hDataVol = 0;
int lastVol=9999999,volDif=0,volID=0;
TBits* bitTop=0;


const int    kNRow   = 650*2; 
const int    kNCol   = 750*2;

const int    kNBinClVol = 6; // number of module multiplicity bins
const double kBinsClVol[kNBinClVol+1] = {0.5, 5.5, 15.5, 31.5, 65.5, 125.5, 500.5};
//
void ProcessModule(int nc);


void anatop(const char* topTreeFile="clusterTopology.root")
{
  // analyze tree with cluster topologies
  AliGeomManager::LoadGeometry("geometry.root");
  AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
  Int_t nLayers = gm->GetNLayers();
  //
  fl = TFile::Open(topTreeFile);
  tr = (TTree*)fl->Get("clTop");
  UShort_t mnCol,mnRow;
  tr->SetBranchAddress("tp",&bitTop);
  tr->SetBranchAddress("volID",&volID);
  tr->SetBranchAddress("mnCol",&mnCol);
  tr->SetBranchAddress("mnRow",&mnRow);
  //
  int nCl = tr->GetEntries();
  //
  hDataVol     = new TH1F("hDataVol","Data Volume", kNBinClVol,kBinsClVol);
  hVolDiffFreq = new TH1F("hVolDiffFreq","Delta Vol ID",gm->GetNChips(),0,gm->GetNChips());
  hColDiffFreq = new TH2F("hColDiffFreq","Delta Col",kNBinClVol,kBinsClVol,2*kNCol,-kNCol,kNCol);
  hColFreq     = new TH2F("hColFreq","Col"          ,kNBinClVol,kBinsClVol,kNCol,0,kNCol);
  hRowDiffFreq = new TH2F("hRowDiffFreq","Delta Row",kNBinClVol,kBinsClVol,kNRow,0,kNRow);
  //
  printf("%d clusters in the tree\n",nCl);
  modClStore = new TClonesArray("TBits");
  modClStore->ExpandCreate(5000);
  sortIndex.Set(5000);
  sortColRow.Set(5000);
  //
  int lastEv=-1,nClVol=0;
  Bool_t newEv;
  nPatterns  = 0;
  //
  if (nCl<1) {
    printf("No clusters\n");
    return;
  }
  for (int icl=0;icl<nCl;icl++) {
    if (icl%(nCl/10)==0) {
      printf("Done %d clusters of %d, seen %d patterns\n",icl,nCl,nPatterns);
    }
    newEv = kFALSE;
    tr->GetEntry(icl);
    //
    if (volID != lastVol) {
      if (volID<lastVol) {
	newEv = kTRUE;
	lastVol = 0;
	lastEv++;
      }
      // store previous module info in the frequency histos
      if (nClVol>0) ProcessModule(nClVol);
      nClVol = 0;
      //
      //
    }
    //-------------------------
    //
    *((TBits*)modClStore->At(nClVol)) = *bitTop;
    sortColRow[nClVol] = (mnRow<<16) + mnCol;
    nClVol++;

    //-------------------------
  }
  //
  printf("Sorting %d patterns\n",nPatterns);
  // sort patterns
  sortIndex.Set(nPatterns);
  TMath::Sort(nPatterns,pattCount.GetArray(),sortIndex.GetArray());
  //
  hPattFreq = new TH1F("hPattFreq","Patterns Frequency",nPatterns,0,nPatterns);
  for (int i=0;i<nPatterns;i++) {
    int ind = sortIndex[i];
    TBits* patt = (TBits*)patterns[ind];
    pattCount[ind] /= nCl;
    Double_t freq = pattCount[ind];
    hPattFreq->SetBinContent(i+1,freq);
    //
    int nrow = (patt->GetUniqueID()>>8) & 0xff;
    int ncol = patt->GetUniqueID() & 0xff;
    printf("\nPattern %5d, Freq: %e Ncol=%d Nrow=%d\n",i,freq,ncol,nrow);
    for (int ir=0;ir<nrow;ir++) {
      for (int ic=0;ic<ncol;ic++) {
	int bt = ir*ncol + ic;
	printf("%c",patt->TestBitNumber(bt) ? '*':'.');
      }
      printf("\n");
    }
  }

}


//___________________________________________
void ProcessModule(int nc)
{
  // sort clusters
  int lastCol = 0;
  int lastRow = 0;
  TMath::Sort(nc,sortColRow.GetArray(),sortIndex.GetArray(),kFALSE);
  hVolDiffFreq->Fill(volID - lastVol);
  for (int iclv=0;iclv<nc;iclv++) {
    hDataVol->Fill(nc);
    int ind = sortIndex[iclv];
    int col = sortColRow[ind]&0xffff;
    int row = (sortColRow[ind]>>16)&0xffff;
    TBits *refPat=(TBits*)modClStore->At(ind);
    //    printf("%5d %5d | %4d %4d  | %4d %4d\n",volID,lastVol, col,row, lastCol,lastRow);
    hColDiffFreq->Fill(nc, col - lastCol);
    hColFreq->Fill(nc, col);
    hRowDiffFreq->Fill(nc, row - lastRow);
    lastCol = col;
    lastRow = row;
    //
    // check if the pattern is already encountered
    Bool_t newPatt = kTRUE;
    for (int ip=0;ip<nPatterns;ip++) {
      TBits* patt = (TBits*)patterns.At(ip);
      if ( *patt == *refPat && patt->GetUniqueID()==refPat->GetUniqueID()) {
	// require that the matrix size and the bit patterns are equal
	newPatt = kFALSE;
	pattCount[ip]++;
	break;
      }
    }
    //
    if (newPatt) {
      TBits* pt = new TBits(*refPat);
      patterns.AddLast(pt);
      if (pattCount.GetSize()<nPatterns+100) {
	pattCount.Set(nPatterns+100);
      }
      pattCount[nPatterns++] = 1;
    }
    //-------------------------

  }  
  lastVol = volID;
  //
}
 anatop.C:1
 anatop.C:2
 anatop.C:3
 anatop.C:4
 anatop.C:5
 anatop.C:6
 anatop.C:7
 anatop.C:8
 anatop.C:9
 anatop.C:10
 anatop.C:11
 anatop.C:12
 anatop.C:13
 anatop.C:14
 anatop.C:15
 anatop.C:16
 anatop.C:17
 anatop.C:18
 anatop.C:19
 anatop.C:20
 anatop.C:21
 anatop.C:22
 anatop.C:23
 anatop.C:24
 anatop.C:25
 anatop.C:26
 anatop.C:27
 anatop.C:28
 anatop.C:29
 anatop.C:30
 anatop.C:31
 anatop.C:32
 anatop.C:33
 anatop.C:34
 anatop.C:35
 anatop.C:36
 anatop.C:37
 anatop.C:38
 anatop.C:39
 anatop.C:40
 anatop.C:41
 anatop.C:42
 anatop.C:43
 anatop.C:44
 anatop.C:45
 anatop.C:46
 anatop.C:47
 anatop.C:48
 anatop.C:49
 anatop.C:50
 anatop.C:51
 anatop.C:52
 anatop.C:53
 anatop.C:54
 anatop.C:55
 anatop.C:56
 anatop.C:57
 anatop.C:58
 anatop.C:59
 anatop.C:60
 anatop.C:61
 anatop.C:62
 anatop.C:63
 anatop.C:64
 anatop.C:65
 anatop.C:66
 anatop.C:67
 anatop.C:68
 anatop.C:69
 anatop.C:70
 anatop.C:71
 anatop.C:72
 anatop.C:73
 anatop.C:74
 anatop.C:75
 anatop.C:76
 anatop.C:77
 anatop.C:78
 anatop.C:79
 anatop.C:80
 anatop.C:81
 anatop.C:82
 anatop.C:83
 anatop.C:84
 anatop.C:85
 anatop.C:86
 anatop.C:87
 anatop.C:88
 anatop.C:89
 anatop.C:90
 anatop.C:91
 anatop.C:92
 anatop.C:93
 anatop.C:94
 anatop.C:95
 anatop.C:96
 anatop.C:97
 anatop.C:98
 anatop.C:99
 anatop.C:100
 anatop.C:101
 anatop.C:102
 anatop.C:103
 anatop.C:104
 anatop.C:105
 anatop.C:106
 anatop.C:107
 anatop.C:108
 anatop.C:109
 anatop.C:110
 anatop.C:111
 anatop.C:112
 anatop.C:113
 anatop.C:114
 anatop.C:115
 anatop.C:116
 anatop.C:117
 anatop.C:118
 anatop.C:119
 anatop.C:120
 anatop.C:121
 anatop.C:122
 anatop.C:123
 anatop.C:124
 anatop.C:125
 anatop.C:126
 anatop.C:127
 anatop.C:128
 anatop.C:129
 anatop.C:130
 anatop.C:131
 anatop.C:132
 anatop.C:133
 anatop.C:134
 anatop.C:135
 anatop.C:136
 anatop.C:137
 anatop.C:138
 anatop.C:139
 anatop.C:140
 anatop.C:141
 anatop.C:142
 anatop.C:143
 anatop.C:144
 anatop.C:145
 anatop.C:146
 anatop.C:147
 anatop.C:148
 anatop.C:149
 anatop.C:150
 anatop.C:151
 anatop.C:152
 anatop.C:153
 anatop.C:154
 anatop.C:155
 anatop.C:156
 anatop.C:157
 anatop.C:158
 anatop.C:159
 anatop.C:160
 anatop.C:161
 anatop.C:162
 anatop.C:163
 anatop.C:164
 anatop.C:165
 anatop.C:166
 anatop.C:167
 anatop.C:168
 anatop.C:169
 anatop.C:170
 anatop.C:171
 anatop.C:172
 anatop.C:173
 anatop.C:174
 anatop.C:175
 anatop.C:176
 anatop.C:177
 anatop.C:178
 anatop.C:179
 anatop.C:180
 anatop.C:181
 anatop.C:182
 anatop.C:183
 anatop.C:184
 anatop.C:185
 anatop.C:186
 anatop.C:187