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 // TODO:
  19 //
  20 // - Introduce a rough pt estimation from the difference in phi ?
  21 // - Allow for a more refined selection criterium in case of multiple
  22 //   candidates (for instance by introducing weights for the difference
  23 //   in Phi and Zeta).
  24 //
  25 //____________________________________________________________________
  26
  27 #include "AliITSMultReconstructor.h"
  28
  29 #include "TTree.h"
  30 #include "TH1F.h"
  31 #include "TH2F.h"
  32
  33
  34 #include "AliITSRecPoint.h"
  35 #include "AliITSgeom.h"
  36 #include "AliLog.h"
  37
  38 //____________________________________________________________________
  39 ClassImp(AliITSMultReconstructor)
  40
  41 //____________________________________________________________________
  42 AliITSMultReconstructor::AliITSMultReconstructor() {
  43
  44   fGeometry =0;
  45
  46   SetHistOn();
  47   SetPhiWindow();
  48   SetZetaWindow();
  49   SetOnlyOneTrackletPerC2();
  50
  51   fClustersLay1       = new Float_t*[300000];
  52   fClustersLay2       = new Float_t*[300000];
  53   fTracklets          = new Float_t*[300000];
  54   fAssociationFlag    = new Bool_t[300000];
  55
  56   for(Int_t i=0; i<300000; i++) {
  57     fClustersLay1[i]       = new Float_t[3];
  58     fClustersLay2[i]       = new Float_t[3];
  59     fTracklets[i]          = new Float_t[3];
  60     fAssociationFlag[i]    = kFALSE;
  61   }
  62
  63   // definition of histograms
  64   fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,-0.1,0.1);
  65   fhClustersDPhiAcc->SetDirectory(0);
  66   fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
  67   fhClustersDThetaAcc->SetDirectory(0);
  68   fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
  69   fhClustersDZetaAcc->SetDirectory(0);
  70
  71   fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,-0.1,0.1);
  72   fhDPhiVsDZetaAcc->SetDirectory(0);
  73   fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
  74   fhDPhiVsDThetaAcc->SetDirectory(0);
  75
  76   fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,-0.5,0.5);
  77   fhClustersDPhiAll->SetDirectory(0);
  78   fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
  79   fhClustersDThetaAll->SetDirectory(0);
  80   fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
  81   fhClustersDZetaAll->SetDirectory(0);
  82
  83   fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,-0.5,0.5);
  84   fhDPhiVsDZetaAll->SetDirectory(0);
  85   fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,-0.5,0.5);
  86   fhDPhiVsDThetaAll->SetDirectory(0);
  87
  88   fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
  89   fhetaTracklets->SetDirectory(0);
  90   fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100,-3.14159,3.14159);
  91   fhphiTracklets->SetDirectory(0);
  92   fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
  93   fhetaClustersLay1->SetDirectory(0);
  94   fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100,-3.141,3.141);
  95   fhphiClustersLay1->SetDirectory(0);
  96 }
  97 //____________________________________________________________________
  98 AliITSMultReconstructor::~AliITSMultReconstructor() {
  99   // Destructor
 100
 101   fGeometry = 0x0;
 102   for(Int_t i=0; i<300000; i++) {
 103     delete [] fClustersLay1[i];
 104     delete [] fClustersLay2[i];
 105     delete [] fTracklets[i];
 106   }
 107
 108   delete [] fClustersLay1;
 109   delete [] fClustersLay2;
 110   delete [] fTracklets;
 111   delete [] fAssociationFlag;
 112
 113   delete fhClustersDPhiAcc;
 114   delete fhClustersDThetaAcc;
 115   delete fhClustersDZetaAcc;
 116   delete fhClustersDPhiAll;
 117   delete fhClustersDThetaAll;
 118   delete fhClustersDZetaAll;
 119
 120   delete fhDPhiVsDThetaAll;
 121   delete fhDPhiVsDThetaAcc;
 122   delete fhDPhiVsDZetaAll;
 123   delete fhDPhiVsDZetaAcc;
 124
 125   delete fhetaTracklets;
 126   delete fhphiTracklets;
 127   delete fhetaClustersLay1;
 128   delete fhphiClustersLay1;
 129 }
 130
 131 //____________________________________________________________________
 132 void
 133 AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
 134   //
 135   // - calls LoadClusterArray that finds the position of the clusters
 136   //   (in global coord)
 137   // - convert the cluster coordinates to theta, phi (seen from the
 138   //   interaction vertex). The third coordinate is used for ....
 139   // - makes an array of tracklets
 140   //
 141   // After this method has been called, the clusters of the two layers
 142   // and the tracklets can be retrieved by calling the Get'er methods.
 143
 144   // reset counters
 145   fNClustersLay1 = 0;
 146   fNClustersLay2 = 0;
 147   fNTracklets = 0;
 148
 149   // loading the clusters
 150   LoadClusterArrays(clusterTree);
 151
 152   // find the tracklets
 153   AliDebug(1,"Looking for tracklets... ");
 154
 155   //###########################################################
 156   // Loop on layer 1 : finding theta, phi and z
 157   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 158     Float_t x = fClustersLay1[iC1][0] - vtx[0];
 159     Float_t y = fClustersLay1[iC1][1] - vtx[1];
 160     Float_t z = fClustersLay1[iC1][2] - vtx[2];
 161
 162     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
 163                                TMath::Power(y,2) +
 164                                TMath::Power(z,2));
 165
 166     fClustersLay1[iC1][0] = TMath::ACos(z/r);  // Store Theta
 167     fClustersLay1[iC1][1] = TMath::ATan2(x,y);  // Store Phi
 168     fClustersLay1[iC1][2] = z/r;               // Store scaled z
 169     if (fHistOn) {
 170       Float_t eta=fClustersLay1[iC1][0];
 171       eta= TMath::Tan(eta/2.);
 172       eta=-TMath::Log(eta);
 173       fhetaClustersLay1->Fill(eta);
 174       fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
 175     }
 176 }
 177
 178   // Loop on layer 2 : finding theta, phi and r
 179   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 180     Float_t x = fClustersLay2[iC2][0] - vtx[0];
 181     Float_t y = fClustersLay2[iC2][1] - vtx[1];
 182     Float_t z = fClustersLay2[iC2][2] - vtx[2];
 183
 184     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
 185                                TMath::Power(y,2) +
 186                                TMath::Power(z,2));
 187
 188     fClustersLay2[iC2][0] = TMath::ACos(z/r);  // Store Theta
 189     fClustersLay2[iC2][1] = TMath::ATan2(x,y);  // Store Phi
 190     fClustersLay2[iC2][2] = z;                 // Store z
 191
 192  // this only needs to be initialized for the fNClustersLay2 first associations
 193     fAssociationFlag[iC2] = kFALSE;
 194   }
 195
 196   //###########################################################
 197   // Loop on layer 1
 198   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 199
 200     // reset of variables for multiple candidates
 201     Int_t  iC2WithBestDist = 0;     // reset
 202     Float_t Distmin        = 100.;  // just to put a huge number!
 203     Float_t dPhimin        = 0.;  // Used for histograms only!
 204     Float_t dThetamin      = 0.;  // Used for histograms only!
 205     Float_t dZetamin       = 0.;  // Used for histograms only!
 206
 207     // Loop on layer 2
 208     for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 209
 210       // The following excludes double associations
 211       if (!fAssociationFlag[iC2]) {
 212
 213         // find the difference in angles
 214         Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
 215         Float_t dPhi   = fClustersLay2[iC2][1] - fClustersLay1[iC1][1];
 216
 217         // find the difference in z (between linear projection from layer 1
 218         // and the actual point: Dzeta= z1/r1*r2 -z2)
 219         Float_t r2   = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
 220         Float_t dZeta  = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
 221
 222         if (fHistOn) {
 223           fhClustersDPhiAll->Fill(dPhi);
 224           fhClustersDThetaAll->Fill(dTheta);
 225           fhClustersDZetaAll->Fill(dZeta);
 226           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
 227           fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
 228         }
 229         // make "elliptical" cut in Phi and Zeta!
 230         Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
 231         if (d>1) continue;
 232
 233         //look for the minimum distance: the minimum is in iC2WithBestDist
 234         if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < Distmin ) {
 235           Distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
 236           dPhimin = dPhi;
 237           dThetamin = dTheta;
 238           dZetamin = dZeta;
 239           iC2WithBestDist = iC2;
 240         }
 241       }
 242     } // end of loop over clusters in layer 2
 243
 244     if (Distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
 245
 246         if (fHistOn) {
 247           fhClustersDPhiAcc->Fill(dPhimin);
 248           fhClustersDThetaAcc->Fill(dThetamin);
 249           fhClustersDZetaAcc->Fill(dZetamin);
 250           fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
 251           fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
 252         }
 253
 254       if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
 255
 256       // store the tracklet
 257
 258       // use the theta from the clusters in the first layer
 259       fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
 260       // use the phi from the clusters in the first layer
 261       fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
 262       // Store the difference between phi1 and phi2
 263       fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
 264
 265         if (fHistOn) {
 266           Float_t eta=fTracklets[fNTracklets][0];
 267           eta= TMath::Tan(eta/2.);
 268           eta=-TMath::Log(eta);
 269           fhetaTracklets->Fill(eta);
 270           fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
 271         }
 272       fNTracklets++;
 273
 274       AliDebug(1,Form(" Adding tracklet candidate %d (cluster %d  of layer 1 and %d  of layer 2)", fNTracklets, iC1));
 275     }
 276   } // end of loop over clusters in layer 1
 277
 278   AliDebug(1,Form("%d tracklets found", fNTracklets));
 279 }
 280
 281 //____________________________________________________________________
 282 void
 283 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
 284   // This method
 285   // - gets the clusters from the cluster tree
 286   // - convert them into global coordinates
 287   // - store them in the internal arrays
 288
 289   AliDebug(1,"Loading clusters ...");
 290
 291   fNClustersLay1 = 0;
 292   fNClustersLay2 = 0;
 293
 294   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
 295   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
 296
 297   itsClusterBranch->SetAddress(&itsClusters);
 298
 299   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
 300
 301   // loop over the its subdetectors
 302   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
 303
 304     if (!itsClusterTree->GetEvent(iIts))
 305       continue;
 306
 307     Int_t nClusters = itsClusters->GetEntriesFast();
 308
 309     // stuff needed to get the global coordinates
 310     Double_t rot[9];   fGeometry->GetRotMatrix(iIts,rot);
 311     Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det);
 312     Float_t tx,ty,tz;  fGeometry->GetTrans(lay,lad,det,tx,ty,tz);
 313
 314     // Below:
 315     // "alpha" is the angle from the global X-axis to the
 316     //         local GEANT X'-axis  ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
 317     // "phi" is the angle from the global X-axis to the
 318     //       local cluster X"-axis
 319
 320     Double_t alpha   = TMath::ATan2(rot[1],rot[0])+TMath::Pi();
 321     Double_t itsPhi = TMath::Pi()/2+alpha;
 322
 323     if (lay==1) itsPhi+=TMath::Pi();
 324     Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi);
 325     Double_t r=tx*cp+ty*sp;
 326
 327     // loop over clusters
 328     while(nClusters--) {
 329       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 330
 331       if (cluster->GetLayer()>1)
 332         continue;
 333
 334       Float_t x = r*cp - cluster->GetY()*sp;
 335       Float_t y = r*sp + cluster->GetY()*cp;
 336       Float_t z = cluster->GetZ();
 337
 338       if (cluster->GetLayer()==0) {
 339         fClustersLay1[fNClustersLay1][0] = x;
 340         fClustersLay1[fNClustersLay1][1] = y;
 341         fClustersLay1[fNClustersLay1][2] = z;
 342         fNClustersLay1++;
 343       }
 344       if (cluster->GetLayer()==1) {
 345         fClustersLay2[fNClustersLay2][0] = x;
 346         fClustersLay2[fNClustersLay2][1] = y;
 347         fClustersLay2[fNClustersLay2][2] = z;
 348         fNClustersLay2++;
 349       }
 350
 351     }// end of cluster loop
 352   } // end of its "subdetector" loop
 353
 354   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
 355 }
 356 //____________________________________________________________________
 357 void
 358 AliITSMultReconstructor::SaveHists() {
 359
 360   if (!fHistOn)
 361     return;
 362
 363   fhClustersDPhiAll->Write();
 364   fhClustersDThetaAll->Write();
 365   fhClustersDZetaAll->Write();
 366   fhDPhiVsDThetaAll->Write();
 367   fhDPhiVsDZetaAll->Write();
 368
 369   fhClustersDPhiAcc->Write();
 370   fhClustersDThetaAcc->Write();
 371   fhClustersDZetaAcc->Write();
 372   fhDPhiVsDThetaAcc->Write();
 373   fhDPhiVsDZetaAcc->Write();
 374
 375   fhetaTracklets->Write();
 376   fhphiTracklets->Write();
 377   fhetaClustersLay1->Write();
 378   fhphiClustersLay1->Write();
 379 }