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 // Last revision: 13-05-09 Enrico Fragiacomo //
23 // enrico.fragiacomo@ts.infn.it //
25 ///////////////////////////////////////////////////////////////////////////
27 #include "AliITSClusterFinderV2SSD.h"
29 #include <Riostream.h>
30 #include <TGeoGlobalMagField.h>
34 #include "AliITSRecPoint.h"
35 #include "AliITSRecPointContainer.h"
36 #include "AliITSgeomTGeo.h"
37 #include "AliRawReader.h"
38 #include "AliITSRawStreamSSD.h"
39 #include <TClonesArray.h>
40 #include <TCollection.h>
41 #include "AliITSdigitSSD.h"
42 #include "AliITSReconstructor.h"
43 #include "AliITSCalibrationSSD.h"
44 #include "AliITSsegmentationSSD.h"
46 Short_t *AliITSClusterFinderV2SSD::fgPairs = 0x0;
47 Int_t AliITSClusterFinderV2SSD::fgPairsSize = 0;
48 const Float_t AliITSClusterFinderV2SSD::fgkThreshold = 5.;
50 const Float_t AliITSClusterFinderV2SSD::fgkCosmic2008StripShifts[16][9] =
51 {{-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 512
52 {-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35,-0.35}, // DDL 513
53 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 514
54 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 515
55 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 516
56 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 517
57 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 518
58 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 519
59 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.25,-0.15}, // DDL 520
60 {-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15,-0.15}, // DDL 521
61 {-0.10,-0.10,-0.10,-0.40,-0.40,-0.40,-0.10,-0.10,-0.45}, // DDL 522
62 {-0.10,-0.10,-0.10,-0.35,-0.35,-0.35,-0.10,-0.35,-0.50}, // DDL 523
63 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 524
64 { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}, // DDL 525
65 { 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35, 0.35}, // DDL 526
66 { 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45}}; // DDL 527
68 ClassImp(AliITSClusterFinderV2SSD)
71 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1), fLorentzShiftP(0), fLorentzShiftN(0)
74 static AliITSRecoParam *repa = NULL;
76 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
78 repa = AliITSRecoParam::GetHighFluxParam();
79 AliWarning("Using default AliITSRecoParam class");
83 if (repa->GetCorrectLorentzAngleSSD()) {
84 AliMagF* field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
86 AliError("Cannot get magnetic field from TGeoGlobalMagField");
89 Float_t bField = field->SolenoidField();
90 // NB: spatial shift has opposite sign for lay 5 and 6, but strip numbering also changes direction, so no sign-change
91 // Shift due to ExB on drift N-side, units: strip width
92 fLorentzShiftP = -repa->GetTanLorentzAngleElectronsSSD() * 150.e-4/95.e-4 * bField / 5.0;
93 // Shift due to ExB on drift P-side, units: strip width
94 fLorentzShiftN = -repa->GetTanLorentzAngleHolesSSD() * 150.e-4/95.e-4 * bField / 5.0;
95 AliDebug(1,Form("bField %f Lorentz Shift P-side %f N-side %f",bField,fLorentzShiftN,fLorentzShiftP));
100 //______________________________________________________________________
101 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinder(cf), fLastSSD1(cf.fLastSSD1), fLorentzShiftP(cf.fLorentzShiftP), fLorentzShiftN(cf.fLorentzShiftN)
106 //______________________________________________________________________
107 AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD& cf ){
108 // Assignment operator
110 this->~AliITSClusterFinderV2SSD();
111 new(this) AliITSClusterFinderV2SSD(cf);
116 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
120 FindClustersSSD(fDigits);
124 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
125 //------------------------------------------------------------
126 // Actual SSD cluster finder
127 //------------------------------------------------------------
128 Int_t smaxall=alldigits->GetEntriesFast();
129 if (smaxall==0) return;
132 //---------------------------------------
133 // load recoparam and calibration
135 static AliITSRecoParam *repa = NULL;
137 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
139 repa = AliITSRecoParam::GetHighFluxParam();
140 AliWarning("Using default AliITSRecoParam class");
144 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
147 //---------------------------------------
150 //------------------------------------
151 // fill the digits array with zero-suppression condition
152 // Signal is converted in KeV
155 for (Int_t i=0;i<smaxall; i++){
156 AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
158 if(d->IsSideP()) noise = cal->GetNoiseP(d->GetStripNumber());
159 else noise = cal->GetNoiseN(d->GetStripNumber());
160 if (d->GetSignal()<3.*noise) continue;
162 if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());
163 else gain = cal->GetGainN(d->GetStripNumber());
165 Float_t q=gain*d->GetSignal(); //
166 q=cal->ADCToKeV(q); // converts the charge in KeV from ADC units
167 d->SetSignal(Int_t(q));
171 Int_t smax = digits.GetEntriesFast();
173 //------------------------------------
176 const Int_t kMax=1000;
178 Ali1Dcluster pos[kMax], neg[kMax];
179 Float_t y=0., q=0., qmax=0.;
180 Int_t lab[4]={-2,-2,-2,-2};
184 cout<<"-----------------------------"<<endl;
185 cout<<"this is module "<<fModule;
190 if (fModule>fLastSSD1)
193 //--------------------------------------------------------
194 // start 1D-clustering from the first digit in the digits array
196 AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
198 y += d->GetCoord2()*d->GetSignal();
200 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
203 noise = cal->GetNoiseP(d->GetStripNumber());
204 gain = cal->GetGainP(d->GetStripNumber());
207 noise = cal->GetNoiseN(d->GetStripNumber());
208 gain = cal->GetGainN(d->GetStripNumber());
211 noise=cal->ADCToKeV(noise); // converts noise in KeV from ADC units
213 if(qmax>fgkThreshold*noise) flag5=1; // seed for the cluster
216 cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
217 d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
220 Int_t curr=d->GetCoord2();
221 Int_t flag=d->GetCoord1();
223 // Note: the first side which will be processed is supposed to be the
224 // P-side which is neg
227 if(flag) {n=&np; c=pos;} // in case we have only Nstrips (P was bad!)
231 for (Int_t ilab=0;ilab<10;ilab++){
234 milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
237 //----------------------------------------------------------
238 // search for neighboring digits
240 for (Int_t s=1; s<smax; s++) {
241 d=(AliITSdigitSSD*)digits.UncheckedAt(s);
242 Int_t strip=d->GetCoord2();
244 // if digits is not a neighbour or side did not change
245 // and at least one of the previous digits met the seed condition
246 // then creates a new 1D cluster
247 if ( ( ((strip-curr) > 1) || (flag!=d->GetCoord1()) ) ) {
250 //cout<<"here1"<<endl;
251 Float_t dLorentz = 0;
252 if (!flag) { // P-side is neg clust
253 dLorentz = fLorentzShiftN;
255 else { // N-side is p clust
256 dLorentz = fLorentzShiftP;
258 c[*n].SetY(y/q+dLorentz);
262 c[*n].SetLabels(milab);
264 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
265 // Note: fUseUnfoldingInClusterFinderSSD=kFALSE by default in RecoParam
267 //Split suspiciously big cluster
269 c[*n].SetY(y/q-0.25*nd+dLorentz);
273 Error("FindClustersSSD","Too many 1D clusters !");
276 c[*n].SetY(y/q+0.25*nd+dLorentz);
279 c[*n].SetLabels(milab);
286 Error("FindClustersSSD","Too many 1D clusters !");
296 lab[0]=lab[1]=lab[2]=-2;
297 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
299 // if side changed from P to N, switch to pos 1D clusters
300 // (if for some reason the side changed from N to P then do the opposite)
301 if (flag!=d->GetCoord1())
302 { if(!flag) {n=&np; c=pos;} else {n=&nn; c=neg;} }
304 } // end create new 1D cluster from previous neighboring digits
306 // continues adding digits to the previous cluster
307 // or start a new one
310 y += d->GetCoord2()*d->GetSignal();
314 noise = cal->GetNoiseP(d->GetStripNumber());
315 gain = cal->GetGainP(d->GetStripNumber());
318 noise = cal->GetNoiseN(d->GetStripNumber());
319 gain = cal->GetGainN(d->GetStripNumber());
322 noise=cal->ADCToKeV(noise); // converts the charge in KeV from ADC units
324 if(d->GetSignal()>fgkThreshold*noise) flag5=1;
327 cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
328 d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
331 if (d->GetSignal()>qmax) {
333 lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
335 for (Int_t ilab=0;ilab<10;ilab++) {
336 if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab)));
341 } // loop over digits, no more digits in the digits array
344 // add the last 1D cluster
347 // cout<<"here2"<<endl;
348 Float_t dLorentz = 0;
349 if (!flag) { // P-side is neg clust
350 dLorentz = fLorentzShiftN;
352 else { // N-side is p clust
353 dLorentz = fLorentzShiftP;
356 c[*n].SetY(y/q + dLorentz);
359 c[*n].SetLabels(lab);
361 if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
363 //Split suspiciously big cluster
365 c[*n].SetY(y/q-0.25*nd + dLorentz);
369 Error("FindClustersSSD","Too many 1D clusters !");
372 c[*n].SetY(y/q+0.25*nd + dLorentz);
375 c[*n].SetLabels(lab);
381 Error("FindClustersSSD","Too many 1D clusters !");
385 } // if flag5 last 1D cluster added
388 //------------------------------------------------------
389 // call FindClustersSSD to pair neg and pos 1D clusters
390 // and create recpoints from the crosses
391 // Note1: neg are Pside and pos are Nside!!
392 // Note2: if there are no Pside digits nn=0 (bad strips!!) (same for Nside)
394 // cout<<nn<<" Pside and "<<np<<" Nside clusters"<<endl;
396 AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();
398 TClonesArray* clusters = rpc->UncheckedGetClusters(fModule);
400 FindClustersSSD(neg, nn, pos, np, clusters);
403 while ((irp = (AliITSRecPoint*)itr.Next())) fDetTypeRec->AddRecPoint(*irp);
405 //-----------------------------------------------------
409 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader){
411 //------------------------------------------------------------
412 // This function creates ITS clusters from raw data
413 //------------------------------------------------------------
416 AliITSRawStreamSSD inputSSD(rawReader);
417 FindClustersSSD(&inputSSD);
422 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input)
424 //------------------------------------------------------------
425 // Actual SSD cluster finder for raw data
426 //------------------------------------------------------------
428 AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();
429 static AliITSRecoParam *repa = NULL;
431 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
433 repa = AliITSRecoParam::GetHighFluxParam();
434 AliWarning("Using default AliITSRecoParam class");
437 if (fRawID2ClusID) { // RS: reset references from 1D clusters to rawID's
438 fRawIDRef[0].Reset();
439 fRawIDRef[1].Reset();
441 Int_t nClustersSSD = 0;
442 const Int_t kNADC = 12;
443 const Int_t kMaxADCClusters = 1000;
445 Int_t strips[kNADC][2][kMaxADCClusters][3]; // [ADC],[side],[istrip], [0]=istrip [1]=signal [2]=rawID (for embedding, RS)
446 Int_t nStrips[kNADC][2];
448 for( int i=0; i<kNADC; i++ ){
457 //* Loop over modules DDL+AD
459 int countRW = 0; //RS
460 if (fRawID2ClusID) fRawID2ClusID->Reset(); //RS if array was provided, we shall store the rawID -> ClusterID
464 bool next = input->Next();
467 //* Continue if corrupted input
470 if( (!next)&&(input->flag) ){
471 AliWarning(Form("ITSClustersFinderSSD: Corrupted data: warning from RawReader"));
475 Int_t newDDL = input->GetDDL();
476 Int_t newAD = input->GetAD();
479 if( newDDL<0 || newDDL>15 ){
480 AliWarning(Form("ITSClustersFinderSSD: Corrupted data: wrong DDL number (%d)",newDDL));
484 if( newAD<1 || newAD>9 ){
485 AliWarning(Form("ITSClustersFinderSSD: Corrupted data: wrong AD number (%d)",newAD));
490 bool newModule = ( !next || ddl!= newDDL || ad!=newAD );
492 if( newModule && ddl>=0 && ad>=0 ){
495 //* Reconstruct the previous block of 12 modules --- actual clusterfinder
498 for( int adc = 0; adc<kNADC; adc++ ){
502 Ali1Dcluster clusters1D[2][kMaxADCClusters]; // per ADC, per side
503 Int_t nClusters1D[2] = {0,0};
504 //int nstat[2] = {0,0};
505 fModule = AliITSRawStreamSSD::GetModuleNumber(ddl, (ad - 1) * 12 + adc );
508 // AliWarning(Form("ITSClustersFinderSSD: Corrupted data: module (ddl %d ad %d adc %d) not found in the map",ddl,ad,adc));
509 //CM channels are always present even everything is suppressed
514 if (fModule>fLastSSD1)
517 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)fDetTypeRec->GetCalibrationModel(fModule);
519 AliWarning(Form("ITSClustersFinderSSD: No calibration found for module (ddl %d ad %d adc %d)",ddl,ad,adc));
525 if( repa->GetUseCosmicRunShiftsSSD()) { // Special condition for 2007/2008 cosmic data
526 dStrip = fgkCosmic2008StripShifts[ddl][ad-1];
527 if (TMath::Abs(dStrip) > 1.5){
528 AliWarning(Form("Indexing error in Cosmic calibration: ddl = %d, dStrip %f\n",ddl,dStrip));
533 for( int side=0; side<=1; side++ ){
535 Int_t lab[3]={-2,-2,-2};
542 Float_t dLorentz = 0;
543 if (side==0) { // P-side is neg clust
544 dLorentz = fLorentzShiftN;
546 else { // N-side is pos clust
547 dLorentz = fLorentzShiftP;
550 Int_t n = nStrips[adc][side];
551 for( int istr = 0; istr<n+1; istr++ ){
554 Int_t strip=0, rwID = 0;
555 Float_t signal=0.0, noise=0.0, gain=0.0;
558 strip = strips[adc][side][istr][0];
559 signal = strips[adc][side][istr][1];
560 rwID = strips[adc][side][istr][2]; // RS
561 //cout<<"strip "<<adc<<" / "<<side<<": "<<strip<<endl;
564 noise = side ?cal->GetNoiseN(strip) :cal->GetNoiseP(strip);
565 gain = side ?cal->GetGainN(strip) :cal->GetGainP(strip);
566 stripOK = ( noise>=1. && signal>=3.0*noise
567 //&& !cal->IsPChannelBad(strip)
570 } else stripOK = 0; // end of data
572 bool newCluster = ( (abs(strip-ostrip)!=1) || !stripOK );
576 //* Store the previous cluster
578 if( nDigits>0 && q>0 && snFlag ){
580 if (nClusters1D[side] >= kMaxADCClusters-1 ) {
581 AliWarning("HLT ClustersFinderSSD: Too many 1D clusters !");
584 Ali1Dcluster &cluster = clusters1D[side][nClusters1D[side]++];
585 cluster.SetY( y / q + dStrip + dLorentz);
587 cluster.SetNd(nDigits);
588 cluster.SetLabels(lab);
589 //cout<<"cluster 1D side "<<side<<": y= "<<y<<" q= "<<q<<" d="<<dStrip<<" Y="<<cluster.GetY()<<endl;
590 //Split suspiciously big cluster
592 if( repa->GetUseUnfoldingInClusterFinderSSD()
593 && nDigits > 4 && nDigits < 25
595 cluster.SetY(y/q + dStrip - 0.25*nDigits + dLorentz);
597 Ali1Dcluster& cluster2 = clusters1D[side][nClusters1D[side]++];
598 cluster2.SetY(y/q + dStrip + 0.25*nDigits + dLorentz);
599 cluster2.SetQ(0.5*q);
600 cluster2.SetNd(nDigits);
601 cluster2.SetLabels(lab);
609 } //* End store the previous cluster
611 if( stripOK ){ // add new signal to the cluster
612 // signal = (Int_t) (signal * gain); // signal is corrected for gain
613 if( signal>fgkThreshold*noise) snFlag = 1;
614 signal = signal * gain; // signal is corrected for gain
615 // if( cal ) signal = (Int_t) cal->ADCToKeV( signal ); // signal is converted in KeV
616 if( cal ) signal = cal->ADCToKeV( signal ); // signal is converted in KeV
617 q += signal; // add digit to current cluster
622 if (fRawID2ClusID) fRawIDRef[side].AddReference(nClusters1D[side],rwID);
625 } //* end loop over strips
627 } //* end loop over ADC sides
631 if( nClusters1D[0] && nClusters1D[1] && fModule>=0 ){
632 TClonesArray* clusters = rpc->UncheckedGetClusters(fModule);
633 FindClustersSSD( clusters1D[0], nClusters1D[0], clusters1D[1], nClusters1D[1], clusters);
634 Int_t nClustersn = clusters->GetEntriesFast();
635 nClustersSSD += nClustersn;
638 //cout<<"SG: "<<ddl<<" "<<ad<<" "<<adc<<": strips "<<nstat[0]<<"+"<<nstat[1]<<", clusters 1D= "<<nClusters1D[0]<<" + "<<nClusters1D[1]<<", 2D= "<<clusters.size()<<endl;
640 }//* end loop over adc
642 }//* end of reconstruction of previous block of 12 modules
647 //* Clean up arrays and set new module
650 for( int i=0; i<kNADC; i++ ){
660 //* Exit main loop when there is no more input
666 //* Fill the current strip information
669 Int_t adc = input->GetADC();
670 if( adc<0 || adc>=kNADC+2 || (adc>5&&adc<8) ){
671 AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong adc number (%d)", adc));
675 if( adc>7 ) adc-= 2; // shift ADC numbers 8-13 to 6-11
677 Bool_t side = input->GetSideFlag();
678 Int_t strip = input->GetStrip();
679 Int_t signal = input->GetSignal();
682 //cout<<"SSD: "<<ddl<<" "<<ad<<" "<<adc<<" "<<side<<" "<<strip<<" : "<<signal<<endl;
685 AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong strip number (ddl %d ad %d adc %d side %d, strip %d",
686 ddl, ad, adc, side,strip) );
689 if (strip < 0) continue;
691 int &n = nStrips[adc][side];
693 Int_t oldStrip = strips[adc][side][n-1][0];
695 if( strip==oldStrip ){
696 AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: duplicated signal: ddl %d ad %d adc %d, side %d, strip %d",
697 ddl, ad, adc, side, strip ));
701 strips[adc][side][n][0] = strip;
702 strips[adc][side][n][1] = signal;
703 strips[adc][side][n][2] = countRW;
706 //cout<<"SSD: "<<input->GetDDL()<<" "<<input->GetAD()<<" "
707 //<<input->GetADC()<<" "<<input->GetSideFlag()<<" "<<((int)input->GetStrip())<<" "<<strip<<" : "<<input->GetSignal()<<endl;
710 } //* End main loop over the input
712 AliDebug(1,Form("found clusters in ITS SSD: %d", nClustersSSD));
716 void AliITSClusterFinderV2SSD::
717 FindClustersSSD(const Ali1Dcluster* neg, Int_t nn,
718 const Ali1Dcluster* pos, Int_t np,
719 TClonesArray *clusters) {
720 //------------------------------------------------------------
721 // Actual SSD cluster finder
722 //------------------------------------------------------------
724 const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
726 //---------------------------------------
729 static AliITSRecoParam *repa = NULL;
731 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
733 repa = AliITSRecoParam::GetHighFluxParam();
734 AliWarning("Using default AliITSRecoParam class");
738 // TClonesArray &cl=*clusters;
740 AliITSsegmentationSSD *seg = static_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
741 if (fModule>fLastSSD1)
746 Float_t hwSSD = seg->Dx()*1e-4/2;
747 Float_t hlSSD = seg->Dz()*1e-4/2;
749 Int_t idet=fNdet[fModule];
753 Int_t *cnegative = new Int_t[np];
754 Int_t *cused1 = new Int_t[np];
755 Int_t *negativepair = new Int_t[10*np];
756 Int_t *cpositive = new Int_t[nn];
757 Int_t *cused2 = new Int_t[nn];
758 Int_t *positivepair = new Int_t[10*nn];
759 for (Int_t i=0;i<np;i++) {cnegative[i]=0; cused1[i]=0;}
760 for (Int_t i=0;i<nn;i++) {cpositive[i]=0; cused2[i]=0;}
761 for (Int_t i=0;i<10*np;i++) {negativepair[i]=0;}
762 for (Int_t i=0;i<10*nn;i++) {positivepair[i]=0;}
764 if ((np*nn) > fgPairsSize) {
767 fgPairsSize = 2*np*nn;
768 fgPairs = new Short_t[fgPairsSize];
770 memset(fgPairs,0,sizeof(Short_t)*np*nn);
773 // find available pairs
776 for (Int_t i=0; i<np; i++) {
777 Float_t yp=pos[i].GetY();
778 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
779 for (Int_t j=0; j<nn; j++) {
780 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
781 Float_t yn=neg[j].GetY();
784 seg->GetPadCxz(yn, yp, xt, zt);
785 //cout<<yn<<" "<<yp<<" "<<xt<<" "<<zt<<endl;
787 if (TMath::Abs(xt)<hwSSD)
788 if (TMath::Abs(zt)<hlSSD) {
789 Int_t in = i*10+cnegative[i];
790 Int_t ip = j*10+cpositive[j];
791 if ((in < 10*np) && (ip < 10*nn)) {
792 negativepair[in] =j; //index
794 cnegative[i]++; //counters
800 AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
808 delete [] negativepair;
811 delete [] positivepair;
814 //why not to allocate memorey here? if(!clusters) clusters = new TClonesArray("AliITSRecPoint", ncross);
817 // try to recover points out of but close to the module boundaries
819 for (Int_t i=0; i<np; i++) {
820 Float_t yp=pos[i].GetY();
821 if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
822 for (Int_t j=0; j<nn; j++) {
823 if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
824 // if both 1Dclusters have an other cross continue
825 if (cpositive[j]&&cnegative[i]) continue;
826 Float_t yn=neg[j].GetY();
829 seg->GetPadCxz(yn, yp, xt, zt);
831 if (TMath::Abs(xt)<hwSSD+0.1)
832 if (TMath::Abs(zt)<hlSSD+0.15) {
833 // tag 1Dcluster (eventually will produce low quality recpoint)
834 if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
835 if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
836 Int_t in = i*10+cnegative[i];
837 Int_t ip = j*10+cpositive[j];
838 if ((in < 10*np) && (ip < 10*nn)) {
839 negativepair[in] =j; //index
841 cnegative[i]++; //counters
846 AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
858 if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
864 for (Int_t ip=0;ip<np;ip++){
865 Float_t xbest=1000,zbest=1000,qbest=0;
867 // select gold clusters
868 if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
869 Float_t yp=pos[ip].GetY();
870 Int_t j = negativepair[10*ip];
872 if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
873 // both bad, hence continue;
874 // mark both as used (to avoid recover at the end)
880 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
881 //cout<<"ratio="<<ratio<<endl;
883 // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
884 if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
887 Float_t yn=neg[j].GetY();
890 seg->GetPadCxz(yn, yp, xt, zt);
895 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
896 if( (pos[ip].GetQ()==0)||(neg[j].GetQ()==0)) qbest*=2; // in case of bad strips on one side keep all charge from the other one
899 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
900 mT2L->MasterToLocal(loc,trk);
905 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
906 for (Int_t ilab=0;ilab<3;ilab++){
907 milab[ilab] = pos[ip].GetLabel(ilab);
908 milab[ilab+3] = neg[j].GetLabel(ilab);
912 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
913 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
915 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
916 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
917 // out-of-diagonal element of covariance matrix
918 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
919 else if ( (info[0]>1) && (info[1]>1) ) {
920 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
921 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
925 lp[2]=4.80e-06; // 0.00219*0.00219
926 lp[3]=0.0093; // 0.0964*0.0964;
931 lp[2]=2.79e-06; // 0.0017*0.0017;
932 lp[3]=0.00935; // 0.967*0.967;
937 AliITSRecPoint * cl2;
938 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
940 cl2->SetChargeRatio(ratio);
944 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
949 if(pos[ip].GetQ()==0) cl2->SetType(3);
950 if(neg[j].GetQ()==0) cl2->SetType(4);
959 for (Int_t ip=0;ip<np;ip++){
960 Float_t xbest=1000,zbest=1000,qbest=0;
963 // select "silber" cluster
964 if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
965 Int_t in = negativepair[10*ip];
966 Int_t ip2 = positivepair[10*in];
967 if (ip2==ip) ip2 = positivepair[10*in+1];
968 Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
972 ratio = (pcharge-neg[in].GetQ())/(pcharge+neg[in].GetQ());
973 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
974 //if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
978 if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) { //
980 Float_t yp=pos[ip].GetY();
981 Float_t yn=neg[in].GetY();
984 seg->GetPadCxz(yn, yp, xt, zt);
988 qbest =pos[ip].GetQ();
989 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
990 mT2L->MasterToLocal(loc,trk);
995 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
996 for (Int_t ilab=0;ilab<3;ilab++){
997 milab[ilab] = pos[ip].GetLabel(ilab);
998 milab[ilab+3] = neg[in].GetLabel(ilab);
1001 CheckLabels2(milab);
1002 ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
1003 milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
1004 Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1006 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1007 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1008 // out-of-diagonal element of covariance matrix
1009 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1010 else if ( (info[0]>1) && (info[1]>1) ) {
1011 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1012 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1016 lp[2]=4.80e-06; // 0.00219*0.00219
1017 lp[3]=0.0093; // 0.0964*0.0964;
1022 lp[2]=2.79e-06; // 0.0017*0.0017;
1023 lp[3]=0.00935; // 0.967*0.967;
1028 AliITSRecPoint * cl2;
1029 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1030 cl2->SetChargeRatio(ratio);
1032 fgPairs[ip*nn+in] = 5;
1033 if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1035 fgPairs[ip*nn+in] = 6;
1044 // if (!(cused1[ip2] || cused2[in])){ //
1045 if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
1047 Float_t yp=pos[ip2].GetY();
1048 Float_t yn=neg[in].GetY();
1051 seg->GetPadCxz(yn, yp, xt, zt);
1055 qbest =pos[ip2].GetQ();
1057 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1058 mT2L->MasterToLocal(loc,trk);
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[ip2].GetLabel(ilab);
1066 milab[ilab+3] = neg[in].GetLabel(ilab);
1069 CheckLabels2(milab);
1070 ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
1071 milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
1072 Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1074 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1075 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1076 // out-of-diagonal element of covariance matrix
1077 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1078 else if ( (info[0]>1) && (info[1]>1) ) {
1079 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1080 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1084 lp[2]=4.80e-06; // 0.00219*0.00219
1085 lp[3]=0.0093; // 0.0964*0.0964;
1090 lp[2]=2.79e-06; // 0.0017*0.0017;
1091 lp[3]=0.00935; // 0.967*0.967;
1096 AliITSRecPoint * cl2;
1097 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1099 cl2->SetChargeRatio(ratio);
1101 fgPairs[ip2*nn+in] =5;
1102 if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1104 fgPairs[ip2*nn+in] =6;
1113 } // charge matching condition
1115 } // 2 Pside cross 1 Nside
1116 } // loop over Pside clusters
1121 for (Int_t jn=0;jn<nn;jn++){
1122 if (cused2[jn]) continue;
1123 Float_t xbest=1000,zbest=1000,qbest=0;
1124 // select "silber" cluster
1125 if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
1126 Int_t ip = positivepair[10*jn];
1127 Int_t jn2 = negativepair[10*ip];
1128 if (jn2==jn) jn2 = negativepair[10*ip+1];
1129 Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
1133 ratio = (pcharge-pos[ip].GetQ())/(pcharge+pos[ip].GetQ());
1134 if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
1137 if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) && // charge matching
1138 (pcharge!=0) ) { // reject combinations of bad strips
1144 // if (!(cused1[ip]||cused2[jn])){
1145 if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) { //
1147 Float_t yn=neg[jn].GetY();
1148 Float_t yp=pos[ip].GetY();
1151 seg->GetPadCxz(yn, yp, xt, zt);
1155 qbest =neg[jn].GetQ();
1158 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1159 mT2L->MasterToLocal(loc,trk);
1165 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1166 for (Int_t ilab=0;ilab<3;ilab++){
1167 milab[ilab] = pos[ip].GetLabel(ilab);
1168 milab[ilab+3] = neg[jn].GetLabel(ilab);
1171 CheckLabels2(milab);
1172 ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1173 milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1174 Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1176 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1177 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1178 // out-of-diagonal element of covariance matrix
1179 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1180 else if ( (info[0]>1) && (info[1]>1) ) {
1181 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1182 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1186 lp[2]=4.80e-06; // 0.00219*0.00219
1187 lp[3]=0.0093; // 0.0964*0.0964;
1192 lp[2]=2.79e-06; // 0.0017*0.0017;
1193 lp[3]=0.00935; // 0.967*0.967;
1198 AliITSRecPoint * cl2;
1199 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1201 cl2->SetChargeRatio(ratio);
1203 fgPairs[ip*nn+jn] =7;
1204 if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1206 fgPairs[ip*nn+jn]=8;
1212 // if (!(cused1[ip]||cused2[jn2])){
1213 if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) { //
1215 Float_t yn=neg[jn2].GetY();
1216 Double_t yp=pos[ip].GetY();
1219 seg->GetPadCxz(yn, yp, xt, zt);
1223 qbest =neg[jn2].GetQ();
1226 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1227 mT2L->MasterToLocal(loc,trk);
1233 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1234 for (Int_t ilab=0;ilab<3;ilab++){
1235 milab[ilab] = pos[ip].GetLabel(ilab);
1236 milab[ilab+3] = neg[jn2].GetLabel(ilab);
1239 CheckLabels2(milab);
1240 ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1241 milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1242 Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
1244 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1245 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1246 // out-of-diagonal element of covariance matrix
1247 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1248 else if ( (info[0]>1) && (info[1]>1) ) {
1249 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1250 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1254 lp[2]=4.80e-06; // 0.00219*0.00219
1255 lp[3]=0.0093; // 0.0964*0.0964;
1260 lp[2]=2.79e-06; // 0.0017*0.0017;
1261 lp[3]=0.00935; // 0.967*0.967;
1266 AliITSRecPoint * cl2;
1267 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1270 cl2->SetChargeRatio(ratio);
1271 fgPairs[ip*nn+jn2]=7;
1273 if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1275 fgPairs[ip*nn+jn2]=8;
1283 } // charge matching condition
1285 } // 2 Nside cross 1 Pside
1286 } // loop over Pside clusters
1290 for (Int_t ip=0;ip<np;ip++){
1292 if(cused1[ip]) continue;
1295 Float_t xbest=1000,zbest=1000,qbest=0;
1299 if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){
1300 Float_t minchargediff =4.;
1301 Float_t minchargeratio =0.2;
1304 for (Int_t di=0;di<cnegative[ip];di++){
1305 Int_t jc = negativepair[ip*10+di];
1306 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1307 ratio = (pos[ip].GetQ()-neg[jc].GetQ())/(pos[ip].GetQ()+neg[jc].GetQ());
1308 //if (TMath::Abs(chargedif)<minchargediff){
1309 if (TMath::Abs(ratio)<0.2){
1311 minchargediff = TMath::Abs(chargedif);
1312 minchargeratio = TMath::Abs(ratio);
1315 if (j<0) continue; // not proper cluster
1319 for (Int_t di=0;di<cnegative[ip];di++){
1320 Int_t jc = negativepair[ip*10+di];
1321 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1322 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1324 if (count>1) continue; // more than one "proper" cluster for positive
1328 for (Int_t dj=0;dj<cpositive[j];dj++){
1329 Int_t ic = positivepair[j*10+dj];
1330 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1331 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1333 if (count>1) continue; // more than one "proper" cluster for negative
1338 for (Int_t dj=0;dj<cnegative[jp];dj++){
1339 Int_t ic = positivepair[jp*10+dj];
1340 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1341 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1343 if (count>1) continue;
1344 if (fgPairs[ip*nn+j]<100) continue;
1349 //almost gold clusters
1350 Float_t yp=pos[ip].GetY();
1351 Float_t yn=neg[j].GetY();
1353 seg->GetPadCxz(yn, yp, xt, zt);
1355 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1357 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1358 mT2L->MasterToLocal(loc,trk);
1363 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1364 for (Int_t ilab=0;ilab<3;ilab++){
1365 milab[ilab] = pos[ip].GetLabel(ilab);
1366 milab[ilab+3] = neg[j].GetLabel(ilab);
1369 CheckLabels2(milab);
1370 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1371 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1372 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1373 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1375 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1376 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1377 // out-of-diagonal element of covariance matrix
1378 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1379 else if ( (info[0]>1) && (info[1]>1) ) {
1380 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1381 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1385 lp[2]=4.80e-06; // 0.00219*0.00219
1386 lp[3]=0.0093; // 0.0964*0.0964;
1391 lp[2]=2.79e-06; // 0.0017*0.0017;
1392 lp[3]=0.00935; // 0.967*0.967;
1397 AliITSRecPoint * cl2;
1398 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1400 cl2->SetChargeRatio(ratio);
1402 fgPairs[ip*nn+j]=10;
1403 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1405 fgPairs[ip*nn+j]=11;
1412 } // loop over Pside 1Dclusters
1416 for (Int_t ip=0;ip<np;ip++){
1418 if(cused1[ip]) continue;
1421 Float_t xbest=1000,zbest=1000,qbest=0;
1423 // manyxmany clusters
1425 if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1426 Float_t minchargediff =4.;
1428 for (Int_t di=0;di<cnegative[ip];di++){
1429 Int_t jc = negativepair[ip*10+di];
1430 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1431 if (TMath::Abs(chargedif)<minchargediff){
1433 minchargediff = TMath::Abs(chargedif);
1436 if (j<0) continue; // not proper cluster
1439 for (Int_t di=0;di<cnegative[ip];di++){
1440 Int_t jc = negativepair[ip*10+di];
1441 Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1442 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1444 if (count>1) continue; // more than one "proper" cluster for positive
1448 for (Int_t dj=0;dj<cpositive[j];dj++){
1449 Int_t ic = positivepair[j*10+dj];
1450 Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1451 if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1453 if (count>1) continue; // more than one "proper" cluster for negative
1458 for (Int_t dj=0;dj<cnegative[jp];dj++){
1459 Int_t ic = positivepair[jp*10+dj];
1460 Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1461 if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1463 if (count>1) continue;
1464 if (fgPairs[ip*nn+j]<100) continue;
1467 //almost gold clusters
1468 Float_t yp=pos[ip].GetY();
1469 Float_t yn=neg[j].GetY();
1473 seg->GetPadCxz(yn, yp, xt, zt);
1477 qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1480 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1481 mT2L->MasterToLocal(loc,trk);
1486 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1487 for (Int_t ilab=0;ilab<3;ilab++){
1488 milab[ilab] = pos[ip].GetLabel(ilab);
1489 milab[ilab+3] = neg[j].GetLabel(ilab);
1492 CheckLabels2(milab);
1493 if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1494 ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1495 milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1496 Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1498 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1499 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1500 // out-of-diagonal element of covariance matrix
1501 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1502 else if ( (info[0]>1) && (info[1]>1) ) {
1503 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1504 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1508 lp[2]=4.80e-06; // 0.00219*0.00219
1509 lp[3]=0.0093; // 0.0964*0.0964;
1514 lp[2]=2.79e-06; // 0.0017*0.0017;
1515 lp[3]=0.00935; // 0.967*0.967;
1520 if (fRawID2ClusID) { // set rawID <-> clusterID correspondence for embedding
1521 const int kMaxRefRW = 200;
1522 UInt_t nrefsRW,refsRW[kMaxRefRW];
1523 nrefsRW = fRawIDRef[0].GetReferences(j,refsRW,kMaxRefRW); // n-side
1524 for (int ir=nrefsRW;ir--;) {
1525 int rwid = (int)refsRW[ir];
1526 if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
1527 (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
1530 nrefsRW = fRawIDRef[1].GetReferences(ip,refsRW,kMaxRefRW); // p-side
1531 for (int ir=nrefsRW;ir--;) {
1532 int rwid = (int)refsRW[ir];
1533 if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
1534 (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
1537 milab[0] = fNClusters+1; // RS: assign id as cluster label
1540 AliITSRecPoint * cl2;
1541 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1543 cl2->SetChargeRatio(ratio);
1545 fgPairs[ip*nn+j]=12;
1546 if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1548 fgPairs[ip*nn+j]=13;
1556 } // loop over Pside 1Dclusters
1558 } // use charge matching
1561 // recover all the other crosses
1563 for (Int_t i=0; i<np; i++) {
1564 Float_t xbest=1000,zbest=1000,qbest=0;
1565 Float_t yp=pos[i].GetY();
1566 if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
1567 for (Int_t j=0; j<nn; j++) {
1568 // for (Int_t di = 0;di<cpositive[i];di++){
1569 // Int_t j = negativepair[10*i+di];
1570 if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1572 if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1574 if (cused2[j]||cused1[i]) continue;
1575 if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1576 ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1577 Float_t yn=neg[j].GetY();
1580 seg->GetPadCxz(yn, yp, xt, zt);
1582 if (TMath::Abs(xt)<hwSSD)
1583 if (TMath::Abs(zt)<hlSSD) {
1586 qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1589 Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1590 mT2L->MasterToLocal(loc,trk);
1595 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1596 for (Int_t ilab=0;ilab<3;ilab++){
1597 milab[ilab] = pos[i].GetLabel(ilab);
1598 milab[ilab+3] = neg[j].GetLabel(ilab);
1601 CheckLabels2(milab);
1602 milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1603 Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1605 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1606 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1607 // out-of-diagonal element of covariance matrix
1608 if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1609 else if ( (info[0]>1) && (info[1]>1) ) {
1610 lp[2]=2.63e-06; // 0.0016*0.0016; //SigmaY2
1611 lp[3]=0.0065; // 0.08*0.08; //SigmaZ2
1615 lp[2]=4.80e-06; // 0.00219*0.00219
1616 lp[3]=0.0093; // 0.0964*0.0964;
1621 lp[2]=2.79e-06; // 0.0017*0.0017;
1622 lp[3]=0.00935; // 0.967*0.967;
1627 AliITSRecPoint * cl2;
1628 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1630 cl2->SetChargeRatio(ratio);
1631 cl2->SetType(100+cpositive[j]+cnegative[i]);
1633 if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1634 if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1642 if(repa->GetUseBadChannelsInClusterFinderSSD()==kTRUE) {
1644 //---------------------------------------------------------
1645 // recover crosses of good 1D clusters with bad strips on the other side
1646 // Note1: at first iteration skip modules with a bad side (or almost), (would produce too many fake!)
1647 // Note2: for modules with a bad side see below
1649 AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
1650 Int_t countPbad=0, countNbad=0;
1651 for(Int_t ib=0; ib<768; ib++) {
1652 if(cal->IsPChannelBad(ib)) countPbad++;
1653 if(cal->IsNChannelBad(ib)) countNbad++;
1655 // AliInfo(Form("module %d has %d P- and %d N-bad strips",fModule,countPbad,countNbad));
1657 if( (countPbad<100) && (countNbad<100) ) { // no bad side!!
1659 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1660 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1662 // for(Int_t ib=0; ib<768; ib++) { // loop over all Pstrips
1663 for(Int_t ib=15; ib<753; ib++) { // loop over all Pstrips
1665 if(cal->IsPChannelBad(ib)) { // check if strips is bad
1666 Float_t yN=pos[i].GetY();
1668 seg->GetPadCxz(1.*ib, yN, xt, zt);
1671 // bad Pstrip is crossing the Nside 1Dcluster -> create recpoint
1673 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1674 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1675 mT2L->MasterToLocal(loc,trk);
1678 lp[4]=pos[i].GetQ(); //Q
1679 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1680 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1681 CheckLabels2(milab);
1682 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1683 Int_t info[3] = {pos[i].GetNd(),0,fNlayer[fModule]};
1685 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1686 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1687 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1694 AliITSRecPoint * cl2;
1695 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1696 cl2->SetChargeRatio(1.);
1699 } // cross is within the detector
1705 } // end loop over Pstrips
1707 } // end loop over Nside 1D clusters
1709 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1710 if(cpositive[j]) continue;
1712 // for(Int_t ib=0; ib<768; ib++) { // loop over all Nside strips
1713 for(Int_t ib=15; ib<753; ib++) { // loop over all Nside strips
1715 if(cal->IsNChannelBad(ib)) { // check if strip is bad
1716 Float_t yP=neg[j].GetY();
1718 seg->GetPadCxz(yP, 1.*ib, xt, zt);
1721 // bad Nstrip is crossing the Pside 1Dcluster -> create recpoint
1723 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1724 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1725 mT2L->MasterToLocal(loc,trk);
1728 lp[4]=neg[j].GetQ(); //Q
1729 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1730 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1731 CheckLabels2(milab);
1732 milab[3]=( j << 10 ) + idet; // pos|neg|det
1733 Int_t info[3]={0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1735 lp[2]=4.968e-06; // 0.00223*0.00223; //SigmaY2
1736 lp[3]=0.012; // 0.110*0.110; //SigmaZ2
1737 lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1744 AliITSRecPoint * cl2;
1745 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1746 cl2->SetChargeRatio(1.);
1749 } // cross is within the detector
1754 } // end loop over Nstrips
1755 } // end loop over Pside 1D clusters
1759 //---------------------------------------------------------
1761 else if( (countPbad>700) && (countNbad<100) ) { // bad Pside!!
1763 for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1764 if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1767 Float_t yN=pos[i].GetY();
1769 if (seg->GetLayer()==5) yP = yN + (7.6/1.9);
1770 else yP = yN - (7.6/1.9);
1771 seg->GetPadCxz(yP, yN, xt, zt);
1773 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1774 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1775 mT2L->MasterToLocal(loc,trk);
1778 lp[4]=pos[i].GetQ(); //Q
1779 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1780 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1781 CheckLabels2(milab);
1782 milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1783 Int_t info[3] = {(Int_t)pos[i].GetNd(),0,fNlayer[fModule]};
1785 lp[2]=0.00098; // 0.031*0.031; //SigmaY2
1786 lp[3]=1.329; // 1.15*1.15; //SigmaZ2
1788 if(info[0]>1) lp[2]=0.00097;
1790 AliITSRecPoint * cl2;
1791 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1792 cl2->SetChargeRatio(1.);
1795 } // cross is within the detector
1799 } // end loop over Nside 1D clusters
1801 } // bad Pside module
1803 else if( (countNbad>700) && (countPbad<100) ) { // bad Nside!!
1805 for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1806 if(cpositive[j]) continue;
1809 Float_t yP=neg[j].GetY();
1811 if (seg->GetLayer()==5) yN = yP - (7.6/1.9);
1812 else yN = yP + (7.6/1.9);
1813 seg->GetPadCxz(yP, yN, xt, zt);
1815 if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1816 Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1817 mT2L->MasterToLocal(loc,trk);
1820 lp[4]=neg[j].GetQ(); //Q
1821 for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1822 for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1823 CheckLabels2(milab);
1824 milab[3]=( j << 10 ) + idet; // pos|neg|det
1825 Int_t info[3] = {0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1827 lp[2]=7.27e-05; // 0.0085*0.0085; //SigmaY2
1828 lp[3]=1.33; // 1.15*1.15; //SigmaZ2
1830 if(info[1]>1) lp[2]=6.91e-05;
1832 AliITSRecPoint * cl2;
1833 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1834 cl2->SetChargeRatio(1.);
1837 } // cross is within the detector
1841 } // end loop over Pside 1D clusters
1843 } // bad Nside module
1845 //---------------------------------------------------------
1847 } // use bad channels
1849 //cout<<ncl<<" clusters for this module"<<endl;
1851 delete [] cnegative;
1853 delete [] negativepair;
1854 delete [] cpositive;
1856 delete [] positivepair;