Fixing problems in the generation of HTML documentation
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterErrors.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
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 "TF1.h"
20 #include "AliTrackReference.h"
21 #endif    
22 #include "AliTRDclusterErrors.h"
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 }