]>
Commit | Line | Data |
---|---|---|
ae7d73d2 | 1 | #if !defined(__CINT__) || defined(__MAKECINT__) |
2 | #include <Riostream.h> | |
3 | ||
4 | #include "AliITStrackerMI.h" | |
5 | #include "AliITSclusterV2.h" | |
6 | ||
7 | #include "AliITShit.h" | |
8 | #include "AliITS.h" | |
9 | #include "AliITSgeom.h" | |
10 | //#include "AliITSgeometryDetail.h" | |
11 | //#include "AliITSparameter.h" | |
12 | #include "alles.h" | |
13 | #include "TFile.h" | |
14 | #include "TStopwatch.h" | |
15 | #include "Rtypes.h" | |
16 | #include "TTree.h" | |
17 | ||
18 | ||
19 | #include "AliRunLoader.h" | |
20 | #include "AliStack.h" | |
21 | #include "TF1.h" | |
22 | #include "AliTrackReference.h" | |
23 | #endif | |
24 | #include "AliITSclusterComparison.h" | |
25 | ||
26 | ||
27 | ClassImp(AliITSCI) | |
28 | ClassImp(AliITSClusterErrAnal) | |
29 | ||
30 | ||
31 | ||
32 | AliITSClusterErrAnal::AliITSClusterErrAnal(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 | fITSLoader = fRunLoader->GetLoader("ITSLoader"); | |
47 | if (fITSLoader == 0x0){ | |
48 | cerr<<"Can not get ITS Loader"<<endl; | |
49 | return ; | |
50 | } | |
51 | if (fRunLoader->LoadgAlice()){ | |
52 | cerr<<"Error occured while l"<<endl; | |
53 | return; | |
54 | } | |
55 | fRunLoader->CdGAFile(); | |
56 | fITS = (AliITS*)gAlice->GetDetector("ITS"); | |
57 | if (!fITS) { | |
58 | cerr<<"AliITSclusterComparison.C : Can not find the ITS detector !"<<endl; | |
59 | } | |
60 | AliITSgeom *geom = fITS->GetITSgeom(); | |
61 | //An instance of the ITS tracker | |
62 | fTracker = new AliITStrackerMI(geom); | |
63 | // | |
64 | // | |
65 | AliITSCI * clinfo = new AliITSCI(); | |
66 | fFileA = new TFile("itsclusteranal.root","recreate"); | |
67 | fFileA->cd(); | |
68 | fTreeA = new TTree("itscl","itscl"); | |
69 | fTreeA->Branch("itscl","AliITSCI",&clinfo); | |
70 | ||
71 | AliITSclusterV2 * cl = new AliITSclusterV2; | |
72 | fTreeB = new TTree("Clusters","Clusters"); | |
73 | fTreeB->Branch("cl","AliITSclusterV2",&cl); | |
74 | ||
75 | fClusters = 0; | |
76 | delete clinfo; | |
77 | } | |
78 | ||
79 | void AliITSClusterErrAnal::SetIO(Int_t event) | |
80 | { | |
81 | // | |
82 | //set input output for given event | |
83 | fRunLoader->SetEventNumber(event); | |
84 | fRunLoader->LoadHeader(); | |
85 | fRunLoader->LoadKinematics(); | |
86 | fRunLoader->LoadTrackRefs(); | |
87 | fITSLoader->LoadHits(); | |
88 | fITSLoader->LoadRecPoints("read"); | |
89 | // | |
90 | fStack = fRunLoader->Stack(); | |
91 | fHitTree = fITSLoader->TreeH(); | |
92 | fClusterTree = fITSLoader->TreeR(); | |
93 | fReferenceTree = fRunLoader->TreeTR(); | |
94 | // | |
95 | } | |
96 | ||
97 | void AliITSClusterErrAnal::LoadClusters() | |
98 | { | |
99 | // | |
100 | // | |
101 | // Load clusters | |
102 | if (fClusters) { | |
103 | fClusters->Delete(); | |
104 | delete fClusters; | |
105 | } | |
106 | Int_t nparticles = fStack->GetNtrack(); | |
107 | fClusters = new TObjArray(nparticles); | |
108 | // | |
109 | TObjArray *ClusterArray = new TClonesArray("AliITSclusterV2"); | |
110 | TObjArray carray(2000); | |
111 | TBranch *branch=fClusterTree->GetBranch("Clusters"); | |
112 | if (!branch) { | |
113 | Error("ReadClusters","Can't get the branch !"); | |
114 | return; | |
115 | } | |
116 | branch->SetAddress(&ClusterArray); | |
117 | TBranch *brcl = fTreeB->GetBranch("cl"); | |
118 | AliITSclusterV2* cluster= new AliITSclusterV2; | |
119 | brcl->SetAddress(&cluster); | |
120 | Int_t nentries = (Int_t)fClusterTree->GetEntries(); | |
121 | for (Int_t i=0;i<nentries;i++){ | |
122 | fClusterTree->GetEvent(i); | |
123 | Int_t nCluster = ClusterArray->GetEntriesFast(); | |
124 | SignDeltas(ClusterArray,-3.8); | |
125 | for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { | |
126 | cluster = (AliITSclusterV2*)ClusterArray->UncheckedAt(iCluster); | |
127 | carray.AddLast(new AliITSclusterV2(*cluster)); | |
128 | fTreeB->Fill(); | |
129 | } | |
130 | } | |
131 | // | |
132 | Int_t nClusters = carray.GetEntriesFast(); | |
133 | printf("Total number of clusters %d \n", nClusters); | |
134 | fTracker->LoadClusters(fClusterTree); | |
135 | // | |
136 | // | |
137 | //SORT clusters | |
138 | // | |
139 | Int_t all=0; | |
140 | Int_t noisecl =0; | |
141 | for (Int_t i=0;i<nClusters;i++){ | |
142 | AliITSclusterV2 *cl = (AliITSclusterV2 *) carray.UncheckedAt(i); | |
143 | if (cl->GetLabel(0)<0) noisecl++; | |
144 | // | |
145 | for (Int_t itrack=0; itrack<3;itrack++){ | |
146 | Int_t lab = cl->GetLabel(itrack); | |
147 | if (lab>=0){ | |
148 | TObjArray * array = (TObjArray*)fClusters->At(lab); | |
149 | if (array==0){ | |
150 | array = new TObjArray(20); | |
151 | fClusters->AddAt(array,lab); | |
152 | } | |
153 | array->AddLast(cl); | |
154 | all++; | |
155 | } | |
156 | } | |
157 | } | |
158 | printf("Total number of track clusters %d \n", all); | |
159 | printf("Total number of noise clusters %d \n", noisecl); | |
160 | } | |
161 | ||
162 | void AliITSClusterErrAnal::SignDeltas( TObjArray *ClusterArray, Float_t vz) | |
163 | { | |
164 | // | |
165 | // | |
166 | Int_t entries = ClusterArray->GetEntriesFast(); | |
167 | if (entries<4) return; | |
168 | AliITSclusterV2* cluster = (AliITSclusterV2*)ClusterArray->At(0); | |
169 | Int_t layer = cluster->GetLayer(); | |
170 | if (layer>1) return; | |
171 | Int_t index[10000]; | |
172 | Int_t ncandidates=0; | |
173 | Float_t r = (layer>0)? 7:4; | |
174 | // | |
175 | for (Int_t i=0;i<entries;i++){ | |
176 | AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(i); | |
177 | Float_t nz = 1+(cl0->GetZ()-vz)/r; | |
178 | if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue; | |
179 | index[ncandidates] = i; //candidate to belong to delta electron track | |
180 | ncandidates++; | |
181 | if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) { | |
182 | cl0->SetDeltaProbability(1); | |
183 | } | |
184 | } | |
185 | // | |
186 | // | |
187 | // | |
188 | for (Int_t i=0;i<ncandidates;i++){ | |
189 | AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(index[i]); | |
190 | if (cl0->GetDeltaProbability()>0.8) continue; | |
191 | // | |
192 | Int_t ncl = 0; | |
193 | Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw; | |
194 | sumy=sumz=sumy2=sumyz=sumw=0.0; | |
195 | for (Int_t j=0;j<ncandidates;j++){ | |
196 | if (i==j) continue; | |
197 | AliITSclusterV2* cl1 = (AliITSclusterV2*)ClusterArray->At(index[j]); | |
198 | // | |
199 | Float_t dz = cl0->GetZ()-cl1->GetZ(); | |
200 | Float_t dy = cl0->GetY()-cl1->GetY(); | |
201 | if (TMath::Sqrt(dz*dz+dy*dy)<0.2){ | |
202 | Float_t weight = cl1->GetNy()+cl1->GetNz()-2; | |
203 | y[ncl] = cl1->GetY(); | |
204 | z[ncl] = cl1->GetZ(); | |
205 | sumy+= y[ncl]*weight; | |
206 | sumz+= z[ncl]*weight; | |
207 | sumy2+=y[ncl]*y[ncl]*weight; | |
208 | sumyz+=y[ncl]*z[ncl]*weight; | |
209 | sumw+=weight; | |
210 | ncl++; | |
211 | } | |
212 | } | |
213 | if (ncl<4) continue; | |
214 | Float_t det = sumw*sumy2 - sumy*sumy; | |
215 | Float_t delta=1000; | |
216 | if (TMath::Abs(det)>0.01){ | |
217 | Float_t z0 = (sumy2*sumz - sumy*sumyz)/det; | |
218 | Float_t k = (sumyz*sumw - sumy*sumz)/det; | |
219 | delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY())); | |
220 | } | |
221 | else{ | |
222 | Float_t z0 = sumyz/sumy; | |
223 | delta = TMath::Abs(cl0->GetZ()-z0); | |
224 | } | |
225 | if ( delta<0.05) { | |
226 | cl0->SetDeltaProbability(1-20.*delta); | |
227 | } | |
228 | } | |
229 | } | |
230 | ||
231 | void AliITSClusterErrAnal::LoadParticles() | |
232 | { | |
233 | // | |
234 | // | |
235 | // fReferences = new TObjArray; | |
236 | } | |
237 | ||
238 | void AliITSClusterErrAnal::SortReferences() | |
239 | { | |
240 | // | |
241 | // | |
242 | // | |
243 | printf("Sorting references\n"); | |
244 | fReferences = new TObjArray; | |
245 | Int_t ntracks = fStack->GetNtrack(); | |
246 | fReferences->Expand(ntracks); | |
247 | Int_t nentries = (Int_t)fReferenceTree->GetEntries(); | |
248 | TClonesArray * arr = new TClonesArray("AliTrackReference"); | |
249 | TBranch * br = fReferenceTree->GetBranch("ITS"); | |
250 | br->SetAddress(&arr); | |
251 | // | |
252 | TClonesArray *labarr=0; | |
253 | Int_t nreferences=0; | |
254 | Int_t nreftracks=0; | |
255 | for (Int_t iprim=0;iprim<nentries;iprim++){ | |
256 | if (br->GetEntry(iprim)){ | |
257 | for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){ | |
258 | AliTrackReference *ref =(AliTrackReference*)arr->At(iref); | |
259 | if (!ref) continue; | |
260 | Int_t lab = ref->GetTrack(); | |
261 | //if ( (lab<0) || (lab>ntracks)) continue; | |
262 | // | |
263 | if (fReferences->At(lab)==0) { | |
264 | labarr = new TClonesArray("AliTrackReference"); | |
265 | fReferences->AddAt(labarr,lab); | |
266 | nreftracks++; | |
267 | } | |
268 | TClonesArray &larr = *labarr; | |
269 | new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref); | |
270 | nreferences++; | |
271 | } | |
272 | } | |
273 | } | |
274 | printf("Total number of references = \t%d\n", nreferences); | |
275 | printf("Total number of tracks with references = \t%d\n", nreftracks); | |
276 | printf("End - Sorting references\n"); | |
277 | ||
278 | } | |
279 | ||
280 | ||
281 | void AliITSClusterErrAnal::GetNTeor(Int_t layer, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz) | |
282 | { | |
283 | // | |
284 | //get "mean shape" | |
285 | if (layer==0){ | |
286 | ny = 1.+TMath::Abs(phi)*3.2; | |
287 | nz = 1.+TMath::Abs(theta)*0.34; | |
288 | return; | |
289 | } | |
290 | if (layer==1){ | |
291 | ny = 1.+TMath::Abs(phi)*3.2; | |
292 | nz = 1.+TMath::Abs(theta)*0.28; | |
293 | return; | |
294 | } | |
295 | ||
296 | if (layer>3){ | |
297 | ny = 2.02+TMath::Abs(phi)*1.95; | |
298 | nz = 2.02+TMath::Abs(phi)*2.35; | |
299 | return; | |
300 | } | |
301 | ny = 6.6-2.7*abs(phi); | |
302 | nz = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta); | |
303 | } | |
304 | ||
305 | ||
306 | Int_t AliITSClusterErrAnal::GetError(Int_t layer, const AliITSclusterV2*cl, Float_t theta, Float_t phi,Float_t expQ, Float_t &erry, Float_t &errz) | |
307 | { | |
308 | //calculate cluster position error | |
309 | // | |
310 | Float_t nz,ny; | |
311 | GetNTeor(layer, theta,phi,ny,nz); | |
312 | erry = TMath::Sqrt(cl->GetSigmaY2()); | |
313 | errz = TMath::Sqrt(cl->GetSigmaZ2()); | |
314 | // | |
315 | // PIXELS | |
316 | if (layer<2){ | |
317 | ||
318 | if (TMath::Abs(ny-cl->GetNy())>0.6) { | |
319 | if (ny<cl->GetNy()){ | |
320 | erry*=0.4+TMath::Abs(ny-cl->GetNy()); | |
321 | errz*=0.4+TMath::Abs(ny-cl->GetNy()); | |
322 | }else{ | |
323 | erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy()); | |
324 | errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy()); | |
325 | } | |
326 | } | |
327 | if (TMath::Abs(nz-cl->GetNz())>1.) { | |
328 | erry*=TMath::Abs(nz-cl->GetNz()); | |
329 | errz*=TMath::Abs(nz-cl->GetNz()); | |
330 | } | |
331 | erry*=0.85; | |
332 | errz*=0.85; | |
333 | erry= TMath::Min(erry,float(0.005)); | |
334 | errz= TMath::Min(errz,float(0.03)); | |
335 | return 10; | |
336 | } | |
337 | //STRIPS | |
338 | if (layer>3){ | |
339 | if (cl->GetNy()==100||cl->GetNz()==100){ | |
340 | erry = 0.004; | |
341 | errz = 0.2; | |
342 | return 100; | |
343 | } | |
344 | if (cl->GetNy()+cl->GetNz()>12){ | |
345 | erry = 0.06; | |
346 | errz = 0.57; | |
347 | return 100; | |
348 | } | |
349 | Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi)); | |
350 | Float_t chargematch = TMath::Max(double(normq/expQ),2.); | |
351 | // | |
352 | if (cl->GetType()==1 || cl->GetType()==10 ){ | |
353 | if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){ | |
354 | errz = 0.043; | |
355 | erry = 0.00094; | |
356 | return 101; | |
357 | } | |
358 | if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){ | |
359 | errz = 0.06; | |
360 | erry =0.0013; | |
361 | return 102; | |
362 | } | |
363 | erry = 0.0027; | |
364 | errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15); | |
365 | return 103; | |
366 | } | |
367 | if (cl->GetType()==2 || cl->GetType()==11 ){ | |
368 | erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05); | |
369 | errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5); | |
370 | return 104; | |
371 | } | |
372 | ||
373 | if (cl->GetType()>100 ){ | |
374 | if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){ | |
375 | errz = 0.05; | |
376 | erry = 0.00096; | |
377 | return 105; | |
378 | } | |
379 | if (cl->GetNy()+cl->GetNz()-nz-ny<1){ | |
380 | errz = 0.10; | |
381 | erry = 0.0025; | |
382 | return 106; | |
383 | } | |
384 | ||
385 | errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4); | |
386 | erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05); | |
387 | return 107; | |
388 | } | |
389 | Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz; | |
390 | if (diff<1) diff=1; | |
391 | if (diff>4) diff=4; | |
392 | ||
393 | if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){ | |
394 | errz = 0.14*diff; | |
395 | erry = 0.003*diff; | |
396 | return 108; | |
397 | } | |
398 | erry = 0.04*diff; | |
399 | errz = 0.06*diff; | |
400 | return 109; | |
401 | } | |
402 | ||
403 | //DRIFTS | |
404 | Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi)); | |
405 | Float_t chargematch = normq/expQ; | |
406 | Float_t factorz=1; | |
407 | Int_t cnz = cl->GetNz()%10; | |
408 | //charge match | |
409 | if (cl->GetType()==1){ | |
410 | if (chargematch<1.25){ | |
411 | erry = 0.0028*(1.+6./cl->GetQ()); // gold clusters | |
412 | } | |
413 | else{ | |
414 | erry = 0.003*chargematch; | |
415 | if (cl->GetNz()==3) erry*=1.5; | |
416 | } | |
417 | if (chargematch<1.0){ | |
418 | errz = 0.0011*(1.+6./cl->GetQ()); | |
419 | } | |
420 | else{ | |
421 | errz = 0.002*(1+2*(chargematch-1.)); | |
422 | } | |
423 | if (cnz>nz+0.6) { | |
424 | erry*=(cnz-nz+0.5); | |
425 | errz*=1.4*(cnz-nz+0.5); | |
426 | } | |
427 | } | |
428 | if (cl->GetType()>1){ | |
429 | if (chargematch<1){ | |
430 | erry = 0.00385*(1.+6./cl->GetQ()); // gold clusters | |
431 | errz = 0.0016*(1.+6./cl->GetQ()); | |
432 | } | |
433 | else{ | |
434 | errz = 0.0014*(1+3*(chargematch-1.)); | |
435 | erry = 0.003*(1+3*(chargematch-1.)); | |
436 | } | |
437 | if (cnz>nz+0.6) { | |
438 | erry*=(cnz-nz+0.5); | |
439 | errz*=1.4*(cnz-nz+0.5); | |
440 | } | |
441 | } | |
442 | ||
443 | if (TMath::Abs(cl->GetY())>2.5){ | |
444 | factorz*=1+2*(TMath::Abs(cl->GetY())-2.5); | |
445 | } | |
446 | if (TMath::Abs(cl->GetY())<1){ | |
447 | factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.); | |
448 | } | |
449 | factorz= TMath::Min(factorz,float(4.)); | |
450 | errz*=factorz; | |
451 | ||
452 | erry= TMath::Min(erry,float(0.05)); | |
453 | errz= TMath::Min(errz,float(0.05)); | |
454 | return 200; | |
455 | } | |
456 | ||
457 | ||
458 | AliITSclusterV2 * AliITSClusterErrAnal::FindNearestCluster(AliITSCI * clinfo, Float_t dmax) | |
459 | { | |
460 | // | |
461 | // | |
462 | AliTrackReference * ref = &(clinfo->fRef); | |
463 | TObjArray * clusters = (TObjArray*)(fClusters->At(ref->GetTrack())); | |
464 | if (!clusters) return 0; | |
465 | Int_t nclusters = clusters->GetEntriesFast(); | |
466 | Float_t maxd = dmax; | |
467 | AliITSclusterV2 * cluster=0; | |
468 | Float_t radius = ref->R(); | |
469 | clinfo->fNClusters=0; | |
470 | for (Int_t icluster=0;icluster<nclusters;icluster++){ | |
471 | AliITSclusterV2 * cl = (AliITSclusterV2*)clusters->At(icluster); | |
472 | if (!cl) continue; | |
473 | Float_t dz = cl->GetZ()-ref->Z(); | |
474 | // | |
475 | if (TMath::Abs(dz)>dmax) continue; | |
476 | Int_t idet=cl->GetDetectorIndex(); | |
477 | Double_t r=-1; | |
478 | Float_t dy =0; | |
479 | // for (Int_t i=cl->fLayer;i<=cl->fLayer-1;i++){ | |
480 | ||
481 | AliITStrackerMI::AliITSdetector & det = fTracker->GetDetector(cl->fLayer,idet); | |
482 | if (TMath::Abs(radius-det.GetR())<1.5) { | |
483 | Bool_t isOK=kFALSE; | |
484 | for (Int_t itrack=0;itrack<3;itrack++){ | |
485 | if (cl->fTracks[itrack]== clinfo->fRef.GetTrack()) { | |
486 | clinfo->fNClusters++; | |
487 | isOK= kTRUE; | |
488 | } | |
489 | } | |
490 | if (!isOK) continue; | |
491 | r = det.GetR(); | |
492 | Double_t phicl= det.GetPhi(); | |
493 | Double_t ca=TMath::Cos(phicl), sa=TMath::Sin(phicl); | |
494 | Double_t x = double(ref->X())*ca +double(ref->Y())*sa; | |
495 | Double_t y = -double(ref->X())*sa +double(ref->Y())*ca; | |
496 | if (TMath::Abs(x-r-0.012)>0.1) continue; | |
497 | // | |
498 | Double_t px = ref->Px(), py = ref->Py(); | |
499 | Double_t lpx = double(px)*ca +double(py)*sa; | |
500 | Double_t lpy = -double(px)*sa +double(py)*ca; | |
501 | Double_t lphi = TMath::ATan2(lpy,lpx); | |
502 | ||
503 | Double_t theta = ref->Pz()/ref->Pt(); | |
504 | Double_t dx = x-r; | |
505 | dy = y -cl->GetY()-dx*lphi; | |
506 | dz = ref->Z()-cl->GetZ()-dx*theta; | |
507 | //dx-=0.012; | |
508 | Double_t d = TMath::Sqrt(dy*dy+dz*dz); | |
509 | if (d<maxd){ | |
510 | cluster=cl; | |
511 | maxd=d; | |
512 | clinfo->Update(fTracker); | |
513 | clinfo->fCl = *cluster; | |
514 | clinfo->fLx = x; | |
515 | clinfo->fLy = y; | |
516 | clinfo->fLayer = cl->fLayer; | |
517 | clinfo->fPhiCl = phicl; | |
518 | clinfo->fRCl = det.GetR(); | |
519 | clinfo->fPhiL = lphi; | |
520 | clinfo->fDz = dz; | |
521 | clinfo->fDy = dy; | |
522 | cl->SetSigmaZ2(TMath::Abs(cl->GetSigmaZ2())); | |
523 | cl->SetSigmaY2(TMath::Abs(cl->GetSigmaY2())); | |
524 | clinfo->fPoolZ = dz/TMath::Sqrt(cl->GetSigmaZ2()); | |
525 | clinfo->fPoolY = dy/TMath::Sqrt(cl->GetSigmaY2()); | |
526 | clinfo->fPoolY2 = dy/TMath::Sqrt(cl->GetSigmaY2()); | |
527 | clinfo->fErrY = TMath::Sqrt(cl->GetSigmaY2()); | |
528 | clinfo->fErrZ = TMath::Sqrt(cl->GetSigmaZ2()); | |
529 | Float_t erry,errz,ny,nz; | |
530 | clinfo->fErrType = GetError(cl->fLayer,cl,theta,lphi,clinfo->fExpectedQ,erry,errz); | |
531 | GetNTeor(cl->fLayer, theta,lphi,ny,nz); | |
532 | clinfo->fNormQ = cl->GetQ()/TMath::Sqrt(1+theta*theta+lpy*lpy/(lpx*lpx)); | |
533 | clinfo->fTNy = ny; | |
534 | clinfo->fTNz = nz; | |
535 | clinfo->fErrY = erry; | |
536 | clinfo->fErrZ = errz; | |
537 | clinfo->fPoolY2 = dy/clinfo->fErrY; | |
538 | clinfo->fPoolZ2 = dz/clinfo->fErrZ; | |
539 | clinfo->fLabPos =-1; | |
540 | Double_t xyz[3]; | |
541 | det.GetGlobalXYZ(cl,xyz); | |
542 | clinfo->fGx = xyz[0]; | |
543 | clinfo->fGy = xyz[1]; | |
544 | for (Int_t i=0;i<3;i++) if (cl->GetLabel(i)==clinfo->fRef.GetTrack()) clinfo->fLabPos = i; | |
545 | } | |
546 | } | |
547 | } | |
548 | // } | |
549 | return cluster; | |
550 | } | |
551 | ||
552 | ||
553 | ||
554 | Int_t AliITSClusterErrAnal::Analyze(Int_t trackmax) { | |
555 | // | |
556 | // | |
557 | // dummy cluster to be fill if not cluster info | |
558 | // | |
559 | AliITSCI * clinfo = new AliITSCI; | |
560 | AliITSclusterV2 dummy; | |
561 | TBranch * branch = fTreeA->GetBranch("itscl"); | |
562 | branch->SetAddress(&clinfo); | |
563 | Int_t nall =0; | |
564 | Int_t withcl =0; | |
565 | Int_t ntracks = fReferences->GetEntriesFast(); | |
566 | ntracks = TMath::Min(trackmax,ntracks); | |
567 | for (Int_t itrack=0;itrack<ntracks;itrack++){ | |
568 | TClonesArray * references = (TClonesArray*) fReferences->At(itrack); | |
569 | if (!references) continue; | |
570 | Int_t nreferences = references->GetEntriesFast(); | |
571 | if (nreferences<4) continue; | |
572 | TObjArray * clarray = (TObjArray*)fClusters->At(itrack); | |
573 | for (Int_t iref=0;iref<nreferences;iref++){ | |
574 | AliTrackReference* trackref = (AliTrackReference*)references->At(iref); | |
575 | if (trackref&& trackref->P()>0.01){ | |
576 | TParticle * p = fStack->Particle(trackref->GetTrack()); | |
577 | if (p) clinfo->fP = *p; | |
578 | clinfo->fRef = *trackref; | |
579 | clinfo->fStatus = 0; | |
580 | clinfo->Update(fTracker); | |
581 | nall++; | |
582 | AliITSclusterV2 * cluster = FindNearestCluster(clinfo,3.0); | |
583 | if (cluster){ | |
584 | //clinfo->fCl = *cluster; | |
585 | clinfo->fStatus=1; | |
586 | withcl++; | |
587 | }else{ | |
588 | clinfo->fCl= dummy; | |
589 | } | |
590 | clinfo->Update(fTracker); | |
591 | fTreeA->Fill(); | |
592 | } | |
593 | } | |
594 | } | |
595 | ||
596 | fFileA->cd(); | |
597 | fTreeA->Write(); | |
598 | fTreeB->Write(); | |
599 | fFileA->Close(); | |
600 | return 0; | |
601 | } | |
602 | ||
603 | ||
604 | ||
605 | void AliITSClusterErrAnal::MakeSeeds(Double_t zv) | |
606 | { | |
607 | // | |
608 | // | |
609 | // | |
610 | ||
611 | AliITStrackerMI::AliITSlayer & layer3 = fTracker->GetLayer(3); | |
612 | AliITStrackerMI::AliITSlayer & layer2 = fTracker->GetLayer(2); | |
613 | AliITStrackerMI::AliITSlayer & layer1 = fTracker->GetLayer(1); | |
614 | AliITStrackerMI::AliITSlayer & layer0 = fTracker->GetLayer(0); | |
615 | Int_t ncl3 = layer3.GetNumberOfClusters(); | |
616 | Int_t ncl2 = layer2.GetNumberOfClusters(); | |
617 | Int_t ncl1 = layer1.GetNumberOfClusters(); | |
618 | Int_t ncl0 = layer0.GetNumberOfClusters(); | |
619 | TH1F *hc = new TH1F("hc","hc",100,-0.01,0.01); | |
620 | TH1F *hrc = new TH1F("hrc","hrc",100,-1,1); | |
621 | TH1F *hl = new TH1F("hl","hl",100,-0.05,0.05); | |
622 | TH1F *hlr = new TH1F("hlr","hlr",1000,-30.,30.); | |
623 | TH2F *hcl = new TH2F("hcl","hcl",1000,-0.01,0.01,1000,-0.03,0.03); | |
624 | TH2F *hcc = new TH2F("hcc","hcc",1000,-0.01,0.01,1000,-0.03,0.03); | |
625 | TH1F *hr = new TH1F("hr","hr",100,0,2000.); | |
626 | // | |
627 | TH1F *hpoolc = new TH1F("hpoolc","hpoolc",100,-5,5); | |
628 | TH1F *hpooll = new TH1F("hpooll","hpooll",100,-5,5); | |
629 | // | |
630 | Double_t ratio23 = layer2.GetR()/layer3.GetR(); | |
631 | Double_t ratio12 = layer1.GetR()/layer2.GetR(); | |
632 | Double_t referenceR = (layer0.GetR()+layer1.GetR()+layer2.GetR()+layer3.GetR())*0.25; | |
633 | // | |
634 | Double_t *clusterxyzr3 = new Double_t[ncl3*6]; | |
635 | Double_t *clusterxyzr2 = new Double_t[ncl2*6]; | |
636 | Double_t *clusterxyzr1 = new Double_t[ncl1*6]; | |
637 | Double_t *clusterxyzr0 = new Double_t[ncl0*6]; | |
638 | // | |
639 | // Get global coordinates 3 | |
640 | for (Int_t icl3=0;icl3<ncl3;icl3++){ | |
641 | AliITSclusterV2* cl3 = layer3.GetCluster(icl3); | |
642 | Double_t *xyz3 = &clusterxyzr3[6*icl3]; | |
643 | AliITStrackerMI::AliITSdetector & det = layer3.GetDetector(cl3->GetDetectorIndex()); | |
644 | det.GetGlobalXYZ(cl3,xyz3); | |
645 | xyz3[3] = TMath::Sqrt(xyz3[0]*xyz3[0]+xyz3[1]*xyz3[1]); | |
646 | Double_t fi3 = TMath::ATan2(xyz3[1],xyz3[0]); | |
647 | xyz3[4] = TMath::Cos(fi3); | |
648 | xyz3[5] = TMath::Sin(fi3); | |
649 | } | |
650 | // | |
651 | // Get global coordinates 2 | |
652 | for (Int_t icl2=0;icl2<ncl2;icl2++){ | |
653 | AliITSclusterV2* cl2 = layer2.GetCluster(icl2); | |
654 | Double_t *xyz2 = &clusterxyzr2[6*icl2]; | |
655 | AliITStrackerMI::AliITSdetector & det = layer2.GetDetector(cl2->GetDetectorIndex()); | |
656 | det.GetGlobalXYZ(cl2,xyz2); | |
657 | xyz2[3] = TMath::Sqrt(xyz2[0]*xyz2[0]+xyz2[1]*xyz2[1]); | |
658 | Double_t fi2 = TMath::ATan2(xyz2[1],xyz2[0]); | |
659 | xyz2[4] = TMath::Cos(fi2); | |
660 | xyz2[5] = TMath::Sin(fi2); | |
661 | } | |
662 | // | |
663 | // Get global coordinates 1 | |
664 | for (Int_t icl1=0;icl1<ncl1;icl1++){ | |
665 | AliITSclusterV2* cl1 = layer1.GetCluster(icl1); | |
666 | Double_t *xyz1 = &clusterxyzr1[6*icl1]; | |
667 | AliITStrackerMI::AliITSdetector & det = layer1.GetDetector(cl1->GetDetectorIndex()); | |
668 | det.GetGlobalXYZ(cl1,xyz1); | |
669 | xyz1[3] = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]); | |
670 | Double_t fi1 = TMath::ATan2(xyz1[1],xyz1[0]); | |
671 | xyz1[4] = TMath::Cos(fi1); | |
672 | xyz1[5] = TMath::Sin(fi1); | |
673 | } | |
674 | // | |
675 | // Get global coordinates 0 | |
676 | for (Int_t icl0=0;icl0<ncl0;icl0++){ | |
677 | AliITSclusterV2* cl0 = layer0.GetCluster(icl0); | |
678 | Double_t *xyz0 = &clusterxyzr0[6*icl0]; | |
679 | AliITStrackerMI::AliITSdetector & det = layer0.GetDetector(cl0->GetDetectorIndex()); | |
680 | det.GetGlobalXYZ(cl0,xyz0); | |
681 | xyz0[3] = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]); | |
682 | Double_t fi0 = TMath::ATan2(xyz0[1],xyz0[0]); | |
683 | xyz0[4] = TMath::Cos(fi0); | |
684 | xyz0[5] = TMath::Sin(fi0); | |
685 | ||
686 | } | |
687 | // | |
688 | // | |
689 | Int_t ntracks=0; | |
690 | Int_t ntracks0=0; | |
691 | Int_t good0=0; | |
692 | Int_t accepted0=0; | |
693 | // | |
694 | Int_t nfoundtracks=0; | |
695 | Int_t good1 =0;; | |
696 | Int_t accepted1=0; | |
697 | Int_t accepted2=0; | |
698 | // | |
699 | Int_t *trackindex = new Int_t[ncl3*4*300]; | |
700 | const AliITSclusterV2 **trackclusters = new const AliITSclusterV2*[ncl3*4*300]; | |
701 | // | |
702 | for (Int_t icl3=0;icl3<ncl3;icl3++){ | |
703 | AliITSclusterV2* cl3 = layer3.GetCluster(icl3); | |
704 | Double_t *xyz3 = &clusterxyzr3[icl3*6]; | |
705 | Double_t &r3 = clusterxyzr3[icl3*6+3]; | |
706 | ntracks=0; //number of possible tracks corresponding to given cluster | |
707 | // | |
708 | // | |
709 | // find pairs in layer 2 | |
710 | Float_t zc = (xyz3[2]-zv)*ratio23+zv; | |
711 | Float_t lphi = TMath::ATan2(xyz3[1],xyz3[0]); | |
712 | layer2.SelectClusters(zc-1,zc+1, layer2.GetR()*lphi-1, layer2.GetR()*lphi+1.); | |
713 | const AliITSclusterV2 * cl2=0; | |
714 | Int_t npairs2=0; | |
715 | Int_t ci2[100]; | |
716 | while ( (cl2=layer2.GetNextCluster(ci2[npairs2]))!=0){ | |
717 | Int_t index2 = ci2[npairs2]; | |
718 | Double_t *xyz2 = &clusterxyzr2[index2*6]; | |
719 | Double_t &r2 = clusterxyzr2[index2*6+3]; | |
720 | // | |
721 | if (TMath::Abs((xyz3[2]-zv)*r2/r3-(xyz2[2]-zv))>0.3) continue; | |
722 | npairs2++; | |
723 | } | |
724 | // | |
725 | // | |
726 | // find pairs in layer 2 | |
727 | for (Int_t jcl2 =0;jcl2<npairs2; jcl2++){ | |
728 | Int_t index2 = ci2[jcl2]; | |
729 | Double_t *xyz2 = &clusterxyzr2[index2*6]; | |
730 | Float_t zc = (xyz2[2]-zv)*ratio12+zv; | |
731 | Float_t lphi = TMath::ATan2(xyz2[1],xyz2[0]); | |
732 | cl2 = layer2.GetCluster(index2); | |
733 | // | |
734 | // | |
735 | layer1.SelectClusters(zc-1,zc+1, layer1.GetR()*lphi-1, layer1.GetR()*lphi+1.); | |
736 | const AliITSclusterV2 * cl1=0; | |
737 | Int_t ci1; | |
738 | // | |
739 | Double_t lr2 = (xyz3[0]-xyz2[0])*(xyz3[0]-xyz2[0])+(xyz3[1]-xyz2[1])*(xyz3[1]-xyz2[1]); | |
740 | Double_t lr = TMath::Sqrt(lr2); | |
741 | Double_t sinphi = (xyz3[1]-xyz2[1])/lr; | |
742 | Double_t cosphi = (xyz3[0]-xyz2[0])/lr; | |
743 | Double_t t1 = (xyz3[2]-xyz2[2])/lr; | |
744 | // | |
745 | while ( (cl1=layer1.GetNextCluster(ci1))!=0){ | |
746 | if (cl1->GetQ()==0) continue; | |
747 | Int_t index1 = ci1; | |
748 | Double_t *xyz1 = &clusterxyzr1[index1*6]; | |
749 | AliITSclusterV2* cl1 = layer1.GetCluster(index1); | |
750 | ||
751 | if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) good0++; | |
752 | Double_t loc1[2]; | |
753 | // | |
754 | loc1[0] = (xyz1[0]-xyz2[0])*cosphi+(xyz1[1]-xyz2[1])*sinphi; | |
755 | loc1[1] = -(xyz1[0]-xyz2[0])*sinphi+(xyz1[1]-xyz2[1])*cosphi; | |
756 | // | |
757 | Double_t c1 = 2.*loc1[1]/(loc1[0]*(loc1[0]-lr)); //curvature approximated using 2 derivation | |
758 | Double_t c=c1; | |
759 | Double_t t2 = (xyz2[2]-xyz1[2])/TMath::Sqrt(loc1[0]*loc1[0]+loc1[1]*loc1[1]); | |
760 | if (TMath::Abs(c)>0.008) continue; | |
761 | if (TMath::Abs(t2-t1)>0.02) continue; | |
762 | // | |
763 | trackclusters[ntracks*4] = 0; | |
764 | trackclusters[ntracks*4+1] = cl1; | |
765 | trackclusters[ntracks*4+2] = cl2; | |
766 | trackclusters[ntracks*4+3] = cl3; | |
767 | trackindex[ntracks*4+0] = 0; | |
768 | trackindex[ntracks*4+1] = index1; | |
769 | trackindex[ntracks*4+2] = index2; | |
770 | trackindex[ntracks*4+3] = icl3; | |
771 | if (ntracks>=300*ncl3){ | |
772 | break; | |
773 | } | |
774 | if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) { | |
775 | //printf("%f\n",c); | |
776 | ||
777 | //hc->Fill(c); | |
778 | //hl->Fill(t2-t1); | |
779 | //hr->Fill(1/c); | |
780 | //hcc->Fill(c,TMath::Abs(c1)); | |
781 | //hcl->Fill(c1,t2-t1); | |
782 | ||
783 | accepted0++; | |
784 | } | |
785 | ntracks++; | |
786 | } | |
787 | } // end - find pairs in layer 2 | |
788 | // | |
789 | ntracks0+=ntracks; | |
790 | //hcc->Draw(); | |
791 | Int_t ntracks2=0; | |
792 | Int_t cindex2[10000*4]; | |
793 | Float_t pools[10000]; | |
794 | ||
795 | for (Int_t itrack=0;itrack<ntracks;itrack++){ | |
796 | const AliITSclusterV2 *cl0, *cl1,*cl2,*cl3; | |
797 | cl3 = trackclusters[itrack*4+3]; | |
798 | cl2 = trackclusters[itrack*4+2]; | |
799 | cl1 = trackclusters[itrack*4+1]; | |
800 | Int_t index1 = trackindex[itrack*4+1]; | |
801 | Int_t index2 = trackindex[itrack*4+2]; | |
802 | Int_t index3 = trackindex[itrack*4+3]; | |
803 | // | |
804 | Double_t *xyz1 = &clusterxyzr1[index1*6],*xyz2 = &clusterxyzr2[index2*6],*xyz3=&clusterxyzr3[index3*6]; | |
805 | Double_t &r1 = clusterxyzr1[index1*6+3],&r2 = clusterxyzr2[index2*6+3], &r3 = clusterxyzr3[index3*6+3]; | |
806 | ||
807 | Double_t r0 = layer0.GetR(); | |
808 | Float_t zc = xyz1[2] + (xyz2[2]-xyz1[2])*(r0-r1)/(r2-r1); | |
809 | Float_t lphi = TMath::ATan2(xyz1[1],xyz1[0]); | |
810 | ||
811 | layer0.SelectClusters(zc-1,zc+1, layer0.GetR()*lphi-0.3, layer0.GetR()*lphi+0.3); | |
812 | Int_t ci0; | |
813 | // | |
814 | // | |
815 | Double_t lr2 = (xyz2[0]-xyz1[0])*(xyz2[0]-xyz1[0])+(xyz2[1]-xyz1[1])*(xyz2[1]-xyz1[1]); | |
816 | Double_t lr = TMath::Sqrt(lr2); | |
817 | Double_t sinphi = (xyz2[1]-xyz1[1])/lr; | |
818 | Double_t cosphi = (xyz2[0]-xyz1[0])/lr; | |
819 | // | |
820 | Double_t loc3[2]; // cluster3 in local coordindate | |
821 | loc3[0] = (xyz3[0]-xyz1[0])*cosphi+(xyz3[1]-xyz1[1])*sinphi; | |
822 | loc3[1] = -(xyz3[0]-xyz1[0])*sinphi+(xyz3[1]-xyz1[1])*cosphi; | |
823 | ||
824 | Double_t tan1 = (xyz3[2]-xyz2[2])/TMath::Sqrt((xyz3[0]-xyz2[0])*(xyz3[0]-xyz2[0])+(xyz3[1]-xyz2[1])*(xyz3[1]-xyz2[1])); | |
825 | Double_t c3 = 2.*loc3[1]/(loc3[0]*(loc3[0]-lr)); //curvature approximated using 2 derivation | |
826 | // | |
827 | while ( (cl0=layer0.GetNextCluster(ci0))!=0){ | |
828 | if (cl0->GetQ()==0) continue; | |
829 | Double_t *xyz0 = &clusterxyzr0[ci0*6]; | |
830 | if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0) | |
831 | && cl3->GetLabel(0)==cl1->GetLabel(0)&&cl3->GetLabel(0)==cl0->GetLabel(0)){ | |
832 | good1++; | |
833 | } | |
834 | Double_t loc0[2]; // cluster3 in local coordindate | |
835 | loc0[0] = (xyz0[0]-xyz1[0])*cosphi+(xyz0[1]-xyz1[1])*sinphi; | |
836 | loc0[1] = -(xyz0[0]-xyz1[0])*sinphi+(xyz0[1]-xyz1[1])*cosphi; | |
837 | // | |
838 | Double_t c0 = 2.*loc0[1]/(loc0[0]*(loc0[0]-lr)); //curvature approximated using 2 derivation | |
839 | Double_t tan2 = (xyz1[2]-xyz0[2])/TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])+(xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1])); | |
840 | //if (TMath::Abs(tan1-tan2)>0.05) continue; | |
841 | //if (TMath::Abs(c3-c0)>0.006) continue; | |
842 | Double_t poolphi = (c3-c0)/0.0011; | |
843 | poolphi*=poolphi; | |
844 | Double_t poolr = ((c3-c0)/(TMath::Abs(c3+c0)+0.001))/0.14; | |
845 | Double_t poolth = (tan1-tan2)/0.0094; | |
846 | poolth*=poolth; | |
847 | Double_t poolthr = ((tan1-tan2)/(TMath::Abs(c3+c0)+0.005))/0.83; | |
848 | if (poolr*poolr+poolthr*poolthr>40) continue; | |
849 | //if ( TMath::Abs((c3-c0)/(c3+c0))>0.6) continue; | |
850 | cindex2[ntracks2*4+0] = ci0; | |
851 | cindex2[ntracks2*4+1] = index1; | |
852 | cindex2[ntracks2*4+2] = index2; | |
853 | cindex2[ntracks2*4+3] = index3; | |
854 | pools[ntracks2] = poolr*poolr+poolthr*poolthr; | |
855 | // | |
856 | ntracks2++; | |
857 | if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0) | |
858 | && cl3->GetLabel(0)==cl1->GetLabel(0)){ | |
859 | hl->Fill(tan1-tan2); | |
860 | hrc->Fill(poolr); | |
861 | hc->Fill((c3-c0)); | |
862 | hlr->Fill(poolthr); | |
863 | hpoolc->Fill(poolr); | |
864 | hpooll->Fill(poolthr); | |
865 | accepted1++; | |
866 | } | |
867 | } | |
868 | } | |
869 | Int_t sindex[10000]; | |
870 | TMath::Sort(ntracks2,pools,sindex,kFALSE); | |
871 | Float_t minpool = pools[sindex[0]]; | |
872 | Int_t max = TMath::Min(ntracks2,20); | |
873 | for (Int_t itrack=0;itrack<ntracks2;itrack++){ | |
874 | if (itrack>=max) continue; | |
875 | //if (pools[sindex[itrack]]>minpool+5) continue; | |
876 | const AliITSclusterV2 *cl0, *cl1,*cl2,*cl3; | |
877 | cl3 = layer3.GetCluster(cindex2[sindex[itrack]*4+3]); | |
878 | cl2 = layer2.GetCluster(cindex2[sindex[itrack]*4+2]); | |
879 | cl1 = layer1.GetCluster(cindex2[sindex[itrack]*4+1]); | |
880 | cl0 = layer0.GetCluster(cindex2[sindex[itrack]*4+0]); | |
881 | if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0) | |
882 | && cl3->GetLabel(0)==cl1->GetLabel(0)){ | |
883 | accepted2++; | |
884 | } | |
885 | nfoundtracks++; | |
886 | } | |
887 | } | |
888 | hc->Draw(); | |
889 | printf("First pass \t%d\t%d\t%d\n",ntracks0,good0,accepted0); | |
890 | printf("Found tracks\t%d\t%d\t%d\t%d\n",nfoundtracks,good1, accepted1,accepted2); | |
891 | ||
892 | } |