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