1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //____________________________________________________________________
20 // AliITSMultReconstructor - find clusters in the pixels (theta and
21 // phi) and tracklets.
23 // These can be used to extract charged particles multiplcicity from the ITS.
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)
34 // Two methods return the number of traklets and the number of clusters
35 // in the first SPD layer (GetNTracklets GetNSingleClusters)
37 // -----------------------------------------------------------------
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.
43 // Author : Tiziano Virgili
47 //____________________________________________________________________
49 #include <TClonesArray.h>
54 #include "AliITSMultReconstructor.h"
55 #include "AliITSRecPoint.h"
56 #include "AliITSgeom.h"
59 //____________________________________________________________________
60 ClassImp(AliITSMultReconstructor)
63 //____________________________________________________________________
64 AliITSMultReconstructor::AliITSMultReconstructor():
77 fOnlyOneTrackletPerC2(0),
80 fhClustersDThetaAcc(0),
81 fhClustersDZetaAcc(0),
83 fhClustersDThetaAll(0),
84 fhClustersDZetaAll(0),
93 // Method to reconstruct the charged particles multiplicity with the
101 SetOnlyOneTrackletPerC2();
103 fClustersLay1 = new Float_t*[300000];
104 fClustersLay2 = new Float_t*[300000];
105 fTracklets = new Float_t*[300000];
106 fSClusters = new Float_t*[300000];
107 fAssociationFlag = new Bool_t[300000];
109 for(Int_t i=0; i<300000; i++) {
110 fClustersLay1[i] = new Float_t[6];
111 fClustersLay2[i] = new Float_t[6];
112 fTracklets[i] = new Float_t[4];
113 fSClusters[i] = new Float_t[2];
114 fAssociationFlag[i] = kFALSE;
117 // definition of histograms
118 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
119 fhClustersDPhiAcc->SetDirectory(0);
120 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
121 fhClustersDThetaAcc->SetDirectory(0);
122 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
123 fhClustersDZetaAcc->SetDirectory(0);
125 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
126 fhDPhiVsDZetaAcc->SetDirectory(0);
127 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
128 fhDPhiVsDThetaAcc->SetDirectory(0);
130 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
131 fhClustersDPhiAll->SetDirectory(0);
132 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
133 fhClustersDThetaAll->SetDirectory(0);
134 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
135 fhClustersDZetaAll->SetDirectory(0);
137 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
138 fhDPhiVsDZetaAll->SetDirectory(0);
139 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
140 fhDPhiVsDThetaAll->SetDirectory(0);
142 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
143 fhphiTracklets = new TH1F("phiTracklets", "phi", 100,-3.14159,3.14159);
144 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
145 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100,-3.141,3.141);
149 //______________________________________________________________________
150 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
151 fGeometry(mr.fGeometry),
152 fClustersLay1(mr.fClustersLay1),
153 fClustersLay2(mr.fClustersLay2),
154 fTracklets(mr.fTracklets),
155 fSClusters(mr.fSClusters),
156 fAssociationFlag(mr.fAssociationFlag),
157 fNClustersLay1(mr.fNClustersLay1),
158 fNClustersLay2(mr.fNClustersLay2),
159 fNTracklets(mr.fNTracklets),
160 fNSingleCluster(mr.fNSingleCluster),
161 fPhiWindow(mr.fPhiWindow),
162 fZetaWindow(mr.fZetaWindow),
163 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
165 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
166 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
167 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
168 fhClustersDPhiAll(mr.fhClustersDPhiAll),
169 fhClustersDThetaAll(mr.fhClustersDThetaAll),
170 fhClustersDZetaAll(mr.fhClustersDZetaAll),
171 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
172 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
173 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
174 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
175 fhetaTracklets(mr.fhetaTracklets),
176 fhphiTracklets(mr.fhphiTracklets),
177 fhetaClustersLay1(mr.fhetaClustersLay1),
178 fhphiClustersLay1(mr.fhphiClustersLay1) {
183 //______________________________________________________________________
184 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
185 // Assignment operator
186 this->~AliITSMultReconstructor();
187 new(this) AliITSMultReconstructor(mr);
191 //______________________________________________________________________
192 AliITSMultReconstructor::~AliITSMultReconstructor(){
196 delete fhClustersDPhiAcc;
197 delete fhClustersDThetaAcc;
198 delete fhClustersDZetaAcc;
199 delete fhClustersDPhiAll;
200 delete fhClustersDThetaAll;
201 delete fhClustersDZetaAll;
202 delete fhDPhiVsDThetaAll;
203 delete fhDPhiVsDThetaAcc;
204 delete fhDPhiVsDZetaAll;
205 delete fhDPhiVsDZetaAcc;
206 delete fhetaTracklets;
207 delete fhphiTracklets;
208 delete fhetaClustersLay1;
209 delete fhphiClustersLay1;
212 for(Int_t i=0; i<300000; i++) {
213 delete [] fClustersLay1[i];
214 delete [] fClustersLay2[i];
215 delete [] fTracklets[i];
216 delete [] fSClusters[i];
218 delete [] fClustersLay1;
219 delete [] fClustersLay2;
220 delete [] fTracklets;
221 delete [] fSClusters;
223 delete [] fAssociationFlag;
226 //____________________________________________________________________
228 AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
230 // - calls LoadClusterArray that finds the position of the clusters
232 // - convert the cluster coordinates to theta, phi (seen from the
233 // interaction vertex). The third coordinate is used for ....
234 // - makes an array of tracklets
236 // After this method has been called, the clusters of the two layers
237 // and the tracklets can be retrieved by calling the Get'er methods.
244 // loading the clusters
245 LoadClusterArrays(clusterTree);
247 // find the tracklets
248 AliDebug(1,"Looking for tracklets... ");
250 //###########################################################
251 // Loop on layer 1 : finding theta, phi and z
252 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
253 Float_t x = fClustersLay1[iC1][0] - vtx[0];
254 Float_t y = fClustersLay1[iC1][1] - vtx[1];
255 Float_t z = fClustersLay1[iC1][2] - vtx[2];
257 Float_t r = TMath::Sqrt(TMath::Power(x,2) +
261 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
262 fClustersLay1[iC1][1] = TMath::ATan2(y,x); // Store Phi
263 fClustersLay1[iC1][2] = z/r; // Store scaled z
265 Float_t eta=fClustersLay1[iC1][0];
266 eta= TMath::Tan(eta/2.);
267 eta=-TMath::Log(eta);
268 fhetaClustersLay1->Fill(eta);
269 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
273 // Loop on layer 2 : finding theta, phi and r
274 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
275 Float_t x = fClustersLay2[iC2][0] - vtx[0];
276 Float_t y = fClustersLay2[iC2][1] - vtx[1];
277 Float_t z = fClustersLay2[iC2][2] - vtx[2];
279 Float_t r = TMath::Sqrt(TMath::Power(x,2) +
283 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
284 fClustersLay2[iC2][1] = TMath::ATan2(y,x); // Store Phi
285 fClustersLay2[iC2][2] = z; // Store z
287 // this only needs to be initialized for the fNClustersLay2 first associations
288 fAssociationFlag[iC2] = kFALSE;
291 //###########################################################
293 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
295 // reset of variables for multiple candidates
296 Int_t iC2WithBestDist = 0; // reset
297 Float_t distmin = 100.; // just to put a huge number!
298 Float_t dPhimin = 0.; // Used for histograms only!
299 Float_t dThetamin = 0.; // Used for histograms only!
300 Float_t dZetamin = 0.; // Used for histograms only!
303 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
305 // The following excludes double associations
306 if (!fAssociationFlag[iC2]) {
308 // find the difference in angles
309 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
310 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
311 // take into account boundary condition
312 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
314 // find the difference in z (between linear projection from layer 1
315 // and the actual point: Dzeta= z1/r1*r2 -z2)
316 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
317 Float_t dZeta = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
320 fhClustersDPhiAll->Fill(dPhi);
321 fhClustersDThetaAll->Fill(dTheta);
322 fhClustersDZetaAll->Fill(dZeta);
323 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
324 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
326 // make "elliptical" cut in Phi and Zeta!
327 Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
331 //look for the minimum distance: the minimum is in iC2WithBestDist
332 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
333 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
337 iC2WithBestDist = iC2;
340 } // end of loop over clusters in layer 2
342 if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
345 fhClustersDPhiAcc->Fill(dPhimin);
346 fhClustersDThetaAcc->Fill(dThetamin);
347 fhClustersDZetaAcc->Fill(dZetamin);
348 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
349 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
352 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
354 // store the tracklet
356 // use the theta from the clusters in the first layer
357 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
358 // use the phi from the clusters in the first layer
359 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
360 // Store the difference between phi1 and phi2
361 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
368 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
381 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));
382 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
386 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));
387 fTracklets[fNTracklets][3] = -2;
391 Float_t eta=fTracklets[fNTracklets][0];
392 eta= TMath::Tan(eta/2.);
393 eta=-TMath::Log(eta);
394 fhetaTracklets->Fill(eta);
395 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
398 AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
399 AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1,
404 // Delete the following else if you do not want to save Clusters!
406 else { // This means that the cluster has not been associated
410 fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
411 fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
412 AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
413 fNSingleCluster, iC1));
417 } // end of loop over clusters in layer 1
419 AliDebug(1,Form("%d tracklets found", fNTracklets));
422 //____________________________________________________________________
424 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
426 // - gets the clusters from the cluster tree
427 // - convert them into global coordinates
428 // - store them in the internal arrays
430 AliDebug(1,"Loading clusters ...");
435 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
436 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
438 itsClusterBranch->SetAddress(&itsClusters);
440 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
442 // loop over the its subdetectors
443 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
445 if (!itsClusterTree->GetEvent(iIts))
448 Int_t nClusters = itsClusters->GetEntriesFast();
450 // stuff needed to get the global coordinates
451 Double_t rot[9]; fGeometry->GetRotMatrix(iIts,rot);
452 Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det);
453 Float_t tx,ty,tz; fGeometry->GetTrans(lay,lad,det,tx,ty,tz);
456 // "alpha" is the angle from the global X-axis to the
457 // local GEANT X'-axis ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
458 // "phi" is the angle from the global X-axis to the
459 // local cluster X"-axis
461 Double_t alpha = TMath::ATan2(rot[1],rot[0])+TMath::Pi();
462 Double_t itsPhi = TMath::Pi()/2+alpha;
464 if (lay==1) itsPhi+=TMath::Pi();
465 Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi);
466 Double_t r=tx*cp+ty*sp;
468 // loop over clusters
470 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
472 if (cluster->GetLayer()>1)
475 Float_t x = r*cp - cluster->GetY()*sp;
476 Float_t y = r*sp + cluster->GetY()*cp;
477 Float_t z = cluster->GetZ();
479 if (cluster->GetLayer()==0) {
480 fClustersLay1[fNClustersLay1][0] = x;
481 fClustersLay1[fNClustersLay1][1] = y;
482 fClustersLay1[fNClustersLay1][2] = z;
483 for (Int_t i=0; i<3; i++)
484 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
487 if (cluster->GetLayer()==1) {
488 fClustersLay2[fNClustersLay2][0] = x;
489 fClustersLay2[fNClustersLay2][1] = y;
490 fClustersLay2[fNClustersLay2][2] = z;
491 for (Int_t i=0; i<3; i++)
492 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
496 }// end of cluster loop
497 } // end of its "subdetector" loop
500 itsClusters->Delete();
504 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
506 //____________________________________________________________________
508 AliITSMultReconstructor::SaveHists() {
509 // This method save the histograms on the output file
510 // (only if fHistOn is TRUE).
515 fhClustersDPhiAll->Write();
516 fhClustersDThetaAll->Write();
517 fhClustersDZetaAll->Write();
518 fhDPhiVsDThetaAll->Write();
519 fhDPhiVsDZetaAll->Write();
521 fhClustersDPhiAcc->Write();
522 fhClustersDThetaAcc->Write();
523 fhClustersDZetaAcc->Write();
524 fhDPhiVsDThetaAcc->Write();
525 fhDPhiVsDZetaAcc->Write();
527 fhetaTracklets->Write();
528 fhphiTracklets->Write();
529 fhetaClustersLay1->Write();
530 fhphiClustersLay1->Write();