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