]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSMultReconstructor.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New classes for finding multiple vertices (in case of pile-up). They will be used...
[u/mrichter/AliRoot.git] / ITS / AliITSMultReconstructor.cxx
CommitLineData
7ca4655f 1/**************************************************************************
eefb3acc 2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
7ca4655f 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
ac903f1b 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.
de4c520e 30// The parameter AssociationChoice allows to control if two clusters
ac903f1b 31// in layer 2 can be associated to the same cluster in layer 1 or not.
02a95988 32// (TRUE means double associations exluded; default = TRUE)
ac903f1b 33//
968e8539 34// Two methods return the number of traklets and the number of clusters
35// in the first SPD layer (GetNTracklets GetNSingleClusters)
36//
ac903f1b 37// -----------------------------------------------------------------
38//
02a95988 39// NOTE: The cuts on phi and zeta depend on the interacting system (p-p
3ef75756 40// or Pb-Pb). Please, check the file AliITSMultReconstructor.h and be
41// sure that SetPhiWindow and SetZetaWindow are defined accordingly.
ac903f1b 42//
968e8539 43// Author : Tiziano Virgili
3ef75756 44//
f606f16a 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()
9b373e9a 49// - fGeometry removed
50// - number of fired chips on the two layers
7b116aa1 51// - option to avoid duplicates in the overlaps (fRemoveClustersFromOverlaps)
52// - options and fiducial cuts via AliITSRecoParam
ac903f1b 53//
54//____________________________________________________________________
55
7ca4655f 56#include <TClonesArray.h>
57#include <TH1F.h>
58#include <TH2F.h>
59#include <TTree.h>
ac903f1b 60
7ca4655f 61#include "AliITSMultReconstructor.h"
7b116aa1 62#include "AliITSReconstructor.h"
9b373e9a 63#include "AliITSsegmentationSPD.h"
b51872de 64#include "AliITSRecPoint.h"
ac903f1b 65#include "AliITSgeom.h"
66#include "AliLog.h"
67
68//____________________________________________________________________
0762f3a8 69ClassImp(AliITSMultReconstructor)
ac903f1b 70
3ef75756 71
ac903f1b 72//____________________________________________________________________
7537d03c 73AliITSMultReconstructor::AliITSMultReconstructor():
7537d03c 74fClustersLay1(0),
75fClustersLay2(0),
7b116aa1 76fDetectorIndexClustersLay1(0),
77fDetectorIndexClustersLay2(0),
78fOverlapFlagClustersLay1(0),
79fOverlapFlagClustersLay2(0),
7537d03c 80fTracklets(0),
968e8539 81fSClusters(0),
7537d03c 82fAssociationFlag(0),
83fNClustersLay1(0),
84fNClustersLay2(0),
85fNTracklets(0),
968e8539 86fNSingleCluster(0),
7b116aa1 87fOnlyOneTrackletPerC2(0),
7537d03c 88fPhiWindow(0),
89fZetaWindow(0),
7b116aa1 90fRemoveClustersFromOverlaps(0),
91fPhiOverlapCut(0),
92fZetaOverlapCut(0),
7537d03c 93fHistOn(0),
94fhClustersDPhiAcc(0),
95fhClustersDThetaAcc(0),
96fhClustersDZetaAcc(0),
97fhClustersDPhiAll(0),
98fhClustersDThetaAll(0),
99fhClustersDZetaAll(0),
100fhDPhiVsDThetaAll(0),
101fhDPhiVsDThetaAcc(0),
102fhDPhiVsDZetaAll(0),
103fhDPhiVsDZetaAcc(0),
104fhetaTracklets(0),
105fhphiTracklets(0),
106fhetaClustersLay1(0),
107fhphiClustersLay1(0){
9b373e9a 108
109 fNFiredChips[0] = 0;
110 fNFiredChips[1] = 0;
111
3ef75756 112 // Method to reconstruct the charged particles multiplicity with the
113 // SPD (tracklets).
ac903f1b 114
ac903f1b 115
116 SetHistOn();
ac903f1b 117
7b116aa1 118 if(AliITSReconstructor::GetRecoParam()) {
119 SetOnlyOneTrackletPerC2(AliITSReconstructor::GetRecoParam()->GetTrackleterOnlyOneTrackletPerC2());
120 SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow());
121 SetZetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaWindow());
122 SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps());
123 SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut());
124 SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut());
125 } else {
126 SetOnlyOneTrackletPerC2();
127 SetPhiWindow();
128 SetZetaWindow();
129 SetRemoveClustersFromOverlaps();
130 SetPhiOverlapCut();
131 SetZetaOverlapCut();
132 }
133
134
135 fClustersLay1 = new Float_t*[300000];
136 fClustersLay2 = new Float_t*[300000];
137 fDetectorIndexClustersLay1 = new Int_t[300000];
138 fDetectorIndexClustersLay2 = new Int_t[300000];
139 fOverlapFlagClustersLay1 = new Bool_t[300000];
140 fOverlapFlagClustersLay2 = new Bool_t[300000];
141 fTracklets = new Float_t*[300000];
142 fSClusters = new Float_t*[300000];
143 fAssociationFlag = new Bool_t[300000];
ac903f1b 144
145 for(Int_t i=0; i<300000; i++) {
de4c520e 146 fClustersLay1[i] = new Float_t[6];
147 fClustersLay2[i] = new Float_t[6];
0939e22a 148 fTracklets[i] = new Float_t[5];
968e8539 149 fSClusters[i] = new Float_t[2];
7b116aa1 150 fOverlapFlagClustersLay1[i] = kFALSE;
151 fOverlapFlagClustersLay2[i] = kFALSE;
ac903f1b 152 fAssociationFlag[i] = kFALSE;
153 }
154
155 // definition of histograms
02a95988 156 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
ddced3c8 157 fhClustersDPhiAcc->SetDirectory(0);
158 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
159 fhClustersDThetaAcc->SetDirectory(0);
160 fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
161 fhClustersDZetaAcc->SetDirectory(0);
162
02a95988 163 fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
ddced3c8 164 fhDPhiVsDZetaAcc->SetDirectory(0);
02a95988 165 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
ac903f1b 166 fhDPhiVsDThetaAcc->SetDirectory(0);
167
02a95988 168 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
ddced3c8 169 fhClustersDPhiAll->SetDirectory(0);
170 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
171 fhClustersDThetaAll->SetDirectory(0);
172 fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
173 fhClustersDZetaAll->SetDirectory(0);
174
02a95988 175 fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
ddced3c8 176 fhDPhiVsDZetaAll->SetDirectory(0);
02a95988 177 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
ddced3c8 178 fhDPhiVsDThetaAll->SetDirectory(0);
179
180 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
96c2c35d 181 fhetaTracklets->SetDirectory(0);
f606f16a 182 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
96c2c35d 183 fhphiTracklets->SetDirectory(0);
ddced3c8 184 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
96c2c35d 185 fhetaClustersLay1->SetDirectory(0);
f606f16a 186 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
96c2c35d 187 fhphiClustersLay1->SetDirectory(0);
ac903f1b 188}
ddced3c8 189
3ef75756 190//______________________________________________________________________
7537d03c 191AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
7537d03c 192fClustersLay1(mr.fClustersLay1),
193fClustersLay2(mr.fClustersLay2),
7b116aa1 194fDetectorIndexClustersLay1(mr.fDetectorIndexClustersLay1),
195fDetectorIndexClustersLay2(mr.fDetectorIndexClustersLay2),
196fOverlapFlagClustersLay1(mr.fOverlapFlagClustersLay1),
197fOverlapFlagClustersLay2(mr.fOverlapFlagClustersLay2),
7537d03c 198fTracklets(mr.fTracklets),
968e8539 199fSClusters(mr.fSClusters),
7537d03c 200fAssociationFlag(mr.fAssociationFlag),
201fNClustersLay1(mr.fNClustersLay1),
202fNClustersLay2(mr.fNClustersLay2),
203fNTracklets(mr.fNTracklets),
968e8539 204fNSingleCluster(mr.fNSingleCluster),
7b116aa1 205fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
7537d03c 206fPhiWindow(mr.fPhiWindow),
207fZetaWindow(mr.fZetaWindow),
7b116aa1 208fRemoveClustersFromOverlaps(mr.fRemoveClustersFromOverlaps),
209fPhiOverlapCut(mr.fPhiOverlapCut),
210fZetaOverlapCut(mr.fZetaOverlapCut),
7537d03c 211fHistOn(mr.fHistOn),
212fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
213fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
214fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
215fhClustersDPhiAll(mr.fhClustersDPhiAll),
216fhClustersDThetaAll(mr.fhClustersDThetaAll),
217fhClustersDZetaAll(mr.fhClustersDZetaAll),
218fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
219fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
220fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
221fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
222fhetaTracklets(mr.fhetaTracklets),
223fhphiTracklets(mr.fhphiTracklets),
224fhetaClustersLay1(mr.fhetaClustersLay1),
225fhphiClustersLay1(mr.fhphiClustersLay1) {
3ef75756 226 // Copy constructor
7537d03c 227
3ef75756 228}
229
230//______________________________________________________________________
7537d03c 231AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
3ef75756 232 // Assignment operator
7537d03c 233 this->~AliITSMultReconstructor();
234 new(this) AliITSMultReconstructor(mr);
3ef75756 235 return *this;
236}
237
238//______________________________________________________________________
239AliITSMultReconstructor::~AliITSMultReconstructor(){
240 // Destructor
1ba5b31c 241
242 // delete histograms
243 delete fhClustersDPhiAcc;
244 delete fhClustersDThetaAcc;
245 delete fhClustersDZetaAcc;
246 delete fhClustersDPhiAll;
247 delete fhClustersDThetaAll;
248 delete fhClustersDZetaAll;
249 delete fhDPhiVsDThetaAll;
250 delete fhDPhiVsDThetaAcc;
251 delete fhDPhiVsDZetaAll;
252 delete fhDPhiVsDZetaAcc;
253 delete fhetaTracklets;
254 delete fhphiTracklets;
255 delete fhetaClustersLay1;
256 delete fhphiClustersLay1;
257
258 // delete arrays
259 for(Int_t i=0; i<300000; i++) {
260 delete [] fClustersLay1[i];
261 delete [] fClustersLay2[i];
262 delete [] fTracklets[i];
968e8539 263 delete [] fSClusters[i];
ddced3c8 264 }
1ba5b31c 265 delete [] fClustersLay1;
266 delete [] fClustersLay2;
7b116aa1 267 delete [] fDetectorIndexClustersLay1;
268 delete [] fDetectorIndexClustersLay2;
269 delete [] fOverlapFlagClustersLay1;
270 delete [] fOverlapFlagClustersLay2;
1ba5b31c 271 delete [] fTracklets;
968e8539 272 delete [] fSClusters;
1ba5b31c 273
274 delete [] fAssociationFlag;
ddced3c8 275}
ac903f1b 276
277//____________________________________________________________________
278void
279AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
280 //
281 // - calls LoadClusterArray that finds the position of the clusters
282 // (in global coord)
283 // - convert the cluster coordinates to theta, phi (seen from the
284 // interaction vertex). The third coordinate is used for ....
285 // - makes an array of tracklets
286 //
287 // After this method has been called, the clusters of the two layers
288 // and the tracklets can be retrieved by calling the Get'er methods.
289
ac903f1b 290 // reset counters
291 fNClustersLay1 = 0;
292 fNClustersLay2 = 0;
293 fNTracklets = 0;
968e8539 294 fNSingleCluster = 0;
ac903f1b 295 // loading the clusters
296 LoadClusterArrays(clusterTree);
3ef75756 297
ac903f1b 298 // find the tracklets
299 AliDebug(1,"Looking for tracklets... ");
300
301 //###########################################################
302 // Loop on layer 1 : finding theta, phi and z
303 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
304 Float_t x = fClustersLay1[iC1][0] - vtx[0];
305 Float_t y = fClustersLay1[iC1][1] - vtx[1];
306 Float_t z = fClustersLay1[iC1][2] - vtx[2];
ddced3c8 307
ac903f1b 308 Float_t r = TMath::Sqrt(TMath::Power(x,2) +
309 TMath::Power(y,2) +
310 TMath::Power(z,2));
311
eefb3acc 312 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
313 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
314 fClustersLay1[iC1][2] = z/r; // Store scaled z
ddced3c8 315 if (fHistOn) {
316 Float_t eta=fClustersLay1[iC1][0];
317 eta= TMath::Tan(eta/2.);
318 eta=-TMath::Log(eta);
319 fhetaClustersLay1->Fill(eta);
de4c520e 320 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
ddced3c8 321 }
96c2c35d 322 }
ac903f1b 323
324 // Loop on layer 2 : finding theta, phi and r
325 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
326 Float_t x = fClustersLay2[iC2][0] - vtx[0];
327 Float_t y = fClustersLay2[iC2][1] - vtx[1];
328 Float_t z = fClustersLay2[iC2][2] - vtx[2];
ddced3c8 329
ac903f1b 330 Float_t r = TMath::Sqrt(TMath::Power(x,2) +
331 TMath::Power(y,2) +
332 TMath::Power(z,2));
333
eefb3acc 334 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
335 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
336 fClustersLay2[iC2][2] = z; // Store z
ac903f1b 337
ddced3c8 338 // this only needs to be initialized for the fNClustersLay2 first associations
ac903f1b 339 fAssociationFlag[iC2] = kFALSE;
340 }
341
342 //###########################################################
343 // Loop on layer 1
344 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
345
346 // reset of variables for multiple candidates
ddced3c8 347 Int_t iC2WithBestDist = 0; // reset
3ef75756 348 Float_t distmin = 100.; // just to put a huge number!
ddced3c8 349 Float_t dPhimin = 0.; // Used for histograms only!
350 Float_t dThetamin = 0.; // Used for histograms only!
351 Float_t dZetamin = 0.; // Used for histograms only!
7b116aa1 352
353 if (fOverlapFlagClustersLay1[iC1]) continue;
354
ac903f1b 355 // Loop on layer 2
356 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
7b116aa1 357 if (fOverlapFlagClustersLay2[iC2]) continue;
ac903f1b 358 // The following excludes double associations
359 if (!fAssociationFlag[iC2]) {
360
361 // find the difference in angles
362 Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
02a95988 363 Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
364 // take into account boundary condition
365 if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
366
ac903f1b 367 // find the difference in z (between linear projection from layer 1
368 // and the actual point: Dzeta= z1/r1*r2 -z2)
ddced3c8 369 Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
de4c520e 370 Float_t dZeta = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
ddced3c8 371
372 if (fHistOn) {
373 fhClustersDPhiAll->Fill(dPhi);
374 fhClustersDThetaAll->Fill(dTheta);
375 fhClustersDZetaAll->Fill(dZeta);
ac903f1b 376 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
ddced3c8 377 fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
ac903f1b 378 }
379 // make "elliptical" cut in Phi and Zeta!
380 Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
3ef75756 381
ac903f1b 382 if (d>1) continue;
383
ddced3c8 384 //look for the minimum distance: the minimum is in iC2WithBestDist
3ef75756 385 if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
386 distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
ddced3c8 387 dPhimin = dPhi;
388 dThetamin = dTheta;
389 dZetamin = dZeta;
390 iC2WithBestDist = iC2;
ac903f1b 391 }
392 }
393 } // end of loop over clusters in layer 2
394
3ef75756 395 if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
396
397 if (fHistOn) {
de4c520e 398 fhClustersDPhiAcc->Fill(dPhimin);
3ef75756 399 fhClustersDThetaAcc->Fill(dThetamin);
400 fhClustersDZetaAcc->Fill(dZetamin);
401 fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
402 fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
403 }
ac903f1b 404
ddced3c8 405 if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
ac903f1b 406
407 // store the tracklet
408
de4c520e 409 // use the theta from the clusters in the first layer
ddced3c8 410 fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
de4c520e 411 // use the phi from the clusters in the first layer
ac903f1b 412 fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
35e2e4eb 413 // store the difference between phi1 and phi2
de4c520e 414 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
415
35e2e4eb 416 // define dphi in the range [0,pi] with proper sign (track charge correlated)
417 if (fTracklets[fNTracklets][2] > TMath::Pi())
418 fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]-2.*TMath::Pi();
419 if (fTracklets[fNTracklets][2] < -TMath::Pi())
420 fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]+2.*TMath::Pi();
421
de4c520e 422 // find label
0939e22a 423 // if equal label in both clusters found this label is assigned
424 // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
de4c520e 425 Int_t label1 = 0;
426 Int_t label2 = 0;
427 while (label2 < 3)
428 {
429 if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
430 break;
de4c520e 431 label1++;
432 if (label1 == 3)
433 {
434 label1 = 0;
435 label2++;
436 }
437 }
438
439 if (label2 < 3)
440 {
441 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));
442 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
0939e22a 443 fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3+label2];
de4c520e 444 }
445 else
446 {
447 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));
0939e22a 448 fTracklets[fNTracklets][3] = fClustersLay1[iC1][3];
449 fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3];
de4c520e 450 }
451
3ef75756 452 if (fHistOn) {
453 Float_t eta=fTracklets[fNTracklets][0];
454 eta= TMath::Tan(eta/2.);
455 eta=-TMath::Log(eta);
456 fhetaTracklets->Fill(eta);
457 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
458 }
ac903f1b 459
3ef75756 460 AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
461 AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1,
462 iC2WithBestDist));
463 fNTracklets++;
7b116aa1 464
465 if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (iC1,iC2WithBestDist);
466
ac903f1b 467 }
3ef75756 468
469 // Delete the following else if you do not want to save Clusters!
470
de4c520e 471 else { // This means that the cluster has not been associated
3ef75756 472
473 // store the cluster
474
968e8539 475 fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
476 fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
de4c520e 477 AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
968e8539 478 fNSingleCluster, iC1));
479 fNSingleCluster++;
3ef75756 480 }
481
ac903f1b 482 } // end of loop over clusters in layer 1
483
484 AliDebug(1,Form("%d tracklets found", fNTracklets));
485}
486
487//____________________________________________________________________
488void
489AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
490 // This method
491 // - gets the clusters from the cluster tree
492 // - convert them into global coordinates
493 // - store them in the internal arrays
9b373e9a 494 // - count the number of cluster-fired chips
ac903f1b 495
9b373e9a 496 AliDebug(1,"Loading clusters and cluster-fired chips ...");
ac903f1b 497
498 fNClustersLay1 = 0;
499 fNClustersLay2 = 0;
9b373e9a 500 fNFiredChips[0] = 0;
501 fNFiredChips[1] = 0;
ac903f1b 502
9b373e9a 503 AliITSsegmentationSPD *seg = new AliITSsegmentationSPD();
504
b51872de 505 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
506 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
ddced3c8 507
ac903f1b 508 itsClusterBranch->SetAddress(&itsClusters);
ddced3c8 509
ac903f1b 510 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
f606f16a 511 Float_t cluGlo[3]={0.,0.,0.};
ddced3c8 512
ac903f1b 513 // loop over the its subdetectors
514 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
515
516 if (!itsClusterTree->GetEvent(iIts))
517 continue;
518
519 Int_t nClusters = itsClusters->GetEntriesFast();
9b373e9a 520
521 // number of clusters in each chip of the current module
522 Int_t nClustersInChip[5] = {0,0,0,0,0};
523 Int_t layer = 0;
ac903f1b 524
ac903f1b 525 // loop over clusters
526 while(nClusters--) {
de4c520e 527 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
ac903f1b 528
9b373e9a 529 layer = cluster->GetLayer();
530 if (layer>1) continue;
ac903f1b 531
f606f16a 532 cluster->GetGlobalXYZ(cluGlo);
533 Float_t x = cluGlo[0];
534 Float_t y = cluGlo[1];
535 Float_t z = cluGlo[2];
9b373e9a 536
537 // find the chip for the current cluster
538 Float_t locz = cluster->GetDetLocalZ();
539 Int_t iChip = seg->GetChipFromLocal(0,locz);
540 nClustersInChip[iChip]++;
ac903f1b 541
9b373e9a 542 if (layer==0) {
ac903f1b 543 fClustersLay1[fNClustersLay1][0] = x;
544 fClustersLay1[fNClustersLay1][1] = y;
545 fClustersLay1[fNClustersLay1][2] = z;
7b116aa1 546
547 fDetectorIndexClustersLay1[fNClustersLay1]=iIts;
548
de4c520e 549 for (Int_t i=0; i<3; i++)
550 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
ac903f1b 551 fNClustersLay1++;
552 }
9b373e9a 553 if (layer==1) {
ac903f1b 554 fClustersLay2[fNClustersLay2][0] = x;
555 fClustersLay2[fNClustersLay2][1] = y;
556 fClustersLay2[fNClustersLay2][2] = z;
7b116aa1 557
558 fDetectorIndexClustersLay2[fNClustersLay2]=iIts;
559
de4c520e 560 for (Int_t i=0; i<3; i++)
561 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
ac903f1b 562 fNClustersLay2++;
563 }
564
565 }// end of cluster loop
9b373e9a 566
567 // get number of fired chips in the current module
568 if(layer<2)
569 for(Int_t ifChip=0; ifChip<5; ifChip++) {
570 if(nClustersInChip[ifChip] >= 1) fNFiredChips[layer]++;
571 }
572
ac903f1b 573 } // end of its "subdetector" loop
9b373e9a 574
cedf398d 575 if (itsClusters) {
576 itsClusters->Delete();
577 delete itsClusters;
9b373e9a 578 delete seg;
cedf398d 579 itsClusters = 0;
580 }
ac903f1b 581 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
9b373e9a 582 AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
583}
584//____________________________________________________________________
585void
586AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
587 // This method
588 // - gets the clusters from the cluster tree
589 // - counts the number of (cluster)fired chips
590
591 AliDebug(1,"Loading cluster-fired chips ...");
592
593 fNFiredChips[0] = 0;
594 fNFiredChips[1] = 0;
595
596 AliITSsegmentationSPD *seg = new AliITSsegmentationSPD();
597
598 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
599 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
600
601 itsClusterBranch->SetAddress(&itsClusters);
602
603 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
604
605 // loop over the its subdetectors
606 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
607
608 if (!itsClusterTree->GetEvent(iIts))
609 continue;
610
611 Int_t nClusters = itsClusters->GetEntriesFast();
612
613 // number of clusters in each chip of the current module
614 Int_t nClustersInChip[5] = {0,0,0,0,0};
615 Int_t layer = 0;
616
617 // loop over clusters
618 while(nClusters--) {
619 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
620
621 layer = cluster->GetLayer();
622 if (layer>1) continue;
623
624 // find the chip for the current cluster
625 Float_t locz = cluster->GetDetLocalZ();
626 Int_t iChip = seg->GetChipFromLocal(0,locz);
627 nClustersInChip[iChip]++;
628
629 }// end of cluster loop
630
631 // get number of fired chips in the current module
632 if(layer<2)
633 for(Int_t ifChip=0; ifChip<5; ifChip++) {
634 if(nClustersInChip[ifChip] >= 1) fNFiredChips[layer]++;
635 }
636
637 } // end of its "subdetector" loop
638
639 if (itsClusters) {
640 itsClusters->Delete();
641 delete itsClusters;
642 delete seg;
643 itsClusters = 0;
644 }
645 AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
ac903f1b 646}
647//____________________________________________________________________
648void
649AliITSMultReconstructor::SaveHists() {
3ef75756 650 // This method save the histograms on the output file
651 // (only if fHistOn is TRUE).
ac903f1b 652
653 if (!fHistOn)
654 return;
655
ddced3c8 656 fhClustersDPhiAll->Write();
657 fhClustersDThetaAll->Write();
658 fhClustersDZetaAll->Write();
ac903f1b 659 fhDPhiVsDThetaAll->Write();
ddced3c8 660 fhDPhiVsDZetaAll->Write();
661
662 fhClustersDPhiAcc->Write();
663 fhClustersDThetaAcc->Write();
664 fhClustersDZetaAcc->Write();
ac903f1b 665 fhDPhiVsDThetaAcc->Write();
ddced3c8 666 fhDPhiVsDZetaAcc->Write();
667
668 fhetaTracklets->Write();
669 fhphiTracklets->Write();
670 fhetaClustersLay1->Write();
671 fhphiClustersLay1->Write();
ac903f1b 672}
7b116aa1 673
674//____________________________________________________________________
675void
676AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
677
678 Float_t distClSameMod=0.;
679 Float_t distClSameModMin=0.;
680 Int_t iClOverlap =0;
681 Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
682 Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
683
684 Float_t zproj1=0.;
685 Float_t zproj2=0.;
686 Float_t deZproj=0.;
687
688 // Loop on inner layer clusters
689 for (Int_t iiC1=0; iiC1<fNClustersLay1; iiC1++) {
690 if (!fOverlapFlagClustersLay1[iiC1]) {
691 // only for adjacent modules
692 if ((TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==4)||
693 (TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==76)) {
694 Float_t dePhi=TMath::Abs(fClustersLay1[iiC1][1]-fClustersLay1[iC1][1]);
695 if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
696
697 zproj1=meanRadiusLay1/TMath::Tan(fClustersLay1[iC1][0]);
698 zproj2=meanRadiusLay1/TMath::Tan(fClustersLay1[iiC1][0]);
699
700 deZproj=TMath::Abs(zproj1-zproj2);
701
702 distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
703 if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
704
705// if (distClSameMod<=1.) {
706// if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
707// distClSameModMin=distClSameMod;
708// iClOverlap=iiC1;
709// }
710// }
711
712
713 } // end adjacent modules
714 }
715 } // end Loop on inner layer clusters
716
717// if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
718
719 distClSameMod=0.;
720 distClSameModMin=0.;
721 iClOverlap =0;
722 // Loop on outer layer clusters
723 for (Int_t iiC2=0; iiC2<fNClustersLay2; iiC2++) {
724 if (!fOverlapFlagClustersLay2[iiC2]) {
725 // only for adjacent modules
726 if ((TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==4) ||
727 (TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==156)) {
728 Float_t dePhi=TMath::Abs(fClustersLay2[iiC2][1]-fClustersLay2[iC2WithBestDist][1]);
729 if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
730
731 zproj1=meanRadiusLay2/TMath::Tan(fClustersLay2[iC2WithBestDist][0]);
732 zproj2=meanRadiusLay2/TMath::Tan(fClustersLay2[iiC2][0]);
733
734 deZproj=TMath::Abs(zproj1-zproj2);
735 distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
736 if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE;
737
738// if (distClSameMod<=1.) {
739// if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
740// distClSameModMin=distClSameMod;
741// iClOverlap=iiC2;
742// }
743// }
744
745 } // end adjacent modules
746 }
747 } // end Loop on outer layer clusters
748
749// if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;
750
751}