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 (fOverlapFlagClustersLay1[partners[iC2]] || fOverlapFlagClustersLay2[iC2]) continue;
 446     if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (partners[iC2],iC2);
 447
 448     fTracklets[fNTracklets] = new Float_t[6];
 449
 450     // use the theta from the clusters in the first layer
 451     fTracklets[fNTracklets][0] = fClustersLay1[partners[iC2]][0];
 452     // use the phi from the clusters in the first layer
 453     fTracklets[fNTracklets][1] = fClustersLay1[partners[iC2]][1];
 454     // store the difference between phi1 and phi2
 455     fTracklets[fNTracklets][2] = fClustersLay1[partners[iC2]][1] - fClustersLay2[iC2][1];
 456
 457     // define dphi in the range [0,pi] with proper sign (track charge correlated)
 458     if (fTracklets[fNTracklets][2] > TMath::Pi())
 459       fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]-2.*TMath::Pi();
 460     if (fTracklets[fNTracklets][2] < -TMath::Pi())
 461       fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]+2.*TMath::Pi();
 462
 463     // store the difference between theta1 and theta2
 464     fTracklets[fNTracklets][3] = fClustersLay1[partners[iC2]][0] - fClustersLay2[iC2][0];
 465
 466     if (fHistOn) {
 467       fhClustersDPhiAcc->Fill(fTracklets[fNTracklets][2]);
 468       fhClustersDThetaAcc->Fill(fTracklets[fNTracklets][3]);
 469       fhDPhiVsDThetaAcc->Fill(fTracklets[fNTracklets][3],fTracklets[fNTracklets][2]);
 470     }
 471
 472     // find label
 473     // if equal label in both clusters found this label is assigned
 474     // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
 475     Int_t label1 = 0;
 476     Int_t label2 = 0;
 477     while (label2 < 3) {
 478       if ((Int_t) fClustersLay1[partners[iC2]][3+label1] != -2 && (Int_t) fClustersLay1[partners[iC2]][3+label1] == (Int_t) fClustersLay2[iC2][3+label2])
 479         break;
 480       label1++;
 481       if (label1 == 3) {
 482         label1 = 0;
 483         label2++;
 484       }
 485     }
 486     if (label2 < 3) {
 487       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));
 488       fTracklets[fNTracklets][4] = fClustersLay1[partners[iC2]][3+label1];
 489       fTracklets[fNTracklets][5] = fClustersLay2[iC2][3+label2];
 490     } else {
 491       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));
 492       fTracklets[fNTracklets][4] = fClustersLay1[partners[iC2]][3];
 493       fTracklets[fNTracklets][5] = fClustersLay2[iC2][3];
 494     }
 495
 496     if (fHistOn) {
 497       Float_t eta=fTracklets[fNTracklets][0];
 498       eta= TMath::Tan(eta/2.);
 499       eta=-TMath::Log(eta);
 500       fhetaTracklets->Fill(eta);
 501       fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
 502     }
 503
 504     AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
 505     AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", partners[iC2], iC2));
 506     fNTracklets++;
 507
 508     associatedLay1[partners[iC2]] = 1;
 509   }
 510
 511   // Delete the following else if you do not want to save Clusters!
 512   // store the cluster
 513   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 514     if (associatedLay1[iC1]==2||associatedLay1[iC1]==0) {
 515       fSClusters[fNSingleCluster] = new Float_t[2];
 516       fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
 517       fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
 518       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
 519                 fNSingleCluster, iC1));
 520       fNSingleCluster++;
 521     }
 522   }
 523
 524   delete[] partners;
 525   delete[] minDists;
 526
 527   for (Int_t i=0; i<fNClustersLay1; i++)
 528     if (blacklist[i])
 529       delete blacklist[i];
 530   delete[] blacklist;
 531
 532   AliDebug(1,Form("%d tracklets found", fNTracklets));
 533 }
 534
 535 //____________________________________________________________________
 536 void
 537 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
 538   // This method
 539   // - gets the clusters from the cluster tree
 540   // - convert them into global coordinates
 541   // - store them in the internal arrays
 542   // - count the number of cluster-fired chips
 543
 544   AliDebug(1,"Loading clusters and cluster-fired chips ...");
 545
 546   fNClustersLay1 = 0;
 547   fNClustersLay2 = 0;
 548   fNFiredChips[0] = 0;
 549   fNFiredChips[1] = 0;
 550
 551   AliITSsegmentationSPD seg;
 552
 553   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 554   TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
 555   if(!rpcont->IsSPDActive()){
 556     AliWarning("No SPD rec points found, multiplicity not calculated");
 557     return;
 558   }
 559   Float_t cluGlo[3]={0.,0.,0.};
 560
 561
 562   // count clusters
 563   // loop over the SPD subdetectors
 564   Int_t nSPDL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
 565   for (Int_t iIts=0; iIts < nSPDL1; iIts++) {
 566     itsClusters=rpcont->UncheckedGetClusters(iIts);
 567     fNClustersLay1 += itsClusters->GetEntriesFast();
 568   }
 569   Int_t nSPDL2=AliITSgeomTGeo::GetModuleIndex(3,1,1);
 570   for (Int_t iIts=nSPDL1; iIts < nSPDL2; iIts++) {
 571     itsClusters=rpcont->UncheckedGetClusters(iIts);
 572     fNClustersLay2 += itsClusters->GetEntriesFast();
 573   }
 574
 575   // create arrays
 576   fClustersLay1              = new Float_t*[fNClustersLay1];
 577   fDetectorIndexClustersLay1 = new Int_t[fNClustersLay1];
 578   fOverlapFlagClustersLay1   = new Bool_t[fNClustersLay1];
 579
 580   fClustersLay2              = new Float_t*[fNClustersLay2];
 581   fDetectorIndexClustersLay2 = new Int_t[fNClustersLay2];
 582   fOverlapFlagClustersLay2   = new Bool_t[fNClustersLay2];
 583
 584   // no double association allowed
 585   fTracklets                 = new Float_t*[TMath::Min(fNClustersLay1, fNClustersLay2)];
 586   fSClusters                 = new Float_t*[fNClustersLay1];
 587
 588   for (Int_t i=0; i<fNClustersLay1; i++) {
 589     fClustersLay1[i]       = new Float_t[6];
 590     fOverlapFlagClustersLay1[i]   = kFALSE;
 591     fSClusters[i] = 0;
 592   }
 593
 594   for (Int_t i=0; i<fNClustersLay2; i++) {
 595     fClustersLay2[i]       = new Float_t[6];
 596     fOverlapFlagClustersLay2[i]   = kFALSE;
 597   }
 598
 599   for (Int_t i=0; i<TMath::Min(fNClustersLay1, fNClustersLay2); i++)
 600     fTracklets[i] = 0;
 601
 602   // fill clusters
 603   // loop over the its subdetectors
 604   fNClustersLay1 = 0; // reset to 0
 605   fNClustersLay2 = 0;
 606   for (Int_t iIts=0; iIts < nSPDL2; iIts++) {
 607
 608     itsClusters=rpcont->UncheckedGetClusters(iIts);
 609
 610     Int_t nClusters = itsClusters->GetEntriesFast();
 611
 612     // number of clusters in each chip of the current module
 613     Int_t nClustersInChip[5] = {0,0,0,0,0};
 614     Int_t layer = 0;
 615
 616     // loop over clusters
 617     while(nClusters--) {
 618       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 619
 620       layer = cluster->GetLayer();
 621       if (layer>1) continue;
 622
 623       cluster->GetGlobalXYZ(cluGlo);
 624       Float_t x = cluGlo[0];
 625       Float_t y = cluGlo[1];
 626       Float_t z = cluGlo[2];
 627
 628       // find the chip for the current cluster
 629       Float_t locz = cluster->GetDetLocalZ();
 630       Int_t iChip = seg.GetChipFromLocal(0,locz);
 631       nClustersInChip[iChip]++;
 632
 633       if (layer==0) {
 634         fClustersLay1[fNClustersLay1][0] = x;
 635         fClustersLay1[fNClustersLay1][1] = y;
 636         fClustersLay1[fNClustersLay1][2] = z;
 637
 638         fDetectorIndexClustersLay1[fNClustersLay1]=iIts;
 639
 640         for (Int_t i=0; i<3; i++)
 641                 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
 642         fNClustersLay1++;
 643       }
 644       if (layer==1) {
 645         fClustersLay2[fNClustersLay2][0] = x;
 646         fClustersLay2[fNClustersLay2][1] = y;
 647         fClustersLay2[fNClustersLay2][2] = z;
 648
 649         fDetectorIndexClustersLay2[fNClustersLay2]=iIts;
 650
 651         for (Int_t i=0; i<3; i++)
 652                 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
 653         fNClustersLay2++;
 654       }
 655
 656     }// end of cluster loop
 657
 658     // get number of fired chips in the current module
 659
 660     for(Int_t ifChip=0; ifChip<5; ifChip++) {
 661       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
 662     }
 663
 664   } // end of its "subdetector" loop
 665
 666   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
 667   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 668 }
 669 //____________________________________________________________________
 670 void
 671 AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
 672   // This method
 673   // - gets the clusters from the cluster tree
 674   // - counts the number of (cluster)fired chips
 675
 676   AliDebug(1,"Loading cluster-fired chips ...");
 677
 678   fNFiredChips[0] = 0;
 679   fNFiredChips[1] = 0;
 680
 681   AliITSsegmentationSPD seg;
 682   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 683   TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
 684   if(!rpcont->IsSPDActive()){
 685     AliWarning("No SPD rec points found, multiplicity not calculated");
 686     return;
 687   }
 688
 689   // loop over the its subdetectors
 690   Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1);
 691   for (Int_t iIts=0; iIts < nSPDmodules; iIts++) {
 692     itsClusters=rpcont->UncheckedGetClusters(iIts);
 693     Int_t nClusters = itsClusters->GetEntriesFast();
 694
 695     // number of clusters in each chip of the current module
 696     Int_t nClustersInChip[5] = {0,0,0,0,0};
 697     Int_t layer = 0;
 698
 699     // loop over clusters
 700     while(nClusters--) {
 701       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 702
 703       layer = cluster->GetLayer();
 704       if (layer>1) continue;
 705
 706       // find the chip for the current cluster
 707       Float_t locz = cluster->GetDetLocalZ();
 708       Int_t iChip = seg.GetChipFromLocal(0,locz);
 709       nClustersInChip[iChip]++;
 710
 711     }// end of cluster loop
 712
 713     // get number of fired chips in the current module
 714     for(Int_t ifChip=0; ifChip<5; ifChip++) {
 715       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
 716     }
 717
 718   } // end of its "subdetector" loop
 719
 720
 721   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 722 }
 723 //____________________________________________________________________
 724 void
 725 AliITSMultReconstructor::SaveHists() {
 726   // This method save the histograms on the output file
 727   // (only if fHistOn is TRUE).
 728
 729   if (!fHistOn)
 730     return;
 731
 732   fhClustersDPhiAll->Write();
 733   fhClustersDThetaAll->Write();
 734   fhDPhiVsDThetaAll->Write();
 735
 736   fhClustersDPhiAcc->Write();
 737   fhClustersDThetaAcc->Write();
 738   fhDPhiVsDThetaAcc->Write();
 739
 740   fhetaTracklets->Write();
 741   fhphiTracklets->Write();
 742   fhetaClustersLay1->Write();
 743   fhphiClustersLay1->Write();
 744 }
 745
 746 //____________________________________________________________________
 747 void
 748 AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
 749
 750   Float_t distClSameMod=0.;
 751   Float_t distClSameModMin=0.;
 752   Int_t   iClOverlap =0;
 753   Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
 754   Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
 755
 756   Float_t zproj1=0.;
 757   Float_t zproj2=0.;
 758   Float_t deZproj=0.;
 759
 760   // Loop on inner layer clusters
 761   for (Int_t iiC1=0; iiC1<fNClustersLay1; iiC1++) {
 762     if (!fOverlapFlagClustersLay1[iiC1]) {
 763       // only for adjacent modules
 764       if ((TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==4)||
 765          (TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==76)) {
 766         Float_t dePhi=TMath::Abs(fClustersLay1[iiC1][1]-fClustersLay1[iC1][1]);
 767         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 768
 769         zproj1=meanRadiusLay1/TMath::Tan(fClustersLay1[iC1][0]);
 770         zproj2=meanRadiusLay1/TMath::Tan(fClustersLay1[iiC1][0]);
 771
 772         deZproj=TMath::Abs(zproj1-zproj2);
 773
 774         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 775         if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
 776
 777 //        if (distClSameMod<=1.) {
 778 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 779 //            distClSameModMin=distClSameMod;
 780 //            iClOverlap=iiC1;
 781 //          }
 782 //        }
 783
 784
 785       } // end adjacent modules
 786     }
 787   } // end Loop on inner layer clusters
 788
 789 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
 790
 791   distClSameMod=0.;
 792   distClSameModMin=0.;
 793   iClOverlap =0;
 794   // Loop on outer layer clusters
 795   for (Int_t iiC2=0; iiC2<fNClustersLay2; iiC2++) {
 796     if (!fOverlapFlagClustersLay2[iiC2]) {
 797       // only for adjacent modules
 798       if ((TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==4) ||
 799          (TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==156)) {
 800         Float_t dePhi=TMath::Abs(fClustersLay2[iiC2][1]-fClustersLay2[iC2WithBestDist][1]);
 801         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 802
 803         zproj1=meanRadiusLay2/TMath::Tan(fClustersLay2[iC2WithBestDist][0]);
 804         zproj2=meanRadiusLay2/TMath::Tan(fClustersLay2[iiC2][0]);
 805
 806         deZproj=TMath::Abs(zproj1-zproj2);
 807         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 808         if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE;
 809
 810 //        if (distClSameMod<=1.) {
 811 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 812 //            distClSameModMin=distClSameMod;
 813 //            iClOverlap=iiC2;
 814 //          }
 815 //        }
 816
 817       } // end adjacent modules
 818     }
 819   } // end Loop on outer layer clusters
 820
 821 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;
 822
 823 }