]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliITSclusterComparison.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added function returning number in truncated gaussian (user passes mean, sigma and...
[u/mrichter/AliRoot.git] / STEER / AliITSclusterComparison.C
CommitLineData
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
27ClassImp(AliITSCI)
28ClassImp(AliITSClusterErrAnal)
29
30
31
32AliITSClusterErrAnal::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
79void 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
97void 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
162void 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
231void AliITSClusterErrAnal::LoadParticles()
232{
233 //
234 //
235 // fReferences = new TObjArray;
236}
237
238void 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
281void 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
306Int_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
458AliITSclusterV2 * 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
554Int_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
605void 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}