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
 518         AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)fDetTypeRec->GetCalibrationModel(fModule);
 519         if( !cal ){
 520           AliWarning(Form("ITSClustersFinderSSD: No calibration found for module (ddl %d ad %d adc %d)",ddl,ad,adc));
 521           continue;
 522         }
 523
 524         Float_t dStrip = 0;
 525
 526         if( repa->GetUseCosmicRunShiftsSSD()) {  // Special condition for 2007/2008 cosmic data
 527           dStrip = fgkCosmic2008StripShifts[ddl][ad-1];
 528           if (TMath::Abs(dStrip) > 1.5){
 529             AliWarning(Form("Indexing error in Cosmic calibration: ddl = %d, dStrip %f\n",ddl,dStrip));
 530             dStrip = 0;
 531           }
 532         }
 533
 534         for( int side=0; side<=1; side++ ){
 535
 536           Int_t lab[3]={-2,-2,-2};
 537           Float_t q = 0.;
 538           Float_t y = 0.;
 539           Int_t nDigits = 0;
 540           Int_t ostrip = -2;
 541           Bool_t snFlag = 0;
 542
 543           Float_t dLorentz = 0;
 544           if (side==0) { // P-side is neg clust
 545             dLorentz = fLorentzShiftN;
 546           }
 547           else { // N-side is pos clust
 548             dLorentz = fLorentzShiftP;
 549           }
 550
 551           Int_t n = nStrips[adc][side];
 552           for( int istr = 0; istr<n+1; istr++ ){
 553
 554             bool stripOK = 1;
 555             Int_t strip=0, rwID = 0;
 556             Float_t signal=0.0, noise=0.0, gain=0.0;
 557
 558             if( istr<n ){
 559               strip = strips[adc][side][istr][0];
 560               signal = strips[adc][side][istr][1];
 561               rwID   = strips[adc][side][istr][2]; // RS
 562               //cout<<"strip "<<adc<<" / "<<side<<": "<<strip<<endl;
 563
 564               if( cal ){
 565                 noise = side ?cal->GetNoiseN(strip) :cal->GetNoiseP(strip);
 566                 gain = side ?cal->GetGainN(strip) :cal->GetGainP(strip);
 567                 stripOK = ( noise>=1. && signal>=3.0*noise
 568                             //&& !cal->IsPChannelBad(strip)
 569                             );
 570               }
 571             } else stripOK = 0; // end of data
 572
 573             bool newCluster = ( (abs(strip-ostrip)!=1) || !stripOK );
 574
 575             if( newCluster ){
 576
 577               //* Store the previous cluster
 578
 579               if( nDigits>0 && q>0 && snFlag ){
 580
 581                 if (nClusters1D[side] >= kMaxADCClusters-1 ) {
 582                   AliWarning("HLT ClustersFinderSSD: Too many 1D clusters !");
 583                 }else {
 584
 585                   Ali1Dcluster &cluster = clusters1D[side][nClusters1D[side]++];
 586                   cluster.SetY( y / q + dStrip + dLorentz);
 587                   cluster.SetQ(q);
 588                   cluster.SetNd(nDigits);
 589                   cluster.SetLabels(lab);
 590                   //cout<<"cluster 1D side "<<side<<": y= "<<y<<" q= "<<q<<" d="<<dStrip<<" Y="<<cluster.GetY()<<endl;
 591                   //Split suspiciously big cluster
 592
 593                   if( repa->GetUseUnfoldingInClusterFinderSSD()
 594                       && nDigits > 4 && nDigits < 25
 595                       ){
 596                     cluster.SetY(y/q + dStrip - 0.25*nDigits + dLorentz);
 597                     cluster.SetQ(0.5*q);
 598                     Ali1Dcluster& cluster2 = clusters1D[side][nClusters1D[side]++];
 599                     cluster2.SetY(y/q + dStrip + 0.25*nDigits + dLorentz);
 600                     cluster2.SetQ(0.5*q);
 601                     cluster2.SetNd(nDigits);
 602                     cluster2.SetLabels(lab);
 603                   } // unfolding is on
 604                 }
 605               }
 606               y = q = 0.;
 607               nDigits = 0;
 608               snFlag = 0;
 609
 610             } //* End store the previous cluster
 611
 612             if( stripOK ){ // add new signal to the cluster
 613 //            signal = (Int_t) (signal * gain); // signal is corrected for gain
 614               if( signal>fgkThreshold*noise) snFlag = 1;
 615               signal = signal * gain; // signal is corrected for gain
 616 //            if( cal ) signal = (Int_t) cal->ADCToKeV( signal ); // signal is  converted in KeV
 617               if( cal ) signal = cal->ADCToKeV( signal ); // signal is  converted in KeV
 618               q += signal;        // add digit to current cluster
 619               y += strip * signal;
 620               nDigits++;
 621               //nstat[side]++;
 622               ostrip = strip;
 623               if (fRawID2ClusID) fRawIDRef[side].AddReference(nClusters1D[side],rwID);
 624
 625             }
 626           } //* end loop over strips
 627
 628         } //* end loop over ADC sides
 629
 630
 631         //* 2D clusterfinder
 632         if( nClusters1D[0] && nClusters1D[1] && fModule>=0 ){
 633            TClonesArray* clusters = rpc->UncheckedGetClusters(fModule);
 634                FindClustersSSD( clusters1D[0], nClusters1D[0], clusters1D[1], nClusters1D[1], clusters);
 635            Int_t nClustersn = clusters->GetEntriesFast();
 636                nClustersSSD += nClustersn;
 637         }
 638
 639         //cout<<"SG: "<<ddl<<" "<<ad<<" "<<adc<<": strips "<<nstat[0]<<"+"<<nstat[1]<<", clusters 1D= "<<nClusters1D[0]<<" + "<<nClusters1D[1]<<", 2D= "<<clusters.size()<<endl;
 640
 641       }//* end loop over adc
 642
 643     }//* end of reconstruction of previous block of 12 modules
 644
 645     if( newModule ){
 646
 647       //*
 648       //* Clean up arrays and set new module
 649       //*
 650
 651       for( int i=0; i<kNADC; i++ ){
 652         nStrips[i][0] = 0;
 653         nStrips[i][1] = 0;
 654       }
 655       ddl = newDDL;
 656       ad = newAD;
 657     }
 658
 659
 660     //*
 661     //* Exit main loop when there is no more input
 662     //*
 663
 664     if( !next ) break;
 665
 666     //*
 667     //* Fill the current strip information
 668     //*
 669
 670     Int_t adc = input->GetADC();
 671     if( adc<0 || adc>=kNADC+2 || (adc>5&&adc<8) ){
 672       AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong adc number (%d)", adc));
 673       continue;
 674     }
 675
 676     if( adc>7 ) adc-= 2; // shift ADC numbers 8-13 to 6-11
 677
 678     Bool_t side = input->GetSideFlag();
 679     Int_t strip = input->GetStrip();
 680     Int_t signal = input->GetSignal();
 681
 682
 683     //cout<<"SSD: "<<ddl<<" "<<ad<<" "<<adc<<" "<<side<<" "<<strip<<" : "<<signal<<endl;
 684
 685     if( strip>767 ){
 686       AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: wrong strip number (ddl %d ad %d adc %d side %d, strip %d",
 687                       ddl, ad, adc, side,strip) );
 688       continue;
 689     }
 690     if (strip < 0) continue;
 691
 692     int &n = nStrips[adc][side];
 693     if( n >0 ){
 694       Int_t oldStrip = strips[adc][side][n-1][0];
 695
 696       if( strip==oldStrip ){
 697         AliWarning(Form("HLT ClustersFinderSSD: Corrupted data: duplicated signal: ddl %d ad %d adc %d, side %d, strip %d",
 698                         ddl, ad, adc, side, strip ));
 699         continue;
 700       }
 701     }
 702     strips[adc][side][n][0] = strip;
 703     strips[adc][side][n][1] = signal;
 704     strips[adc][side][n][2] = countRW;
 705     n++;
 706
 707     //cout<<"SSD: "<<input->GetDDL()<<" "<<input->GetAD()<<" "
 708     //<<input->GetADC()<<" "<<input->GetSideFlag()<<" "<<((int)input->GetStrip())<<" "<<strip<<" : "<<input->GetSignal()<<endl;
 709     //
 710     countRW++; //RS
 711   } //* End main loop over the input
 712
 713   AliDebug(1,Form("found clusters in ITS SSD: %d", nClustersSSD));
 714 }
 715
 716
 717 void AliITSClusterFinderV2SSD::
 718 FindClustersSSD(const Ali1Dcluster* neg, Int_t nn,
 719                 const Ali1Dcluster* pos, Int_t np,
 720                 TClonesArray *clusters) {
 721   //------------------------------------------------------------
 722   // Actual SSD cluster finder
 723   //------------------------------------------------------------
 724
 725   const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
 726
 727   //---------------------------------------
 728   // load recoparam
 729   //
 730   static AliITSRecoParam *repa = NULL;
 731   if(!repa){
 732     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
 733     if(!repa){
 734       repa = AliITSRecoParam::GetHighFluxParam();
 735       AliWarning("Using default AliITSRecoParam class");
 736     }
 737   }
 738
 739 //  TClonesArray &cl=*clusters;
 740
 741   AliITSsegmentationSSD *seg = static_cast<AliITSsegmentationSSD*>(fDetTypeRec->GetSegmentationModel(2));
 742   if (fModule>fLastSSD1)
 743     seg->SetLayer(6);
 744   else
 745     seg->SetLayer(5);
 746
 747   Float_t hwSSD = seg->Dx()*1e-4/2;
 748   Float_t hlSSD = seg->Dz()*1e-4/2;
 749
 750   Int_t idet=fNdet[fModule];
 751   Int_t ncl=0;
 752
 753   //
 754   Int_t *cnegative = new Int_t[np];
 755   Int_t *cused1 = new Int_t[np];
 756   Int_t *negativepair = new Int_t[10*np];
 757   Int_t *cpositive = new Int_t[nn];
 758   Int_t *cused2 = new Int_t[nn];
 759   Int_t *positivepair = new Int_t[10*nn];
 760   for (Int_t i=0;i<np;i++) {cnegative[i]=0; cused1[i]=0;}
 761   for (Int_t i=0;i<nn;i++) {cpositive[i]=0; cused2[i]=0;}
 762   for (Int_t i=0;i<10*np;i++) {negativepair[i]=0;}
 763   for (Int_t i=0;i<10*nn;i++) {positivepair[i]=0;}
 764
 765   if ((np*nn) > fgPairsSize) {
 766
 767     delete [] fgPairs;
 768     fgPairsSize = 2*np*nn;
 769     fgPairs = new Short_t[fgPairsSize];
 770   }
 771   memset(fgPairs,0,sizeof(Short_t)*np*nn);
 772
 773   //
 774   // find available pairs
 775   //
 776   Int_t ncross = 0;
 777   for (Int_t i=0; i<np; i++) {
 778     Float_t yp=pos[i].GetY();
 779     if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
 780     for (Int_t j=0; j<nn; j++) {
 781       if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
 782       Float_t yn=neg[j].GetY();
 783
 784       Float_t xt, zt;
 785       seg->GetPadCxz(yn, yp, xt, zt);
 786       //cout<<yn<<" "<<yp<<" "<<xt<<" "<<zt<<endl;
 787
 788       if (TMath::Abs(xt)<hwSSD)
 789       if (TMath::Abs(zt)<hlSSD) {
 790         Int_t in = i*10+cnegative[i];
 791         Int_t ip = j*10+cpositive[j];
 792         if ((in < 10*np) && (ip < 10*nn)) {
 793           negativepair[in] =j;  //index
 794           positivepair[ip] =i;
 795           cnegative[i]++;  //counters
 796           cpositive[j]++;
 797           ncross++;
 798           fgPairs[i*nn+j]=100;
 799         }
 800         else
 801           AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
 802       }
 803     }
 804   }
 805
 806   if (!ncross) {
 807     delete [] cnegative;
 808     delete [] cused1;
 809     delete [] negativepair;
 810     delete [] cpositive;
 811     delete [] cused2;
 812     delete [] positivepair;
 813     return;
 814   }
 815 //why not to allocate memorey here?  if(!clusters) clusters = new TClonesArray("AliITSRecPoint", ncross);
 816
 817   /* //
 818   // try to recover points out of but close to the module boundaries
 819   //
 820   for (Int_t i=0; i<np; i++) {
 821     Float_t yp=pos[i].GetY();
 822     if ( (pos[i].GetQ()>0) && (pos[i].GetQ()<3) ) continue;
 823     for (Int_t j=0; j<nn; j++) {
 824       if ( (neg[j].GetQ()>0) && (neg[j].GetQ()<3) ) continue;
 825       // if both 1Dclusters have an other cross continue
 826       if (cpositive[j]&&cnegative[i]) continue;
 827       Float_t yn=neg[j].GetY();
 828
 829       Float_t xt, zt;
 830       seg->GetPadCxz(yn, yp, xt, zt);
 831
 832       if (TMath::Abs(xt)<hwSSD+0.1)
 833       if (TMath::Abs(zt)<hlSSD+0.15) {
 834         // tag 1Dcluster (eventually will produce low quality recpoint)
 835         if (cnegative[i]==0) pos[i].SetNd(100);  // not available pair
 836         if (cpositive[j]==0) neg[j].SetNd(100);  // not available pair
 837         Int_t in = i*10+cnegative[i];
 838         Int_t ip = j*10+cpositive[j];
 839         if ((in < 10*np) && (ip < 10*nn)) {
 840           negativepair[in] =j;  //index
 841           positivepair[ip] =i;
 842           cnegative[i]++;  //counters
 843           cpositive[j]++;
 844           fgPairs[i*nn+j]=100;
 845         }
 846         else
 847           AliError(Form("Index out of range: ip=%d, in=%d",ip,in));
 848       }
 849     }
 850   }
 851   */
 852
 853   //
 854   Float_t lp[6];
 855   Int_t milab[10];
 856   Double_t ratio;
 857
 858
 859   if(repa->GetUseChargeMatchingInClusterFinderSSD()==kTRUE) {
 860
 861
 862     //
 863     // sign gold tracks
 864     //
 865     for (Int_t ip=0;ip<np;ip++){
 866       Float_t xbest=1000,zbest=1000,qbest=0;
 867       //
 868       // select gold clusters
 869       if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
 870         Float_t yp=pos[ip].GetY();
 871         Int_t j = negativepair[10*ip];
 872
 873         if( (pos[ip].GetQ()==0) && (neg[j].GetQ() ==0) ) {
 874           // both bad, hence continue;
 875           // mark both as used (to avoid recover at the end)
 876           cused1[ip]++;
 877           cused2[j]++;
 878           continue;
 879         }
 880
 881         ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
 882         //cout<<"ratio="<<ratio<<endl;
 883
 884         // charge matching (note that if posQ or negQ is 0 -> ratio=1 and the following condition is met
 885         if (TMath::Abs(ratio)>0.2) continue; // note: 0.2=3xsigma_ratio calculated in cosmics tests
 886
 887         //
 888         Float_t yn=neg[j].GetY();
 889
 890         Float_t xt, zt;
 891         seg->GetPadCxz(yn, yp, xt, zt);
 892
 893         xbest=xt; zbest=zt;
 894
 895
 896         qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
 897         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
 898
 899         {
 900           Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
 901           mT2L->MasterToLocal(loc,trk);
 902           lp[0]=trk[1];
 903           lp[1]=trk[2];
 904         }
 905         lp[4]=qbest;        //Q
 906         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 907         for (Int_t ilab=0;ilab<3;ilab++){
 908           milab[ilab] = pos[ip].GetLabel(ilab);
 909           milab[ilab+3] = neg[j].GetLabel(ilab);
 910         }
 911         //
 912         CheckLabels2(milab);
 913         milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
 914         Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
 915
 916         lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
 917         lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
 918         // out-of-diagonal element of covariance matrix
 919         if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
 920         else if ( (info[0]>1) && (info[1]>1) ) {
 921           lp[2]=2.63e-06;    // 0.0016*0.0016;  //SigmaY2
 922           lp[3]=0.0065;      // 0.08*0.08;   //SigmaZ2
 923           lp[5]=-6.48e-05;
 924         }
 925         else {
 926           lp[2]=4.80e-06;      // 0.00219*0.00219
 927           lp[3]=0.0093;        // 0.0964*0.0964;
 928           if (info[0]==1) {
 929             lp[5]=-0.00014;
 930           }
 931           else {
 932             lp[2]=2.79e-06;    // 0.0017*0.0017;
 933             lp[3]=0.00935;     // 0.967*0.967;
 934             lp[5]=-4.32e-05;
 935           }
 936         }
 937
 938         AliITSRecPoint * cl2;
 939         cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
 940
 941     cl2->SetChargeRatio(ratio);
 942     cl2->SetType(1);
 943     fgPairs[ip*nn+j]=1;
 944
 945     if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
 946       cl2->SetType(2);
 947       fgPairs[ip*nn+j]=2;
 948     }
 949
 950     if(pos[ip].GetQ()==0) cl2->SetType(3);
 951     if(neg[j].GetQ()==0) cl2->SetType(4);
 952
 953     cused1[ip]++;
 954     cused2[j]++;
 955
 956         ncl++;
 957       }
 958     }
 959
 960     for (Int_t ip=0;ip<np;ip++){
 961       Float_t xbest=1000,zbest=1000,qbest=0;
 962       //
 963       //
 964       // select "silber" cluster
 965       if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
 966         Int_t in  = negativepair[10*ip];
 967         Int_t ip2 = positivepair[10*in];
 968         if (ip2==ip) ip2 =  positivepair[10*in+1];
 969         Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
 970
 971
 972
 973         ratio = (pcharge-neg[in].GetQ())/(pcharge+neg[in].GetQ());
 974         if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
 975           //if ( (TMath::Abs(pcharge-neg[in].GetQ())<30) && (pcharge!=0) ) { //
 976
 977           //
 978           // add first pair
 979           if ( (fgPairs[ip*nn+in]==100)&&(pos[ip].GetQ() ) ) {  //
 980
 981             Float_t yp=pos[ip].GetY();
 982             Float_t yn=neg[in].GetY();
 983
 984             Float_t xt, zt;
 985             seg->GetPadCxz(yn, yp, xt, zt);
 986
 987             xbest=xt; zbest=zt;
 988
 989             qbest =pos[ip].GetQ();
 990             Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
 991             mT2L->MasterToLocal(loc,trk);
 992             lp[0]=trk[1];
 993             lp[1]=trk[2];
 994
 995             lp[4]=qbest;        //Q
 996             for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
 997             for (Int_t ilab=0;ilab<3;ilab++){
 998               milab[ilab] = pos[ip].GetLabel(ilab);
 999               milab[ilab+3] = neg[in].GetLabel(ilab);
1000             }
1001             //
1002             CheckLabels2(milab);
1003             ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
1004             milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
1005             Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1006
1007         lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
1008         lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
1009         // out-of-diagonal element of covariance matrix
1010         if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1011         else if ( (info[0]>1) && (info[1]>1) ) {
1012           lp[2]=2.63e-06;    // 0.0016*0.0016;  //SigmaY2
1013           lp[3]=0.0065;      // 0.08*0.08;   //SigmaZ2
1014           lp[5]=-6.48e-05;
1015         }
1016         else {
1017           lp[2]=4.80e-06;      // 0.00219*0.00219
1018           lp[3]=0.0093;        // 0.0964*0.0964;
1019           if (info[0]==1) {
1020             lp[5]=-0.00014;
1021           }
1022           else {
1023             lp[2]=2.79e-06;    // 0.0017*0.0017;
1024             lp[3]=0.00935;     // 0.967*0.967;
1025             lp[5]=-4.32e-05;
1026           }
1027         }
1028
1029         AliITSRecPoint * cl2;
1030               cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1031               cl2->SetChargeRatio(ratio);
1032               cl2->SetType(5);
1033               fgPairs[ip*nn+in] = 5;
1034               if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
1035                 cl2->SetType(6);
1036                 fgPairs[ip*nn+in] = 6;
1037               }
1038             ncl++;
1039           }
1040
1041
1042           //
1043           // add second pair
1044
1045           //    if (!(cused1[ip2] || cused2[in])){  //
1046           if ( (fgPairs[ip2*nn+in]==100) && (pos[ip2].GetQ()) ) {
1047
1048             Float_t yp=pos[ip2].GetY();
1049             Float_t yn=neg[in].GetY();
1050
1051             Float_t xt, zt;
1052             seg->GetPadCxz(yn, yp, xt, zt);
1053
1054             xbest=xt; zbest=zt;
1055
1056             qbest =pos[ip2].GetQ();
1057
1058             Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1059             mT2L->MasterToLocal(loc,trk);
1060             lp[0]=trk[1];
1061             lp[1]=trk[2];
1062
1063             lp[4]=qbest;        //Q
1064             for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1065             for (Int_t ilab=0;ilab<3;ilab++){
1066               milab[ilab] = pos[ip2].GetLabel(ilab);
1067               milab[ilab+3] = neg[in].GetLabel(ilab);
1068             }
1069             //
1070             CheckLabels2(milab);
1071             ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
1072             milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
1073             Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fModule]};
1074
1075         lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
1076         lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
1077         // out-of-diagonal element of covariance matrix
1078         if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1079         else if ( (info[0]>1) && (info[1]>1) ) {
1080           lp[2]=2.63e-06;    // 0.0016*0.0016;  //SigmaY2
1081           lp[3]=0.0065;      // 0.08*0.08;   //SigmaZ2
1082           lp[5]=-6.48e-05;
1083         }
1084         else {
1085           lp[2]=4.80e-06;      // 0.00219*0.00219
1086           lp[3]=0.0093;        // 0.0964*0.0964;
1087           if (info[0]==1) {
1088             lp[5]=-0.00014;
1089           }
1090           else {
1091             lp[2]=2.79e-06;    // 0.0017*0.0017;
1092             lp[3]=0.00935;     // 0.967*0.967;
1093             lp[5]=-4.32e-05;
1094           }
1095         }
1096
1097             AliITSRecPoint * cl2;
1098               cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1099
1100               cl2->SetChargeRatio(ratio);
1101               cl2->SetType(5);
1102               fgPairs[ip2*nn+in] =5;
1103               if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
1104                 cl2->SetType(6);
1105                 fgPairs[ip2*nn+in] =6;
1106               }
1107             ncl++;
1108           }
1109
1110           cused1[ip]++;
1111           cused1[ip2]++;
1112           cused2[in]++;
1113
1114         } // charge matching condition
1115
1116       } // 2 Pside cross 1 Nside
1117     } // loop over Pside clusters
1118
1119
1120
1121       //
1122     for (Int_t jn=0;jn<nn;jn++){
1123       if (cused2[jn]) continue;
1124       Float_t xbest=1000,zbest=1000,qbest=0;
1125       // select "silber" cluster
1126       if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
1127         Int_t ip  = positivepair[10*jn];
1128         Int_t jn2 = negativepair[10*ip];
1129         if (jn2==jn) jn2 =  negativepair[10*ip+1];
1130         Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
1131         //
1132
1133
1134         ratio = (pcharge-pos[ip].GetQ())/(pcharge+pos[ip].GetQ());
1135         if ( (TMath::Abs(ratio)<0.2) && (pcharge!=0) ) {
1136
1137           /*
1138         if ( (TMath::Abs(pcharge-pos[ip].GetQ())<30) &&  // charge matching
1139              (pcharge!=0) ) { // reject combinations of bad strips
1140           */
1141
1142
1143           //
1144           // add first pair
1145           //    if (!(cused1[ip]||cused2[jn])){
1146           if ( (fgPairs[ip*nn+jn]==100) && (neg[jn].GetQ()) ) {  //
1147
1148             Float_t yn=neg[jn].GetY();
1149             Float_t yp=pos[ip].GetY();
1150
1151             Float_t xt, zt;
1152             seg->GetPadCxz(yn, yp, xt, zt);
1153
1154             xbest=xt; zbest=zt;
1155
1156             qbest =neg[jn].GetQ();
1157
1158             {
1159               Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1160               mT2L->MasterToLocal(loc,trk);
1161               lp[0]=trk[1];
1162               lp[1]=trk[2];
1163           }
1164
1165           lp[4]=qbest;        //Q
1166           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1167           for (Int_t ilab=0;ilab<3;ilab++){
1168             milab[ilab] = pos[ip].GetLabel(ilab);
1169             milab[ilab+3] = neg[jn].GetLabel(ilab);
1170           }
1171           //
1172           CheckLabels2(milab);
1173           ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
1174           milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
1175           Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fModule]};
1176
1177         lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
1178         lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
1179         // out-of-diagonal element of covariance matrix
1180         if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1181         else if ( (info[0]>1) && (info[1]>1) ) {
1182           lp[2]=2.63e-06;    // 0.0016*0.0016;  //SigmaY2
1183           lp[3]=0.0065;      // 0.08*0.08;   //SigmaZ2
1184           lp[5]=-6.48e-05;
1185         }
1186         else {
1187           lp[2]=4.80e-06;      // 0.00219*0.00219
1188           lp[3]=0.0093;        // 0.0964*0.0964;
1189           if (info[0]==1) {
1190             lp[5]=-0.00014;
1191           }
1192           else {
1193             lp[2]=2.79e-06;    // 0.0017*0.0017;
1194             lp[3]=0.00935;     // 0.967*0.967;
1195             lp[5]=-4.32e-05;
1196           }
1197         }
1198
1199           AliITSRecPoint * cl2;
1200             cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1201
1202             cl2->SetChargeRatio(ratio);
1203             cl2->SetType(7);
1204             fgPairs[ip*nn+jn] =7;
1205             if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
1206               cl2->SetType(8);
1207               fgPairs[ip*nn+jn]=8;
1208             }
1209           ncl++;
1210         }
1211         //
1212         // add second pair
1213         //      if (!(cused1[ip]||cused2[jn2])){
1214         if ( (fgPairs[ip*nn+jn2]==100)&&(neg[jn2].GetQ() ) ) {  //
1215
1216           Float_t yn=neg[jn2].GetY();
1217           Double_t yp=pos[ip].GetY();
1218
1219           Float_t xt, zt;
1220           seg->GetPadCxz(yn, yp, xt, zt);
1221
1222           xbest=xt; zbest=zt;
1223
1224           qbest =neg[jn2].GetQ();
1225
1226           {
1227           Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1228           mT2L->MasterToLocal(loc,trk);
1229           lp[0]=trk[1];
1230           lp[1]=trk[2];
1231           }
1232
1233           lp[4]=qbest;        //Q
1234           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1235           for (Int_t ilab=0;ilab<3;ilab++){
1236             milab[ilab] = pos[ip].GetLabel(ilab);
1237             milab[ilab+3] = neg[jn2].GetLabel(ilab);
1238           }
1239           //
1240           CheckLabels2(milab);
1241           ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
1242           milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
1243           Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fModule]};
1244
1245         lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
1246         lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
1247         // out-of-diagonal element of covariance matrix
1248         if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1249         else if ( (info[0]>1) && (info[1]>1) ) {
1250           lp[2]=2.63e-06;    // 0.0016*0.0016;  //SigmaY2
1251           lp[3]=0.0065;      // 0.08*0.08;   //SigmaZ2
1252           lp[5]=-6.48e-05;
1253         }
1254         else {
1255           lp[2]=4.80e-06;      // 0.00219*0.00219
1256           lp[3]=0.0093;        // 0.0964*0.0964;
1257           if (info[0]==1) {
1258             lp[5]=-0.00014;
1259           }
1260           else {
1261             lp[2]=2.79e-06;    // 0.0017*0.0017;
1262             lp[3]=0.00935;     // 0.967*0.967;
1263             lp[5]=-4.32e-05;
1264           }
1265         }
1266
1267           AliITSRecPoint * cl2;
1268             cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1269
1270
1271             cl2->SetChargeRatio(ratio);
1272             fgPairs[ip*nn+jn2]=7;
1273             cl2->SetType(7);
1274             if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
1275               cl2->SetType(8);
1276               fgPairs[ip*nn+jn2]=8;
1277             }
1278           ncl++;
1279         }
1280         cused1[ip]++;
1281         cused2[jn]++;
1282         cused2[jn2]++;
1283
1284         } // charge matching condition
1285
1286       } // 2 Nside cross 1 Pside
1287     } // loop over Pside clusters
1288
1289
1290
1291     for (Int_t ip=0;ip<np;ip++){
1292
1293       if(cused1[ip]) continue;
1294
1295
1296       Float_t xbest=1000,zbest=1000,qbest=0;
1297       //
1298       // 2x2 clusters
1299       //
1300       if ( (cnegative[ip]==2) && cpositive[negativepair[10*ip]]==2){
1301         Float_t minchargediff =4.;
1302         //      Float_t minchargeratio =0.2;
1303
1304         Int_t j=-1;
1305         for (Int_t di=0;di<cnegative[ip];di++){
1306           Int_t   jc = negativepair[ip*10+di];
1307           Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1308           ratio = (pos[ip].GetQ()-neg[jc].GetQ())/(pos[ip].GetQ()+neg[jc].GetQ());
1309           //if (TMath::Abs(chargedif)<minchargediff){
1310           if (TMath::Abs(ratio)<0.2){
1311             j =jc;
1312             minchargediff = TMath::Abs(chargedif);
1313             //      minchargeratio = TMath::Abs(ratio);
1314           }
1315         }
1316         if (j<0) continue;  // not proper cluster
1317
1318
1319         Int_t count =0;
1320         for (Int_t di=0;di<cnegative[ip];di++){
1321           Int_t   jc = negativepair[ip*10+di];
1322           Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1323           if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1324         }
1325         if (count>1) continue;  // more than one "proper" cluster for positive
1326         //
1327
1328         count =0;
1329         for (Int_t dj=0;dj<cpositive[j];dj++){
1330           Int_t   ic  = positivepair[j*10+dj];
1331           Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1332           if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1333         }
1334         if (count>1) continue;  // more than one "proper" cluster for negative
1335
1336         Int_t jp = 0;
1337
1338         count =0;
1339         for (Int_t dj=0;dj<cnegative[jp];dj++){
1340           Int_t   ic = positivepair[jp*10+dj];
1341           Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1342           if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1343         }
1344         if (count>1) continue;
1345         if (fgPairs[ip*nn+j]<100) continue;
1346         //
1347
1348
1349
1350         //almost gold clusters
1351         Float_t yp=pos[ip].GetY();
1352         Float_t yn=neg[j].GetY();
1353         Float_t xt, zt;
1354         seg->GetPadCxz(yn, yp, xt, zt);
1355         xbest=xt; zbest=zt;
1356         qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1357         {
1358           Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1359           mT2L->MasterToLocal(loc,trk);
1360           lp[0]=trk[1];
1361           lp[1]=trk[2];
1362         }
1363         lp[4]=qbest;        //Q
1364         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1365         for (Int_t ilab=0;ilab<3;ilab++){
1366           milab[ilab] = pos[ip].GetLabel(ilab);
1367           milab[ilab+3] = neg[j].GetLabel(ilab);
1368         }
1369         //
1370         CheckLabels2(milab);
1371         if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1372         ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1373         milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1374         Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1375
1376         lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
1377         lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
1378         // out-of-diagonal element of covariance matrix
1379         if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1380         else if ( (info[0]>1) && (info[1]>1) ) {
1381           lp[2]=2.63e-06;    // 0.0016*0.0016;  //SigmaY2
1382           lp[3]=0.0065;      // 0.08*0.08;   //SigmaZ2
1383           lp[5]=-6.48e-05;
1384         }
1385         else {
1386           lp[2]=4.80e-06;      // 0.00219*0.00219
1387           lp[3]=0.0093;        // 0.0964*0.0964;
1388           if (info[0]==1) {
1389             lp[5]=-0.00014;
1390           }
1391           else {
1392             lp[2]=2.79e-06;    // 0.0017*0.0017;
1393             lp[3]=0.00935;     // 0.967*0.967;
1394             lp[5]=-4.32e-05;
1395           }
1396         }
1397
1398         AliITSRecPoint * cl2;
1399           cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1400
1401           cl2->SetChargeRatio(ratio);
1402           cl2->SetType(10);
1403           fgPairs[ip*nn+j]=10;
1404           if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1405             cl2->SetType(11);
1406             fgPairs[ip*nn+j]=11;
1407           }
1408           cused1[ip]++;
1409           cused2[j]++;
1410         ncl++;
1411
1412       } // 2X2
1413     } // loop over Pside 1Dclusters
1414
1415
1416
1417     for (Int_t ip=0;ip<np;ip++){
1418
1419       if(cused1[ip]) continue;
1420
1421
1422       Float_t xbest=1000,zbest=1000,qbest=0;
1423       //
1424       // manyxmany clusters
1425       //
1426       if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
1427         Float_t minchargediff =4.;
1428         Int_t j=-1;
1429         for (Int_t di=0;di<cnegative[ip];di++){
1430           Int_t   jc = negativepair[ip*10+di];
1431           Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1432           if (TMath::Abs(chargedif)<minchargediff){
1433             j =jc;
1434             minchargediff = TMath::Abs(chargedif);
1435           }
1436         }
1437         if (j<0) continue;  // not proper cluster
1438
1439         Int_t count =0;
1440         for (Int_t di=0;di<cnegative[ip];di++){
1441           Int_t   jc = negativepair[ip*10+di];
1442           Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
1443           if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1444         }
1445         if (count>1) continue;  // more than one "proper" cluster for positive
1446         //
1447
1448         count =0;
1449         for (Int_t dj=0;dj<cpositive[j];dj++){
1450           Int_t   ic  = positivepair[j*10+dj];
1451           Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
1452           if (TMath::Abs(chargedif)<minchargediff+3.) count++;
1453         }
1454         if (count>1) continue;  // more than one "proper" cluster for negative
1455
1456         Int_t jp = 0;
1457
1458         count =0;
1459         for (Int_t dj=0;dj<cnegative[jp];dj++){
1460           Int_t   ic = positivepair[jp*10+dj];
1461           Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
1462           if (TMath::Abs(chargedif)<minchargediff+4.) count++;
1463         }
1464         if (count>1) continue;
1465         if (fgPairs[ip*nn+j]<100) continue;
1466         //
1467
1468         //almost gold clusters
1469         Float_t yp=pos[ip].GetY();
1470         Float_t yn=neg[j].GetY();
1471
1472
1473         Float_t xt, zt;
1474         seg->GetPadCxz(yn, yp, xt, zt);
1475
1476         xbest=xt; zbest=zt;
1477
1478         qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
1479
1480         {
1481           Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1482           mT2L->MasterToLocal(loc,trk);
1483           lp[0]=trk[1];
1484           lp[1]=trk[2];
1485         }
1486         lp[4]=qbest;        //Q
1487         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1488         for (Int_t ilab=0;ilab<3;ilab++){
1489           milab[ilab] = pos[ip].GetLabel(ilab);
1490           milab[ilab+3] = neg[j].GetLabel(ilab);
1491         }
1492         //
1493         CheckLabels2(milab);
1494         if ((neg[j].GetQ()==0)&&(pos[ip].GetQ()==0)) continue; // reject crosses of bad strips!!
1495         ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
1496         milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
1497         Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1498
1499         lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
1500         lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
1501         // out-of-diagonal element of covariance matrix
1502         if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1503         else if ( (info[0]>1) && (info[1]>1) ) {
1504           lp[2]=2.63e-06;    // 0.0016*0.0016;  //SigmaY2
1505           lp[3]=0.0065;      // 0.08*0.08;   //SigmaZ2
1506           lp[5]=-6.48e-05;
1507         }
1508         else {
1509           lp[2]=4.80e-06;      // 0.00219*0.00219
1510           lp[3]=0.0093;        // 0.0964*0.0964;
1511           if (info[0]==1) {
1512             lp[5]=-0.00014;
1513           }
1514           else {
1515             lp[2]=2.79e-06;    // 0.0017*0.0017;
1516             lp[3]=0.00935;     // 0.967*0.967;
1517             lp[5]=-4.32e-05;
1518           }
1519         }
1520         //
1521         if (fRawID2ClusID) { // set rawID <-> clusterID correspondence for embedding
1522           const int kMaxRefRW = 200;
1523           UInt_t nrefsRW,refsRW[kMaxRefRW];
1524           nrefsRW = fRawIDRef[0].GetReferences(j,refsRW,kMaxRefRW); // n-side
1525           for (int ir=nrefsRW;ir--;) {
1526             int rwid = (int)refsRW[ir];
1527             if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
1528             (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
1529           }
1530           //
1531           nrefsRW = fRawIDRef[1].GetReferences(ip,refsRW,kMaxRefRW); // p-side
1532           for (int ir=nrefsRW;ir--;) {
1533             int rwid = (int)refsRW[ir];
1534             if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
1535             (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
1536           }
1537           //
1538           milab[0] = fNClusters+1;  // RS: assign id as cluster label
1539         }
1540
1541         AliITSRecPoint * cl2;
1542           cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1543
1544           cl2->SetChargeRatio(ratio);
1545           cl2->SetType(12);
1546           fgPairs[ip*nn+j]=12;
1547           if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
1548             cl2->SetType(13);
1549             fgPairs[ip*nn+j]=13;
1550           }
1551           cused1[ip]++;
1552           cused2[j]++;
1553           ncl++;
1554           fNClusters++;
1555
1556       } // manyXmany
1557     } // loop over Pside 1Dclusters
1558
1559   } // use charge matching
1560
1561
1562   // recover all the other crosses
1563   //
1564   for (Int_t i=0; i<np; i++) {
1565     Float_t xbest=1000,zbest=1000,qbest=0;
1566     Float_t yp=pos[i].GetY();
1567     if ((pos[i].GetQ()>0)&&(pos[i].GetQ()<3)) continue;
1568     for (Int_t j=0; j<nn; j++) {
1569     //    for (Int_t di = 0;di<cpositive[i];di++){
1570     //  Int_t j = negativepair[10*i+di];
1571       if ((neg[j].GetQ()>0)&&(neg[j].GetQ()<3)) continue;
1572
1573       if ((neg[j].GetQ()==0)&&(pos[i].GetQ()==0)) continue; // reject crosses of bad strips!!
1574
1575       if (cused2[j]||cused1[i]) continue;
1576       if (fgPairs[i*nn+j]>0 &&fgPairs[i*nn+j]<100) continue;
1577       ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ());
1578       Float_t yn=neg[j].GetY();
1579
1580       Float_t xt, zt;
1581       seg->GetPadCxz(yn, yp, xt, zt);
1582
1583       if (TMath::Abs(xt)<hwSSD)
1584       if (TMath::Abs(zt)<hlSSD) {
1585         xbest=xt; zbest=zt;
1586
1587         qbest=0.5*(pos[i].GetQ()+neg[j].GetQ());
1588
1589         {
1590         Double_t loc[3]={xbest,0.,zbest},trk[3]={0.,0.,0.};
1591         mT2L->MasterToLocal(loc,trk);
1592         lp[0]=trk[1];
1593         lp[1]=trk[2];
1594         }
1595         lp[4]=qbest;        //Q
1596         for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1597         for (Int_t ilab=0;ilab<3;ilab++){
1598           milab[ilab] = pos[i].GetLabel(ilab);
1599           milab[ilab+3] = neg[j].GetLabel(ilab);
1600         }
1601         //
1602         CheckLabels2(milab);
1603         milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det
1604         Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fModule]};
1605
1606         lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
1607         lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
1608         // out-of-diagonal element of covariance matrix
1609         if( (info[0]==1) && (info[1]==1) ) lp[5]=-0.00012;
1610         else if ( (info[0]>1) && (info[1]>1) ) {
1611           lp[2]=2.63e-06;    // 0.0016*0.0016;  //SigmaY2
1612           lp[3]=0.0065;      // 0.08*0.08;   //SigmaZ2
1613           lp[5]=-6.48e-05;
1614         }
1615         else {
1616           lp[2]=4.80e-06;      // 0.00219*0.00219
1617           lp[3]=0.0093;        // 0.0964*0.0964;
1618           if (info[0]==1) {
1619             lp[5]=-0.00014;
1620           }
1621           else {
1622             lp[2]=2.79e-06;    // 0.0017*0.0017;
1623             lp[3]=0.00935;     // 0.967*0.967;
1624             lp[5]=-4.32e-05;
1625           }
1626         }
1627
1628         AliITSRecPoint * cl2;
1629           cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1630
1631           cl2->SetChargeRatio(ratio);
1632           cl2->SetType(100+cpositive[j]+cnegative[i]);
1633
1634           if(pos[i].GetQ()==0) cl2->SetType(200+cpositive[j]+cnegative[i]);
1635           if(neg[j].GetQ()==0) cl2->SetType(300+cpositive[j]+cnegative[i]);
1636         ncl++;
1637       }
1638     }
1639   }
1640
1641
1642
1643   if(repa->GetUseBadChannelsInClusterFinderSSD()==kTRUE) {
1644
1645     //---------------------------------------------------------
1646     // recover crosses of good 1D clusters with bad strips on the other side
1647     // Note1: at first iteration skip modules with a bad side (or almost), (would produce too many fake!)
1648     // Note2: for modules with a bad side see below
1649
1650     AliITSCalibrationSSD* cal = (AliITSCalibrationSSD*)GetResp(fModule);
1651     Int_t countPbad=0, countNbad=0;
1652     for(Int_t ib=0; ib<768; ib++) {
1653       if(cal->IsPChannelBad(ib)) countPbad++;
1654       if(cal->IsNChannelBad(ib)) countNbad++;
1655     }
1656     //  AliInfo(Form("module %d has %d P- and %d N-bad strips",fModule,countPbad,countNbad));
1657
1658     if( (countPbad<100) && (countNbad<100) ) { // no bad side!!
1659
1660       for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1661         if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1662
1663         //      for(Int_t ib=0; ib<768; ib++) { // loop over all Pstrips
1664         for(Int_t ib=15; ib<753; ib++) { // loop over all Pstrips
1665
1666           if(cal->IsPChannelBad(ib)) { // check if strips is bad
1667             Float_t yN=pos[i].GetY();
1668             Float_t xt, zt;
1669             seg->GetPadCxz(1.*ib, yN, xt, zt);
1670
1671             //----------
1672             // bad Pstrip is crossing the Nside 1Dcluster -> create recpoint
1673             //
1674             if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1675               Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1676               mT2L->MasterToLocal(loc,trk);
1677               lp[0]=trk[1];
1678               lp[1]=trk[2];
1679               lp[4]=pos[i].GetQ(); //Q
1680               for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1681               for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1682               CheckLabels2(milab);
1683               milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1684               Int_t info[3] = {pos[i].GetNd(),0,fNlayer[fModule]};
1685
1686               lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
1687               lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
1688               lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1689               if (info[0]>1) {
1690                 lp[2]=4.80e-06;
1691                 lp[3]=0.0093;
1692                 lp[5]=0.00014;
1693               }
1694
1695               AliITSRecPoint * cl2;
1696                 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1697                 cl2->SetChargeRatio(1.);
1698               cl2->SetType(50);
1699               ncl++;
1700             } // cross is within the detector
1701             //
1702             //--------------
1703
1704           } // bad Pstrip
1705
1706         } // end loop over Pstrips
1707
1708       } // end loop over Nside 1D clusters
1709
1710       for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1711         if(cpositive[j]) continue;
1712
1713         //      for(Int_t ib=0; ib<768; ib++) { // loop over all Nside strips
1714         for(Int_t ib=15; ib<753; ib++) { // loop over all Nside strips
1715
1716           if(cal->IsNChannelBad(ib)) { // check if strip is bad
1717             Float_t yP=neg[j].GetY();
1718             Float_t xt, zt;
1719             seg->GetPadCxz(yP, 1.*ib, xt, zt);
1720
1721             //----------
1722             // bad Nstrip is crossing the Pside 1Dcluster -> create recpoint
1723             //
1724             if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1725               Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1726               mT2L->MasterToLocal(loc,trk);
1727               lp[0]=trk[1];
1728               lp[1]=trk[2];
1729               lp[4]=neg[j].GetQ(); //Q
1730               for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1731               for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1732               CheckLabels2(milab);
1733               milab[3]=( j << 10 ) + idet; // pos|neg|det
1734               Int_t info[3]={0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1735
1736               lp[2]=4.968e-06;     // 0.00223*0.00223;  //SigmaY2
1737               lp[3]=0.012;         // 0.110*0.110;  //SigmaZ2
1738               lp[5]=-0.00012; // out-of-diagonal element of covariance matrix
1739               if (info[0]>1) {
1740                 lp[2]=2.79e-06;
1741                 lp[3]=0.00935;
1742                 lp[5]=-4.32e-05;
1743               }
1744
1745               AliITSRecPoint * cl2;
1746                 cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1747                 cl2->SetChargeRatio(1.);
1748                 cl2->SetType(60);
1749               ncl++;
1750             } // cross is within the detector
1751             //
1752             //--------------
1753
1754           } // bad Nstrip
1755         } // end loop over Nstrips
1756       } // end loop over Pside 1D clusters
1757
1758     } // no bad sides
1759
1760     //---------------------------------------------------------
1761
1762     else if( (countPbad>700) && (countNbad<100) ) { // bad Pside!!
1763
1764       for (Int_t i=0; i<np; i++) { // loop over Nside 1Dclusters with no crosses
1765         if(cnegative[i]) continue; // if intersecting Pside clusters continue;
1766
1767         Float_t xt, zt;
1768         Float_t yN=pos[i].GetY();
1769         Float_t yP=0.;
1770         if (seg->GetLayer()==5) yP = yN + (7.6/1.9);
1771         else yP = yN - (7.6/1.9);
1772         seg->GetPadCxz(yP, yN, xt, zt);
1773
1774         if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1775           Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1776           mT2L->MasterToLocal(loc,trk);
1777           lp[0]=trk[1];
1778           lp[1]=trk[2];
1779           lp[4]=pos[i].GetQ(); //Q
1780           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1781           for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = pos[i].GetLabel(ilab);
1782           CheckLabels2(milab);
1783           milab[3]=( (i<<10) << 10 ) + idet; // pos|neg|det
1784           Int_t info[3] = {(Int_t)pos[i].GetNd(),0,fNlayer[fModule]};
1785
1786           lp[2]=0.00098;    // 0.031*0.031;  //SigmaY2
1787           lp[3]=1.329;      // 1.15*1.15;  //SigmaZ2
1788           lp[5]=-0.0359;
1789           if(info[0]>1) lp[2]=0.00097;
1790
1791           AliITSRecPoint * cl2;
1792             cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1793             cl2->SetChargeRatio(1.);
1794             cl2->SetType(70);
1795           ncl++;
1796         } // cross is within the detector
1797         //
1798         //--------------
1799
1800       } // end loop over Nside 1D clusters
1801
1802     } // bad Pside module
1803
1804     else if( (countNbad>700) && (countPbad<100) ) { // bad Nside!!
1805
1806       for (Int_t j=0; j<nn; j++) { // loop over Pside 1D clusters with no crosses
1807         if(cpositive[j]) continue;
1808
1809         Float_t xt, zt;
1810         Float_t yP=neg[j].GetY();
1811         Float_t yN=0.;
1812         if (seg->GetLayer()==5) yN = yP - (7.6/1.9);
1813         else yN = yP + (7.6/1.9);
1814         seg->GetPadCxz(yP, yN, xt, zt);
1815
1816         if ( (TMath::Abs(xt)<hwSSD) && (TMath::Abs(zt)<hlSSD) ) {
1817           Double_t loc[3]={xt,0.,zt},trk[3]={0.,0.,0.};
1818           mT2L->MasterToLocal(loc,trk);
1819           lp[0]=trk[1];
1820           lp[1]=trk[2];
1821           lp[4]=neg[j].GetQ(); //Q
1822           for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
1823           for (Int_t ilab=0;ilab<3;ilab++) milab[ilab] = neg[j].GetLabel(ilab);
1824           CheckLabels2(milab);
1825           milab[3]=( j << 10 ) + idet; // pos|neg|det
1826           Int_t info[3] = {0,(Int_t)neg[j].GetNd(),fNlayer[fModule]};
1827
1828           lp[2]=7.27e-05;   // 0.0085*0.0085;  //SigmaY2
1829           lp[3]=1.33;       // 1.15*1.15;  //SigmaZ2
1830           lp[5]=0.00931;
1831           if(info[1]>1) lp[2]=6.91e-05;
1832
1833           AliITSRecPoint * cl2;
1834             cl2 = new ((*clusters)[ncl]) AliITSRecPoint(milab,lp,info);
1835             cl2->SetChargeRatio(1.);
1836             cl2->SetType(80);
1837           ncl++;
1838         } // cross is within the detector
1839         //
1840         //--------------
1841
1842       } // end loop over Pside 1D clusters
1843
1844     } // bad Nside module
1845
1846     //---------------------------------------------------------
1847
1848   } // use bad channels
1849
1850   //cout<<ncl<<" clusters for this module"<<endl;
1851
1852   delete [] cnegative;
1853   delete [] cused1;
1854   delete [] negativepair;
1855   delete [] cpositive;
1856   delete [] cused2;
1857   delete [] positivepair;
1858
1859 }