1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 ////////////////////////////////////////////////////////////////////////////
19 // Implementation of the ITS clusterer V2 class //
21 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
22 // Revised: Enrico Fragiacomo, enrico.fragiacomo@ts.infn.it //
24 ///////////////////////////////////////////////////////////////////////////
26 #include <Riostream.h>
29 #include "AliITSClusterFinderV2SSD.h"
30 #include "AliITSRecPoint.h"
31 #include "AliITSgeomTGeo.h"
32 #include "AliITSDetTypeRec.h"
33 #include "AliRawReader.h"
34 #include "AliITSRawStreamSSD.h"
35 #include <TClonesArray.h>
36 #include "AliITSdigitSSD.h"
37 #include "AliITSCalibrationSSD.h"
39 Short_t *AliITSClusterFinderV2SSD::fgPairs = 0x0;
40 Int_t AliITSClusterFinderV2SSD::fgPairsSize = 0;
42 ClassImp(AliITSClusterFinderV2SSD)
45 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinderV2(dettyp),
46 fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1),
57 //______________________________________________________________________
58 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinderV2(cf), fLastSSD1(cf.fLastSSD1),
59 fYpitchSSD(cf.fYpitchSSD),
69 //______________________________________________________________________
70 AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD& cf ){
71 // Assignment operator
73 this->~AliITSClusterFinderV2SSD();
74 new(this) AliITSClusterFinderV2SSD(cf);
79 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
83 FindClustersSSD(fDigits);
87 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
88 //------------------------------------------------------------
89 // Actual SSD cluster finder
90 //------------------------------------------------------------
91 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
94 Int_t smaxall=alldigits->GetEntriesFast();
95 if (smaxall==0) return;
96 // TObjArray *digits = new TObjArray;
98 for (Int_t i=0;i<smaxall; i++){
99 AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
101 if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());
102 else gain = cal->GetGainN(d->GetStripNumber());
104 Float_t q=gain*d->GetSignal(); // calibration brings mip peaks around 120 (in ADC units)
105 q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
106 //Float_t q=d->GetSignal()/4.29;// temp. fix (for PID purposed - normalis. to be checked)
107 d->SetSignal(Int_t(q));
109 if (d->GetSignal()<3) continue;
112 Int_t smax = digits.GetEntriesFast();
115 const Int_t kMax=1000;
117 Ali1Dcluster pos[kMax], neg[kMax];
118 Float_t y=0., q=0., qmax=0.;
119 Int_t lab[4]={-2,-2,-2,-2};
121 AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
123 y += d->GetCoord2()*d->GetSignal();
125 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
126 Int_t curr=d->GetCoord2();
127 Int_t flag=d->GetCoord1();
132 for (Int_t ilab=0;ilab<10;ilab++){
135 milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
137 for (Int_t s=1; s<smax; s++) {
138 d=(AliITSdigitSSD*)digits.UncheckedAt(s);
139 Int_t strip=d->GetCoord2();
140 if ((strip-curr) > 1 || flag!=d->GetCoord1()) {
145 c[*n].SetLabels(milab);
146 //Split suspiciously big cluster
148 c[*n].SetY(y/q-0.25*nd);
152 Error("FindClustersSSD","Too many 1D clusters !");
155 c[*n].SetY(y/q+0.25*nd);
158 c[*n].SetLabels(milab);
162 Error("FindClustersSSD","Too many 1D clusters !");
167 lab[0]=lab[1]=lab[2]=-2;
169 for (Int_t ilab=0;ilab<10;ilab++){
173 if (flag!=d->GetCoord1()) { n=&np; c=pos; }
177 y += d->GetCoord2()*d->GetSignal();
179 if (d->GetSignal()>qmax) {
181 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
183 for (Int_t ilab=0;ilab<10;ilab++) {
184 if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
191 c[*n].SetLabels(lab);
192 //Split suspiciously big cluster
194 c[*n].SetY(y/q-0.25*nd);
198 Error("FindClustersSSD","Too many 1D clusters !");
201 c[*n].SetY(y/q+0.25*nd);
204 c[*n].SetLabels(lab);
208 Error("FindClustersSSD","Too many 1D clusters !");
212 FindClustersSSD(neg, nn, pos, np);
216 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader,TClonesArray** clusters){
218 //------------------------------------------------------------
219 // This function creates ITS clusters from raw data
220 //------------------------------------------------------------
223 const UInt_t *evid; evid = rawReader->GetEventId();
224 cout<<"Event="<<evid[0]<<endl;
226 AliITSRawStreamSSD inputSSD(rawReader);
227 // rawReader->SelectEquipment(-1,0,15);
228 FindClustersSSD(&inputSSD,clusters);
232 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input,
233 TClonesArray** clusters)
235 //------------------------------------------------------------
236 // Actual SSD cluster finder for raw data
237 //------------------------------------------------------------
238 Int_t nClustersSSD = 0;
239 const Int_t kMax = 1000;
240 Ali1Dcluster clusters1D[2][kMax];
241 Int_t nClusters[2] = {0, 0};
242 Int_t lab[3]={-2,-2,-2};
250 AliITSCalibrationSSD* cal=NULL;
252 Int_t matrix[12][1536];
259 Int_t osignal = 65535;
263 // read raw data input stream
266 // reset signal matrix
267 for(Int_t i=0; i<12; i++) { for(Int_t j=0; j<1536; j++) { matrix[i][j] = 65535;} }
271 matrix[oadc][ostrip] = osignal; // recover data from previous occurence of input->Next()
274 // buffer data for ddl=iddl and ad=iad
277 next = input->Next();
278 if((!next)&&(input->flag)) continue;
279 Int_t ddl=input->GetDDL();
280 Int_t ad=input->GetAD();
281 Int_t adc = input->GetADC(); adc = (adc<6)? adc : adc - 2;
282 Int_t strip = input->GetStrip();
283 if(input->GetSideFlag()) strip=1535-strip;
284 Int_t signal = input->GetSignal();
285 //cout<<ddl<<" "<<ad<<" "<<adc<<" "<<strip<<" "<<signal<<endl;
287 if((ddl==iddl)&&(ad==iad)) {n++; matrix[adc][strip] = signal;}
288 else {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
290 if(!next) {oddl=iddl; oad=iad; oadc = adc; ostrip = strip; osignal=signal; iddl=ddl; iad=ad; break;}
295 if(!next && oddl<0) break;
297 if(n==0) continue; // first occurence
301 for(Int_t iadc=0; iadc<12; iadc++) { // loop over ADC index for ddl=oddl and ad=oad
303 Int_t iimod = (oad - 1) * 12 + iadc;
304 Int_t iModule = AliITSRawStreamSSD::GetModuleNumber(oddl,iimod);
305 if(iModule==-1) continue;
306 // cout<<"ddl="<<oddl<<" ad"<<oad<<" module="<<iModule<<endl;
307 cal = (AliITSCalibrationSSD*)GetResp(iModule);
311 for(Int_t istrip=0; istrip<768; istrip++) { // P-side
312 Int_t signal = matrix[iadc][istrip];
313 pedestal = cal->GetPedestalP(istrip);
314 matrix[iadc][istrip]=signal-(Int_t)pedestal;
318 for(Int_t l=0; l<6; l++) {
320 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
322 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
326 for(Int_t istrip=0; istrip<768; istrip++) { // P-side
328 Int_t signal = TMath::Abs(matrix[iadc][istrip]);
331 noise = cal->GetNoiseP(istrip); if(noise<1.) signal = 65535;
332 if(signal<5*noise) signal = 65535; // in case ZS was not done in hw do it now
333 if( (signal<30.) || (istrip<10) || (istrip>758) ) signal=65535;
336 gain = cal->GetGainP(istrip);
337 signal = (Int_t) ( signal * gain ); // signal is corrected for gain
338 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
340 q += signal; // add digit to current cluster
341 y += istrip * signal;
348 if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
350 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
353 cluster.SetNd(nDigits);
354 cluster.SetLabels(lab);
356 //Split suspiciously big cluster
357 if (nDigits > 4&&nDigits < 25) {
358 cluster.SetY(y/q - 0.25*nDigits);
360 if (nClusters[0] == kMax) {
361 Error("FindClustersSSD", "Too many 1D clusters !");
364 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
365 cluster2.SetY(y/q + 0.25*nDigits);
366 cluster2.SetQ(0.5*q);
367 cluster2.SetNd(nDigits);
368 cluster2.SetLabels(lab);
378 } // loop over strip on P-side
380 // if last strip does have signal
383 if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
385 Ali1Dcluster& cluster = clusters1D[0][nClusters[0]++];
388 cluster.SetNd(nDigits);
389 cluster.SetLabels(lab);
391 //Split suspiciously big cluster
392 if (nDigits > 4&&nDigits < 25) {
393 cluster.SetY(y/q - 0.25*nDigits);
395 if (nClusters[0] == kMax) {
396 Error("FindClustersSSD", "Too many 1D clusters !");
399 Ali1Dcluster& cluster2 = clusters1D[0][nClusters[0]++];
400 cluster2.SetY(y/q + 0.25*nDigits);
401 cluster2.SetQ(0.5*q);
402 cluster2.SetNd(nDigits);
403 cluster2.SetLabels(lab);
412 for(Int_t istrip=768; istrip<1536; istrip++) { // P-side
413 Int_t signal = matrix[iadc][istrip];
414 pedestal = cal->GetPedestalN(1535-istrip);
415 matrix[iadc][istrip]=signal-(Int_t)pedestal;
418 for(Int_t l=6; l<12; l++) {
420 for(Int_t n=20; n<108; n++) cmode+=matrix[iadc][l*128+n];
422 for(Int_t n=0; n<128; n++) matrix[iadc][l*128+n]-=(Int_t)cmode;
427 for(Int_t istrip=768; istrip<1536; istrip++) { // N-side
429 Int_t signal = TMath::Abs(matrix[iadc][istrip]);
430 //cout<<"####"<<" "<<oddl<<" "<<oad<<" "<<iadc<<" "<<istrip<<" "<<signal<<endl;
432 Int_t strip = 1535-istrip;
435 noise = cal->GetNoiseN(strip); if(noise<1.) signal=65535;
436 if(signal<5*noise) signal = 65535; // in case ZS was not done in hw do it now
437 if( (signal<30.) || (istrip<778) || (istrip>1526) ) signal=65535;
440 // cout<<"ddl="<<oddl<<" ad"<<oad<<" module="<<iModule<<" strip= "<<istrip<<
441 // " sig="<<signal<<" "<<cal->GetPedestalN(strip)<<endl;
442 gain = cal->GetGainN(strip);
443 signal = (Int_t) ( signal * gain); // signal is corrected for gain
444 signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
446 // add digit to current cluster
455 if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
457 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
460 cluster.SetNd(nDigits);
461 cluster.SetLabels(lab);
463 //Split suspiciously big cluster
464 if (nDigits > 4&&nDigits < 25) {
465 cluster.SetY(y/q - 0.25*nDigits);
467 if (nClusters[1] == kMax) {
468 Error("FindClustersSSD", "Too many 1D clusters !");
471 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
472 cluster2.SetY(y/q + 0.25*nDigits);
473 cluster2.SetQ(0.5*q);
474 cluster2.SetNd(nDigits);
475 cluster2.SetLabels(lab);
484 } // loop over strips on N-side
488 if ( ((nDigits==1)&&(q>5*oldnoise)) || (nDigits>1) ) {
490 Ali1Dcluster& cluster = clusters1D[1][nClusters[1]++];
493 cluster.SetNd(nDigits);
494 cluster.SetLabels(lab);
496 //Split suspiciously big cluster
497 if (nDigits > 4&&nDigits < 25) {
498 cluster.SetY(y/q - 0.25*nDigits);
500 if (nClusters[1] == kMax) {
501 Error("FindClustersSSD", "Too many 1D clusters !");
504 Ali1Dcluster& cluster2 = clusters1D[1][nClusters[1]++];
505 cluster2.SetY(y/q + 0.25*nDigits);
506 cluster2.SetQ(0.5*q);
507 cluster2.SetNd(nDigits);
508 cluster2.SetLabels(lab);
518 if((nClusters[0])&&(nClusters[1])) {
520 //cout<<"creating recpoint for module="<<iModule<<" "<<nClusters[0]<<" "<<nClusters[1]<<endl;
521 clusters[iModule] = new TClonesArray("AliITSRecPoint");
524 FindClustersSSD(&clusters1D[0][0], nClusters[0],
525 &clusters1D[1][0], nClusters[1], clusters[iModule]);
526 Int_t nClusters = clusters[iModule]->GetEntriesFast();
527 nClustersSSD += nClusters;
530 nClusters[0] = nClusters[1] = 0;
539 Info("FindClustersSSD", "found clusters in ITS SSD: %d", nClustersSSD);
542 void AliITSClusterFinderV2SSD::
543 FindClustersSSD(Ali1Dcluster* neg, Int_t nn,
544 Ali1Dcluster* pos, Int_t np,
545 TClonesArray *clusters) {
546 //------------------------------------------------------------
547 // Actual SSD cluster finder
548 //------------------------------------------------------------
552 const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
554 TClonesArray &cl=*clusters;
556 Float_t tanp=fTanP, tann=fTanN;
557 if (fModule>fLastSSD1) {tann=fTanP; tanp=fTanN;}
558 Int_t idet=fNdet[fModule];
561 Int_t negativepair[30000];
562 Int_t cnegative[3000];
564 Int_t positivepair[30000];
565 Int_t cpositive[3000];
567 for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
568 for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
569 for (Int_t i=0;i<30000;i++) {negativepair[i]=0; positivepair[i]=0;}
571 if ((np*nn) > fgPairsSize) {
572 if (fgPairs) delete [] fgPairs;
573 fgPairsSize = 4*np*nn;
574 fgPairs = new Short_t[fgPairsSize];
576 memset(fgPairs,0,sizeof(Short_t)*np*nn);
579 // find available pairs
581 for (Int_t i=0; i<np; i++) {
582 Float_t yp=pos[i].GetY()*fYpitchSSD;
583 if (pos[i].GetQ()<3) continue;
584 for (Int_t j=0; j<nn; j++) {
585 if (neg[j].GetQ()<3) continue;
586 Float_t yn=neg[j].GetY()*fYpitchSSD;
587 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
588 Float_t yt=yn + tann*zt;
589 zt-=fHlSSD; yt-=fHwSSD;
590 if (TMath::Abs(yt)<fHwSSD+0.01)
591 if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
592 negativepair[i*10+cnegative[i]] =j; //index
593 positivepair[j*10+cpositive[j]] =i;
594 cnegative[i]++; //counters
602 // try to recover points out of but close to the module boundaries
604 for (Int_t i=0; i<np; i++) {
605 Float_t yp=pos[i].GetY()*fYpitchSSD;
606 if (pos[i].GetQ()<3) continue;
607 for (Int_t j=0; j<nn; j++) {
608 if (neg[j].GetQ()<3) continue;
609 // if both 1Dclusters have an other cross continue
610 if (cpositive[j]&&cnegative[i]) continue;
611 Float_t yn=neg[j].GetY()*fYpitchSSD;
612 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
613 Float_t yt=yn + tann*zt;
614 zt-=fHlSSD; yt-=fHwSSD;
615 if (TMath::Abs(yt)<fHwSSD+0.1)
616 if (TMath::Abs(zt)<fHlSSD+0.15) {
617 // tag 1Dcluster (eventually will produce low quality recpoint)
618 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
619 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
620 negativepair[i*10+cnegative[i]] =j; //index
621 positivepair[j*10+cpositive[j]] =i;
622 cnegative[i]++; //counters
637 for (Int_t ip=0;ip<np;ip++){
638 Float_t ybest=1000,zbest=1000,qbest=0;
640 // select gold clusters
641 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
642 Float_t yp=pos[ip].GetY()*fYpitchSSD;
643 Int_t j = negativepair[10*ip];
644 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
646 Float_t yn=neg[j].GetY()*fYpitchSSD;
647 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
648 Float_t yt=yn + tann*zt;
649 zt-=fHlSSD; yt-=fHwSSD;
651 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
653 //cout<<yt<<" "<<zt<<" "<<qbest<<endl;
656 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
657 mT2L->MasterToLocal(loc,trk);
661 lp[2]=0.0025*0.0025; //SigmaY2
662 lp[3]=0.110*0.110; //SigmaZ2
665 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
666 for (Int_t ilab=0;ilab<3;ilab++){
667 milab[ilab] = pos[ip].GetLabel(ilab);
668 milab[ilab+3] = neg[j].GetLabel(ilab);
672 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
673 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
674 AliITSRecPoint * cl2;
676 if(clusters){ // Note clusters != 0 when method is called for rawdata
679 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
681 // cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
683 cl2->SetChargeRatio(ratio);
686 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
694 else{ // Note clusters == 0 when method is called for digits
696 cl2 = new AliITSRecPoint(milab,lp,info);
698 // cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
700 cl2->SetChargeRatio(ratio);
703 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
709 //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" gold"<<endl;
710 fDetTypeRec->AddRecPoint(*cl2);
716 for (Int_t ip=0;ip<np;ip++){
717 Float_t ybest=1000,zbest=1000,qbest=0;
720 // select "silber" cluster
721 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
722 Int_t in = negativepair[10*ip];
723 Int_t ip2 = positivepair[10*in];
724 if (ip2==ip) ip2 = positivepair[10*in+1];
725 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
726 if (TMath::Abs(pcharge-neg[in].GetQ())<10){
729 if (fgPairs[ip*nn+in]==100){ //
730 Float_t yp=pos[ip].GetY()*fYpitchSSD;
731 Float_t yn=neg[in].GetY()*fYpitchSSD;
732 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
733 Float_t yt=yn + tann*zt;
734 zt-=fHlSSD; yt-=fHwSSD;
736 qbest =pos[ip].GetQ();
738 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
739 mT2L->MasterToLocal(loc,trk);
743 lp[2]=0.0025*0.0025; //SigmaY2
744 lp[3]=0.110*0.110; //SigmaZ2
747 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
748 for (Int_t ilab=0;ilab<3;ilab++){
749 milab[ilab] = pos[ip].GetLabel(ilab);
750 milab[ilab+3] = neg[in].GetLabel(ilab);
754 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
755 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
756 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
758 AliITSRecPoint * cl2;
761 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
763 // cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
765 cl2->SetChargeRatio(ratio);
767 fgPairs[ip*nn+in] = 5;
768 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
770 fgPairs[ip*nn+in] = 6;
774 cl2 = new AliITSRecPoint(milab,lp,info);
775 cl2->SetChargeRatio(ratio);
777 fgPairs[ip*nn+in] = 5;
778 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
780 fgPairs[ip*nn+in] = 6;
782 //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silver1"<<endl;
784 fDetTypeRec->AddRecPoint(*cl2);
792 // if (!(cused1[ip2] || cused2[in])){ //
793 if (fgPairs[ip2*nn+in]==100){
794 Float_t yp=pos[ip2].GetY()*fYpitchSSD;
795 Float_t yn=neg[in].GetY()*fYpitchSSD;
796 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
797 Float_t yt=yn + tann*zt;
798 zt-=fHlSSD; yt-=fHwSSD;
800 qbest =pos[ip2].GetQ();
802 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
803 mT2L->MasterToLocal(loc,trk);
807 lp[2]=0.0025*0.0025; //SigmaY2
808 lp[3]=0.110*0.110; //SigmaZ2
811 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
812 for (Int_t ilab=0;ilab<3;ilab++){
813 milab[ilab] = pos[ip2].GetLabel(ilab);
814 milab[ilab+3] = neg[in].GetLabel(ilab);
818 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
819 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
820 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
822 AliITSRecPoint * cl2;
824 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
826 // cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
828 cl2->SetChargeRatio(ratio);
830 fgPairs[ip2*nn+in] =5;
831 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
833 fgPairs[ip2*nn+in] =6;
837 cl2 = new AliITSRecPoint(milab,lp,info);
838 cl2->SetChargeRatio(ratio);
840 fgPairs[ip2*nn+in] =5;
841 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
843 fgPairs[ip2*nn+in] =6;
846 // cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silver2"<<endl;
847 fDetTypeRec->AddRecPoint(*cl2);
859 for (Int_t jn=0;jn<nn;jn++){
860 if (cused2[jn]) continue;
861 Float_t ybest=1000,zbest=1000,qbest=0;
862 // select "silber" cluster
863 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
864 Int_t ip = positivepair[10*jn];
865 Int_t jn2 = negativepair[10*ip];
866 if (jn2==jn) jn2 = negativepair[10*ip+1];
867 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
869 if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
872 // if (!(cused1[ip]||cused2[jn])){
873 if (fgPairs[ip*nn+jn]==100){
874 Float_t yn=neg[jn].GetY()*fYpitchSSD;
875 Float_t yp=pos[ip].GetY()*fYpitchSSD;
876 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
877 Float_t yt=yn + tann*zt;
878 zt-=fHlSSD; yt-=fHwSSD;
880 qbest =neg[jn].GetQ();
882 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
883 mT2L->MasterToLocal(loc,trk);
887 lp[2]=0.0025*0.0025; //SigmaY2
888 lp[3]=0.110*0.110; //SigmaZ2
891 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
892 for (Int_t ilab=0;ilab<3;ilab++){
893 milab[ilab] = pos[ip].GetLabel(ilab);
894 milab[ilab+3] = neg[jn].GetLabel(ilab);
898 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
899 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
900 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
902 AliITSRecPoint * cl2;
904 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
906 // cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
908 cl2->SetChargeRatio(ratio);
910 fgPairs[ip*nn+jn] =7;
911 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
918 cl2 = new AliITSRecPoint(milab,lp,info);
919 cl2->SetChargeRatio(ratio);
921 fgPairs[ip*nn+jn] =7;
922 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
926 //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silverN1"<<endl;
928 fDetTypeRec->AddRecPoint(*cl2);
934 // if (!(cused1[ip]||cused2[jn2])){
935 if (fgPairs[ip*nn+jn2]==100){
936 Float_t yn=neg[jn2].GetY()*fYpitchSSD;
937 Double_t yp=pos[ip].GetY()*fYpitchSSD;
938 Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
939 Double_t yt=yn + tann*zt;
940 zt-=fHlSSD; yt-=fHwSSD;
942 qbest =neg[jn2].GetQ();
944 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
945 mT2L->MasterToLocal(loc,trk);
949 lp[2]=0.0025*0.0025; //SigmaY2
950 lp[3]=0.110*0.110; //SigmaZ2
953 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
954 for (Int_t ilab=0;ilab<3;ilab++){
955 milab[ilab] = pos[ip].GetLabel(ilab);
956 milab[ilab+3] = neg[jn2].GetLabel(ilab);
960 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
961 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
962 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
963 AliITSRecPoint * cl2;
965 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
967 // cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
969 cl2->SetChargeRatio(ratio);
970 fgPairs[ip*nn+jn2]=7;
972 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
974 fgPairs[ip*nn+jn2]=8;
979 cl2 = new AliITSRecPoint(milab,lp,info);
980 cl2->SetChargeRatio(ratio);
981 fgPairs[ip*nn+jn2]=7;
983 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
985 fgPairs[ip*nn+jn2]=8;
987 //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" silverN2"<<endl;
989 fDetTypeRec->AddRecPoint(*cl2);
1001 for (Int_t ip=0;ip<np;ip++){
1002 Float_t ybest=1000,zbest=1000,qbest=0;
1006 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1007 Float_t minchargediff =4.;
1009 for (Int_t di=0;di<cnegative[ip];di++){
1010 Int_t jc = negativepair[ip*10+di];
1011 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1012 if (TMath::Abs(chargedif)<minchargediff){
1014 minchargediff = TMath::Abs(chargedif);
1017 if (j<0) continue; // not proper cluster
1020 for (Int_t di=0;di<cnegative[ip];di++){
1021 Int_t jc = negativepair[ip*10+di];
1022 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1023 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1025 if (count>1) continue; // more than one "proper" cluster for positive
1028 for (Int_t dj=0;dj<cpositive[j];dj++){
1029 Int_t ic = positivepair[j*10+dj];
1030 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1031 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1033 if (count>1) continue; // more than one "proper" cluster for negative
1038 for (Int_t dj=0;dj<cnegative[jp];dj++){
1039 Int_t ic = positivepair[jp*10+dj];
1040 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1041 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1043 if (count>1) continue;
1044 if (fgPairs[ip*nn+j]<100) continue;
1046 //almost gold clusters
1047 Float_t yp=pos[ip].GetY()*fYpitchSSD;
1048 Float_t yn=neg[j].GetY()*fYpitchSSD;
1049 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1050 Float_t yt=yn + tann*zt;
1051 zt-=fHlSSD; yt-=fHwSSD;
1053 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1055 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1056 mT2L->MasterToLocal(loc,trk);
1060 lp[2]=0.0025*0.0025; //SigmaY2
1061 lp[3]=0.110*0.110; //SigmaZ2
1063 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1064 for (Int_t ilab=0;ilab<3;ilab++){
1065 milab[ilab] = pos[ip].GetLabel(ilab);
1066 milab[ilab+3] = neg[j].GetLabel(ilab);
1069 CheckLabels2(milab);
1070 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1071 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1072 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1073 AliITSRecPoint * cl2;
1075 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1077 // cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
1079 cl2->SetChargeRatio(ratio);
1081 fgPairs[ip*nn+j]=10;
1082 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1084 fgPairs[ip*nn+j]=11;
1090 cl2 = new AliITSRecPoint(milab,lp,info);
1091 cl2->SetChargeRatio(ratio);
1093 fgPairs[ip*nn+j]=10;
1094 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1096 fgPairs[ip*nn+j]=11;
1101 //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" 2x2"<<endl;
1103 fDetTypeRec->AddRecPoint(*cl2);
1111 for (Int_t i=0; i<np; i++) {
1112 Float_t ybest=1000,zbest=1000,qbest=0;
1113 Float_t yp=pos[i].GetY()*fYpitchSSD;
1114 if (pos[i].GetQ()<3) continue;
1115 for (Int_t j=0; j<nn; j++) {
1116 // for (Int_t di = 0;di<cpositive[i];di++){
1117 // Int_t j = negativepair[10*i+di];
1118 if (neg[j].GetQ()<3) continue;
1119 if (cused2[j]||cused1[i]) continue;
1120 if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1121 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1122 Float_t yn=neg[j].GetY()*fYpitchSSD;
1123 Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
1124 Float_t yt=yn + tann*zt;
1125 zt-=fHlSSD; yt-=fHwSSD;
1126 if (TMath::Abs(yt)<fHwSSD+0.01)
1127 if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
1129 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1131 Double_t loc[3]={ybest,0.,zbest},trk[3]={0.,0.,0.};
1132 mT2L->MasterToLocal(loc,trk);
1136 lp[2]=0.0025*0.0025; //SigmaY2
1137 lp[3]=0.110*0.110; //SigmaZ2
1140 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1141 for (Int_t ilab=0;ilab<3;ilab++){
1142 milab[ilab] = pos[i].GetLabel(ilab);
1143 milab[ilab+3] = neg[j].GetLabel(ilab);
1146 CheckLabels2(milab);
1147 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1148 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1149 AliITSRecPoint * cl2;
1151 cl2 = new (cl[ncl]) AliITSRecPoint(milab,lp,info);
1153 // cl2-> GetGlobalXYZ(xyz); cout<<"rec "<<xyz[0]<<" "<<xyz[1]<<" "<<xyz[2]<<endl;
1155 cl2->SetChargeRatio(ratio);
1156 cl2->SetType(100+cpositive[j]+cnegative[i]);
1159 cl2 = new AliITSRecPoint(milab,lp,info);
1160 cl2->SetChargeRatio(ratio);
1161 cl2->SetType(100+cpositive[j]+cnegative[i]);
1163 //cout<<"AliITSClusterFinderV2SSD "<<fModule<<" other"<<endl;
1165 fDetTypeRec->AddRecPoint(*cl2);
1170 if (fgPairs[i*nn+j]<100){
1171 printf("problem:- %d\n", fgPairs[i*nn+j]);
1173 if (cnegative[i]<2&&cpositive[j]<2){
1174 printf("problem:- %d\n", fgPairs[i*nn+j]);