ITS/UPGRADE/compression/PatternAnalysis/SecondcompClusHits.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include "TObjArray.h"
   3 #include "TFile.h"
   4 #include "TTree.h"
   5 #include "TH1F.h"
   6 #include "TBits.h"
   7 #include "../ITS/UPGRADE/AliITSUClusterPix.h"
   8 #include "../ITS/UPGRADE/AliITSURecoLayer.h"
   9 #include "../ITS/UPGRADE/AliITSURecoDet.h"
  10 #include "../ITS/UPGRADE/AliITSUHit.h"
  11 #include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
  12 #include "AliITSsegmentation.h"
  13 #include "AliGeomManager.h"
  14 #include "AliStack.h"
  15 #include "AliLoader.h"
  16 #include "AliCDBManager.h"
  17
  18 #include "TROOT.h"
  19 #include "TStyle.h"
  20 #include "TGeoMatrix.h"
  21 #include "TParticle.h"
  22 #include "TCanvas.h"
  23 #include "TPaveStats.h"
  24 #include "TClonesArray.h"
  25
  26 #endif
  27
  28 TObjArray *pattDB=0;
  29 TVectorF*  pattFR=0;
  30 TVectorF*  xCentrePix=0;
  31 TVectorF*  zCentrePix=0;
  32 TVectorF*  xCentreShift=0;
  33 TVectorF*  zCentreShift=0;
  34 TObjArray histoArr;
  35 enum {kNPixAll=0,kNPixSPL=1,kDR=0,kDTXodd,kDTXeven,kDTZ, kDTXoddSPL,kDTXevenSPL,kDTZSPL};
  36
  37 TPaveStats* GetStPad(TH1* hst);
  38 TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl=-1,Int_t col=-1);
  39 TCanvas* DrawNP(int np, TObjArray* harr=0, TCanvas* cnv=0);
  40 TH1* GetHistoClSize(int npix,int id,TObjArray* harr=0);
  41 void DrawReport(const char* psname, TObjArray* harr=0);
  42 void LoadDB(const char* fname);
  43
  44
  45 typedef struct {
  46   Int_t evID;
  47   Int_t volID;
  48   Int_t lrID;
  49   Int_t clID;
  50   Int_t nPix;
  51   Int_t nX;
  52   Int_t nZ;
  53   Int_t q;
  54   Float_t pt;
  55   Float_t eta;
  56   Float_t phi;
  57   Float_t xyz[3];
  58   Float_t dX;
  59   Float_t dY;
  60   Float_t dZ;
  61   Bool_t split;
  62   Bool_t prim;
  63   Int_t  pdg;
  64   Int_t  ntr;
  65   Float_t alpha; // alpha is the angle in y-radius plane in local frame
  66   Float_t beta;  // beta is the angle in xz plane, taken from z axis, growing counterclockwise
  67   Int_t nRowPatt;
  68   Int_t nColPatt;
  69   Int_t pattID;
  70   Float_t freq;
  71   Float_t xCen;
  72   Float_t zCen;
  73   Float_t zShift;
  74   Float_t xShift;
  75 } clSumm;
  76
  77 TBits clTop;
  78
  79 void SecondcompClusHits(int nev=-1)
  80 {
  81   const int kSplit=0x1<<22;
  82   const int kSplCheck=0x1<<23;
  83   //
  84   gSystem->Load("libITSUpgradeBase");
  85   gSystem->Load("libITSUpgradeSim");
  86   gSystem->Load("libITSUpgradeRec");
  87   gROOT->SetStyle("Plain");
  88
  89   AliCDBManager* man = AliCDBManager::Instance();
  90   man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  91   man->SetSpecificStorage("GRP/GRP/Data",
  92                           Form("local://%s",gSystem->pwd()));
  93   man->SetSpecificStorage("ITS/Align/Data",
  94                           Form("local://%s",gSystem->pwd()));
  95   man->SetSpecificStorage("ITS/Calib/RecoParam",
  96                           Form("local://%s",gSystem->pwd()));
  97   man->SetRun(0);
  98
  99
 100   gAlice=NULL;
 101   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
 102   runLoader->LoadgAlice();
 103
 104   gAlice = runLoader->GetAliRun();
 105
 106   runLoader->LoadHeader();
 107   runLoader->LoadKinematics();
 108   runLoader->LoadRecPoints();
 109   runLoader->LoadSDigits();
 110   runLoader->LoadHits();
 111
 112   AliLoader *dl = runLoader->GetDetectorLoader("ITS");
 113
 114   AliGeomManager::LoadGeometry("geometry.root");
 115   TObjArray algITS;
 116   AliGeomManager::LoadAlignObjsFromCDBSingleDet("ITS",algITS);
 117   AliGeomManager::ApplyAlignObjsToGeom(algITS);
 118   //
 119   AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
 120   AliITSUClusterPix::SetGeom(gm);
 121   //
 122   AliITSURecoDet *its = new AliITSURecoDet(gm, "ITSinterface");
 123   its->CreateClusterArrays();
 124   //
 125   Double_t xg1,yg1,zg1=0.,xg0,yg0,zg0=0.,tg0;
 126   Double_t xExit,yExit,zExit,xEnt,yEnt,zEnt,tof1;
 127   //
 128   TTree * cluTree = 0x0;
 129   TTree *hitTree = 0x0;
 130   TClonesArray *hitList=new TClonesArray("AliITSUHit");
 131   //
 132   TObjArray arrMCTracks; // array of hit arrays for each particle
 133   //
 134   Float_t xyzClGloF[3];
 135   Double_t xyzClGlo[3],xyzClTr[3];
 136   Int_t labels[3];
 137   int nLab = 0;
 138   int nlr=its->GetNLayersActive();
 139   int ntotev = (Int_t)runLoader->GetNumberOfEvents();
 140   printf("N Events : %i \n",ntotev);
 141   if (nev>0) ntotev = TMath::Min(nev,ntotev);
 142   //
 143   // load topology database
 144   LoadDB("clusterTopology.root");
 145
 146   Int_t nPatterns = pattDB->GetEntriesFast();
 147
 148   // output tree
 149   TFile* flOut = TFile::Open("clInfo.root","recreate");
 150   TTree* trOut = new TTree("clitsu","clitsu");
 151   clSumm cSum;
 152   trOut->Branch("evID", &cSum.evID ,"evID/I");
 153   trOut->Branch("volID",&cSum.volID,"volID/I");
 154   trOut->Branch("lrID", &cSum.lrID ,"lrID/I");
 155   trOut->Branch("clID", &cSum.clID ,"clID/I");
 156   trOut->Branch("nPix", &cSum.nPix ,"nPix/I");
 157   trOut->Branch("nX"  , &cSum.nX   ,"nX/I");
 158   trOut->Branch("nZ"  , &cSum.nZ   ,"nZ/I");
 159   trOut->Branch("q"   , &cSum.q    ,"q/I");
 160   trOut->Branch("pt"  , &cSum.pt   ,"pt/F");
 161   trOut->Branch("eta"  ,&cSum.eta  ,"eta/F");
 162   trOut->Branch("phi"  , &cSum.phi  ,"phi/F");
 163   trOut->Branch("xyz",   cSum.xyz,  "xyz[3]/F");
 164   trOut->Branch("dX"  , &cSum.dX   ,"dX/F");
 165   trOut->Branch("dY"  , &cSum.dY   ,"dY/F");
 166   trOut->Branch("dZ"  , &cSum.dZ   ,"dZ/F");
 167   trOut->Branch("split",&cSum.split,"split/O");
 168   trOut->Branch("prim", &cSum.prim, "prim/O");
 169   trOut->Branch("pdg",  &cSum.pdg,  "pdg/I");
 170   trOut->Branch("ntr",  &cSum.ntr,  "ntr/I");
 171   trOut->Branch("alpha", &cSum.alpha, "alpha/F");
 172   trOut->Branch("beta", &cSum.beta, "beta/F");
 173   trOut->Branch("nRowPatt", &cSum.nRowPatt, "nRowPatt/I");
 174   trOut->Branch("nColPatt", &cSum.nColPatt, "nColPatt/I");
 175   trOut->Branch("pattID", &cSum.pattID,"pattID/I");
 176   trOut->Branch("freq", &cSum.freq, "freq/F");
 177   trOut->Branch("xCen", &cSum.xCen, "xCen/F");
 178   trOut->Branch("zCen", &cSum.zCen, "zCen/F");
 179   trOut->Branch("xShift", &cSum.xShift, "xShift/F");
 180   trOut->Branch("zShift", &cSum.zShift, "zShift/F");
 181   //
 182   for (Int_t iEvent = 0; iEvent < ntotev; iEvent++) {
 183     printf("\n Event %i \n",iEvent);
 184     runLoader->GetEvent(iEvent);
 185     AliStack *stack = runLoader->Stack();
 186     cluTree=dl->TreeR();
 187     hitTree=dl->TreeH();
 188     hitTree->SetBranchAddress("ITS",&hitList);
 189     //
 190     // read clusters
 191     for (int ilr=nlr;ilr--;) {
 192       TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",ilr));
 193       if (!br) {printf("Did not find cluster branch for lr %d\n",ilr); exit(1);}
 194       br->SetAddress(its->GetLayerActive(ilr)->GetClustersAddress());
 195     }
 196     cluTree->GetEntry(0);
 197     its->ProcessClusters();
 198     //
 199     // read hits
 200     for(Int_t iEnt=0;iEnt<hitTree->GetEntries();iEnt++){//entries loop degli hits
 201       hitTree->GetEntry(iEnt);
 202       int nh = hitList->GetEntries();
 203       for(Int_t iHit=0; iHit<nh;iHit++){
 204              AliITSUHit *pHit = (AliITSUHit*)hitList->At(iHit);
 205              int mcID = pHit->GetTrack();
 206              TClonesArray* harr = arrMCTracks.GetEntriesFast()>mcID ? (TClonesArray*)arrMCTracks.At(mcID) : 0;
 207              if (!harr) {
 208                harr = new TClonesArray("AliITSUHit"); // 1st encounter of the MC track
 209                arrMCTracks.AddAtAndExpand(harr,mcID);
 210              }
 211                //
 212              new ( (*harr)[harr->GetEntriesFast()] ) AliITSUHit(*pHit);
 213       }
 214     }
 215     //
 216     // compare clusters and hits
 217     //
 218     printf(" tree entries: %lld\n",cluTree->GetEntries());
 219     //
 220     for (int ilr=0;ilr<nlr;ilr++) {
 221       AliITSURecoLayer* lr = its->GetLayerActive(ilr);
 222       TClonesArray* clr = lr->GetClusters();
 223       int nClu = clr->GetEntries();
 224       printf("Layer %d : %d clusters\n",ilr,nClu);
 225       //
 226       for (int icl=0;icl<nClu;icl++) {
 227               AliITSUClusterPix *cl = (AliITSUClusterPix*)clr->At(icl);
 228               int modID = cl->GetVolumeId();
 229
 230         //Compare clusters with DB
 231
 232         clTop.Clear();
 233         Int_t tempPattID=-1;
 234         Float_t tempFreq=0;
 235         Float_t tempxCen=0;
 236         Float_t tempzCen=0;
 237         Float_t tempxShift=0;
 238         Float_t tempzShift=0;
 239         UShort_t rs = cl->GetPatternRowSpan();
 240         UShort_t cs = cl->GetPatternColSpan();
 241
 242         for(int ir=0; ir<rs; ir++)
 243           for(int ic=0; ic<cs; ic++)
 244             if(cl->TestPixel(ir,ic)) clTop.SetBitNumber(ir*cs+ic);//Set cluster pattern
 245         clTop.SetUniqueID((rs<<16)+cs);//Set cluster ID
 246
 247         for(int iDB=0; iDB<nPatterns; iDB++){
 248           TBits* refPatt = (TBits*) pattDB->At(iDB);
 249           if (*refPatt == clTop && refPatt->GetUniqueID()== clTop.GetUniqueID()){
 250             tempPattID=iDB;
 251             tempFreq=(*pattFR)[iDB];
 252             tempxCen=(*xCentrePix)[iDB];
 253             tempzCen=(*zCentrePix)[iDB];
 254             tempxShift=(*xCentreShift)[iDB];
 255             tempzShift=(*zCentreShift)[iDB];
 256             break;
 257           }
 258         }
 259
 260              //------------ check if this is a split cluster
 261              if (!cl->TestBit(kSplCheck)) {
 262                cl->SetBit(kSplCheck);
 263                // check if there is no other cluster with same label on this module
 264                AliITSURecoSens* sens = lr->GetSensorFromID(modID);
 265                int nclSn = sens->GetNClusters();
 266                int offs = sens->GetFirstClusterId();
 267                //       printf("To check for %d (mod:%d) N=%d from %d\n",icl,modID,nclSn,offs);
 268                for (int ics=0;ics<nclSn;ics++) {
 269                  AliITSUClusterPix* clusT = (AliITSUClusterPix*)lr->GetCluster(offs+ics); // access to clusters
 270                  if (clusT==cl) continue;
 271                  for (int ilb0=0;ilb0<3;ilb0++) {
 272                    int lb0 = cl->GetLabel(ilb0); if (lb0<=-1) break;
 273                    for (int ilb1=0;ilb1<3;ilb1++) {
 274                             int lb1 = clusT->GetLabel(ilb1); if (lb1<=-1) break;
 275                             if (lb1==lb0) {
 276                               cl->SetBit(kSplit);
 277                               clusT->SetBit(kSplit);
 278                               /*
 279                               printf("Discard clusters of module %d:\n",modID);
 280                               cl->Print();
 281                               clusT->Print();
 282                               */
 283                               break;
 284                             }
 285                     }
 286                   }
 287                 }
 288               }
 289              //------------
 290              const AliITSsegmentation* segm = gm->GetSegmentation(ilr);
 291              //
 292              cl->GetGlobalXYZ(xyzClGloF);
 293              int clsize = cl->GetNPix();
 294              for (int i=3;i--;) xyzClGlo[i] = xyzClGloF[i];
 295              const TGeoHMatrix* mat = gm->GetMatrixSens(modID);
 296              if (!mat) {printf("failed to get matrix for module %d\n",cl->GetVolumeId());}
 297              mat->MasterToLocal(xyzClGlo,xyzClTr);
 298              //
 299              int col,row;
 300              segm->LocalToDet(xyzClTr[0],xyzClTr[2],row,col); // effective col/row
 301              nLab = 0;
 302              for (int il=0;il<3;il++) {
 303                if (cl->GetLabel(il)>=0) labels[nLab++] = cl->GetLabel(il);
 304                else break;
 305               }
 306               // find hit info
 307          for (int il=0;il<nLab;il++) {
 308            TClonesArray* htArr = (TClonesArray*)arrMCTracks.At(labels[il]);
 309            if (!htArr) {printf("did not find MChits for label %d ",labels[il]); cl->Print(); continue;}
 310            //
 311            int nh = htArr->GetEntriesFast();
 312            AliITSUHit *pHit=0;
 313            double dst2Max = 1e33;
 314            for (int ih=nh;ih--;) {
 315              AliITSUHit* tHit = (AliITSUHit*)htArr->At(ih);
 316              if (tHit->GetChip()!=modID) continue;
 317              tHit->GetPositionG(xg1,yg1,zg1);
 318              tHit->GetPositionG0(xg0,yg0,zg0,tg0);
 319              double gxyzHDif[3] = { (xg1+xg0)/2 - xyzClGlo[0], (yg1+yg0)/2 - xyzClGlo[1], (zg1+zg0)/2 - xyzClGlo[2] };
 320              double dst2 = gxyzHDif[0]*gxyzHDif[0] + gxyzHDif[1]*gxyzHDif[1] + gxyzHDif[2]*gxyzHDif[2];
 321              if (dst2<dst2Max) {
 322                pHit = tHit;
 323                dst2Max = dst2;
 324               }
 325             }
 326                 if (!pHit) {
 327                  printf("did not find MChit for label %d on module %d ",il,modID);
 328                  cl->Print();
 329                  htArr->Print();
 330                  continue;
 331                 }
 332                 //
 333                 pHit->GetPositionG(xg1,yg1,zg1);
 334                 pHit->GetPositionG0(xg0,yg0,zg0,tg0);
 335                 //
 336                 double txyzH[3],gxyzH[3] = { (xg1+xg0)/2, (yg1+yg0)/2, (zg1+zg0)/2 };
 337                 mat->MasterToLocal(gxyzH,txyzH);
 338                 double rcl = TMath::Sqrt(xyzClTr[0]*xyzClTr[0]+xyzClTr[1]*xyzClTr[1]);
 339                 double rht = TMath::Sqrt(txyzH[0]*txyzH[0]+txyzH[1]*txyzH[1]);
 340                 //
 341           //
 342           //Angles determination
 343
 344           pHit->GetPositionL(xExit,yExit,zExit);
 345           pHit->GetPositionL0(xEnt,yEnt,zEnt,tof1);
 346
 347           Double_t dirHit[3]={(xExit-xEnt),(yExit-yEnt),(zExit-zEnt)};
 348
 349           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
 350           Float_t alpha2 = (Float_t) alpha1; //convert to float
 351           cSum.alpha = alpha2;
 352
 353           Double_t beta1;
 354           beta1 = TMath::ATan2(TMath::Abs(dirHit[0]),TMath::Abs(dirHit[2])); //Azimuthal angle, values from -Pi to Pi
 355           Float_t beta2 = (Float_t) beta1;
 356           cSum.beta = beta2;
 357
 358                 GetHistoClSize(clsize,kDR,&histoArr)->Fill((rht-rcl)*1e4);
 359                 if (cl->TestBit(kSplit)) {
 360                   if (col%2) GetHistoClSize(clsize,kDTXoddSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
 361                   else       GetHistoClSize(clsize,kDTXevenSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
 362                   GetHistoClSize(clsize,kDTZSPL,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
 363                   GetHistoClSize(0,kNPixSPL,&histoArr)->Fill(clsize);
 364                 }
 365                 if (col%2) GetHistoClSize(clsize,kDTXodd,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
 366                 else       GetHistoClSize(clsize,kDTXeven,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
 367                 GetHistoClSize(clsize,kDTZ,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
 368                 GetHistoClSize(0,kNPixAll,&histoArr)->Fill(clsize);
 369                 //
 370                 cSum.evID = iEvent;
 371                 cSum.volID = cl->GetVolumeId();
 372           cSum.lrID = ilr;
 373           cSum.clID = icl;
 374           cSum.nPix = cl->GetNPix();
 375           cSum.nX   = cl->GetNx();
 376           cSum.nZ   = cl->GetNz();
 377           cSum.q    = cl->GetQ();
 378           cSum.split = cl->TestBit(kSplit);
 379           cSum.dX = (txyzH[0]-xyzClTr[0])*1e4;
 380           cSum.dY = (txyzH[1]-xyzClTr[1])*1e4;
 381           cSum.dZ = (txyzH[2]-xyzClTr[2])*1e4;
 382           cSum.nRowPatt = cl-> GetPatternRowSpan();
 383           cSum.nColPatt = cl-> GetPatternColSpan();
 384           cSum.pattID = tempPattID;
 385           cSum.freq = tempFreq;
 386           cSum.xCen = tempxCen;
 387           cSum.zCen = tempzCen;
 388           cSum.xShift = tempxShift;
 389           cSum.zShift = tempzShift;
 390
 391           int label = cl->GetLabel(0);
 392           TParticle* part = 0;
 393           if (label>=0 && (part=stack->Particle(label)) ) {
 394             cSum.pdg = part->GetPdgCode();
 395             cSum.eta = part->Eta();
 396             cSum.pt  = part->Pt();
 397             cSum.phi = part->Phi();
 398                   cSum.prim = stack->IsPhysicalPrimary(label);
 399                 }
 400                 cSum.ntr = 0;
 401                 for (int ilb=0;ilb<3;ilb++) if (cl->GetLabel(ilb)>=0) cSum.ntr++;
 402                 for (int i=0;i<3;i++) cSum.xyz[i] = xyzClGloF[i];
 403                 //
 404                 trOut->Fill();
 405           /*
 406           if (clsize==5) {
 407             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',
 408                    modID,icl,(txyzH[0]-xyzClTr[0])*1e4,(txyzH[2]-xyzClTr[2])*1e4, row,col,
 409                    gxyzH[0],gxyzH[1],gxyzH[2],xyzClGlo[0],xyzClGlo[1],xyzClGlo[2]);
 410             cl->Print();
 411             pHit->Print();
 412             //
 413             double a0,b0,c0,a1,b1,c1,e0;
 414             pHit->GetPositionL0(a0,b0,c0,e0);
 415             pHit->GetPositionL(a1,b1,c1);
 416             float cloc[3];
 417             cl->GetLocalXYZ(cloc);
 418             printf("LocH: %e %e %e | %e %e %e\n",a0,b0,c0,a1,b1,c1);
 419             printf("LocC: %e %e %e | %e %e %e\n",cloc[0],cloc[1],cloc[2],xyzClTr[0],xyzClTr[1],xyzClTr[2]);
 420           }
 421           */
 422           //
 423         }
 424       }
 425     }
 426
 427     //    layerClus.Clear();
 428     //
 429     arrMCTracks.Delete();
 430   }//event loop
 431   //
 432   flOut->cd();
 433   trOut->Write();
 434   delete trOut;
 435   flOut->Close();
 436   flOut->Delete();
 437   DrawReport("clinfo.ps",&histoArr);
 438   //
 439 }
 440
 441 void DrawReport(const char* psname, TObjArray* harr)
 442 {
 443   gStyle->SetOptFit(1);
 444   if (!harr) harr = &histoArr;
 445   TCanvas* cnv = new TCanvas("cl","cl",900,600);
 446   //
 447   TString psnm1 = psname;
 448   if (psnm1.IsNull()) psnm1 = "clusters.ps";
 449   TString psnm0 = psnm1.Data();
 450   psnm0 += "[";
 451   TString psnm2 = psnm1.Data();
 452   psnm2 += "]";
 453   cnv->Print(psnm0.Data());
 454   //
 455   TH1* clall = GetHistoClSize(0,kNPixAll,harr);
 456   clall->SetLineColor(kRed);
 457   clall->Draw();
 458   TH1* clSpl = GetHistoClSize(0,kNPixSPL,harr);
 459   clSpl->SetLineColor(kBlue);
 460   clSpl->Draw("sames");
 461   gPad->Modified();
 462   gPad->Update();
 463   SetStPadPos(clall,0.75,0.97,0.8,1.,-1,clall->GetLineColor());
 464   SetStPadPos(clSpl,0.75,0.97,0.6,0.8,-1,clSpl->GetLineColor());
 465   gPad->Modified();
 466   gPad->Update();
 467   gPad->SetLogy(1);
 468   //
 469   cnv->cd();
 470   cnv->Print(psnm1.Data());
 471   //
 472   // plot cluster sized from 1 to 10
 473   for (int i=1;i<=10;i++) {
 474     if (clall->GetBinContent(clall->FindBin(i))<100) continue;
 475     DrawNP(i,harr,cnv);
 476     cnv->Print(psnm1.Data());
 477   }
 478   cnv->Print(psnm2.Data());
 479 }
 480
 481 //------------------------------------------
 482 TH1* GetHistoClSize(int npix,int id,TObjArray* harr)
 483 {
 484   // book histos
 485   TH1* h = 0;
 486   if (!harr) harr = &histoArr;
 487   //
 488   if (npix<1) {
 489     if (harr->GetEntriesFast()>=id && (h=(TH1*)harr->At(id))) return h;
 490     h = new TH1F("npixAll","npixAll",150,0.5,54.5);
 491     h->SetDirectory(0);
 492     h->SetLineColor(kRed);
 493     harr->AddAtAndExpand(h, kNPixAll);
 494     //
 495     h = new TH1F("npixSpl","npixSpl",150,0.5,54.5);
 496     h->SetLineColor(kBlue);
 497     h->SetDirectory(0);
 498     harr->AddAtAndExpand(h, kNPixSPL);
 499     //
 500     h = (TH1*)harr->At(id);
 501     if (!h) {printf("Unknown histo id=%d\n",id); exit(1);}
 502     return h;
 503   }
 504   //
 505   int idh = npix*10+id;
 506   if (harr->GetEntriesFast()>=idh && (h=(TH1*)harr->At(idh))) return h;
 507   //
 508   const int nbin=100;
 509   const double kdiff=80;
 510   // need to create set of histos
 511   //
 512   h = new TH1F(Form("dxy_npix%d",npix),Form("dr_npix%d",npix),nbin,-kdiff,kdiff);
 513   h->SetDirectory(0);
 514   harr->AddAtAndExpand(h, npix*10 + kDR);
 515   //
 516   h  = new TH1F(Form("dtxODD_npix%d",npix),Form("dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
 517   h->SetDirectory(0);
 518   h->SetLineColor(kRed);
 519   harr->AddAtAndExpand(h, npix*10 + kDTXodd);
 520   h  = new TH1F(Form("dtxEVN_npix%d",npix),Form("dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
 521   h->SetDirectory(0);
 522   h->SetLineColor(kBlue);
 523   harr->AddAtAndExpand(h, npix*10 + kDTXeven);
 524   //
 525   h  = new TH1F(Form("dtz_npix%d",npix),Form("dtz_npix%d",npix),nbin,-kdiff,kdiff);
 526   h->SetLineColor(kGreen);
 527   h->SetDirectory(0);
 528   harr->AddAtAndExpand(h, npix*10 + kDTZ);
 529   //
 530   //
 531   h  = new TH1F(Form("SPL_dtxODD_npix%d",npix),Form("SPL_dtxODD_npix%d",npix),nbin,-kdiff,kdiff);
 532   h->SetLineColor(kMagenta);
 533   h->SetFillColor(kMagenta);
 534   h->SetFillStyle(3001);
 535   h->SetLineStyle(2);
 536   h->SetDirectory(0);
 537
 538   harr->AddAtAndExpand(h, npix*10 + kDTXoddSPL);
 539   h  = new TH1F(Form("SPL_dtxEVN_npix%d",npix),Form("SPL_dtxEVN_npix%d",npix),nbin,-kdiff,kdiff);
 540   h->SetLineColor(kCyan);
 541   h->SetFillColor(kCyan);
 542   h->SetFillStyle(3006);
 543   h->SetLineStyle(2);
 544   h->SetDirectory(0);
 545   harr->AddAtAndExpand(h, npix*10 + kDTXevenSPL);
 546   //
 547   h  = new TH1F(Form("SPL_dtz_npix%d",npix),Form("SPLdtz_npix%d",npix),nbin,-kdiff,kdiff);
 548   harr->AddAtAndExpand(h, npix*10 + kDTZSPL);
 549   h->SetDirectory(0);
 550   //
 551   h->SetLineColor(kGreen+2);
 552   h->SetFillColor(kGreen+2);
 553   h->SetLineStyle(2);
 554   h->SetFillStyle(3001);
 555   h = (TH1*)harr->At(idh);
 556   if (!h) {printf("Unknown histo id=%d\n",idh); exit(1);}
 557   return h;
 558 }
 559
 560
 561
 562 TCanvas* DrawNP(int np, TObjArray* harr, TCanvas* cnv)
 563 {
 564   if (!harr) harr = &histoArr;
 565   if (!cnv) cnv = new TCanvas(Form("cnv%d",np),Form("cnv%d",np),900,700);
 566   cnv->Clear();
 567   cnv->Divide(2,1);
 568   cnv->cd(1);
 569   //
 570   TH1* dxodd = (TH1*)harr->At(np*10+kDTXodd);
 571   TH1* dxevn = (TH1*)harr->At(np*10+kDTXeven);
 572   TH1* dxoddS =(TH1*)harr->At(np*10+kDTXoddSPL);
 573   TH1* dxevnS =(TH1*)harr->At(np*10+kDTXevenSPL);
 574   double max = TMath::Max(dxodd->GetMaximum(),dxevn->GetMaximum());
 575   dxodd->SetMaximum(1.1*max);
 576   dxodd->GetXaxis()->SetTitle("#DeltaX, #mum");
 577   dxodd->SetTitle(Form("#DeltaX for clSize=%d",np));
 578   dxodd->Fit("gaus","","");
 579   dxevn->Fit("gaus","","sames");
 580   //
 581   dxoddS->Draw("sames");
 582   dxevnS->Draw("sames");
 583   //
 584   gPad->Modified();
 585   gPad->Update();
 586   SetStPadPos(dxodd,0.75,0.97,0.8,1., -1,dxodd->GetLineColor());
 587   SetStPadPos(dxevn,0.75,0.97,0.6,0.8, -1,dxevn->GetLineColor());
 588   SetStPadPos(dxoddS,0.75,0.97,0.4,0.6, -1,dxoddS->GetLineColor());
 589   SetStPadPos(dxevnS,0.75,0.97,0.2,0.4, -1,dxevnS->GetLineColor());
 590   //
 591   cnv->cd(2);
 592   TH1* dz  = (TH1*)harr->At(np*10+kDTZ);
 593   dz->SetTitle(Form("#DeltaZ for clSize=%d",np));
 594   dz->GetXaxis()->SetTitle("#DeltaZ, #mum");
 595   dz->Fit("gaus");
 596   TH1* dzS = (TH1*)harr->At(np*10+kDTZSPL);
 597   dz->Draw("sames");
 598   gPad->Modified();
 599   gPad->Update();
 600   SetStPadPos(dz,0.75,0.97,0.8,1., -1, dz->GetLineColor());
 601   SetStPadPos(dzS,0.75,0.97,0.5,0.7, -1, dzS->GetLineColor());
 602   gPad->Modified();
 603   gPad->Update();
 604   //
 605   cnv->cd();
 606   return cnv;
 607 }
 608
 609
 610 TPaveStats* GetStPad(TH1* hst)
 611 {
 612   TList *lst = hst->GetListOfFunctions();
 613   if (!lst) return 0;
 614   int nf = lst->GetSize();
 615   for (int i=0;i<nf;i++) {
 616     TPaveStats *fnc = (TPaveStats*) lst->At(i);
 617     if (fnc->InheritsFrom("TPaveStats")) return fnc;
 618   }
 619   return 0;
 620   //
 621 }
 622
 623
 624 TPaveStats* SetStPadPos(TH1* hst,float x1,float x2,float y1,float y2, Int_t stl, Int_t col)
 625 {
 626   TPaveStats* pad = GetStPad(hst);
 627   if (!pad) return 0;
 628   pad->SetX1NDC( x1 );
 629   pad->SetX2NDC( x2 );
 630   pad->SetY1NDC( y1 );
 631   pad->SetY2NDC( y2 );
 632   if (stl>=0) pad->SetFillStyle(stl);
 633   if (col>=0) pad->SetTextColor(col);
 634   pad->SetFillColor(0);
 635   //
 636   gPad->Modified();
 637   return pad;
 638 }
 639
 640 void LoadDB(const char* fname)
 641 {
 642   // load database
 643   TFile* fl = TFile::Open(fname);
 644   pattDB = (TObjArray*)fl->Get("TopDB");
 645   pattFR = (TVectorF*)fl->Get("TopFreq");
 646   xCentrePix =(TVectorF*)fl->Get("xCOG");
 647   zCentrePix =(TVectorF*)fl->Get("zCOG");
 648   xCentreShift =(TVectorF*)fl->Get("xShift");
 649   zCentreShift =(TVectorF*)fl->Get("zShift");
 650
 651 }
 652