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