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 }