ITS/AliITSMultReconstructor.cxx

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