]>
Commit | Line | Data |
---|---|---|
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 | ||
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 | } | |
59dfcadd | 69 | |
289478c7 | 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 | // | |
86 | } | |
59dfcadd | 87 | |
289478c7 | 88 | void 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 | 134 | void 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 | 176 | AliTrackReference * 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 | 199 | void 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 | 376 | Int_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 | 477 | AliTRDclusterCorrection* 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 | 508 | TH1F * 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 | 524 | TH1F * 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 | 542 | TH1F * 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 | 557 | void 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 | 564 | TH1F* 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 | 595 | TH1F* 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 | } |