ITS/AliITSMultReconstructor.cxx

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