]>
Commit | Line | Data |
---|---|---|
289478c7 | 1 | //#ifndef __CINT__ |
59dfcadd | 2 | #include <iostream.h> |
3 | ||
4 | #include "AliTRDtracker.h" | |
289478c7 | 5 | #include "AliTRDclusterMI.h" |
59dfcadd | 6 | #include "AliTRDhit.h" |
7 | #include "AliTRDv1.h" | |
8 | #include "AliTRDgeometry.h" | |
289478c7 | 9 | #include "AliTRDgeometryDetail.h" |
59dfcadd | 10 | #include "AliTRDparameter.h" |
11 | ||
12 | #include "alles.h" | |
13 | #include "TFile.h" | |
14 | #include "TStopwatch.h" | |
289478c7 | 15 | #include "Rtypes.h" |
16 | #include "TTree.h" | |
17 | #include "AliRunLoader.h" | |
18 | #include "AliStack.h" | |
19 | #include "AliTRDclusterErrors.h" | |
20 | #include "TF1.h" | |
21 | #include "AliTrackReference.h" | |
22 | //#endif | |
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 | } | |
59dfcadd | 68 | |
289478c7 | 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 | } | |
59dfcadd | 86 | |
289478c7 | 87 | void AliTRDClusterErrAnal::LoadClusters() |
88 | { | |
89 | // | |
90 | // | |
91 | // Load clusters | |
92 | TObjArray *ClusterArray = new TObjArray(400); | |
59dfcadd | 93 | TObjArray carray(2000); |
289478c7 | 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 | } | |
59dfcadd | 110 | Int_t nClusters = carray.GetEntriesFast(); |
59dfcadd | 111 | printf("Total number of clusters %d \n", nClusters); |
289478c7 | 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++; | |
59dfcadd | 130 | } |
289478c7 | 131 | } |
59dfcadd | 132 | |
289478c7 | 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 | } | |
59dfcadd | 167 | |
289478c7 | 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 | } | |
59dfcadd | 190 | |
289478c7 | 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 | } | |
59dfcadd | 355 | |
59dfcadd | 356 | |
59dfcadd | 357 | |
59dfcadd | 358 | |
59dfcadd | 359 | |
59dfcadd | 360 | |
289478c7 | 361 | Int_t AliTRDClusterErrAnal::Analyze(Int_t trackmax) { |
59dfcadd | 362 | |
289478c7 | 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); | |
59dfcadd | 373 | |
289478c7 | 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; | |
59dfcadd | 438 | } |
289478c7 | 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 | } | |
59dfcadd | 461 | |
59dfcadd | 462 | |
289478c7 | 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 | } | |
59dfcadd | 477 | |
59dfcadd | 478 | |
289478c7 | 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 | } | |
59dfcadd | 495 | |
496 | ||
289478c7 | 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); | |
59dfcadd | 502 | |
289478c7 | 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); | |
59dfcadd | 507 | |
289478c7 | 508 | tree->Draw(expression,selectionall); |
509 | return CreateResHisto(&hisdy); | |
510 | } | |
59dfcadd | 511 | |
289478c7 | 512 | void AliTRDClusterErrDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle) |
513 | { | |
514 | histo->GetXaxis()->SetTitle(xAxisTitle); | |
515 | histo->GetYaxis()->SetTitle(yAxisTitle); | |
516 | } | |
59dfcadd | 517 | |
59dfcadd | 518 | |
289478c7 | 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 | } | |
59dfcadd | 544 | |
289478c7 | 545 | return hEff; |
546 | } | |
59dfcadd | 547 | |
548 | ||
59dfcadd | 549 | |
289478c7 | 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()); | |
59dfcadd | 564 | |
289478c7 | 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 | } | |
59dfcadd | 590 | |
289478c7 | 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)); | |
59dfcadd | 605 | } |
289478c7 | 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.); | |
59dfcadd | 621 | } |
289478c7 | 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(); | |
59dfcadd | 654 | } |
655 | ||
289478c7 | 656 | return hRes; |
59dfcadd | 657 | } |