]>
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(); | |
f007cb5f | 85 | fTracker->LoadClusters(fClusterTree); |
289478c7 | 86 | // |
87 | } | |
59dfcadd | 88 | |
289478c7 | 89 | void 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 | 140 | void 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 | 182 | AliTrackReference * 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 | 205 | void 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 | 382 | Int_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 | 485 | AliTRDclusterCorrection* 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 | 516 | TH1F * 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 | 532 | TH1F * 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 | 550 | TH1F * 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 | 565 | void 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 | 572 | TH1F* 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 | 603 | TH1F* 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 | } |