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