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