TRD/AliTRDclusterErrors.C

   1 //#ifndef __CINT__
   2   #include <iostream.h>
   3
   4   #include "AliTRDtracker.h"
   5   #include "AliTRDclusterMI.h"
   6   #include "AliTRDhit.h"
   7   #include "AliTRDv1.h"
   8   #include "AliTRDgeometry.h"
   9   #include "AliTRDgeometryDetail.h"
  10   #include "AliTRDparameter.h"
  11
  12   #include "alles.h"
  13   #include "TFile.h"
  14   #include "TStopwatch.h"
  15   #include "Rtypes.h"
  16   #include "TTree.h"
  17 #include "AliRunLoader.h"
  18 #include "AliStack.h"
  19 #include "AliTRDclusterErrors.h"
  20 #include "TF1.h"
  21 #include "AliTrackReference.h"
  22 //#endif
  23
  24
  25 ClassImp(AliTRDCI)
  26 ClassImp(AliTRDExactPoint)
  27 ClassImp(AliTRDClusterErrAnal)
  28 ClassImp(AliTRDClusterErrDraw)
  29
  30
  31
  32 AliTRDClusterErrAnal::AliTRDClusterErrAnal(Char_t *chloader )
  33 {
  34   //
  35   //SET Input loaders
  36   if (gAlice){
  37     delete gAlice->GetRunLoader();
  38     delete gAlice;
  39     gAlice = 0x0;
  40   }
  41   fRunLoader = AliRunLoader::Open(chloader);
  42   if (fRunLoader == 0x0){
  43     cerr<<"Can not open session"<<endl;
  44     return;
  45   }
  46   fTRDLoader = fRunLoader->GetLoader("TRDLoader");
  47   if (fTRDLoader == 0x0){
  48     cerr<<"Can not get TRD Loader"<<endl;
  49     return ;
  50   }
  51   if (fRunLoader->LoadgAlice()){
  52     cerr<<"Error occured while l"<<endl;
  53     return;
  54   }
  55   fRunLoader->CdGAFile();
  56   fTracker = new AliTRDtracker(gFile);
  57   fParam = (AliTRDparameter*) gFile->Get("TRDparameter");
  58   fGeometry = new AliTRDgeometryDetail();
  59   fTRD      = (AliTRDv1*) gAlice->GetDetector("TRD");
  60   //
  61   AliTRDCI * clinfo = new AliTRDCI();
  62   fFileA  = new TFile("trdclusteranal.root","recreate");
  63   fFileA->cd();
  64   fTreeA  = new TTree("trdcl","trdcl");
  65   fTreeA->Branch("trdcl","AliTRDCI",&clinfo);
  66   delete clinfo;
  67 }
  68
  69 void AliTRDClusterErrAnal::SetIO(Int_t event)
  70 {
  71   //
  72   //set input output for given event
  73   fRunLoader->SetEventNumber(event);
  74   fRunLoader->LoadHeader();
  75   fRunLoader->LoadKinematics();
  76   fRunLoader->LoadTrackRefs();
  77   fTRDLoader->LoadHits();
  78   fTRDLoader->LoadRecPoints("read");
  79   //
  80   fStack = fRunLoader->Stack();
  81   fHitTree = fTRDLoader->TreeH();
  82   fClusterTree = fTRDLoader->TreeR();
  83   fReferenceTree = fRunLoader->TreeTR();
  84   //
  85 }
  86
  87 void AliTRDClusterErrAnal::LoadClusters()
  88 {
  89   //
  90   //
  91   // Load clusters
  92   TObjArray *ClusterArray = new TObjArray(400);
  93   TObjArray carray(2000);
  94   TBranch *branch=fClusterTree->GetBranch("TRDcluster");
  95   if (!branch) {
  96     Error("ReadClusters","Can't get the branch !");
  97     return;
  98   }
  99   branch->SetAddress(&ClusterArray);
 100   Int_t nentries = (Int_t)fClusterTree->GetEntries();
 101   for (Int_t i=0;i<nentries;i++){
 102     fClusterTree->GetEvent(i);
 103     Int_t nCluster = ClusterArray->GetEntriesFast();
 104     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
 105       AliTRDcluster* c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
 106       carray.AddLast(c);
 107       ClusterArray->RemoveAt(iCluster);
 108     }
 109   }
 110   Int_t nClusters = carray.GetEntriesFast();
 111   printf("Total number of clusters %d \n", nClusters);
 112   //
 113   //
 114   //SORT clusters
 115   //
 116   Int_t all=0;
 117   for (Int_t i=0;i<nClusters;i++){
 118     AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i);
 119     Int_t plane = fGeometry->GetPlane(cl->GetDetector());
 120     if (plane>=12) continue;
 121     Int_t time = cl->GetLocalTimeBin();
 122     if (time>=100) continue;
 123     Int_t sector = fGeometry->GetSector(cl->GetDetector());
 124     if (sector>=18){
 125       printf("problem1\n");
 126       continue;
 127     }
 128     fClusters[plane][time][sector].AddLast(cl);
 129     all++;
 130   }
 131 }
 132
 133 void AliTRDClusterErrAnal::SortReferences()
 134 {
 135   //
 136   //
 137   //
 138   printf("Sorting references\n");
 139   fReferences = new TObjArray;
 140   Int_t ntracks = fStack->GetNtrack();
 141   fReferences->Expand(ntracks);
 142   Int_t nentries = (Int_t)fReferenceTree->GetEntries();
 143   TClonesArray * arr = new TClonesArray("AliTrackReference");
 144   TBranch * br = fReferenceTree->GetBranch("TRD");
 145   br->SetAddress(&arr);
 146   //
 147   TClonesArray *labarr=0;
 148   for (Int_t iprim=0;iprim<nentries;iprim++){
 149     if (br->GetEntry(iprim)){
 150       for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){
 151         AliTrackReference *ref =(AliTrackReference*)arr->At(iref);
 152         if (!ref) continue;
 153         Int_t lab = ref->GetTrack();
 154         if ( (lab<0) || (lab>ntracks)) continue;
 155         //
 156         if (fReferences->At(lab)==0) {
 157           labarr = new TClonesArray("AliTrackReference");
 158           fReferences->AddAt(labarr,lab);
 159         }
 160         TClonesArray &larr = *labarr;
 161         new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
 162       }
 163     }
 164   }
 165   printf("end - Sorting references\n");
 166 }
 167
 168 AliTrackReference * AliTRDClusterErrAnal::FindNearestReference(Int_t lab, Float_t pos[3], Float_t dmax)
 169 {
 170   //
 171   //
 172   //
 173   if (fReferences->At(lab)==0) return 0;
 174   AliTrackReference *nearest=0;
 175   TClonesArray * arr = (TClonesArray *)fReferences->At(lab);
 176   for (Int_t iref =0;iref<arr->GetEntriesFast();iref++){
 177     AliTrackReference * ref = ( AliTrackReference *)arr->UncheckedAt(iref);
 178     if (!ref) continue;
 179     Float_t delta = (pos[0]-ref->X())*(pos[0]-ref->X());
 180     delta += (pos[1]-ref->Y())*(pos[1]-ref->Y());
 181     delta += (pos[2]-ref->Z())*(pos[2]-ref->Z());
 182     delta = TMath::Sqrt(delta);
 183     if (delta<dmax){
 184       dmax=delta;
 185       nearest = ref;
 186     }
 187   }
 188   return nearest;
 189 }
 190
 191 void AliTRDClusterErrAnal::MakeExactPoints(Int_t trackmax)
 192 {
 193   //
 194   //make exact points:-)
 195   //
 196   //
 197
 198   fExactPoints.Delete();
 199   fExactPoints.Expand(fStack->GetNtrack());
 200   //
 201   Double_t fSum=0;
 202   Double_t fSumQ =0;
 203   Double_t fSumX=0;
 204   Double_t fSumX2=0;
 205   Double_t fSumXY=0;
 206   Double_t fSumXZ=0;
 207   Double_t fSumY=0;
 208   Double_t fSumZ=0;
 209   //
 210   Int_t entries = Int_t(fHitTree->GetEntries());
 211   printf("Number of primary  entries\t%d\n",entries);
 212   entries = TMath::Min(trackmax,entries);
 213   Int_t nallpoints = 0;
 214
 215   Int_t nalltracks =0;
 216   Int_t pointspertrack =0;
 217
 218   for (Int_t entry=0;entry<entries; entry++){
 219     gAlice->ResetHits();
 220     fHitTree->GetEvent(entry);
 221     Int_t lastlabel    = -1;
 222     Int_t lastdetector = -1;
 223     Int_t lasttimebin   = -1;
 224     Float_t lastpos[3];
 225     //
 226     for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1); hit;
 227         hit = (AliTRDhit *) fTRD->NextHit()) {
 228       //
 229       Int_t label    = hit->Track();
 230       TParticle * particle = fStack->Particle(label);
 231       if (!particle) continue;
 232       if (particle->Pt()<0.05) continue;
 233       Int_t detector = hit->GetDetector();
 234       Int_t plane    = fGeometry->GetPlane(detector);
 235       //
 236       //
 237       if (hit->GetCharge()==0) continue;
 238       Float_t pos[3] = {hit->X(),hit->Y(),hit->Z()};
 239       Int_t indexes[3];
 240       fGeometry->Global2Detector(pos,indexes,fParam);
 241       //
 242       Float_t rot[3];
 243       fGeometry->Rotate(detector,pos,rot);
 244       //
 245       //
 246       Float_t  time0    = fParam->GetTime0(plane);
 247       Int_t    timebin  = Int_t(TMath::Nint(((time0 - rot[0])/fParam->GetTimeBinSize())+ fParam->GetTimeBefore())+0.1);
 248       if (timebin<0) continue;
 249       //
 250       //
 251       if (label!=lastlabel || detector != lastdetector || lasttimebin !=timebin){
 252         //
 253         if (label!=lastlabel){
 254           fExactPoints.AddAt(new TClonesArray("AliTRDExactPoint",0),label);
 255           //printf("new particle\t%d\n",label);
 256           nalltracks++;
 257           //      printf("particle\t%d- hits\t%d\n",lastlabel, pointspertrack);
 258           pointspertrack=0;
 259         }
 260
 261         if ( (fSum>1) && lasttimebin>=0 && lasttimebin<fParam->GetTimeMax() ){
 262           //if we have enough info for given layer time bin - store it
 263           AliTrackReference * ref = FindNearestReference(lastlabel,lastpos,4.);
 264           Float_t rotmom[3];
 265           Float_t rotpos[3];
 266           Float_t refangle=0;
 267           if (ref){
 268             Float_t mom[3] = {ref->Px(),ref->Py(),ref->Pz()};
 269             Float_t refpos[3] = {ref->X(),ref->Y(),ref->Z()};
 270             fGeometry->Rotate(detector,mom,rotmom);
 271             fGeometry->Rotate(detector,refpos,rotpos);
 272             refangle = rotmom[1]/rotmom[0];
 273
 274           }
 275
 276           Double_t ay,by,az,bz;
 277           Double_t det = fSum*fSumX2-fSumX*fSumX;
 278           if (TMath::Abs(det)> 0.000000000000001) {
 279             by = (fSum*fSumXY-fSumX*fSumY)/det;
 280             ay = (fSumX2*fSumY-fSumX*fSumXY)/det;
 281
 282           }else{
 283             ay =fSumXY/fSumX;
 284             by =0;
 285           }
 286           if (TMath::Abs(det)> 0.000000000000001) {
 287             bz = (fSum*fSumXZ-fSumX*fSumZ)/det;
 288             az = (fSumX2*fSumZ-fSumX*fSumXZ)/det;
 289           }else{
 290             az =fSumXZ/fSumX;
 291             bz =0;
 292           }
 293           //
 294           Float_t lastplane = fGeometry->GetPlane(lastdetector);
 295           Float_t time0    = fParam->GetTime0(lastplane);
 296           Float_t xcenter = time0 - (lasttimebin - fParam->GetTimeBefore())*fParam->GetTimeBinSize();
 297           Float_t ty = ay + by * xcenter;
 298           Float_t tz = az + bz * xcenter;
 299           //
 300
 301           TClonesArray * arr = (TClonesArray *) fExactPoints.At(label);
 302           TClonesArray & larr= *arr;
 303           Int_t arrent = arr->GetEntriesFast();
 304           AliTRDExactPoint * point = new (larr[arrent]) AliTRDExactPoint;
 305           nallpoints++;
 306
 307           if (ref){
 308             point->SetReference(ref);
 309             point->fTRefAngleY = rotmom[1]/rotmom[0];
 310           }
 311           point->fTX  = xcenter;
 312           point->fTY  = ty;
 313           point->fTZ  = tz;
 314           point->fTAY = by;
 315           point->fTAZ = bz;
 316           //
 317           point->fGx  = lastpos[0];
 318           point->fGy  = lastpos[1];
 319           point->fGz  = lastpos[2];
 320
 321           //
 322           point->fDetector     = lastdetector;
 323           point->fLocalTimeBin = lasttimebin;
 324           point->fPlane        = fGeometry->GetPlane(lastdetector);
 325           point->fSector      = fGeometry->GetSector(lastdetector);
 326           point->fPlaneMI      = indexes[0];
 327           //
 328           point->fTPrim = fSum;
 329           point->fTQ    = fSumQ;
 330           //
 331         }
 332         lastdetector = detector;
 333         lastlabel    = label;
 334         lasttimebin  = timebin;
 335         fSum=fSumQ=fSumX=fSumX2=fSumXY=fSumXZ=fSumY=fSumZ=0.;
 336       }
 337       //
 338       lastpos[0] = hit->X();
 339       lastpos[1] = hit->Y();
 340       lastpos[2] = hit->Z();
 341       fSum++;
 342       fSumQ = hit->GetCharge();
 343       fSumX  +=rot[0];
 344       fSumX2 +=rot[0]*rot[0];
 345       fSumXY +=rot[0]*rot[1];
 346       fSumXZ +=rot[0]*rot[2];
 347       fSumY  +=rot[1];
 348       fSumZ  +=rot[2];
 349       pointspertrack++;
 350     }
 351   }
 352   //
 353   printf("Found %d exact points\n",nallpoints);
 354 }
 355
 356
 357
 358
 359
 360
 361 Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) {
 362
 363   //
 364   // comparison works with both cluster types MI and old also
 365   //dummy cluster to be fill if not cluster info
 366   AliTRDclusterMI clmi;
 367   TClass * classmi =  clmi.IsA();
 368   //
 369   //SetOutput
 370   AliTRDCI * clinfo = new AliTRDCI();
 371   TBranch * clbr = fTreeA->GetBranch("trdcl");
 372   clbr->SetAddress(&clinfo);
 373
 374   SetIO(0);
 375   SortReferences();
 376   MakeExactPoints(trackmax);
 377   LoadClusters();
 378   //
 379   trackmax =  fStack->GetNtrack();
 380   //
 381   // Get the number of entries in the hit tree
 382   // (Number of primary particles creating a hit somewhere)
 383   Int_t nTrack = (Int_t)fExactPoints.GetEntries();
 384   printf("Found %d charged in TRD in first %d primaries", nTrack, trackmax);
 385   //
 386
 387   for (Int_t itrack = 0; itrack<trackmax; itrack++){
 388     TClonesArray *arrpoints = (TClonesArray*)fExactPoints.At(itrack);
 389
 390     if (!arrpoints) continue;
 391     //printf("new particle\t%d\n",itrack);
 392     TParticle * particle = fStack->Particle(itrack);
 393     if (!particle) continue;
 394     //printf("founded in kine tree \t%d\n",itrack);
 395     Int_t npoints = arrpoints->GetEntriesFast();
 396     if (npoints<10) continue;
 397     //printf("have enough points \t%d\t%d\n",itrack,npoints);
 398
 399     for (Int_t ipoint=0;ipoint<npoints;ipoint++){
 400       AliTRDExactPoint * point = (AliTRDExactPoint*)arrpoints->UncheckedAt(ipoint);
 401       if (!point) continue;
 402       //
 403       Int_t sec = fGeometry->GetSector(point->fDetector);
 404       if (sec>18){
 405         printf("problem2\n");
 406       }
 407       TObjArray & cllocal = fClusters[point->fPlane][point->fLocalTimeBin][sec];
 408       Int_t nclusters = cllocal.GetEntriesFast();
 409       Float_t maxdist = 10;
 410       AliTRDcluster * nearestcluster =0;
 411       //find nearest cluster to hit with given label
 412       for (Int_t icluster =0; icluster<nclusters; icluster++){
 413         AliTRDcluster * cluster = (AliTRDcluster*)cllocal.UncheckedAt(icluster);
 414         if (!cluster) continue;
 415         if ( (cluster->GetLabel(0)!= itrack) &&  (cluster->GetLabel(1)!= itrack)&&(cluster->GetLabel(2)!= itrack))
 416           continue;
 417         Float_t dist = TMath::Abs(cluster->GetY()-point->fTY);
 418         if (TMath::Abs(cluster->GetZ()-point->fTZ)>5.5) continue;
 419         if (dist<maxdist){
 420           maxdist = dist;
 421           nearestcluster = cluster;
 422         }
 423       }
 424       //
 425       clinfo->fEp  = *point;
 426       clinfo->fP   = *particle;
 427       if (!nearestcluster) {
 428         clinfo->fStatus=1;
 429         clinfo->fCl = clmi;
 430       }
 431       else{
 432         clinfo->fStatus=0;
 433         if (nearestcluster->IsA()==classmi){
 434           clinfo->fCl    =*((AliTRDclusterMI*)nearestcluster);
 435         }
 436         else{
 437           clinfo->fCl    = *nearestcluster;
 438         }
 439         //
 440         Float_t dz      = clinfo->fCl.GetZ()-point->fTZ;
 441         Double_t h01    = sin(TMath::Pi() / 180.0 * fParam->GetTiltingAngle());
 442         clinfo->fTDistZ = dz;
 443         clinfo->fDYtilt = h01*dz*((point->fPlane%2)*2.-1.);
 444         //
 445         clinfo->fNTracks =1;
 446         if (nearestcluster->GetLabel(1)>=0)  clinfo->fNTracks++;
 447         if (nearestcluster->GetLabel(2)>=0)  clinfo->fNTracks++;
 448         clinfo->Update();
 449       }
 450       //
 451       fTreeA->Fill();
 452     }
 453   }
 454
 455
 456   fFileA->cd();
 457   fTreeA->Write();
 458   fFileA->Close();
 459   return 0;
 460 }
 461
 462
 463 TH1F * AliTRDClusterErrDraw::ResDyVsAmp(TTree* tree, const char* selection, Float_t t0, Float_t ampmin, Float_t ampmax)
 464 {
 465   //
 466   //
 467   TH2F hisdy("resy","resy",10,ampmin,ampmax,30,-0.3,0.3);
 468   char expression[1000];
 469   sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fQ>>resy",t0);
 470   char selectionall[1000];
 471   sprintf(selectionall,"fStatus==0&&%s",selection);
 472   printf("%s\n",expression);
 473   printf("%s\n",selectionall);
 474   tree->Draw(expression,selectionall);
 475   return CreateResHisto(&hisdy);
 476 }
 477
 478
 479 TH1F * AliTRDClusterErrDraw::ResDyVsRelPos(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
 480 {
 481   //
 482   //
 483   min *=128;
 484   max *=128;
 485   TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
 486   char expression[1000];
 487   sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%.4f*fEp.fTAY:fCl.fRelPos>>resy",t0);
 488   char selectionall[1000];
 489   sprintf(selectionall,"fStatus==0&&%s",selection);
 490   printf("%s\n",expression);
 491   printf("%s\n",selectionall);
 492   tree->Draw(expression,selectionall);
 493   return CreateResHisto(&hisdy);
 494 }
 495
 496
 497 TH1F * AliTRDClusterErrDraw::ResDyVsAngleY(TTree* tree, const char* selection, Float_t t0, Float_t min, Float_t max)
 498 {
 499   //
 500   //
 501   TH2F hisdy("resy","resy",10,min,max,30,-0.3,0.3);
 502
 503   char expression[1000];
 504   sprintf(expression,"fEp.fTY-fCl.fY+fDYtilt+%f*fEp.fTAY:fEp.fTAY>>resy",t0);
 505   char selectionall[1000];
 506   sprintf(selectionall,"fStatus==0&&%s",selection);
 507
 508   tree->Draw(expression,selectionall);
 509   return CreateResHisto(&hisdy);
 510 }
 511
 512 void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
 513 {
 514   histo->GetXaxis()->SetTitle(xAxisTitle);
 515   histo->GetYaxis()->SetTitle(yAxisTitle);
 516 }
 517
 518
 519 TH1F* AliTRDClusterErrDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
 520 {
 521   Int_t nBins = hGen->GetNbinsX();
 522   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
 523   hEff->SetTitle("");
 524   hEff->SetStats(kFALSE);
 525   hEff->SetMinimum(0.);
 526   hEff->SetMaximum(110.);
 527
 528   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
 529     Double_t nGen = hGen->GetBinContent(iBin);
 530     Double_t nRec = hRec->GetBinContent(iBin);
 531     if (nGen > 0) {
 532       Double_t eff = nRec/nGen;
 533       hEff->SetBinContent(iBin, 100. * eff);
 534       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
 535       //      if (error == 0) error = sqrt(0.1/nGen);
 536       //
 537       if (error == 0) error = 0.0001;
 538       hEff->SetBinError(iBin, 100. * error);
 539     } else {
 540       hEff->SetBinContent(iBin, 100. * 0.5);
 541       hEff->SetBinError(iBin, 100. * 0.5);
 542     }
 543   }
 544
 545   return hEff;
 546 }
 547
 548
 549
 550 TH1F* AliTRDClusterErrDraw::CreateResHisto(TH2F* hRes2, Bool_t draw,  Bool_t drawBinFits,
 551                      Bool_t overflowBinFits)
 552 {
 553   TVirtualPad* currentPad = gPad;
 554   TAxis* axis = hRes2->GetXaxis();
 555   Int_t nBins = axis->GetNbins();
 556   TH1F* hRes, *hMean;
 557   if (axis->GetXbins()->GetSize()){
 558     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
 559     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
 560   }
 561   else{
 562     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
 563     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
 564
 565   }
 566   hRes->SetStats(false);
 567   hRes->SetOption("E");
 568   hRes->SetMinimum(0.);
 569   //
 570   hMean->SetStats(false);
 571   hMean->SetOption("E");
 572
 573   // create the fit function
 574   //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
 575   TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 576
 577   fitFunc->SetLineWidth(2);
 578   fitFunc->SetFillStyle(0);
 579   // create canvas for fits
 580   TCanvas* canBinFits = NULL;
 581   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
 582   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
 583   Int_t ny = (nPads-1) / nx + 1;
 584   if (drawBinFits) {
 585     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 586     if (canBinFits) delete canBinFits;
 587     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 588     canBinFits->Divide(nx, ny);
 589   }
 590
 591   // loop over x bins and fit projection
 592   Int_t dBin = ((overflowBinFits) ? 1 : 0);
 593   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 594     if (drawBinFits) canBinFits->cd(bin + dBin);
 595     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
 596     //
 597     if (hBin->GetEntries() > 10) {
 598       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS(),0.02*hBin->GetMaximum());
 599       hBin->Fit(fitFunc,"s");
 600       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 601
 602       if (sigma > 0.){
 603         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 604         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 605       }
 606       else{
 607         hRes->SetBinContent(bin, 0.);
 608         hMean->SetBinContent(bin,0);
 609       }
 610       hRes->SetBinError(bin, fitFunc->GetParError(2));
 611       hMean->SetBinError(bin, fitFunc->GetParError(1));
 612
 613       //
 614       //
 615
 616     } else {
 617       hRes->SetBinContent(bin, 0.);
 618       hRes->SetBinError(bin, 0.);
 619       hMean->SetBinContent(bin, 0.);
 620       hMean->SetBinError(bin, 0.);
 621     }
 622
 623
 624     if (drawBinFits) {
 625       char name[256];
 626       if (bin == 0) {
 627         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
 628       } else if (bin == nBins+1) {
 629         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
 630       } else {
 631         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
 632                 axis->GetTitle(), axis->GetBinUpEdge(bin));
 633       }
 634       canBinFits->cd(bin + dBin);
 635       hBin->SetTitle(name);
 636       hBin->SetStats(kTRUE);
 637       hBin->DrawCopy("E");
 638       canBinFits->Update();
 639       canBinFits->Modified();
 640       canBinFits->Update();
 641     }
 642
 643     delete hBin;
 644   }
 645
 646   delete fitFunc;
 647   currentPad->cd();
 648   if (draw){
 649     currentPad->Divide(1,2);
 650     currentPad->cd(1);
 651     hRes->Draw();
 652     currentPad->cd(2);
 653     hMean->Draw();
 654   }
 655
 656   return hRes;
 657 }