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