1 #if !defined(__CINT__) || defined(__MAKECINT__)
4 #include "AliITStrackerMI.h"
5 #include "AliITSclusterV2.h"
9 #include "AliITSgeom.h"
10 //#include "AliITSgeometryDetail.h"
11 //#include "AliITSparameter.h"
14 #include "TStopwatch.h"
19 #include "AliRunLoader.h"
22 #include "AliTrackReference.h"
24 #include "AliITSclusterComparison.h"
28 ClassImp(AliITSClusterErrAnal)
32 AliITSClusterErrAnal::AliITSClusterErrAnal(Char_t *chloader )
37 delete gAlice->GetRunLoader();
41 fRunLoader = AliRunLoader::Open(chloader);
42 if (fRunLoader == 0x0){
43 cerr<<"Can not open session"<<endl;
46 fITSLoader = fRunLoader->GetLoader("ITSLoader");
47 if (fITSLoader == 0x0){
48 cerr<<"Can not get ITS Loader"<<endl;
51 if (fRunLoader->LoadgAlice()){
52 cerr<<"Error occured while l"<<endl;
55 fRunLoader->CdGAFile();
56 fITS = (AliITS*)gAlice->GetDetector("ITS");
58 cerr<<"AliITSclusterComparison.C : Can not find the ITS detector !"<<endl;
60 AliITSgeom *geom = fITS->GetITSgeom();
61 //An instance of the ITS tracker
62 fTracker = new AliITStrackerMI(geom);
65 AliITSCI * clinfo = new AliITSCI();
66 fFileA = new TFile("itsclusteranal.root","recreate");
68 fTreeA = new TTree("itscl","itscl");
69 fTreeA->Branch("itscl","AliITSCI",&clinfo);
71 AliITSclusterV2 * cl = new AliITSclusterV2;
72 fTreeB = new TTree("Clusters","Clusters");
73 fTreeB->Branch("cl","AliITSclusterV2",&cl);
79 void AliITSClusterErrAnal::SetIO(Int_t event)
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");
90 fStack = fRunLoader->Stack();
91 fHitTree = fITSLoader->TreeH();
92 fClusterTree = fITSLoader->TreeR();
93 fReferenceTree = fRunLoader->TreeTR();
97 void AliITSClusterErrAnal::LoadClusters()
106 Int_t nparticles = fStack->GetNtrack();
107 fClusters = new TObjArray(nparticles);
109 TObjArray *ClusterArray = new TClonesArray("AliITSclusterV2");
110 TObjArray carray(2000);
111 TBranch *branch=fClusterTree->GetBranch("Clusters");
113 Error("ReadClusters","Can't get the branch !");
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));
132 Int_t nClusters = carray.GetEntriesFast();
133 printf("Total number of clusters %d \n", nClusters);
134 fTracker->LoadClusters(fClusterTree);
141 for (Int_t i=0;i<nClusters;i++){
142 AliITSclusterV2 *cl = (AliITSclusterV2 *) carray.UncheckedAt(i);
143 if (cl->GetLabel(0)<0) noisecl++;
145 for (Int_t itrack=0; itrack<3;itrack++){
146 Int_t lab = cl->GetLabel(itrack);
148 TObjArray * array = (TObjArray*)fClusters->At(lab);
150 array = new TObjArray(20);
151 fClusters->AddAt(array,lab);
158 printf("Total number of track clusters %d \n", all);
159 printf("Total number of noise clusters %d \n", noisecl);
162 void AliITSClusterErrAnal::SignDeltas( TObjArray *ClusterArray, Float_t vz)
166 Int_t entries = ClusterArray->GetEntriesFast();
167 if (entries<4) return;
168 AliITSclusterV2* cluster = (AliITSclusterV2*)ClusterArray->At(0);
169 Int_t layer = cluster->GetLayer();
173 Float_t r = (layer>0)? 7:4;
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
181 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
182 cl0->SetDeltaProbability(1);
188 for (Int_t i=0;i<ncandidates;i++){
189 AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(index[i]);
190 if (cl0->GetDeltaProbability()>0.8) continue;
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++){
197 AliITSclusterV2* cl1 = (AliITSclusterV2*)ClusterArray->At(index[j]);
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;
214 Float_t det = sumw*sumy2 - sumy*sumy;
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()));
222 Float_t z0 = sumyz/sumy;
223 delta = TMath::Abs(cl0->GetZ()-z0);
226 cl0->SetDeltaProbability(1-20.*delta);
231 void AliITSClusterErrAnal::LoadParticles()
235 // fReferences = new TObjArray;
238 void AliITSClusterErrAnal::SortReferences()
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);
252 TClonesArray *labarr=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);
260 Int_t lab = ref->GetTrack();
261 //if ( (lab<0) || (lab>ntracks)) continue;
263 if (fReferences->At(lab)==0) {
264 labarr = new TClonesArray("AliTrackReference");
265 fReferences->AddAt(labarr,lab);
268 TClonesArray &larr = *labarr;
269 new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
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");
281 void AliITSClusterErrAnal::GetNTeor(Int_t layer, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz)
286 ny = 1.+TMath::Abs(phi)*3.2;
287 nz = 1.+TMath::Abs(theta)*0.34;
291 ny = 1.+TMath::Abs(phi)*3.2;
292 nz = 1.+TMath::Abs(theta)*0.28;
297 ny = 2.02+TMath::Abs(phi)*1.95;
298 nz = 2.02+TMath::Abs(phi)*2.35;
301 ny = 6.6-2.7*abs(phi);
302 nz = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta);
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)
308 //calculate cluster position error
311 GetNTeor(layer, theta,phi,ny,nz);
312 erry = TMath::Sqrt(cl->GetSigmaY2());
313 errz = TMath::Sqrt(cl->GetSigmaZ2());
318 if (TMath::Abs(ny-cl->GetNy())>0.6) {
320 erry*=0.4+TMath::Abs(ny-cl->GetNy());
321 errz*=0.4+TMath::Abs(ny-cl->GetNy());
323 erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
324 errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
327 if (TMath::Abs(nz-cl->GetNz())>1.) {
328 erry*=TMath::Abs(nz-cl->GetNz());
329 errz*=TMath::Abs(nz-cl->GetNz());
333 erry= TMath::Min(erry,float(0.005));
334 errz= TMath::Min(errz,float(0.03));
339 if (cl->GetNy()==100||cl->GetNz()==100){
344 if (cl->GetNy()+cl->GetNz()>12){
349 Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
350 Float_t chargematch = TMath::Max(double(normq/expQ),2.);
352 if (cl->GetType()==1 || cl->GetType()==10 ){
353 if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
358 if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
364 errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15);
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);
373 if (cl->GetType()>100 ){
374 if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
379 if (cl->GetNy()+cl->GetNz()-nz-ny<1){
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);
389 Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
393 if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
404 Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
405 Float_t chargematch = normq/expQ;
407 Int_t cnz = cl->GetNz()%10;
409 if (cl->GetType()==1){
410 if (chargematch<1.25){
411 erry = 0.0028*(1.+6./cl->GetQ()); // gold clusters
414 erry = 0.003*chargematch;
415 if (cl->GetNz()==3) erry*=1.5;
417 if (chargematch<1.0){
418 errz = 0.0011*(1.+6./cl->GetQ());
421 errz = 0.002*(1+2*(chargematch-1.));
425 errz*=1.4*(cnz-nz+0.5);
428 if (cl->GetType()>1){
430 erry = 0.00385*(1.+6./cl->GetQ()); // gold clusters
431 errz = 0.0016*(1.+6./cl->GetQ());
434 errz = 0.0014*(1+3*(chargematch-1.));
435 erry = 0.003*(1+3*(chargematch-1.));
439 errz*=1.4*(cnz-nz+0.5);
443 if (TMath::Abs(cl->GetY())>2.5){
444 factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
446 if (TMath::Abs(cl->GetY())<1){
447 factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
449 factorz= TMath::Min(factorz,float(4.));
452 erry= TMath::Min(erry,float(0.05));
453 errz= TMath::Min(errz,float(0.05));
458 AliITSclusterV2 * AliITSClusterErrAnal::FindNearestCluster(AliITSCI * clinfo, Float_t dmax)
462 AliTrackReference * ref = &(clinfo->fRef);
463 TObjArray * clusters = (TObjArray*)(fClusters->At(ref->GetTrack()));
464 if (!clusters) return 0;
465 Int_t nclusters = clusters->GetEntriesFast();
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);
473 Float_t dz = cl->GetZ()-ref->Z();
475 if (TMath::Abs(dz)>dmax) continue;
476 Int_t idet=cl->GetDetectorIndex();
479 // for (Int_t i=cl->fLayer;i<=cl->fLayer-1;i++){
481 AliITStrackerMI::AliITSdetector & det = fTracker->GetDetector(cl->fLayer,idet);
482 if (TMath::Abs(radius-det.GetR())<1.5) {
484 for (Int_t itrack=0;itrack<3;itrack++){
485 if (cl->fTracks[itrack]== clinfo->fRef.GetTrack()) {
486 clinfo->fNClusters++;
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;
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);
503 Double_t theta = ref->Pz()/ref->Pt();
505 dy = y -cl->GetY()-dx*lphi;
506 dz = ref->Z()-cl->GetZ()-dx*theta;
508 Double_t d = TMath::Sqrt(dy*dy+dz*dz);
512 clinfo->Update(fTracker);
513 clinfo->fCl = *cluster;
516 clinfo->fLayer = cl->fLayer;
517 clinfo->fPhiCl = phicl;
518 clinfo->fRCl = det.GetR();
519 clinfo->fPhiL = lphi;
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));
535 clinfo->fErrY = erry;
536 clinfo->fErrZ = errz;
537 clinfo->fPoolY2 = dy/clinfo->fErrY;
538 clinfo->fPoolZ2 = dz/clinfo->fErrZ;
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;
554 Int_t AliITSClusterErrAnal::Analyze(Int_t trackmax) {
557 // dummy cluster to be fill if not cluster info
559 AliITSCI * clinfo = new AliITSCI;
560 AliITSclusterV2 dummy;
561 TBranch * branch = fTreeA->GetBranch("itscl");
562 branch->SetAddress(&clinfo);
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;
580 clinfo->Update(fTracker);
582 AliITSclusterV2 * cluster = FindNearestCluster(clinfo,3.0);
584 //clinfo->fCl = *cluster;
590 clinfo->Update(fTracker);
605 void AliITSClusterErrAnal::MakeSeeds(Double_t zv)
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.);
627 TH1F *hpoolc = new TH1F("hpoolc","hpoolc",100,-5,5);
628 TH1F *hpooll = new TH1F("hpooll","hpooll",100,-5,5);
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;
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];
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);
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);
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);
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);
694 Int_t nfoundtracks=0;
699 Int_t *trackindex = new Int_t[ncl3*4*300];
700 const AliITSclusterV2 **trackclusters = new const AliITSclusterV2*[ncl3*4*300];
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
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;
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];
721 if (TMath::Abs((xyz3[2]-zv)*r2/r3-(xyz2[2]-zv))>0.3) continue;
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);
735 layer1.SelectClusters(zc-1,zc+1, layer1.GetR()*lphi-1, layer1.GetR()*lphi+1.);
736 const AliITSclusterV2 * cl1=0;
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;
745 while ( (cl1=layer1.GetNextCluster(ci1))!=0){
746 if (cl1->GetQ()==0) continue;
748 Double_t *xyz1 = &clusterxyzr1[index1*6];
749 AliITSclusterV2* cl1 = layer1.GetCluster(index1);
751 if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) good0++;
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;
757 Double_t c1 = 2.*loc1[1]/(loc1[0]*(loc1[0]-lr)); //curvature approximated using 2 derivation
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;
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){
774 if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) {
780 //hcc->Fill(c,TMath::Abs(c1));
781 //hcl->Fill(c1,t2-t1);
787 } // end - find pairs in layer 2
792 Int_t cindex2[10000*4];
793 Float_t pools[10000];
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];
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];
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]);
811 layer0.SelectClusters(zc-1,zc+1, layer0.GetR()*lphi-0.3, layer0.GetR()*lphi+0.3);
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;
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;
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
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)){
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;
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;
844 Double_t poolr = ((c3-c0)/(TMath::Abs(c3+c0)+0.001))/0.14;
845 Double_t poolth = (tan1-tan2)/0.0094;
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;
857 if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0)
858 && cl3->GetLabel(0)==cl1->GetLabel(0)){
864 hpooll->Fill(poolthr);
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)){
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);