]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDclusterErrors.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Clean-up. (A. Gheata)
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterErrors.C
CommitLineData
d8d3b5b8 1#if !defined(__CINT__) || defined(__MAKECINT__)
aec95541 2#include <Riostream.h>
59dfcadd 3
d8d3b5b8 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"
aec95541 11#include "AliTRDclusterCorrection.h"
d8d3b5b8 12#include "alles.h"
13#include "TFile.h"
14#include "TStopwatch.h"
15#include "Rtypes.h"
16#include "TTree.h"
aec95541 17
289478c7 18#include "AliRunLoader.h"
19#include "AliStack.h"
289478c7 20#include "TF1.h"
21#include "AliTrackReference.h"
d8d3b5b8 22#endif
23#include "AliTRDclusterErrors.h"
289478c7 24
25
26ClassImp(AliTRDCI)
27ClassImp(AliTRDExactPoint)
28ClassImp(AliTRDClusterErrAnal)
29ClassImp(AliTRDClusterErrDraw)
30
31
32
33AliTRDClusterErrAnal::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}
59dfcadd 69
289478c7 70void 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}
59dfcadd 87
289478c7 88void AliTRDClusterErrAnal::LoadClusters()
89{
90 //
91 //
92 // Load clusters
93 TObjArray *ClusterArray = new TObjArray(400);
59dfcadd 94 TObjArray carray(2000);
289478c7 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 }
59dfcadd 111 Int_t nClusters = carray.GetEntriesFast();
59dfcadd 112 printf("Total number of clusters %d \n", nClusters);
289478c7 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++;
59dfcadd 131 }
289478c7 132}
59dfcadd 133
289478c7 134void 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;
aec95541 149 Int_t nreferences=0;
150 Int_t nreftracks=0;
289478c7 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);
aec95541 162 nreftracks++;
289478c7 163 }
164 TClonesArray &larr = *labarr;
165 new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
aec95541 166 nreferences++;
289478c7 167 }
168 }
169 }
aec95541 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
289478c7 174}
59dfcadd 175
289478c7 176AliTrackReference * 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}
59dfcadd 198
289478c7 199void 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);
aec95541 252 //rot[0] *=-1;
253 // rot[1] *=-1;
289478c7 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);
aec95541 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
289478c7 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++;
aec95541 357 fSumQ +=hit->GetCharge();
289478c7 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}
59dfcadd 370
59dfcadd 371
59dfcadd 372
59dfcadd 373
59dfcadd 374
59dfcadd 375
289478c7 376Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) {
59dfcadd 377
289478c7 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);
59dfcadd 388
289478c7 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();
aec95541 399 printf("Found %d charged in TRD in first %d particles", nTrack, trackmax);
289478c7 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;
59dfcadd 453 }
289478c7 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}
59dfcadd 476
aec95541 477AliTRDclusterCorrection* 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}
59dfcadd 507
289478c7 508TH1F * 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}
59dfcadd 522
59dfcadd 523
289478c7 524TH1F * 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}
59dfcadd 540
541
289478c7 542TH1F * 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);
59dfcadd 547
289478c7 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);
59dfcadd 552
289478c7 553 tree->Draw(expression,selectionall);
554 return CreateResHisto(&hisdy);
555}
59dfcadd 556
289478c7 557void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
558{
559 histo->GetXaxis()->SetTitle(xAxisTitle);
560 histo->GetYaxis()->SetTitle(yAxisTitle);
561}
59dfcadd 562
59dfcadd 563
289478c7 564TH1F* 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 }
59dfcadd 589
289478c7 590 return hEff;
591}
59dfcadd 592
593
59dfcadd 594
289478c7 595TH1F* 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());
59dfcadd 609
289478c7 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 }
59dfcadd 635
289478c7 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));
59dfcadd 650 }
289478c7 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.);
59dfcadd 666 }
289478c7 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();
59dfcadd 699 }
700
289478c7 701 return hRes;
59dfcadd 702}