1 /**************************************************************************
2 * Copyright(c) 2007-2009, 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
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()
51 //____________________________________________________________________
53 #include <TClonesArray.h>
58 #include "AliITSMultReconstructor.h"
59 #include "AliITSRecPoint.h"
60 #include "AliITSgeom.h"
63 //____________________________________________________________________
64 ClassImp(AliITSMultReconstructor)
67 //____________________________________________________________________
68 AliITSMultReconstructor::AliITSMultReconstructor():
80 fOnlyOneTrackletPerC2(0),
83 fhClustersDThetaAcc(0),
84 fhClustersDZetaAcc(0),
86 fhClustersDThetaAll(0),
87 fhClustersDZetaAll(0),
96 // Method to reconstruct the charged particles multiplicity with the
103 SetOnlyOneTrackletPerC2();
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];
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;
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);
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);
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);
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);
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);
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),
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) {
187 //______________________________________________________________________
188 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
189 // Assignment operator
190 this->~AliITSMultReconstructor();
191 new(this) AliITSMultReconstructor(mr);
195 //______________________________________________________________________
196 AliITSMultReconstructor::~AliITSMultReconstructor(){
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;
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];
222 delete [] fClustersLay1;
223 delete [] fClustersLay2;
224 delete [] fTracklets;
225 delete [] fSClusters;
227 delete [] fAssociationFlag;
230 //____________________________________________________________________
232 AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
234 // - calls LoadClusterArray that finds the position of the clusters
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
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.
248 // loading the clusters
249 LoadClusterArrays(clusterTree);
251 // find the tracklets
252 AliDebug(1,"Looking for tracklets... ");
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];
261 Float_t r = TMath::Sqrt(TMath::Power(x,2) +
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
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]);
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];
283 Float_t r = TMath::Sqrt(TMath::Power(x,2) +
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
291 // this only needs to be initialized for the fNClustersLay2 first associations
292 fAssociationFlag[iC2] = kFALSE;
295 //###########################################################
297 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
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!
307 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
309 // The following excludes double associations
310 if (!fAssociationFlag[iC2]) {
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;
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];
324 fhClustersDPhiAll->Fill(dPhi);
325 fhClustersDThetaAll->Fill(dTheta);
326 fhClustersDZetaAll->Fill(dZeta);
327 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
328 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
330 // make "elliptical" cut in Phi and Zeta!
331 Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
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));
341 iC2WithBestDist = iC2;
344 } // end of loop over clusters in layer 2
346 if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
349 fhClustersDPhiAcc->Fill(dPhimin);
350 fhClustersDThetaAcc->Fill(dThetamin);
351 fhClustersDZetaAcc->Fill(dZetamin);
352 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
353 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
356 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
358 // store the tracklet
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];
372 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
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];
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;
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]);
401 AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
402 AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1,
407 // Delete the following else if you do not want to save Clusters!
409 else { // This means that the cluster has not been associated
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));
420 } // end of loop over clusters in layer 1
422 AliDebug(1,Form("%d tracklets found", fNTracklets));
425 //____________________________________________________________________
427 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
429 // - gets the clusters from the cluster tree
430 // - convert them into global coordinates
431 // - store them in the internal arrays
433 AliDebug(1,"Loading clusters ...");
438 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
439 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
441 itsClusterBranch->SetAddress(&itsClusters);
443 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
444 Float_t cluGlo[3]={0.,0.,0.};
446 // loop over the its subdetectors
447 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
449 if (!itsClusterTree->GetEvent(iIts))
452 Int_t nClusters = itsClusters->GetEntriesFast();
454 // loop over clusters
456 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
458 if (cluster->GetLayer()>1)
461 cluster->GetGlobalXYZ(cluGlo);
462 Float_t x = cluGlo[0];
463 Float_t y = cluGlo[1];
464 Float_t z = cluGlo[2];
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);
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);
483 }// end of cluster loop
484 } // end of its "subdetector" loop
487 itsClusters->Delete();
491 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
493 //____________________________________________________________________
495 AliITSMultReconstructor::SaveHists() {
496 // This method save the histograms on the output file
497 // (only if fHistOn is TRUE).
502 fhClustersDPhiAll->Write();
503 fhClustersDThetaAll->Write();
504 fhClustersDZetaAll->Write();
505 fhDPhiVsDThetaAll->Write();
506 fhDPhiVsDZetaAll->Write();
508 fhClustersDPhiAcc->Write();
509 fhClustersDThetaAcc->Write();
510 fhClustersDZetaAcc->Write();
511 fhDPhiVsDThetaAcc->Write();
512 fhDPhiVsDZetaAcc->Write();
514 fhetaTracklets->Write();
515 fhphiTracklets->Write();
516 fhetaClustersLay1->Write();
517 fhphiClustersLay1->Write();