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