]> 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
New offline PID code by Prashant
[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();
f007cb5f 85 fTracker->LoadClusters(fClusterTree);
289478c7 86 //
87}
59dfcadd 88
289478c7 89void AliTRDClusterErrAnal::LoadClusters()
90{
91 //
92 //
93 // Load clusters
94 TObjArray *ClusterArray = new TObjArray(400);
59dfcadd 95 TObjArray carray(2000);
289478c7 96 TBranch *branch=fClusterTree->GetBranch("TRDcluster");
97 if (!branch) {
98 Error("ReadClusters","Can't get the branch !");
99 return;
100 }
f007cb5f 101 Int_t over5 =0;
102 Int_t over10=0;
103
289478c7 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);
f007cb5f 113 if (c->GetQ()>5) over5++;
114 if (c->GetQ()>10) over10++;
289478c7 115 }
116 }
59dfcadd 117 Int_t nClusters = carray.GetEntriesFast();
f007cb5f 118 printf("Total number of clusters %d\t%d\t%d\n", nClusters,over5,over10);
289478c7 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++;
59dfcadd 137 }
289478c7 138}
59dfcadd 139
289478c7 140void 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;
aec95541 155 Int_t nreferences=0;
156 Int_t nreftracks=0;
289478c7 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);
aec95541 168 nreftracks++;
289478c7 169 }
170 TClonesArray &larr = *labarr;
171 new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
aec95541 172 nreferences++;
289478c7 173 }
174 }
175 }
aec95541 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
289478c7 180}
59dfcadd 181
289478c7 182AliTrackReference * 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}
59dfcadd 204
289478c7 205void 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);
aec95541 258 //rot[0] *=-1;
259 // rot[1] *=-1;
289478c7 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);
aec95541 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
289478c7 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++;
aec95541 363 fSumQ +=hit->GetCharge();
289478c7 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}
59dfcadd 376
59dfcadd 377
59dfcadd 378
59dfcadd 379
59dfcadd 380
59dfcadd 381
289478c7 382Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) {
59dfcadd 383
289478c7 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);
59dfcadd 394
289478c7 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();
aec95541 405 printf("Found %d charged in TRD in first %d particles", nTrack, trackmax);
289478c7 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;
f007cb5f 432 clinfo->fNClusters=0;
289478c7 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);
f007cb5f 440 if (TMath::Abs(cluster->GetZ()-point->fTZ)>5.5 || dist>3.) continue;
441 clinfo->fNClusters++;
289478c7 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;
59dfcadd 461 }
289478c7 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}
59dfcadd 484
aec95541 485AliTRDclusterCorrection* 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}
59dfcadd 515
289478c7 516TH1F * 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}
59dfcadd 530
59dfcadd 531
289478c7 532TH1F * 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}
59dfcadd 548
549
289478c7 550TH1F * 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);
59dfcadd 555
289478c7 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);
59dfcadd 560
289478c7 561 tree->Draw(expression,selectionall);
562 return CreateResHisto(&hisdy);
563}
59dfcadd 564
289478c7 565void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
566{
567 histo->GetXaxis()->SetTitle(xAxisTitle);
568 histo->GetYaxis()->SetTitle(yAxisTitle);
569}
59dfcadd 570
59dfcadd 571
289478c7 572TH1F* 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 }
59dfcadd 597
289478c7 598 return hEff;
599}
59dfcadd 600
601
59dfcadd 602
289478c7 603TH1F* 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());
59dfcadd 617
289478c7 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 }
59dfcadd 643
289478c7 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));
59dfcadd 658 }
289478c7 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.);
59dfcadd 674 }
289478c7 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();
59dfcadd 707 }
708
289478c7 709 return hRes;
59dfcadd 710}