TRD/AliTRDclusterErrors.C

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