technial fix to suppress the warning
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderV2SSD.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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  **************************************************************************/
15
16 /* $Id$ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //            Implementation of the ITS clusterer V2 class                //
20 //                                                                        //
21 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch            //
22 //          Last revision: 13-05-09 Enrico Fragiacomo                     //
23 //                                  enrico.fragiacomo@ts.infn.it          //
24 //                                                                        //
25 ///////////////////////////////////////////////////////////////////////////
26
27 #include "AliITSClusterFinderV2SSD.h"
28
29 #include <Riostream.h>
30 #include <TGeoGlobalMagField.h>
31
32 #include "AliLog.h"
33 #include "AliMagF.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"
45
46 Short_t *AliITSClusterFinderV2SSD::fgPairs = 0x0;
47 Int_t    AliITSClusterFinderV2SSD::fgPairsSize = 0;
48 const Float_t  AliITSClusterFinderV2SSD::fgkThreshold = 5.;
49
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
67
68 ClassImp(AliITSClusterFinderV2SSD)
69
70
71   AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),fLastSSD1(AliITSgeomTGeo::GetModuleIndex(6,1,1)-1), fLorentzShiftP(0), fLorentzShiftN(0)
72 {
73 //Default constructor
74   static AliITSRecoParam *repa = NULL;  
75   if(!repa){
76     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
77     if(!repa){
78       repa = AliITSRecoParam::GetHighFluxParam();
79       AliWarning("Using default AliITSRecoParam class");
80     }
81   }
82
83   if (repa->GetCorrectLorentzAngleSSD()) {
84     AliMagF* field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField()); 
85     if (field == 0) {
86       AliError("Cannot get magnetic field from TGeoGlobalMagField");
87     }
88     else {
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));
96     }
97   }
98 }
99  
100 //______________________________________________________________________
101 AliITSClusterFinderV2SSD::AliITSClusterFinderV2SSD(const AliITSClusterFinderV2SSD &cf) : AliITSClusterFinder(cf), fLastSSD1(cf.fLastSSD1), fLorentzShiftP(cf.fLorentzShiftP), fLorentzShiftN(cf.fLorentzShiftN)
102 {
103   // Copy constructor
104 }
105
106 //______________________________________________________________________
107 AliITSClusterFinderV2SSD& AliITSClusterFinderV2SSD::operator=(const AliITSClusterFinderV2SSD&  cf ){
108   // Assignment operator
109
110   this->~AliITSClusterFinderV2SSD();
111   new(this) AliITSClusterFinderV2SSD(cf);
112   return *this;
113 }
114
115
116 void AliITSClusterFinderV2SSD::FindRawClusters(Int_t mod){
117
118   //Find clusters V2
119   SetModule(mod);
120   FindClustersSSD(fDigits);
121
122 }
123
124 void AliITSClusterFinderV2SSD::FindClustersSSD(TClonesArray *alldigits) {
125   //------------------------------------------------------------
126   // Actual SSD cluster finder
127   //------------------------------------------------------------
128   Int_t smaxall=alldigits->GetEntriesFast();
129   if (smaxall==0) return;
130
131
132   //---------------------------------------
133   // load recoparam and calibration
134   // 
135   static AliITSRecoParam *repa = NULL;  
136   if(!repa){
137     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
138     if(!repa){
139       repa = AliITSRecoParam::GetHighFluxParam();
140       AliWarning("Using default AliITSRecoParam class");
141     }
142   }
143
144   AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
145   Float_t gain=0;
146   Float_t noise=0;
147   //---------------------------------------
148
149
150   //------------------------------------
151   // fill the digits array with zero-suppression condition
152   // Signal is converted in KeV
153   //
154   TObjArray digits;
155   for (Int_t i=0;i<smaxall; i++){
156     AliITSdigitSSD *d=(AliITSdigitSSD*)alldigits->UncheckedAt(i);
157
158     if(d->IsSideP()) noise = cal->GetNoiseP(d->GetStripNumber());  
159     else noise = cal->GetNoiseN(d->GetStripNumber());
160     if (d->GetSignal()<3.*noise) continue;
161
162     if(d->IsSideP()) gain = cal->GetGainP(d->GetStripNumber());  
163     else gain = cal->GetGainN(d->GetStripNumber());
164
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));
168
169     digits.AddLast(d);
170   }
171   Int_t smax = digits.GetEntriesFast();
172   if (smax==0) return;
173   //------------------------------------
174
175   
176   const Int_t kMax=1000;
177   Int_t np=0, nn=0; 
178   Ali1Dcluster pos[kMax], neg[kMax];
179   Float_t y=0., q=0., qmax=0.; 
180   Int_t lab[4]={-2,-2,-2,-2};  
181   Bool_t flag5 = 0;
182   
183   /*
184   cout<<"-----------------------------"<<endl;
185   cout<<"this is module "<<fModule;
186   cout<<endl;
187   cout<<endl;
188   */
189   Int_t layer = 4;
190   if (fModule>fLastSSD1) 
191     layer = 5;
192
193   //--------------------------------------------------------
194   // start 1D-clustering from the first digit in the digits array
195   //
196   AliITSdigitSSD *d=(AliITSdigitSSD*)digits.UncheckedAt(0);
197   q += d->GetSignal();
198   y += d->GetCoord2()*d->GetSignal();
199   qmax=d->GetSignal();
200   lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
201   
202   if(d->IsSideP()) {
203     noise = cal->GetNoiseP(d->GetStripNumber());  
204     gain = cal->GetGainP(d->GetStripNumber());
205   }
206   else {
207     noise = cal->GetNoiseN(d->GetStripNumber());
208     gain = cal->GetGainN(d->GetStripNumber());
209   }  
210   noise*=gain;
211   noise=cal->ADCToKeV(noise); // converts noise in KeV from ADC units
212
213   if(qmax>fgkThreshold*noise) flag5=1; // seed for the cluster
214
215   /*
216   cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
217     d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
218   */
219
220   Int_t curr=d->GetCoord2();
221   Int_t flag=d->GetCoord1();
222
223   // Note: the first side which will be processed is supposed to be the 
224   // P-side which is neg
225   Int_t *n=&nn;
226   Ali1Dcluster *c=neg;
227   if(flag) {n=&np; c=pos;} // in case we have only Nstrips (P was bad!)
228
229   Int_t nd=1;
230   Int_t milab[10];
231   for (Int_t ilab=0;ilab<10;ilab++){
232     milab[ilab]=-2;
233   }
234   milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2);
235
236
237   //----------------------------------------------------------
238   // search for neighboring digits
239   //
240   for (Int_t s=1; s<smax; s++) {
241       d=(AliITSdigitSSD*)digits.UncheckedAt(s);      
242       Int_t strip=d->GetCoord2();
243
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()) ) ) {
248
249         if(flag5) {
250           //cout<<"here1"<<endl;
251           Float_t dLorentz = 0;
252           if (!flag) { // P-side is neg clust
253             dLorentz = fLorentzShiftN;
254           }
255           else { // N-side is p clust
256             dLorentz = fLorentzShiftP;
257           }
258          c[*n].SetY(y/q+dLorentz);
259          c[*n].SetQ(q);
260          c[*n].SetNd(nd);
261          CheckLabels2(milab);
262          c[*n].SetLabels(milab);
263
264          if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
265            // Note: fUseUnfoldingInClusterFinderSSD=kFALSE by default in RecoParam
266
267            //Split suspiciously big cluster
268            if (nd>4&&nd<25) {
269              c[*n].SetY(y/q-0.25*nd+dLorentz);
270              c[*n].SetQ(0.5*q);
271              (*n)++;
272              if (*n==kMax) {
273                Error("FindClustersSSD","Too many 1D clusters !");
274                return;
275              }
276              c[*n].SetY(y/q+0.25*nd+dLorentz);
277              c[*n].SetQ(0.5*q);
278              c[*n].SetNd(nd);
279              c[*n].SetLabels(milab);
280            }     
281            
282          } // unfolding is on
283
284          (*n)++;
285          if (*n==kMax) {
286           Error("FindClustersSSD","Too many 1D clusters !");
287           return;
288          }
289
290         } // flag5 set
291
292          // reset everything
293          y=q=qmax=0.;
294          nd=0;
295          flag5=0;
296          lab[0]=lab[1]=lab[2]=-2;
297          for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
298
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;} }
303
304       } // end create new 1D cluster from previous neighboring digits
305
306       // continues adding digits to the previous cluster 
307       // or start a new one 
308       flag=d->GetCoord1();
309       q += d->GetSignal();
310       y += d->GetCoord2()*d->GetSignal();
311       nd++;
312
313       if(d->IsSideP()) {
314         noise = cal->GetNoiseP(d->GetStripNumber());  
315         gain = cal->GetGainP(d->GetStripNumber());
316       }
317       else {
318         noise = cal->GetNoiseN(d->GetStripNumber());
319         gain = cal->GetGainN(d->GetStripNumber());
320       }
321       noise*=gain;
322       noise=cal->ADCToKeV(noise); // converts the charge in KeV from ADC units
323
324       if(d->GetSignal()>fgkThreshold*noise) flag5=1;
325
326       /*
327   cout<<d->GetSignal()<<" "<<noise<<" "<<flag5<<" "<<
328     d->GetCoord1()<<" "<<d->GetCoord2()<<endl;
329       */
330
331       if (d->GetSignal()>qmax) {
332          qmax=d->GetSignal();
333          lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2);
334       }
335       for (Int_t ilab=0;ilab<10;ilab++) {
336         if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab))); 
337       }
338       curr=strip;
339
340
341   } // loop over digits, no more digits in the digits array
342
343
344   // add the last 1D cluster 
345   if(flag5) {
346
347     //  cout<<"here2"<<endl;
348     Float_t dLorentz = 0;
349     if (!flag) { // P-side is neg clust
350       dLorentz = fLorentzShiftN;
351     }
352     else { // N-side is p clust
353       dLorentz = fLorentzShiftP;
354     }
355     
356     c[*n].SetY(y/q + dLorentz);
357     c[*n].SetQ(q);
358     c[*n].SetNd(nd);
359     c[*n].SetLabels(lab);
360     
361     if(repa->GetUseUnfoldingInClusterFinderSSD()==kTRUE) {
362       
363       //Split suspiciously big cluster
364       if (nd>4 && nd<25) {
365         c[*n].SetY(y/q-0.25*nd + dLorentz);
366         c[*n].SetQ(0.5*q);
367         (*n)++;
368         if (*n==kMax) {
369           Error("FindClustersSSD","Too many 1D clusters !");
370           return;
371         }
372         c[*n].SetY(y/q+0.25*nd + dLorentz);
373         c[*n].SetQ(0.5*q);
374         c[*n].SetNd(nd);
375         c[*n].SetLabels(lab);
376       }
377     } // unfolding is on
378     
379     (*n)++;
380     if (*n==kMax) {
381       Error("FindClustersSSD","Too many 1D clusters !");
382       return;
383     }
384
385   } // if flag5 last 1D cluster added
386
387
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)
393   //
394   //  cout<<nn<<" Pside and "<<np<<" Nside clusters"<<endl;
395   
396   AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();   
397   if (nn*np > 0) { 
398     TClonesArray* clusters = rpc->UncheckedGetClusters(fModule);
399     clusters->Clear();
400     FindClustersSSD(neg, nn, pos, np, clusters);
401     TIter itr(clusters);
402     AliITSRecPoint *irp;
403     while ((irp = (AliITSRecPoint*)itr.Next()))  fDetTypeRec->AddRecPoint(*irp);
404   }
405   //-----------------------------------------------------
406 }
407
408
409 void AliITSClusterFinderV2SSD::RawdataToClusters(AliRawReader* rawReader){
410
411   //------------------------------------------------------------
412   // This function creates ITS clusters from raw data
413   //------------------------------------------------------------
414   fNClusters = 0;
415   rawReader->Reset();
416   AliITSRawStreamSSD inputSSD(rawReader);
417   FindClustersSSD(&inputSSD);
418   
419 }
420
421
422 void AliITSClusterFinderV2SSD::FindClustersSSD(AliITSRawStreamSSD* input) 
423 {
424   //------------------------------------------------------------
425   // Actual SSD cluster finder for raw data
426   //------------------------------------------------------------
427
428   AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();
429   static AliITSRecoParam *repa = NULL;
430   if(!repa){
431     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
432     if(!repa){
433       repa = AliITSRecoParam::GetHighFluxParam();
434       AliWarning("Using default AliITSRecoParam class");
435     }
436   }
437   if (fRawID2ClusID) { // RS: reset references from 1D clusters to rawID's
438     fRawIDRef[0].Reset();
439     fRawIDRef[1].Reset();
440   }
441   Int_t nClustersSSD = 0;
442   const Int_t kNADC = 12;
443   const Int_t kMaxADCClusters = 1000;
444
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];
447
448   for( int i=0; i<kNADC; i++ ){
449     nStrips[i][0] = 0;
450     nStrips[i][1] = 0;
451   }
452
453   Int_t ddl = -1;
454   Int_t ad = -1;
455   
456   //*
457   //* Loop over modules DDL+AD
458   //*
459   int countRW = 0; //RS
460   if (fRawID2ClusID) fRawID2ClusID->Reset(); //RS if array was provided, we shall store the rawID -> ClusterID
461
462   while (kTRUE) {
463
464     bool next = input->Next();
465     
466     //* 
467     //* Continue if corrupted input
468     //*
469
470     if( (!next)&&(input->flag) ){
471      AliWarning(Form("ITSClustersFinderSSD: Corrupted data: warning from RawReader"));
472       continue; 
473     }
474
475     Int_t newDDL = input->GetDDL(); 
476     Int_t newAD = input->GetAD();
477
478     if( next ){
479       if( newDDL<0 || newDDL>15 ){
480         AliWarning(Form("ITSClustersFinderSSD: Corrupted data: wrong DDL number (%d)",newDDL));
481         continue;
482       }
483       
484       if( newAD<1 || newAD>9 ){
485         AliWarning(Form("ITSClustersFinderSSD: Corrupted data: wrong AD number (%d)",newAD));
486         continue;
487       }
488     }
489
490     bool newModule = ( !next || ddl!= newDDL || ad!=newAD );
491
492     if( newModule && ddl>=0 && ad>=0 ){ 
493
494       //*
495       //* Reconstruct the previous block of 12 modules --- actual clusterfinder
496       //* 
497       //cout<<endl;
498       for( int adc = 0; adc<kNADC; adc++ ){
499         
500         //* 1D clusterfinder
501
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 );
506         
507         if( fModule<0 ){
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 
510           continue;
511         }
512         
513         Int_t layer = 4;
514         if (fModule>fLastSSD1) 
515           layer = 5;
516
517         AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)fDetTypeRec->GetCalibrationModel(fModule);
518         if( !cal ){
519           AliWarning(Form("ITSClustersFinderSSD: No calibration found for module (ddl %d ad %d adc %d)",ddl,ad,adc));       
520           continue;
521         }
522
523         Float_t dStrip = 0;
524
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));
529             dStrip = 0;
530           }     
531         }
532         
533         for( int side=0; side<=1; side++ ){
534
535           Int_t lab[3]={-2,-2,-2};
536           Float_t q = 0.;
537           Float_t y = 0.;
538           Int_t nDigits = 0;
539           Int_t ostrip = -2;
540           Bool_t snFlag = 0;
541
542           Float_t dLorentz = 0;
543           if (side==0) { // P-side is neg clust
544             dLorentz = fLorentzShiftN;
545           }
546           else { // N-side is pos clust
547             dLorentz = fLorentzShiftP;
548           }
549           
550           Int_t n = nStrips[adc][side];
551           for( int istr = 0; istr<n+1; istr++ ){
552             
553             bool stripOK = 1;
554             Int_t strip=0, rwID = 0;
555             Float_t signal=0.0, noise=0.0, gain=0.0;
556             
557             if( istr<n ){
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;
562
563               if( cal ){
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) 
568                             );
569               }
570             } else stripOK = 0; // end of data
571
572             bool newCluster = ( (abs(strip-ostrip)!=1) || !stripOK );     
573                 
574             if( newCluster ){
575
576               //* Store the previous cluster
577
578               if( nDigits>0 && q>0 && snFlag ){
579
580                 if (nClusters1D[side] >= kMaxADCClusters-1 ) {
581                   AliWarning("HLT ClustersFinderSSD: Too many 1D clusters !");
582                 }else {
583                   
584                   Ali1Dcluster &cluster = clusters1D[side][nClusters1D[side]++];
585                   cluster.SetY( y / q + dStrip + dLorentz);
586                   cluster.SetQ(q);
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
591
592                   if( repa->GetUseUnfoldingInClusterFinderSSD()
593                       && nDigits > 4 && nDigits < 25 
594                       ){
595                     cluster.SetY(y/q + dStrip - 0.25*nDigits + dLorentz);           
596                     cluster.SetQ(0.5*q);          
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);      
602                   } // unfolding is on          
603                 }
604               }
605               y = q = 0.;
606               nDigits = 0;
607               snFlag = 0;
608
609             } //* End store the previous cluster
610
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
618               y += strip * signal;        
619               nDigits++;
620               //nstat[side]++;
621               ostrip = strip;
622               if (fRawID2ClusID) fRawIDRef[side].AddReference(nClusters1D[side],rwID);
623
624             }
625           } //* end loop over strips
626
627         } //* end loop over ADC sides
628         
629
630         //* 2D clusterfinder
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;
636         }
637
638         //cout<<"SG: "<<ddl<<" "<<ad<<" "<<adc<<": strips "<<nstat[0]<<"+"<<nstat[1]<<", clusters 1D= "<<nClusters1D[0]<<" + "<<nClusters1D[1]<<", 2D= "<<clusters.size()<<endl;
639
640       }//* end loop over adc
641       
642     }//* end of reconstruction of previous block of 12 modules
643     
644     if( newModule ){
645       
646       //*
647       //* Clean up arrays and set new module
648       //* 
649       
650       for( int i=0; i<kNADC; i++ ){
651         nStrips[i][0] = 0;
652         nStrips[i][1] = 0;
653       }     
654       ddl = newDDL;
655       ad = newAD;
656     } 
657     
658
659     //*
660     //* Exit main loop when there is no more input
661     //* 
662
663     if( !next ) break;
664     
665     //* 
666     //* Fill the current strip information
667     //* 
668
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));
672       continue;
673     }
674
675     if( adc>7 ) adc-= 2; // shift ADC numbers 8-13 to 6-11
676     
677     Bool_t side = input->GetSideFlag();
678     Int_t strip = input->GetStrip();
679     Int_t signal = input->GetSignal();
680     
681
682     //cout<<"SSD: "<<ddl<<" "<<ad<<" "<<adc<<" "<<side<<" "<<strip<<" : "<<signal<<endl;
683
684     if( strip>767 ){    
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) );
687       continue;
688     }
689     if (strip < 0) continue;
690     
691     int &n = nStrips[adc][side];
692     if( n >0 ){
693       Int_t oldStrip = strips[adc][side][n-1][0];
694
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 ));
698         continue;
699       }
700     }
701     strips[adc][side][n][0] = strip;
702     strips[adc][side][n][1] = signal;    
703     strips[adc][side][n][2] = countRW;    
704     n++;
705
706     //cout<<"SSD: "<<input->GetDDL()<<" "<<input->GetAD()<<" "
707     //<<input->GetADC()<<" "<<input->GetSideFlag()<<" "<<((int)input->GetStrip())<<" "<<strip<<" : "<<input->GetSignal()<<endl;
708     //
709     countRW++; //RS
710   } //* End main loop over the input
711   
712   AliDebug(1,Form("found clusters in ITS SSD: %d", nClustersSSD));
713 }
714
715
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   //------------------------------------------------------------
723
724   const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
725
726   //---------------------------------------
727   // load recoparam
728   // 
729   static AliITSRecoParam *repa = NULL;  
730   if(!repa){
731     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
732     if(!repa){
733       repa = AliITSRecoParam::GetHighFluxParam();
734       AliWarning("Using default AliITSRecoParam class");
735     }
736   }
737
738 //  TClonesArray &cl=*clusters;
739   
740   AliITSsegmentationSSD *seg = static_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
741   if (fModule>fLastSSD1) 
742     seg->SetLayer(6);
743   else 
744     seg->SetLayer(5);
745
746   Float_t hwSSD = seg->Dx()*1e-4/2;
747   Float_t hlSSD = seg->Dz()*1e-4/2;
748
749   Int_t idet=fNdet[fModule];
750   Int_t ncl=0;
751
752   //
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;}
763
764   if ((np*nn) > fgPairsSize) {
765
766     delete [] fgPairs;
767     fgPairsSize = 2*np*nn;
768     fgPairs = new Short_t[fgPairsSize];
769   }
770   memset(fgPairs,0,sizeof(Short_t)*np*nn);
771
772   //
773   // find available pairs
774   //
775   Int_t ncross = 0;
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();
782
783       Float_t xt, zt;
784       seg->GetPadCxz(yn, yp, xt, zt);
785       //cout<<yn<<" "<<yp<<" "<<xt<<" "<<zt<<endl;
786       
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
793           positivepair[ip] =i;
794           cnegative[i]++;  //counters
795           cpositive[j]++;
796           ncross++;     
797           fgPairs[i*nn+j]=100;
798         }
799         else
800           AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
801       }
802     }
803   }
804
805   if (!ncross) {
806     delete [] cnegative;
807     delete [] cused1;
808     delete [] negativepair;
809     delete [] cpositive;
810     delete [] cused2;
811     delete [] positivepair;
812     return;
813   }
814 //why not to allocate memorey here?  if(!clusters) clusters = new TClonesArray("AliITSRecPoint", ncross);
815   
816   /* //
817   // try to recover points out of but close to the module boundaries 
818   //
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();
827       
828       Float_t xt, zt;
829       seg->GetPadCxz(yn, yp, xt, zt);
830       
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
840           positivepair[ip] =i;
841           cnegative[i]++;  //counters
842           cpositive[j]++;       
843           fgPairs[i*nn+j]=100;
844         }
845         else
846           AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
847       }
848     }
849   }
850   */
851
852   //
853   Float_t lp[6];
854   Int_t milab[10];
855   Double_t ratio;
856   
857
858   if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
859
860
861     //
862     // sign gold tracks
863     //
864     for (Int_t ip=0;ip<np;ip++){
865       Float_t xbest=1000,zbest=1000,qbest=0;
866       //
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];      
871
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)
875           cused1[ip]++; 
876           cused2[j]++;
877           continue;
878         }
879
880         ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
881         //cout<<"ratio="<<ratio<<endl;
882
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
885
886         //
887         Float_t yn=neg[j].GetY();
888         
889         Float_t xt, zt;
890         seg->GetPadCxz(yn, yp, xt, zt);
891         
892         xbest=xt; zbest=zt; 
893
894         
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
897         
898         {
899           Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
900           mT2L->MasterToLocal(loc,trk);
901           lp[0]=trk[1];
902           lp[1]=trk[2];
903         }
904         lp[4]=qbest;        //Q
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);
909         }
910         //
911         CheckLabels2(milab);
912         milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
913         Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
914
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
922           lp[5]=-6.48e-05;
923         }
924         else {
925           lp[2]=4.80e-06;      // 0.00219*0.00219
926           lp[3]=0.0093;        // 0.0964*0.0964;
927           if (info[0]==1) {
928             lp[5]=-0.00014;
929           }
930           else { 
931             lp[2]=2.79e-06;    // 0.0017*0.0017; 
932             lp[3]=0.00935;     // 0.967*0.967;
933             lp[5]=-4.32e-05;
934           }
935         }
936
937         AliITSRecPoint * cl2;
938         cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
939           
940     cl2->SetChargeRatio(ratio);         
941     cl2->SetType(1);
942     fgPairs[ip*nn+j]=1;
943
944     if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
945       cl2->SetType(2);
946       fgPairs[ip*nn+j]=2;
947     }
948
949     if(pos[ip].GetQ()==0) cl2->SetType(3);
950     if(neg[j].GetQ()==0) cl2->SetType(4);
951
952     cused1[ip]++;
953     cused2[j]++;
954
955         ncl++;
956       }
957     }
958     
959     for (Int_t ip=0;ip<np;ip++){
960       Float_t xbest=1000,zbest=1000,qbest=0;
961       //
962       //
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();
969         
970
971
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) ) { // 
975           
976           //
977           // add first pair
978           if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) {  //
979             
980             Float_t yp=pos[ip].GetY(); 
981             Float_t yn=neg[in].GetY();
982             
983             Float_t xt, zt;
984             seg->GetPadCxz(yn, yp, xt, zt);
985             
986             xbest=xt; zbest=zt; 
987
988             qbest =pos[ip].GetQ();
989             Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
990             mT2L->MasterToLocal(loc,trk);
991             lp[0]=trk[1];
992             lp[1]=trk[2];
993             
994             lp[4]=qbest;        //Q
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);
999             }
1000             //
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]};
1005             
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
1013           lp[5]=-6.48e-05;
1014         }
1015         else {
1016           lp[2]=4.80e-06;      // 0.00219*0.00219
1017           lp[3]=0.0093;        // 0.0964*0.0964;
1018           if (info[0]==1) {
1019             lp[5]=-0.00014;
1020           }
1021           else { 
1022             lp[2]=2.79e-06;    // 0.0017*0.0017; 
1023             lp[3]=0.00935;     // 0.967*0.967;
1024             lp[5]=-4.32e-05;
1025           }
1026         }
1027
1028         AliITSRecPoint * cl2;
1029               cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1030               cl2->SetChargeRatio(ratio);       
1031               cl2->SetType(5);
1032               fgPairs[ip*nn+in] = 5;
1033               if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1034                 cl2->SetType(6);
1035                 fgPairs[ip*nn+in] = 6;
1036               }     
1037             ncl++;
1038           }
1039           
1040           
1041           //
1042           // add second pair
1043           
1044           //    if (!(cused1[ip2] || cused2[in])){  //
1045           if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
1046             
1047             Float_t yp=pos[ip2].GetY();
1048             Float_t yn=neg[in].GetY();
1049             
1050             Float_t xt, zt;
1051             seg->GetPadCxz(yn, yp, xt, zt);
1052             
1053             xbest=xt; zbest=zt; 
1054
1055             qbest =pos[ip2].GetQ();
1056             
1057             Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1058             mT2L->MasterToLocal(loc,trk);
1059             lp[0]=trk[1];
1060             lp[1]=trk[2];
1061             
1062             lp[4]=qbest;        //Q
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);
1067             }
1068             //
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]};
1073             
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
1081           lp[5]=-6.48e-05;
1082         }
1083         else {
1084           lp[2]=4.80e-06;      // 0.00219*0.00219
1085           lp[3]=0.0093;        // 0.0964*0.0964;
1086           if (info[0]==1) {
1087             lp[5]=-0.00014;
1088           }
1089           else { 
1090             lp[2]=2.79e-06;    // 0.0017*0.0017; 
1091             lp[3]=0.00935;     // 0.967*0.967;
1092             lp[5]=-4.32e-05;
1093           }
1094         }
1095
1096             AliITSRecPoint * cl2;
1097               cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1098               
1099               cl2->SetChargeRatio(ratio);       
1100               cl2->SetType(5);
1101               fgPairs[ip2*nn+in] =5;
1102               if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1103                 cl2->SetType(6);
1104                 fgPairs[ip2*nn+in] =6;
1105               }
1106             ncl++;
1107           }
1108           
1109           cused1[ip]++;
1110           cused1[ip2]++;
1111           cused2[in]++;
1112
1113         } // charge matching condition
1114
1115       } // 2 Pside cross 1 Nside
1116     } // loop over Pside clusters
1117     
1118     
1119       
1120       //  
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();
1130         //
1131         
1132
1133         ratio = (pcharge-pos[ip].GetQ())/(pcharge+pos[ip].GetQ());
1134         if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
1135
1136           /*
1137         if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) &&  // charge matching 
1138              (pcharge!=0) ) { // reject combinations of bad strips
1139           */
1140
1141
1142           //
1143           // add first pair
1144           //    if (!(cused1[ip]||cused2[jn])){
1145           if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) {  //
1146             
1147             Float_t yn=neg[jn].GetY(); 
1148             Float_t yp=pos[ip].GetY();
1149
1150             Float_t xt, zt;
1151             seg->GetPadCxz(yn, yp, xt, zt);
1152             
1153             xbest=xt; zbest=zt; 
1154
1155             qbest =neg[jn].GetQ();
1156
1157             {
1158               Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1159               mT2L->MasterToLocal(loc,trk);
1160               lp[0]=trk[1];
1161               lp[1]=trk[2];
1162           }
1163           
1164           lp[4]=qbest;        //Q
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);
1169           }
1170           //
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]};
1175
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
1183           lp[5]=-6.48e-05;
1184         }
1185         else {
1186           lp[2]=4.80e-06;      // 0.00219*0.00219
1187           lp[3]=0.0093;        // 0.0964*0.0964;
1188           if (info[0]==1) {
1189             lp[5]=-0.00014;
1190           }
1191           else { 
1192             lp[2]=2.79e-06;    // 0.0017*0.0017; 
1193             lp[3]=0.00935;     // 0.967*0.967;
1194             lp[5]=-4.32e-05;
1195           }
1196         }
1197
1198           AliITSRecPoint * cl2;
1199             cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1200
1201             cl2->SetChargeRatio(ratio);         
1202             cl2->SetType(7);
1203             fgPairs[ip*nn+jn] =7;
1204             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1205               cl2->SetType(8);
1206               fgPairs[ip*nn+jn]=8;
1207             }
1208           ncl++;
1209         }
1210         //
1211         // add second pair
1212         //      if (!(cused1[ip]||cused2[jn2])){
1213         if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) {  //
1214
1215           Float_t yn=neg[jn2].GetY(); 
1216           Double_t yp=pos[ip].GetY(); 
1217
1218           Float_t xt, zt;
1219           seg->GetPadCxz(yn, yp, xt, zt);
1220           
1221           xbest=xt; zbest=zt; 
1222
1223           qbest =neg[jn2].GetQ();
1224
1225           {
1226           Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1227           mT2L->MasterToLocal(loc,trk);
1228           lp[0]=trk[1];
1229           lp[1]=trk[2];
1230           }
1231
1232           lp[4]=qbest;        //Q
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);
1237           }
1238           //
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]};
1243
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
1251           lp[5]=-6.48e-05;
1252         }
1253         else {
1254           lp[2]=4.80e-06;      // 0.00219*0.00219
1255           lp[3]=0.0093;        // 0.0964*0.0964;
1256           if (info[0]==1) {
1257             lp[5]=-0.00014;
1258           }
1259           else { 
1260             lp[2]=2.79e-06;    // 0.0017*0.0017; 
1261             lp[3]=0.00935;     // 0.967*0.967;
1262             lp[5]=-4.32e-05;
1263           }
1264         }
1265
1266           AliITSRecPoint * cl2;
1267             cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1268
1269
1270             cl2->SetChargeRatio(ratio);         
1271             fgPairs[ip*nn+jn2]=7;
1272             cl2->SetType(7);
1273             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1274               cl2->SetType(8);
1275               fgPairs[ip*nn+jn2]=8;
1276             }
1277           ncl++;
1278         }
1279         cused1[ip]++;
1280         cused2[jn]++;
1281         cused2[jn2]++;
1282
1283         } // charge matching condition
1284
1285       } // 2 Nside cross 1 Pside
1286     } // loop over Pside clusters
1287
1288   
1289     
1290     for (Int_t ip=0;ip<np;ip++){
1291
1292       if(cused1[ip]) continue;
1293
1294
1295       Float_t xbest=1000,zbest=1000,qbest=0;
1296       //
1297       // 2x2 clusters
1298       //
1299       if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){ 
1300         Float_t minchargediff =4.;
1301         Float_t minchargeratio =0.2;
1302
1303         Int_t j=-1;
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){
1310             j =jc;
1311             minchargediff = TMath::Abs(chargedif);
1312             minchargeratio = TMath::Abs(ratio);
1313           }
1314         }
1315         if (j<0) continue;  // not proper cluster      
1316         
1317
1318         Int_t count =0;
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++;
1323         }
1324         if (count>1) continue;  // more than one "proper" cluster for positive
1325         //
1326         
1327         count =0;
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++;
1332         }
1333         if (count>1) continue;  // more than one "proper" cluster for negative
1334         
1335         Int_t jp = 0;
1336         
1337         count =0;
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++;
1342         }
1343         if (count>1) continue;   
1344         if (fgPairs[ip*nn+j]<100) continue;
1345         //
1346         
1347
1348
1349         //almost gold clusters
1350         Float_t yp=pos[ip].GetY(); 
1351         Float_t yn=neg[j].GetY();      
1352         Float_t xt, zt;
1353         seg->GetPadCxz(yn, yp, xt, zt); 
1354         xbest=xt; zbest=zt; 
1355         qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1356         {
1357           Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1358           mT2L->MasterToLocal(loc,trk);
1359           lp[0]=trk[1];
1360           lp[1]=trk[2];
1361         }
1362         lp[4]=qbest;        //Q
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);
1367         }
1368         //
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]};
1374
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
1382           lp[5]=-6.48e-05;
1383         }
1384         else {
1385           lp[2]=4.80e-06;      // 0.00219*0.00219
1386           lp[3]=0.0093;        // 0.0964*0.0964;
1387           if (info[0]==1) {
1388             lp[5]=-0.00014;
1389           }
1390           else { 
1391             lp[2]=2.79e-06;    // 0.0017*0.0017; 
1392             lp[3]=0.00935;     // 0.967*0.967;
1393             lp[5]=-4.32e-05;
1394           }
1395         }
1396
1397         AliITSRecPoint * cl2;
1398           cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1399                   
1400           cl2->SetChargeRatio(ratio);           
1401           cl2->SetType(10);
1402           fgPairs[ip*nn+j]=10;
1403           if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1404             cl2->SetType(11);
1405             fgPairs[ip*nn+j]=11;
1406           }
1407           cused1[ip]++;
1408           cused2[j]++;      
1409         ncl++;
1410         
1411       } // 2X2
1412     } // loop over Pside 1Dclusters
1413
1414
1415
1416     for (Int_t ip=0;ip<np;ip++){
1417
1418       if(cused1[ip]) continue;
1419
1420
1421       Float_t xbest=1000,zbest=1000,qbest=0;
1422       //
1423       // manyxmany clusters
1424       //
1425       if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){ 
1426         Float_t minchargediff =4.;
1427         Int_t j=-1;
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){
1432             j =jc;
1433             minchargediff = TMath::Abs(chargedif);
1434           }
1435         }
1436         if (j<0) continue;  // not proper cluster      
1437         
1438         Int_t count =0;
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++;
1443         }
1444         if (count>1) continue;  // more than one "proper" cluster for positive
1445         //
1446         
1447         count =0;
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++;
1452         }
1453         if (count>1) continue;  // more than one "proper" cluster for negative
1454         
1455         Int_t jp = 0;
1456         
1457         count =0;
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++;
1462         }
1463         if (count>1) continue;   
1464         if (fgPairs[ip*nn+j]<100) continue;
1465         //
1466         
1467         //almost gold clusters
1468         Float_t yp=pos[ip].GetY(); 
1469         Float_t yn=neg[j].GetY();
1470       
1471
1472         Float_t xt, zt;
1473         seg->GetPadCxz(yn, yp, xt, zt);
1474         
1475         xbest=xt; zbest=zt; 
1476
1477         qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1478
1479         {
1480           Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1481           mT2L->MasterToLocal(loc,trk);
1482           lp[0]=trk[1];
1483           lp[1]=trk[2];
1484         }
1485         lp[4]=qbest;        //Q
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);
1490         }
1491         //
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]};
1497
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
1505           lp[5]=-6.48e-05;
1506         }
1507         else {
1508           lp[2]=4.80e-06;      // 0.00219*0.00219
1509           lp[3]=0.0093;        // 0.0964*0.0964;
1510           if (info[0]==1) {
1511             lp[5]=-0.00014;
1512           }
1513           else { 
1514             lp[2]=2.79e-06;    // 0.0017*0.0017; 
1515             lp[3]=0.00935;     // 0.967*0.967;
1516             lp[5]=-4.32e-05;
1517           }
1518         }
1519         // 
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
1528           }
1529           //
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
1535           }
1536           //
1537           milab[0] = fNClusters+1;  // RS: assign id as cluster label
1538         }
1539
1540         AliITSRecPoint * cl2;
1541           cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1542                   
1543           cl2->SetChargeRatio(ratio);           
1544           cl2->SetType(12);
1545           fgPairs[ip*nn+j]=12;
1546           if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1547             cl2->SetType(13);
1548             fgPairs[ip*nn+j]=13;
1549           }
1550           cused1[ip]++;
1551           cused2[j]++;      
1552           ncl++;
1553           fNClusters++;
1554         
1555       } // manyXmany
1556     } // loop over Pside 1Dclusters
1557     
1558   } // use charge matching
1559   
1560   
1561   // recover all the other crosses
1562   //  
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;
1571
1572       if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1573
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();
1578       
1579       Float_t xt, zt;
1580       seg->GetPadCxz(yn, yp, xt, zt);
1581       
1582       if (TMath::Abs(xt)<hwSSD)
1583       if (TMath::Abs(zt)<hlSSD) {
1584         xbest=xt; zbest=zt; 
1585
1586         qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1587
1588         {
1589         Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1590         mT2L->MasterToLocal(loc,trk);
1591         lp[0]=trk[1];
1592         lp[1]=trk[2];
1593         }
1594         lp[4]=qbest;        //Q
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);
1599         }
1600         //
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]};
1604
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
1612           lp[5]=-6.48e-05;
1613         }
1614         else {
1615           lp[2]=4.80e-06;      // 0.00219*0.00219
1616           lp[3]=0.0093;        // 0.0964*0.0964;
1617           if (info[0]==1) {
1618             lp[5]=-0.00014;
1619           }
1620           else { 
1621             lp[2]=2.79e-06;    // 0.0017*0.0017; 
1622             lp[3]=0.00935;     // 0.967*0.967;
1623             lp[5]=-4.32e-05;
1624           }
1625         }
1626
1627         AliITSRecPoint * cl2;
1628           cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1629
1630           cl2->SetChargeRatio(ratio);
1631           cl2->SetType(100+cpositive[j]+cnegative[i]);    
1632
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]);
1635         ncl++;
1636       }
1637     }
1638   }
1639
1640
1641
1642   if(repa->GetUseBadChannelsInClusterFinderSSD()==kTRUE) {
1643     
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 
1648     
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++;
1654     }
1655     //  AliInfo(Form("module %d has %d P- and %d N-bad strips",fModule,countPbad,countNbad));
1656     
1657     if( (countPbad<100) && (countNbad<100) ) { // no bad side!!
1658       
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;
1661         
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
1664           
1665           if(cal->IsPChannelBad(ib)) { // check if strips is bad
1666             Float_t yN=pos[i].GetY();   
1667             Float_t xt, zt;
1668             seg->GetPadCxz(1.*ib, yN, xt, zt);  
1669             
1670             //----------
1671             // bad Pstrip is crossing the Nside 1Dcluster -> create recpoint
1672             // 
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);
1676               lp[0]=trk[1];
1677               lp[1]=trk[2];        
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]};
1684               
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
1688               if (info[0]>1) {
1689                 lp[2]=4.80e-06;
1690                 lp[3]=0.0093;
1691                 lp[5]=0.00014;
1692               }
1693                       
1694               AliITSRecPoint * cl2;
1695                 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);         
1696                 cl2->SetChargeRatio(1.);
1697               cl2->SetType(50);   
1698               ncl++;
1699             } // cross is within the detector
1700             //
1701             //--------------
1702             
1703           } // bad Pstrip
1704           
1705         } // end loop over Pstrips
1706         
1707       } // end loop over Nside 1D clusters
1708       
1709       for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1710         if(cpositive[j]) continue;
1711         
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
1714           
1715           if(cal->IsNChannelBad(ib)) { // check if strip is bad
1716             Float_t yP=neg[j].GetY();   
1717             Float_t xt, zt;
1718             seg->GetPadCxz(yP, 1.*ib, xt, zt);  
1719             
1720             //----------
1721             // bad Nstrip is crossing the Pside 1Dcluster -> create recpoint
1722             // 
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);
1726               lp[0]=trk[1];
1727               lp[1]=trk[2];        
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]};
1734               
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
1738               if (info[0]>1) {
1739                 lp[2]=2.79e-06;
1740                 lp[3]=0.00935;
1741                 lp[5]=-4.32e-05;
1742               }
1743               
1744               AliITSRecPoint * cl2;
1745                 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);         
1746                 cl2->SetChargeRatio(1.);
1747                 cl2->SetType(60);         
1748               ncl++;
1749             } // cross is within the detector
1750             //
1751             //--------------
1752             
1753           } // bad Nstrip
1754         } // end loop over Nstrips
1755       } // end loop over Pside 1D clusters
1756       
1757     } // no bad sides 
1758     
1759     //---------------------------------------------------------
1760     
1761     else if( (countPbad>700) && (countNbad<100) ) { // bad Pside!!
1762       
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;
1765         
1766         Float_t xt, zt;
1767         Float_t yN=pos[i].GetY();       
1768         Float_t yP=0.;
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); 
1772         
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);
1776           lp[0]=trk[1];
1777           lp[1]=trk[2];        
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]};
1784           
1785           lp[2]=0.00098;    // 0.031*0.031;  //SigmaY2
1786           lp[3]=1.329;      // 1.15*1.15;  //SigmaZ2
1787           lp[5]=-0.0359;
1788           if(info[0]>1) lp[2]=0.00097;
1789
1790           AliITSRecPoint * cl2;
1791             cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);     
1792             cl2->SetChargeRatio(1.);
1793             cl2->SetType(70);     
1794           ncl++;
1795         } // cross is within the detector
1796         //
1797         //--------------
1798         
1799       } // end loop over Nside 1D clusters
1800       
1801     } // bad Pside module
1802     
1803     else if( (countNbad>700) && (countPbad<100) ) { // bad Nside!!
1804       
1805       for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1806         if(cpositive[j]) continue;
1807         
1808         Float_t xt, zt;
1809         Float_t yP=neg[j].GetY();       
1810         Float_t yN=0.;
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); 
1814         
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);
1818           lp[0]=trk[1];
1819           lp[1]=trk[2];        
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]};
1826           
1827           lp[2]=7.27e-05;   // 0.0085*0.0085;  //SigmaY2
1828           lp[3]=1.33;       // 1.15*1.15;  //SigmaZ2
1829           lp[5]=0.00931;
1830           if(info[1]>1) lp[2]=6.91e-05;
1831           
1832           AliITSRecPoint * cl2;
1833             cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);     
1834             cl2->SetChargeRatio(1.);
1835             cl2->SetType(80);     
1836           ncl++;
1837         } // cross is within the detector
1838         //
1839         //--------------
1840         
1841       } // end loop over Pside 1D clusters
1842       
1843     } // bad Nside module
1844     
1845     //---------------------------------------------------------
1846     
1847   } // use bad channels
1848     
1849   //cout<<ncl<<" clusters for this module"<<endl;
1850
1851   delete [] cnegative;
1852   delete [] cused1;
1853   delete [] negativepair;
1854   delete [] cpositive;
1855   delete [] cused2;
1856   delete [] positivepair;
1857
1858 }