ITS/AliITSMultReconstructor.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 //
  20 //        Implementation of the ITS-SPD trackleter class
  21 //
  22 // Retrieves clusters in the pixels (theta and phi) and tracklets.
  23 // These can be used to extract charged particle multiplcicity from the ITS.
  24 //
  25 // A tracklet consist of two ITS clusters, one in the first pixel layer and
  26 // one in the second. The clusters are associates if the differencies in
  27 // Phi (azimuth) and Theta (polar) are within fiducial values.
  28 // In case of multiple candidates it is selected the candidate with minimum
  29 // distance.
  30 //
  31 // Two methods return the number of traklets and the number of unassociated
  32 // clusters (i.e. not used in any tracklet) in the first SPD layer
  33 // (GetNTracklets and GetNSingleClusters)
  34 //
  35 // The cuts on phi and theta depend on the interacting system (p-p or Pb-Pb)
  36 // and can be set via AliITSRecoParam class
  37 // (SetPhiWindow and SetThetaWindow)
  38 //
  39 // Origin: Tiziano Virgili
  40 //
  41 // Current support and development:
  42 //         Domenico Elia, Maria Nicassio (INFN Bari)
  43 //         Domenico.Elia@ba.infn.it, Maria.Nicassio@ba.infn.it
  44 //
  45 // Most recent updates:
  46 //     - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE)
  47 //     - phi definition changed to ALICE convention (0,2*TMath::pi())
  48 //     - cluster coordinates taken with GetGlobalXYZ()
  49 //     - fGeometry removed
  50 //     - number of fired chips on the two layers
  51 //     - option to avoid duplicates in the overlaps (fRemoveClustersFromOverlaps)
  52 //     - options and fiducial cuts via AliITSRecoParam
  53 //     - move from DeltaZeta to DeltaTheta cut (M. Nicassio)
  54 //     - update to the new algo (M. Nicassio, J.F.Grosse-Oetringhaus)
  55 //_________________________________________________________________________
  56
  57 #include <TClonesArray.h>
  58 #include <TH1F.h>
  59 #include <TH2F.h>
  60 #include <TTree.h>
  61 #include "TArrayI.h"
  62
  63 #include "AliITSMultReconstructor.h"
  64 #include "AliITSReconstructor.h"
  65 #include "AliITSsegmentationSPD.h"
  66 #include "AliITSRecPoint.h"
  67 #include "AliITSgeom.h"
  68 #include "AliLog.h"
  69 //#include "TGeoGlobalMagField.h"
  70 //#include "AliMagF.h"
  71
  72 //____________________________________________________________________
  73 ClassImp(AliITSMultReconstructor)
  74
  75
  76 //____________________________________________________________________
  77 AliITSMultReconstructor::AliITSMultReconstructor():
  78 TObject(),
  79 fClustersLay1(0),
  80 fClustersLay2(0),
  81 fDetectorIndexClustersLay1(0),
  82 fDetectorIndexClustersLay2(0),
  83 fOverlapFlagClustersLay1(0),
  84 fOverlapFlagClustersLay2(0),
  85 fTracklets(0),
  86 fSClusters(0),
  87 fNClustersLay1(0),
  88 fNClustersLay2(0),
  89 fNTracklets(0),
  90 fNSingleCluster(0),
  91 fPhiWindow(0),
  92 fThetaWindow(0),
  93 fRemoveClustersFromOverlaps(0),
  94 fPhiOverlapCut(0),
  95 fZetaOverlapCut(0),
  96 fHistOn(0),
  97 fhClustersDPhiAcc(0),
  98 fhClustersDThetaAcc(0),
  99 fhClustersDPhiAll(0),
 100 fhClustersDThetaAll(0),
 101 fhDPhiVsDThetaAll(0),
 102 fhDPhiVsDThetaAcc(0),
 103 fhetaTracklets(0),
 104 fhphiTracklets(0),
 105 fhetaClustersLay1(0),
 106 fhphiClustersLay1(0){
 107
 108   fNFiredChips[0] = 0;
 109   fNFiredChips[1] = 0;
 110
 111   // Method to reconstruct the charged particles multiplicity with the
 112   // SPD (tracklets).
 113
 114
 115   SetHistOn();
 116
 117   if(AliITSReconstructor::GetRecoParam()) {
 118     SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow());
 119     SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow());
 120     SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps());
 121     SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut());
 122     SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut());
 123   } else {
 124     SetPhiWindow();
 125     SetThetaWindow();
 126     SetRemoveClustersFromOverlaps();
 127     SetPhiOverlapCut();
 128     SetZetaOverlapCut();
 129   }
 130
 131
 132   fClustersLay1              = new Float_t*[300000];
 133   fClustersLay2              = new Float_t*[300000];
 134   fDetectorIndexClustersLay1 = new Int_t[300000];
 135   fDetectorIndexClustersLay2 = new Int_t[300000];
 136   fOverlapFlagClustersLay1   = new Bool_t[300000];
 137   fOverlapFlagClustersLay2   = new Bool_t[300000];
 138   fTracklets                 = new Float_t*[300000];
 139   fSClusters                 = new Float_t*[300000];
 140
 141   for(Int_t i=0; i<300000; i++) {
 142     fClustersLay1[i]       = new Float_t[6];
 143     fClustersLay2[i]       = new Float_t[6];
 144     fTracklets[i]          = new Float_t[6];
 145     fSClusters[i]           = new Float_t[2];
 146     fOverlapFlagClustersLay1[i]   = kFALSE;
 147     fOverlapFlagClustersLay2[i]   = kFALSE;
 148   }
 149
 150   // definition of histograms
 151   fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,-0.1,0.1);
 152   fhClustersDPhiAcc->SetDirectory(0);
 153   fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
 154   fhClustersDThetaAcc->SetDirectory(0);
 155
 156   fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
 157   fhDPhiVsDThetaAcc->SetDirectory(0);
 158
 159   fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,0.0,0.5);
 160   fhClustersDPhiAll->SetDirectory(0);
 161   fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,0.0,0.5);
 162   fhClustersDThetaAll->SetDirectory(0);
 163
 164   fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,0.,0.5,100,0.,0.5);
 165   fhDPhiVsDThetaAll->SetDirectory(0);
 166
 167   fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
 168   fhetaTracklets->SetDirectory(0);
 169   fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100, 0., 2*TMath::Pi());
 170   fhphiTracklets->SetDirectory(0);
 171   fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
 172   fhetaClustersLay1->SetDirectory(0);
 173   fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
 174   fhphiClustersLay1->SetDirectory(0);
 175 }
 176
 177 //______________________________________________________________________
 178 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
 179 fClustersLay1(mr.fClustersLay1),
 180 fClustersLay2(mr.fClustersLay2),
 181 fDetectorIndexClustersLay1(mr.fDetectorIndexClustersLay1),
 182 fDetectorIndexClustersLay2(mr.fDetectorIndexClustersLay2),
 183 fOverlapFlagClustersLay1(mr.fOverlapFlagClustersLay1),
 184 fOverlapFlagClustersLay2(mr.fOverlapFlagClustersLay2),
 185 fTracklets(mr.fTracklets),
 186 fSClusters(mr.fSClusters),
 187 fNClustersLay1(mr.fNClustersLay1),
 188 fNClustersLay2(mr.fNClustersLay2),
 189 fNTracklets(mr.fNTracklets),
 190 fNSingleCluster(mr.fNSingleCluster),
 191 fPhiWindow(mr.fPhiWindow),
 192 fThetaWindow(mr.fThetaWindow),
 193 fRemoveClustersFromOverlaps(mr.fRemoveClustersFromOverlaps),
 194 fPhiOverlapCut(mr.fPhiOverlapCut),
 195 fZetaOverlapCut(mr.fZetaOverlapCut),
 196 fHistOn(mr.fHistOn),
 197 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
 198 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
 199 fhClustersDPhiAll(mr.fhClustersDPhiAll),
 200 fhClustersDThetaAll(mr.fhClustersDThetaAll),
 201 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
 202 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
 203 fhetaTracklets(mr.fhetaTracklets),
 204 fhphiTracklets(mr.fhphiTracklets),
 205 fhetaClustersLay1(mr.fhetaClustersLay1),
 206 fhphiClustersLay1(mr.fhphiClustersLay1) {
 207   // Copy constructor
 208
 209 }
 210
 211 //______________________________________________________________________
 212 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
 213   // Assignment operator
 214   this->~AliITSMultReconstructor();
 215   new(this) AliITSMultReconstructor(mr);
 216   return *this;
 217 }
 218
 219 //______________________________________________________________________
 220 AliITSMultReconstructor::~AliITSMultReconstructor(){
 221   // Destructor
 222
 223   // delete histograms
 224   delete fhClustersDPhiAcc;
 225   delete fhClustersDThetaAcc;
 226   delete fhClustersDPhiAll;
 227   delete fhClustersDThetaAll;
 228   delete fhDPhiVsDThetaAll;
 229   delete fhDPhiVsDThetaAcc;
 230   delete fhetaTracklets;
 231   delete fhphiTracklets;
 232   delete fhetaClustersLay1;
 233   delete fhphiClustersLay1;
 234
 235   // delete arrays
 236   for(Int_t i=0; i<300000; i++) {
 237     delete [] fClustersLay1[i];
 238     delete [] fClustersLay2[i];
 239     delete [] fTracklets[i];
 240     delete [] fSClusters[i];
 241   }
 242   delete [] fClustersLay1;
 243   delete [] fClustersLay2;
 244   delete [] fDetectorIndexClustersLay1;
 245   delete [] fDetectorIndexClustersLay2;
 246   delete [] fOverlapFlagClustersLay1;
 247   delete [] fOverlapFlagClustersLay2;
 248   delete [] fTracklets;
 249   delete [] fSClusters;
 250 }
 251
 252 //____________________________________________________________________
 253 void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
 254   //
 255   // - calls LoadClusterArray that finds the position of the clusters
 256   //   (in global coord)
 257   // - convert the cluster coordinates to theta, phi (seen from the
 258   //   interaction vertex).
 259   // - makes an array of tracklets
 260   //
 261   // After this method has been called, the clusters of the two layers
 262   // and the tracklets can be retrieved by calling the Get'er methods.
 263
 264   // reset counters
 265   fNClustersLay1 = 0;
 266   fNClustersLay2 = 0;
 267   fNTracklets = 0;
 268   fNSingleCluster = 0;
 269
 270   // loading the clusters
 271   LoadClusterArrays(clusterTree);
 272
 273   const Double_t pi = TMath::Pi();
 274   Int_t* partners = new Int_t[fNClustersLay2];
 275   Float_t* minDists = new Float_t[fNClustersLay2];
 276   Int_t* associatedLay1 = new Int_t[fNClustersLay1];
 277   TArrayI** blacklist = new TArrayI*[fNClustersLay1];
 278
 279   for (Int_t i=0; i<fNClustersLay2; i++) {
 280     partners[i] = -1;
 281     minDists[i] = 2;
 282   }
 283   for (Int_t i=0; i<fNClustersLay1; i++)
 284     associatedLay1[i] = 0;
 285   for (Int_t i=0; i<fNClustersLay1; i++)
 286     blacklist[i] = 0;
 287
 288   // find the tracklets
 289   AliDebug(1,"Looking for tracklets... ");
 290
 291   //###########################################################
 292   // Loop on layer 1 : finding theta, phi and z
 293   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 294     Float_t x = fClustersLay1[iC1][0] - vtx[0];
 295     Float_t y = fClustersLay1[iC1][1] - vtx[1];
 296     Float_t z = fClustersLay1[iC1][2] - vtx[2];
 297
 298     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
 299                                TMath::Power(y,2) +
 300                                TMath::Power(z,2));
 301
 302     fClustersLay1[iC1][0] = TMath::ACos(z/r);                   // Store Theta
 303     fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 304
 305     if (fHistOn) {
 306       Float_t eta=fClustersLay1[iC1][0];
 307       eta= TMath::Tan(eta/2.);
 308       eta=-TMath::Log(eta);
 309       fhetaClustersLay1->Fill(eta);
 310       fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
 311     }
 312   }
 313
 314   // Loop on layer 2 : finding theta, phi and r
 315   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 316     Float_t x = fClustersLay2[iC2][0] - vtx[0];
 317     Float_t y = fClustersLay2[iC2][1] - vtx[1];
 318     Float_t z = fClustersLay2[iC2][2] - vtx[2];
 319
 320     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
 321                                TMath::Power(y,2) +
 322                                TMath::Power(z,2));
 323
 324     fClustersLay2[iC2][0] = TMath::ACos(z/r);                   // Store Theta
 325     fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 326   }
 327
 328   //###########################################################
 329   Int_t found = 1;
 330   while (found > 0) {
 331     //Printf("Iteration...");
 332     found = 0;
 333
 334     // Step1: find all tracklets allowing double assocation
 335     // Loop on layer 1
 336     for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 337
 338       // already used or in the overlap ?
 339       if (associatedLay1[iC1] != 0 || fOverlapFlagClustersLay1[iC1]) continue;
 340
 341       found++;
 342
 343       // reset of variables for multiple candidates
 344       Int_t  iC2WithBestDist = -1;   // reset
 345       Double_t minDist       =  2;   // reset
 346
 347       // Loop on layer 2
 348       for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 349
 350         // in the overlap ?
 351         if (fOverlapFlagClustersLay2[iC2]) continue;
 352
 353         if (blacklist[iC1]) {
 354           Bool_t blacklisted = kFALSE;
 355           for (Int_t i=0; i<blacklist[iC1]->GetSize(); i++) {
 356             if (blacklist[iC1]->At(i) == iC2) {
 357               blacklisted = kTRUE;
 358               break;
 359             }
 360           }
 361           if (blacklisted) continue;
 362         }
 363
 364         // find the difference in angles
 365         Double_t dTheta = TMath::Abs(fClustersLay2[iC2][0] - fClustersLay1[iC1][0]);
 366         Double_t dPhi  = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
 367         // take into account boundary condition
 368         if (dPhi>pi) dPhi=2.*pi-dPhi;
 369
 370         // Account for B-field shifting dphi ? ...
 371         // Double_t pos[3]={0.,0.,0.};
 372         // Double_t B[3]={0.,0.,0.};
 373         // TGeoGlobalMagField::Instance()->Field(pos,B);
 374         // if (B[2]=!0) dPhi-=0.005; // field dependent
 375
 376         if (fHistOn) {
 377           fhClustersDPhiAll->Fill(dPhi);
 378           fhClustersDThetaAll->Fill(dTheta);
 379           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
 380         }
 381         // make "elliptical" cut in Phi and Theta!
 382         Float_t d = TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dTheta/fThetaWindow,2);
 383
 384         // look for the minimum distance: the minimum is in iC2WithBestDist
 385         if (d<1 && d<minDist ) {
 386           minDist=d;
 387           iC2WithBestDist = iC2;
 388         }
 389       } // end of loop over clusters in layer 2
 390
 391       if (minDist<1) { // This means that a cluster in layer 2 was found that matches with iC1
 392
 393         if (minDists[iC2WithBestDist] > minDist) {
 394           Int_t oldPartner = partners[iC2WithBestDist];
 395           partners[iC2WithBestDist] = iC1;
 396           minDists[iC2WithBestDist] = minDist;
 397
 398           // mark as assigned
 399           associatedLay1[iC1] = 1;
 400
 401           if (oldPartner != -1) {
 402             // redo partner search for cluster in L0 (oldPartner), putting this one (best1) on its blacklist
 403             if (blacklist[oldPartner] == 0) {
 404               blacklist[oldPartner] = new TArrayI(1);
 405             } else blacklist[oldPartner]->Set(blacklist[oldPartner]->GetSize()+1);
 406
 407             blacklist[oldPartner]->AddAt(iC2WithBestDist, blacklist[oldPartner]->GetSize()-1);
 408
 409             // mark as free
 410             associatedLay1[iC1] = 0;
 411           }
 412         } else {
 413           // try again to find a cluster without considering iC2WithBestDist
 414           if (blacklist[iC1] == 0) {
 415             blacklist[iC1] = new TArrayI(1);
 416           }
 417           else blacklist[iC1]->Set(blacklist[iC1]->GetSize()+1);
 418
 419           blacklist[iC1]->AddAt(iC2WithBestDist, blacklist[iC1]->GetSize()-1);
 420         }
 421
 422       } else // cluster has no partner; remove
 423         associatedLay1[iC1] = 2;
 424     } // end of loop over clusters in layer 1
 425   }
 426
 427   // Step2: store tracklets; remove used clusters
 428   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 429
 430     if (partners[iC2] == -1) continue;
 431
 432     //Printf("Tracklet7: %d %d %f %d %d %f %f %f", partners[iC2], iC2, minDists[iC2], (Int_t) fClustersLay1[partners[iC2]][3], (Int_t) fClustersLay2[iC2][3], fClustersLay1[partners[iC2]][0], fClustersLay2[iC2][0], fClustersLay1[partners[iC2]][1]);
 433
 434
 435     if (fOverlapFlagClustersLay1[partners[iC2]] || fOverlapFlagClustersLay2[iC2]) continue;
 436     if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (partners[iC2],iC2);
 437
 438     // use the theta from the clusters in the first layer
 439     fTracklets[fNTracklets][0] = fClustersLay1[partners[iC2]][0];
 440     // use the phi from the clusters in the first layer
 441     fTracklets[fNTracklets][1] = fClustersLay1[partners[iC2]][1];
 442     // store the difference between phi1 and phi2
 443     fTracklets[fNTracklets][2] = fClustersLay1[partners[iC2]][1] - fClustersLay2[iC2][1];
 444
 445     // define dphi in the range [0,pi] with proper sign (track charge correlated)
 446     if (fTracklets[fNTracklets][2] > TMath::Pi())
 447       fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]-2.*TMath::Pi();
 448     if (fTracklets[fNTracklets][2] < -TMath::Pi())
 449       fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]+2.*TMath::Pi();
 450
 451     // store the difference between theta1 and theta2
 452     fTracklets[fNTracklets][3] = fClustersLay1[partners[iC2]][0] - fClustersLay2[iC2][0];
 453
 454     if (fHistOn) {
 455       fhClustersDPhiAcc->Fill(fTracklets[fNTracklets][2]);
 456       fhClustersDThetaAcc->Fill(fTracklets[fNTracklets][3]);
 457       fhDPhiVsDThetaAcc->Fill(fTracklets[fNTracklets][3],fTracklets[fNTracklets][2]);
 458     }
 459
 460     // find label
 461     // if equal label in both clusters found this label is assigned
 462     // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
 463     Int_t label1 = 0;
 464     Int_t label2 = 0;
 465     while (label2 < 3) {
 466       if ((Int_t) fClustersLay1[partners[iC2]][3+label1] != -2 && (Int_t) fClustersLay1[partners[iC2]][3+label1] == (Int_t) fClustersLay2[iC2][3+label2])
 467         break;
 468       label1++;
 469       if (label1 == 3) {
 470         label1 = 0;
 471         label2++;
 472       }
 473     }
 474     if (label2 < 3) {
 475       AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) fClustersLay1[partners[iC2]][3+label1], (Int_t) fClustersLay2[iC2][3+label2], fNTracklets));
 476       fTracklets[fNTracklets][4] = fClustersLay1[partners[iC2]][3+label1];
 477       fTracklets[fNTracklets][5] = fClustersLay2[iC2][3+label2];
 478     } else {
 479       AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", (Int_t) fClustersLay1[partners[iC2]][3], (Int_t) fClustersLay1[partners[iC2]][4], (Int_t) fClustersLay1[partners[iC2]][5], (Int_t) fClustersLay2[iC2][3], (Int_t) fClustersLay2[iC2][4], (Int_t) fClustersLay2[iC2][5], fNTracklets));
 480       fTracklets[fNTracklets][4] = fClustersLay1[partners[iC2]][3];
 481       fTracklets[fNTracklets][5] = fClustersLay2[iC2][3];
 482     }
 483
 484     if (fHistOn) {
 485       Float_t eta=fTracklets[fNTracklets][0];
 486       eta= TMath::Tan(eta/2.);
 487       eta=-TMath::Log(eta);
 488       fhetaTracklets->Fill(eta);
 489       fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
 490     }
 491
 492     AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
 493     AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", partners[iC2], iC2));
 494     fNTracklets++;
 495
 496     associatedLay1[partners[iC2]] = 1;
 497   }
 498
 499   // Delete the following else if you do not want to save Clusters!
 500   // store the cluster
 501   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 502     if (associatedLay1[iC1]==2||associatedLay1[iC1]==0) {
 503       fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
 504       fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
 505       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
 506                 fNSingleCluster, iC1));
 507       fNSingleCluster++;
 508     }
 509   }
 510
 511   delete[] partners;
 512   delete[] minDists;
 513
 514   for (Int_t i=0; i<fNClustersLay1; i++)
 515     if (blacklist[i])
 516       delete blacklist[i];
 517   delete[] blacklist;
 518
 519
 520   AliDebug(1,Form("%d tracklets found", fNTracklets));
 521 }
 522
 523 //____________________________________________________________________
 524 void
 525 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
 526   // This method
 527   // - gets the clusters from the cluster tree
 528   // - convert them into global coordinates
 529   // - store them in the internal arrays
 530   // - count the number of cluster-fired chips
 531
 532   AliDebug(1,"Loading clusters and cluster-fired chips ...");
 533
 534   fNClustersLay1 = 0;
 535   fNClustersLay2 = 0;
 536   fNFiredChips[0] = 0;
 537   fNFiredChips[1] = 0;
 538
 539   AliITSsegmentationSPD *seg = new AliITSsegmentationSPD();
 540
 541   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
 542   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
 543
 544   itsClusterBranch->SetAddress(&itsClusters);
 545
 546   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
 547   Float_t cluGlo[3]={0.,0.,0.};
 548
 549   // loop over the its subdetectors
 550   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
 551
 552     if (!itsClusterTree->GetEvent(iIts))
 553       continue;
 554
 555     Int_t nClusters = itsClusters->GetEntriesFast();
 556
 557     // number of clusters in each chip of the current module
 558     Int_t nClustersInChip[5] = {0,0,0,0,0};
 559     Int_t layer = 0;
 560
 561     // loop over clusters
 562     while(nClusters--) {
 563       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 564
 565       layer = cluster->GetLayer();
 566       if (layer>1) continue;
 567
 568       cluster->GetGlobalXYZ(cluGlo);
 569       Float_t x = cluGlo[0];
 570       Float_t y = cluGlo[1];
 571       Float_t z = cluGlo[2];
 572
 573       // find the chip for the current cluster
 574       Float_t locz = cluster->GetDetLocalZ();
 575       Int_t iChip = seg->GetChipFromLocal(0,locz);
 576       nClustersInChip[iChip]++;
 577
 578       if (layer==0) {
 579         fClustersLay1[fNClustersLay1][0] = x;
 580         fClustersLay1[fNClustersLay1][1] = y;
 581         fClustersLay1[fNClustersLay1][2] = z;
 582
 583         fDetectorIndexClustersLay1[fNClustersLay1]=iIts;
 584
 585         for (Int_t i=0; i<3; i++)
 586                 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
 587         fNClustersLay1++;
 588       }
 589       if (layer==1) {
 590         fClustersLay2[fNClustersLay2][0] = x;
 591         fClustersLay2[fNClustersLay2][1] = y;
 592         fClustersLay2[fNClustersLay2][2] = z;
 593
 594         fDetectorIndexClustersLay2[fNClustersLay2]=iIts;
 595
 596         for (Int_t i=0; i<3; i++)
 597                 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
 598         fNClustersLay2++;
 599       }
 600
 601     }// end of cluster loop
 602
 603     // get number of fired chips in the current module
 604     if(layer<2)
 605     for(Int_t ifChip=0; ifChip<5; ifChip++) {
 606       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
 607     }
 608
 609   } // end of its "subdetector" loop
 610
 611   if (itsClusters) {
 612     itsClusters->Delete();
 613     delete itsClusters;
 614     delete seg;
 615     itsClusters = 0;
 616   }
 617   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
 618   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 619 }
 620 //____________________________________________________________________
 621 void
 622 AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
 623   // This method
 624   // - gets the clusters from the cluster tree
 625   // - counts the number of (cluster)fired chips
 626
 627   AliDebug(1,"Loading cluster-fired chips ...");
 628
 629   fNFiredChips[0] = 0;
 630   fNFiredChips[1] = 0;
 631
 632   AliITSsegmentationSPD *seg = new AliITSsegmentationSPD();
 633
 634   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
 635   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
 636
 637   itsClusterBranch->SetAddress(&itsClusters);
 638
 639   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
 640
 641   // loop over the its subdetectors
 642   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
 643
 644     if (!itsClusterTree->GetEvent(iIts))
 645       continue;
 646
 647     Int_t nClusters = itsClusters->GetEntriesFast();
 648
 649     // number of clusters in each chip of the current module
 650     Int_t nClustersInChip[5] = {0,0,0,0,0};
 651     Int_t layer = 0;
 652
 653     // loop over clusters
 654     while(nClusters--) {
 655       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 656
 657       layer = cluster->GetLayer();
 658       if (layer>1) continue;
 659
 660       // find the chip for the current cluster
 661       Float_t locz = cluster->GetDetLocalZ();
 662       Int_t iChip = seg->GetChipFromLocal(0,locz);
 663       nClustersInChip[iChip]++;
 664
 665     }// end of cluster loop
 666
 667     // get number of fired chips in the current module
 668     if(layer<2)
 669     for(Int_t ifChip=0; ifChip<5; ifChip++) {
 670       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
 671     }
 672
 673   } // end of its "subdetector" loop
 674
 675   if (itsClusters) {
 676     itsClusters->Delete();
 677     delete itsClusters;
 678     delete seg;
 679     itsClusters = 0;
 680   }
 681   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 682 }
 683 //____________________________________________________________________
 684 void
 685 AliITSMultReconstructor::SaveHists() {
 686   // This method save the histograms on the output file
 687   // (only if fHistOn is TRUE).
 688
 689   if (!fHistOn)
 690     return;
 691
 692   fhClustersDPhiAll->Write();
 693   fhClustersDThetaAll->Write();
 694   fhDPhiVsDThetaAll->Write();
 695
 696   fhClustersDPhiAcc->Write();
 697   fhClustersDThetaAcc->Write();
 698   fhDPhiVsDThetaAcc->Write();
 699
 700   fhetaTracklets->Write();
 701   fhphiTracklets->Write();
 702   fhetaClustersLay1->Write();
 703   fhphiClustersLay1->Write();
 704 }
 705
 706 //____________________________________________________________________
 707 void
 708 AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
 709
 710   Float_t distClSameMod=0.;
 711   Float_t distClSameModMin=0.;
 712   Int_t   iClOverlap =0;
 713   Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
 714   Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
 715
 716   Float_t zproj1=0.;
 717   Float_t zproj2=0.;
 718   Float_t deZproj=0.;
 719
 720   // Loop on inner layer clusters
 721   for (Int_t iiC1=0; iiC1<fNClustersLay1; iiC1++) {
 722     if (!fOverlapFlagClustersLay1[iiC1]) {
 723       // only for adjacent modules
 724       if ((TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==4)||
 725          (TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==76)) {
 726         Float_t dePhi=TMath::Abs(fClustersLay1[iiC1][1]-fClustersLay1[iC1][1]);
 727         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 728
 729         zproj1=meanRadiusLay1/TMath::Tan(fClustersLay1[iC1][0]);
 730         zproj2=meanRadiusLay1/TMath::Tan(fClustersLay1[iiC1][0]);
 731
 732         deZproj=TMath::Abs(zproj1-zproj2);
 733
 734         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 735         if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
 736
 737 //        if (distClSameMod<=1.) {
 738 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 739 //            distClSameModMin=distClSameMod;
 740 //            iClOverlap=iiC1;
 741 //          }
 742 //        }
 743
 744
 745       } // end adjacent modules
 746     }
 747   } // end Loop on inner layer clusters
 748
 749 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
 750
 751   distClSameMod=0.;
 752   distClSameModMin=0.;
 753   iClOverlap =0;
 754   // Loop on outer layer clusters
 755   for (Int_t iiC2=0; iiC2<fNClustersLay2; iiC2++) {
 756     if (!fOverlapFlagClustersLay2[iiC2]) {
 757       // only for adjacent modules
 758       if ((TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==4) ||
 759          (TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==156)) {
 760         Float_t dePhi=TMath::Abs(fClustersLay2[iiC2][1]-fClustersLay2[iC2WithBestDist][1]);
 761         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 762
 763         zproj1=meanRadiusLay2/TMath::Tan(fClustersLay2[iC2WithBestDist][0]);
 764         zproj2=meanRadiusLay2/TMath::Tan(fClustersLay2[iiC2][0]);
 765
 766         deZproj=TMath::Abs(zproj1-zproj2);
 767         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 768         if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE;
 769
 770 //        if (distClSameMod<=1.) {
 771 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 772 //            distClSameModMin=distClSameMod;
 773 //            iClOverlap=iiC2;
 774 //          }
 775 //        }
 776
 777       } // end adjacent modules
 778     }
 779   } // end Loop on outer layer clusters
 780
 781 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;
 782
 783 }