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 //
  50 //  fGeometry removed
  51 //____________________________________________________________________
  52
  53 #include <TClonesArray.h>
  54 #include <TH1F.h>
  55 #include <TH2F.h>
  56 #include <TTree.h>
  57
  58 #include "AliITSMultReconstructor.h"
  59 #include "AliITSRecPoint.h"
  60 #include "AliITSgeom.h"
  61 #include "AliLog.h"
  62
  63 //____________________________________________________________________
  64 ClassImp(AliITSMultReconstructor)
  65
  66
  67 //____________________________________________________________________
  68 AliITSMultReconstructor::AliITSMultReconstructor():
  69 fClustersLay1(0),
  70 fClustersLay2(0),
  71 fTracklets(0),
  72 fSClusters(0),
  73 fAssociationFlag(0),
  74 fNClustersLay1(0),
  75 fNClustersLay2(0),
  76 fNTracklets(0),
  77 fNSingleCluster(0),
  78 fPhiWindow(0),
  79 fZetaWindow(0),
  80 fOnlyOneTrackletPerC2(0),
  81 fHistOn(0),
  82 fhClustersDPhiAcc(0),
  83 fhClustersDThetaAcc(0),
  84 fhClustersDZetaAcc(0),
  85 fhClustersDPhiAll(0),
  86 fhClustersDThetaAll(0),
  87 fhClustersDZetaAll(0),
  88 fhDPhiVsDThetaAll(0),
  89 fhDPhiVsDThetaAcc(0),
  90 fhDPhiVsDZetaAll(0),
  91 fhDPhiVsDZetaAcc(0),
  92 fhetaTracklets(0),
  93 fhphiTracklets(0),
  94 fhetaClustersLay1(0),
  95 fhphiClustersLay1(0){
  96   // Method to reconstruct the charged particles multiplicity with the
  97   // SPD (tracklets).
  98
  99
 100   SetHistOn();
 101   SetPhiWindow();
 102   SetZetaWindow();
 103   SetOnlyOneTrackletPerC2();
 104
 105   fClustersLay1       = new Float_t*[300000];
 106   fClustersLay2       = new Float_t*[300000];
 107   fTracklets          = new Float_t*[300000];
 108   fSClusters          = new Float_t*[300000];
 109   fAssociationFlag    = new Bool_t[300000];
 110
 111   for(Int_t i=0; i<300000; i++) {
 112     fClustersLay1[i]       = new Float_t[6];
 113     fClustersLay2[i]       = new Float_t[6];
 114     fTracklets[i]          = new Float_t[4];
 115     fSClusters[i]           = new Float_t[2];
 116     fAssociationFlag[i]    = kFALSE;
 117   }
 118
 119   // definition of histograms
 120   fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,0.,0.1);
 121   fhClustersDPhiAcc->SetDirectory(0);
 122   fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
 123   fhClustersDThetaAcc->SetDirectory(0);
 124   fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
 125   fhClustersDZetaAcc->SetDirectory(0);
 126
 127   fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
 128   fhDPhiVsDZetaAcc->SetDirectory(0);
 129   fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
 130   fhDPhiVsDThetaAcc->SetDirectory(0);
 131
 132   fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,0.0,0.5);
 133   fhClustersDPhiAll->SetDirectory(0);
 134   fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
 135   fhClustersDThetaAll->SetDirectory(0);
 136   fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
 137   fhClustersDZetaAll->SetDirectory(0);
 138
 139   fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
 140   fhDPhiVsDZetaAll->SetDirectory(0);
 141   fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
 142   fhDPhiVsDThetaAll->SetDirectory(0);
 143
 144   fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
 145   fhetaTracklets->SetDirectory(0);
 146   fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100, 0., 2*TMath::Pi());
 147   fhphiTracklets->SetDirectory(0);
 148   fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
 149   fhetaClustersLay1->SetDirectory(0);
 150   fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
 151   fhphiClustersLay1->SetDirectory(0);
 152 }
 153
 154 //______________________________________________________________________
 155 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
 156 fClustersLay1(mr.fClustersLay1),
 157 fClustersLay2(mr.fClustersLay2),
 158 fTracklets(mr.fTracklets),
 159 fSClusters(mr.fSClusters),
 160 fAssociationFlag(mr.fAssociationFlag),
 161 fNClustersLay1(mr.fNClustersLay1),
 162 fNClustersLay2(mr.fNClustersLay2),
 163 fNTracklets(mr.fNTracklets),
 164 fNSingleCluster(mr.fNSingleCluster),
 165 fPhiWindow(mr.fPhiWindow),
 166 fZetaWindow(mr.fZetaWindow),
 167 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
 168 fHistOn(mr.fHistOn),
 169 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
 170 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
 171 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
 172 fhClustersDPhiAll(mr.fhClustersDPhiAll),
 173 fhClustersDThetaAll(mr.fhClustersDThetaAll),
 174 fhClustersDZetaAll(mr.fhClustersDZetaAll),
 175 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
 176 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
 177 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
 178 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
 179 fhetaTracklets(mr.fhetaTracklets),
 180 fhphiTracklets(mr.fhphiTracklets),
 181 fhetaClustersLay1(mr.fhetaClustersLay1),
 182 fhphiClustersLay1(mr.fhphiClustersLay1) {
 183   // Copy constructor
 184
 185 }
 186
 187 //______________________________________________________________________
 188 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
 189   // Assignment operator
 190   this->~AliITSMultReconstructor();
 191   new(this) AliITSMultReconstructor(mr);
 192   return *this;
 193 }
 194
 195 //______________________________________________________________________
 196 AliITSMultReconstructor::~AliITSMultReconstructor(){
 197   // Destructor
 198
 199   // delete histograms
 200   delete fhClustersDPhiAcc;
 201   delete fhClustersDThetaAcc;
 202   delete fhClustersDZetaAcc;
 203   delete fhClustersDPhiAll;
 204   delete fhClustersDThetaAll;
 205   delete fhClustersDZetaAll;
 206   delete fhDPhiVsDThetaAll;
 207   delete fhDPhiVsDThetaAcc;
 208   delete fhDPhiVsDZetaAll;
 209   delete fhDPhiVsDZetaAcc;
 210   delete fhetaTracklets;
 211   delete fhphiTracklets;
 212   delete fhetaClustersLay1;
 213   delete fhphiClustersLay1;
 214
 215   // delete arrays
 216   for(Int_t i=0; i<300000; i++) {
 217     delete [] fClustersLay1[i];
 218     delete [] fClustersLay2[i];
 219     delete [] fTracklets[i];
 220     delete [] fSClusters[i];
 221   }
 222   delete [] fClustersLay1;
 223   delete [] fClustersLay2;
 224   delete [] fTracklets;
 225   delete [] fSClusters;
 226
 227   delete [] fAssociationFlag;
 228 }
 229
 230 //____________________________________________________________________
 231 void
 232 AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
 233   //
 234   // - calls LoadClusterArray that finds the position of the clusters
 235   //   (in global coord)
 236   // - convert the cluster coordinates to theta, phi (seen from the
 237   //   interaction vertex). The third coordinate is used for ....
 238   // - makes an array of tracklets
 239   //
 240   // After this method has been called, the clusters of the two layers
 241   // and the tracklets can be retrieved by calling the Get'er methods.
 242
 243   // reset counters
 244   fNClustersLay1 = 0;
 245   fNClustersLay2 = 0;
 246   fNTracklets = 0;
 247   fNSingleCluster = 0;
 248   // loading the clusters
 249   LoadClusterArrays(clusterTree);
 250
 251   // find the tracklets
 252   AliDebug(1,"Looking for tracklets... ");
 253
 254   //###########################################################
 255   // Loop on layer 1 : finding theta, phi and z
 256   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 257     Float_t x = fClustersLay1[iC1][0] - vtx[0];
 258     Float_t y = fClustersLay1[iC1][1] - vtx[1];
 259     Float_t z = fClustersLay1[iC1][2] - vtx[2];
 260
 261     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
 262                                TMath::Power(y,2) +
 263                                TMath::Power(z,2));
 264
 265     fClustersLay1[iC1][0] = TMath::ACos(z/r);                   // Store Theta
 266     fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 267     fClustersLay1[iC1][2] = z/r;                                // Store scaled z
 268     if (fHistOn) {
 269       Float_t eta=fClustersLay1[iC1][0];
 270       eta= TMath::Tan(eta/2.);
 271       eta=-TMath::Log(eta);
 272       fhetaClustersLay1->Fill(eta);
 273       fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
 274     }
 275   }
 276
 277   // Loop on layer 2 : finding theta, phi and r
 278   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 279     Float_t x = fClustersLay2[iC2][0] - vtx[0];
 280     Float_t y = fClustersLay2[iC2][1] - vtx[1];
 281     Float_t z = fClustersLay2[iC2][2] - vtx[2];
 282
 283     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
 284                                TMath::Power(y,2) +
 285                                TMath::Power(z,2));
 286
 287     fClustersLay2[iC2][0] = TMath::ACos(z/r);                   // Store Theta
 288     fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 289     fClustersLay2[iC2][2] = z;                                  // Store z
 290
 291  // this only needs to be initialized for the fNClustersLay2 first associations
 292     fAssociationFlag[iC2] = kFALSE;
 293   }
 294
 295   //###########################################################
 296   // Loop on layer 1
 297   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 298
 299     // reset of variables for multiple candidates
 300     Int_t  iC2WithBestDist = 0;     // reset
 301     Float_t distmin        = 100.;  // just to put a huge number!
 302     Float_t dPhimin        = 0.;  // Used for histograms only!
 303     Float_t dThetamin      = 0.;  // Used for histograms only!
 304     Float_t dZetamin       = 0.;  // Used for histograms only!
 305
 306     // Loop on layer 2
 307     for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 308
 309       // The following excludes double associations
 310       if (!fAssociationFlag[iC2]) {
 311
 312         // find the difference in angles
 313         Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
 314         Float_t dPhi   = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
 315         // take into account boundary condition
 316         if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
 317
 318         // find the difference in z (between linear projection from layer 1
 319         // and the actual point: Dzeta= z1/r1*r2 -z2)
 320         Float_t r2   = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
 321         Float_t dZeta  = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
 322
 323         if (fHistOn) {
 324           fhClustersDPhiAll->Fill(dPhi);
 325           fhClustersDThetaAll->Fill(dTheta);
 326           fhClustersDZetaAll->Fill(dZeta);
 327           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
 328           fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
 329         }
 330         // make "elliptical" cut in Phi and Zeta!
 331         Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
 332
 333         if (d>1) continue;
 334
 335         //look for the minimum distance: the minimum is in iC2WithBestDist
 336         if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
 337           distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
 338           dPhimin = dPhi;
 339           dThetamin = dTheta;
 340           dZetamin = dZeta;
 341           iC2WithBestDist = iC2;
 342         }
 343       }
 344     } // end of loop over clusters in layer 2
 345
 346     if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
 347
 348       if (fHistOn) {
 349         fhClustersDPhiAcc->Fill(dPhimin);
 350         fhClustersDThetaAcc->Fill(dThetamin);
 351         fhClustersDZetaAcc->Fill(dZetamin);
 352         fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
 353         fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
 354       }
 355
 356       if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
 357
 358       // store the tracklet
 359
 360       // use the theta from the clusters in the first layer
 361       fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
 362       // use the phi from the clusters in the first layer
 363       fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
 364       // Store the difference between phi1 and phi2
 365       fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
 366
 367       // find label
 368       Int_t label1 = 0;
 369       Int_t label2 = 0;
 370       while (label2 < 3)
 371       {
 372         if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
 373           break;
 374         label1++;
 375         if (label1 == 3)
 376         {
 377           label1 = 0;
 378           label2++;
 379         }
 380       }
 381
 382       if (label2 < 3)
 383       {
 384         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));
 385         fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
 386       }
 387       else
 388       {
 389         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));
 390         fTracklets[fNTracklets][3] = -2;
 391       }
 392
 393       if (fHistOn) {
 394         Float_t eta=fTracklets[fNTracklets][0];
 395         eta= TMath::Tan(eta/2.);
 396         eta=-TMath::Log(eta);
 397         fhetaTracklets->Fill(eta);
 398         fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
 399       }
 400
 401       AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
 402       AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1,
 403                       iC2WithBestDist));
 404       fNTracklets++;
 405     }
 406
 407     // Delete the following else if you do not want to save Clusters!
 408
 409     else { // This means that the cluster has not been associated
 410
 411       // store the cluster
 412
 413       fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
 414       fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
 415       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
 416                       fNSingleCluster, iC1));
 417       fNSingleCluster++;
 418     }
 419
 420   } // end of loop over clusters in layer 1
 421
 422   AliDebug(1,Form("%d tracklets found", fNTracklets));
 423 }
 424
 425 //____________________________________________________________________
 426 void
 427 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
 428   // This method
 429   // - gets the clusters from the cluster tree
 430   // - convert them into global coordinates
 431   // - store them in the internal arrays
 432
 433   AliDebug(1,"Loading clusters ...");
 434
 435   fNClustersLay1 = 0;
 436   fNClustersLay2 = 0;
 437
 438   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
 439   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
 440
 441   itsClusterBranch->SetAddress(&itsClusters);
 442
 443   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
 444   Float_t cluGlo[3]={0.,0.,0.};
 445
 446   // loop over the its subdetectors
 447   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
 448
 449     if (!itsClusterTree->GetEvent(iIts))
 450       continue;
 451
 452     Int_t nClusters = itsClusters->GetEntriesFast();
 453
 454     // loop over clusters
 455     while(nClusters--) {
 456       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 457
 458       if (cluster->GetLayer()>1)
 459         continue;
 460
 461       cluster->GetGlobalXYZ(cluGlo);
 462       Float_t x = cluGlo[0];
 463       Float_t y = cluGlo[1];
 464       Float_t z = cluGlo[2];
 465
 466       if (cluster->GetLayer()==0) {
 467         fClustersLay1[fNClustersLay1][0] = x;
 468         fClustersLay1[fNClustersLay1][1] = y;
 469         fClustersLay1[fNClustersLay1][2] = z;
 470         for (Int_t i=0; i<3; i++)
 471                 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
 472         fNClustersLay1++;
 473       }
 474       if (cluster->GetLayer()==1) {
 475         fClustersLay2[fNClustersLay2][0] = x;
 476         fClustersLay2[fNClustersLay2][1] = y;
 477         fClustersLay2[fNClustersLay2][2] = z;
 478         for (Int_t i=0; i<3; i++)
 479                 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
 480         fNClustersLay2++;
 481       }
 482
 483     }// end of cluster loop
 484   } // end of its "subdetector" loop
 485
 486   if (itsClusters) {
 487     itsClusters->Delete();
 488     delete itsClusters;
 489     itsClusters = 0;
 490   }
 491   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
 492 }
 493 //____________________________________________________________________
 494 void
 495 AliITSMultReconstructor::SaveHists() {
 496   // This method save the histograms on the output file
 497   // (only if fHistOn is TRUE).
 498
 499   if (!fHistOn)
 500     return;
 501
 502   fhClustersDPhiAll->Write();
 503   fhClustersDThetaAll->Write();
 504   fhClustersDZetaAll->Write();
 505   fhDPhiVsDThetaAll->Write();
 506   fhDPhiVsDZetaAll->Write();
 507
 508   fhClustersDPhiAcc->Write();
 509   fhClustersDThetaAcc->Write();
 510   fhClustersDZetaAcc->Write();
 511   fhDPhiVsDThetaAcc->Write();
 512   fhDPhiVsDZetaAcc->Write();
 513
 514   fhetaTracklets->Write();
 515   fhphiTracklets->Write();
 516   fhetaClustersLay1->Write();
 517   fhphiClustersLay1->Write();
 518 }
 519
 520