]> 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
1. SetDefaultOutputs(kTRUE) used by default. This will automatically archive the...
[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
7284b2b2 18//_________________________________________________________________________
ac903f1b 19//
7284b2b2 20// Implementation of the ITS-SPD trackleter class
ac903f1b 21//
fa9ed8e9 22// It retrieves clusters in the pixels (theta and phi) and finds tracklets.
23// These can be used to extract charged particle multiplicity from the ITS.
ac903f1b 24//
fa9ed8e9 25// A tracklet consists of two ITS clusters, one in the first pixel layer and
26// one in the second. The clusters are associated if the differences in
27// Phi (azimuth) and Theta (polar angle) are within fiducial windows.
28// In case of multiple candidates the candidate with minimum
29// distance is selected.
968e8539 30//
fa9ed8e9 31// Two methods return the number of tracklets and the number of unassociated
7284b2b2 32// clusters (i.e. not used in any tracklet) in the first SPD layer
33// (GetNTracklets and GetNSingleClusters)
34//
35// The cuts on phi and theta depend on the interacting system (p-p or Pb-Pb)
36// and can be set via AliITSRecoParam class
37// (SetPhiWindow and SetThetaWindow)
ac903f1b 38//
7284b2b2 39// Origin: Tiziano Virgili
40//
41// Current support and development:
42// Domenico Elia, Maria Nicassio (INFN Bari)
43// Domenico.Elia@ba.infn.it, Maria.Nicassio@ba.infn.it
44//
45// Most recent updates:
46// - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE)
f606f16a 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
fa9ed8e9 51// - option to cut duplicates in the overlaps
7b116aa1 52// - options and fiducial cuts via AliITSRecoParam
fa9ed8e9 53// - move from DeltaZeta to DeltaTheta cut
54// - update to the new algorithm by Mariella and Jan Fiete
55// - store also DeltaTheta in the ESD
56// - less new and delete calls when creating the needed arrays
7284b2b2 57//_________________________________________________________________________
ac903f1b 58
7ca4655f 59#include <TClonesArray.h>
60#include <TH1F.h>
61#include <TH2F.h>
62#include <TTree.h>
7284b2b2 63#include "TArrayI.h"
ac903f1b 64
7ca4655f 65#include "AliITSMultReconstructor.h"
7b116aa1 66#include "AliITSReconstructor.h"
9b373e9a 67#include "AliITSsegmentationSPD.h"
b51872de 68#include "AliITSRecPoint.h"
b21c1af0 69#include "AliITSRecPointContainer.h"
ac903f1b 70#include "AliITSgeom.h"
b21c1af0 71#include "AliITSgeomTGeo.h"
ac903f1b 72#include "AliLog.h"
fa9ed8e9 73#include "TGeoGlobalMagField.h"
74#include "AliMagF.h"
ac903f1b 75
76//____________________________________________________________________
0762f3a8 77ClassImp(AliITSMultReconstructor)
ac903f1b 78
3ef75756 79
ac903f1b 80//____________________________________________________________________
7537d03c 81AliITSMultReconstructor::AliITSMultReconstructor():
58e8dc31 82TObject(),
7537d03c 83fClustersLay1(0),
84fClustersLay2(0),
7b116aa1 85fDetectorIndexClustersLay1(0),
86fDetectorIndexClustersLay2(0),
87fOverlapFlagClustersLay1(0),
88fOverlapFlagClustersLay2(0),
7537d03c 89fTracklets(0),
968e8539 90fSClusters(0),
7537d03c 91fNClustersLay1(0),
92fNClustersLay2(0),
93fNTracklets(0),
968e8539 94fNSingleCluster(0),
7537d03c 95fPhiWindow(0),
7284b2b2 96fThetaWindow(0),
fa9ed8e9 97fPhiShift(0),
7b116aa1 98fRemoveClustersFromOverlaps(0),
99fPhiOverlapCut(0),
100fZetaOverlapCut(0),
7537d03c 101fHistOn(0),
102fhClustersDPhiAcc(0),
103fhClustersDThetaAcc(0),
7537d03c 104fhClustersDPhiAll(0),
105fhClustersDThetaAll(0),
7537d03c 106fhDPhiVsDThetaAll(0),
107fhDPhiVsDThetaAcc(0),
7537d03c 108fhetaTracklets(0),
109fhphiTracklets(0),
110fhetaClustersLay1(0),
111fhphiClustersLay1(0){
9b373e9a 112
113 fNFiredChips[0] = 0;
114 fNFiredChips[1] = 0;
115
3ef75756 116 // Method to reconstruct the charged particles multiplicity with the
117 // SPD (tracklets).
ac903f1b 118
ac903f1b 119
120 SetHistOn();
ac903f1b 121
7b116aa1 122 if(AliITSReconstructor::GetRecoParam()) {
7b116aa1 123 SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow());
7284b2b2 124 SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow());
fa9ed8e9 125 SetPhiShift(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiShift());
7b116aa1 126 SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps());
127 SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut());
128 SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut());
129 } else {
7b116aa1 130 SetPhiWindow();
7284b2b2 131 SetThetaWindow();
fa9ed8e9 132 SetPhiShift();
7b116aa1 133 SetRemoveClustersFromOverlaps();
134 SetPhiOverlapCut();
135 SetZetaOverlapCut();
136 }
137
138
fa9ed8e9 139 fClustersLay1 = 0;
140 fClustersLay2 = 0;
141 fDetectorIndexClustersLay1 = 0;
142 fDetectorIndexClustersLay2 = 0;
143 fOverlapFlagClustersLay1 = 0;
144 fOverlapFlagClustersLay2 = 0;
145 fTracklets = 0;
146 fSClusters = 0;
ac903f1b 147
148 // definition of histograms
fa9ed8e9 149 Bool_t oldStatus = TH1::AddDirectoryStatus();
150 TH1::AddDirectory(kFALSE);
151
7284b2b2 152 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,-0.1,0.1);
ddced3c8 153 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
ddced3c8 154
7284b2b2 155 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
ac903f1b 156
02a95988 157 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
7284b2b2 158 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,0.0,0.5);
ddced3c8 159
7284b2b2 160 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,0.,0.5,100,0.,0.5);
ddced3c8 161
162 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
f606f16a 163 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
ddced3c8 164 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
f606f16a 165 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
fa9ed8e9 166
167 TH1::AddDirectory(oldStatus);
ac903f1b 168}
ddced3c8 169
3ef75756 170//______________________________________________________________________
7537d03c 171AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
7537d03c 172fClustersLay1(mr.fClustersLay1),
173fClustersLay2(mr.fClustersLay2),
7b116aa1 174fDetectorIndexClustersLay1(mr.fDetectorIndexClustersLay1),
175fDetectorIndexClustersLay2(mr.fDetectorIndexClustersLay2),
176fOverlapFlagClustersLay1(mr.fOverlapFlagClustersLay1),
177fOverlapFlagClustersLay2(mr.fOverlapFlagClustersLay2),
7537d03c 178fTracklets(mr.fTracklets),
968e8539 179fSClusters(mr.fSClusters),
7537d03c 180fNClustersLay1(mr.fNClustersLay1),
181fNClustersLay2(mr.fNClustersLay2),
182fNTracklets(mr.fNTracklets),
968e8539 183fNSingleCluster(mr.fNSingleCluster),
7537d03c 184fPhiWindow(mr.fPhiWindow),
7284b2b2 185fThetaWindow(mr.fThetaWindow),
fa9ed8e9 186fPhiShift(mr.fPhiShift),
7b116aa1 187fRemoveClustersFromOverlaps(mr.fRemoveClustersFromOverlaps),
188fPhiOverlapCut(mr.fPhiOverlapCut),
189fZetaOverlapCut(mr.fZetaOverlapCut),
7537d03c 190fHistOn(mr.fHistOn),
191fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
192fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
7537d03c 193fhClustersDPhiAll(mr.fhClustersDPhiAll),
194fhClustersDThetaAll(mr.fhClustersDThetaAll),
7537d03c 195fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
196fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
7537d03c 197fhetaTracklets(mr.fhetaTracklets),
198fhphiTracklets(mr.fhphiTracklets),
199fhetaClustersLay1(mr.fhetaClustersLay1),
200fhphiClustersLay1(mr.fhphiClustersLay1) {
3ef75756 201 // Copy constructor
7537d03c 202
3ef75756 203}
204
205//______________________________________________________________________
7537d03c 206AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
3ef75756 207 // Assignment operator
7537d03c 208 this->~AliITSMultReconstructor();
209 new(this) AliITSMultReconstructor(mr);
3ef75756 210 return *this;
211}
212
213//______________________________________________________________________
214AliITSMultReconstructor::~AliITSMultReconstructor(){
215 // Destructor
1ba5b31c 216
217 // delete histograms
218 delete fhClustersDPhiAcc;
219 delete fhClustersDThetaAcc;
1ba5b31c 220 delete fhClustersDPhiAll;
221 delete fhClustersDThetaAll;
1ba5b31c 222 delete fhDPhiVsDThetaAll;
223 delete fhDPhiVsDThetaAcc;
1ba5b31c 224 delete fhetaTracklets;
225 delete fhphiTracklets;
226 delete fhetaClustersLay1;
227 delete fhphiClustersLay1;
228
229 // delete arrays
fa9ed8e9 230 for(Int_t i=0; i<fNClustersLay1; i++)
1ba5b31c 231 delete [] fClustersLay1[i];
fa9ed8e9 232
233 for(Int_t i=0; i<fNClustersLay2; i++)
1ba5b31c 234 delete [] fClustersLay2[i];
fa9ed8e9 235
236 for(Int_t i=0; i<fNTracklets; i++)
1ba5b31c 237 delete [] fTracklets[i];
fa9ed8e9 238
239 for(Int_t i=0; i<fNSingleCluster; i++)
968e8539 240 delete [] fSClusters[i];
fa9ed8e9 241
1ba5b31c 242 delete [] fClustersLay1;
243 delete [] fClustersLay2;
7b116aa1 244 delete [] fDetectorIndexClustersLay1;
245 delete [] fDetectorIndexClustersLay2;
246 delete [] fOverlapFlagClustersLay1;
247 delete [] fOverlapFlagClustersLay2;
1ba5b31c 248 delete [] fTracklets;
968e8539 249 delete [] fSClusters;
ddced3c8 250}
ac903f1b 251
252//____________________________________________________________________
7284b2b2 253void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
ac903f1b 254 //
255 // - calls LoadClusterArray that finds the position of the clusters
256 // (in global coord)
257 // - convert the cluster coordinates to theta, phi (seen from the
7284b2b2 258 // interaction vertex).
ac903f1b 259 // - makes an array of tracklets
260 //
261 // After this method has been called, the clusters of the two layers
262 // and the tracklets can be retrieved by calling the Get'er methods.
263
ac903f1b 264 // reset counters
265 fNClustersLay1 = 0;
266 fNClustersLay2 = 0;
267 fNTracklets = 0;
7284b2b2 268 fNSingleCluster = 0;
269
ac903f1b 270 // loading the clusters
271 LoadClusterArrays(clusterTree);
3ef75756 272
7284b2b2 273 const Double_t pi = TMath::Pi();
fa9ed8e9 274
275 // dPhi shift is field dependent
276 // get average magnetic field
277 Float_t bz = 0;
278 AliMagF* field = 0;
279 if (TGeoGlobalMagField::Instance())
280 field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
281 if (!field)
282 {
283 AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor.")
284 }
285 else
286 bz = TMath::Abs(field->SolenoidField());
287
288 const Double_t dPhiShift = fPhiShift / 5 * bz;
289 AliDebug(1, Form("Using phi shift of %f", dPhiShift));
290
291 const Double_t dPhiWindow2 = fPhiWindow * fPhiWindow;
292 const Double_t dThetaWindow2 = fThetaWindow * fThetaWindow;
293
7284b2b2 294 Int_t* partners = new Int_t[fNClustersLay2];
295 Float_t* minDists = new Float_t[fNClustersLay2];
296 Int_t* associatedLay1 = new Int_t[fNClustersLay1];
297 TArrayI** blacklist = new TArrayI*[fNClustersLay1];
298
299 for (Int_t i=0; i<fNClustersLay2; i++) {
300 partners[i] = -1;
301 minDists[i] = 2;
302 }
303 for (Int_t i=0; i<fNClustersLay1; i++)
304 associatedLay1[i] = 0;
305 for (Int_t i=0; i<fNClustersLay1; i++)
306 blacklist[i] = 0;
307
ac903f1b 308 // find the tracklets
309 AliDebug(1,"Looking for tracklets... ");
fa9ed8e9 310
ac903f1b 311 //###########################################################
312 // Loop on layer 1 : finding theta, phi and z
313 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
314 Float_t x = fClustersLay1[iC1][0] - vtx[0];
315 Float_t y = fClustersLay1[iC1][1] - vtx[1];
316 Float_t z = fClustersLay1[iC1][2] - vtx[2];
ddced3c8 317
fa9ed8e9 318 Float_t r = TMath::Sqrt(x*x + y*y + z*z);
ac903f1b 319
eefb3acc 320 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
321 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
fa9ed8e9 322
ddced3c8 323 if (fHistOn) {
324 Float_t eta=fClustersLay1[iC1][0];
325 eta= TMath::Tan(eta/2.);
326 eta=-TMath::Log(eta);
327 fhetaClustersLay1->Fill(eta);
de4c520e 328 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
ddced3c8 329 }
96c2c35d 330 }
ac903f1b 331
332 // Loop on layer 2 : finding theta, phi and r
333 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
334 Float_t x = fClustersLay2[iC2][0] - vtx[0];
335 Float_t y = fClustersLay2[iC2][1] - vtx[1];
336 Float_t z = fClustersLay2[iC2][2] - vtx[2];
ddced3c8 337
fa9ed8e9 338 Float_t r = TMath::Sqrt(x*x + y*y + z*z);
ac903f1b 339
eefb3acc 340 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
341 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
ac903f1b 342 }
343
344 //###########################################################
7284b2b2 345 Int_t found = 1;
346 while (found > 0) {
7284b2b2 347 found = 0;
348
349 // Step1: find all tracklets allowing double assocation
350 // Loop on layer 1
351 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
352
353 // already used or in the overlap ?
354 if (associatedLay1[iC1] != 0 || fOverlapFlagClustersLay1[iC1]) continue;
ac903f1b 355
7284b2b2 356 found++;
357
358 // reset of variables for multiple candidates
359 Int_t iC2WithBestDist = -1; // reset
360 Double_t minDist = 2; // reset
7b116aa1 361
7284b2b2 362 // Loop on layer 2
363 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
364
365 // in the overlap ?
366 if (fOverlapFlagClustersLay2[iC2]) continue;
367
368 if (blacklist[iC1]) {
369 Bool_t blacklisted = kFALSE;
370 for (Int_t i=0; i<blacklist[iC1]->GetSize(); i++) {
371 if (blacklist[iC1]->At(i) == iC2) {
372 blacklisted = kTRUE;
373 break;
374 }
375 }
376 if (blacklisted) continue;
377 }
378
ac903f1b 379 // find the difference in angles
7284b2b2 380 Double_t dTheta = TMath::Abs(fClustersLay2[iC2][0] - fClustersLay1[iC1][0]);
381 Double_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
02a95988 382 // take into account boundary condition
7284b2b2 383 if (dPhi>pi) dPhi=2.*pi-dPhi;
fa9ed8e9 384
ddced3c8 385 if (fHistOn) {
7284b2b2 386 fhClustersDPhiAll->Fill(dPhi);
ddced3c8 387 fhClustersDThetaAll->Fill(dTheta);
ac903f1b 388 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
389 }
fa9ed8e9 390
391 dPhi -= dPhiShift;
392
7284b2b2 393 // make "elliptical" cut in Phi and Theta!
fa9ed8e9 394 Float_t d = dPhi*dPhi/dPhiWindow2 + dTheta*dTheta/dThetaWindow2;
7284b2b2 395
396 // look for the minimum distance: the minimum is in iC2WithBestDist
fa9ed8e9 397 if (d<1 && d<minDist) {
7284b2b2 398 minDist=d;
ddced3c8 399 iC2WithBestDist = iC2;
ac903f1b 400 }
7284b2b2 401 } // end of loop over clusters in layer 2
ac903f1b 402
7284b2b2 403 if (minDist<1) { // This means that a cluster in layer 2 was found that matches with iC1
404
405 if (minDists[iC2WithBestDist] > minDist) {
406 Int_t oldPartner = partners[iC2WithBestDist];
407 partners[iC2WithBestDist] = iC1;
408 minDists[iC2WithBestDist] = minDist;
409
410 // mark as assigned
411 associatedLay1[iC1] = 1;
412
413 if (oldPartner != -1) {
fa9ed8e9 414 // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its blacklist
7284b2b2 415 if (blacklist[oldPartner] == 0) {
416 blacklist[oldPartner] = new TArrayI(1);
417 } else blacklist[oldPartner]->Set(blacklist[oldPartner]->GetSize()+1);
418
419 blacklist[oldPartner]->AddAt(iC2WithBestDist, blacklist[oldPartner]->GetSize()-1);
420
421 // mark as free
fa9ed8e9 422 associatedLay1[oldPartner] = 0;
7284b2b2 423 }
424 } else {
425 // try again to find a cluster without considering iC2WithBestDist
426 if (blacklist[iC1] == 0) {
427 blacklist[iC1] = new TArrayI(1);
428 }
fa9ed8e9 429 else
430 blacklist[iC1]->Set(blacklist[iC1]->GetSize()+1);
7284b2b2 431
432 blacklist[iC1]->AddAt(iC2WithBestDist, blacklist[iC1]->GetSize()-1);
433 }
de4c520e 434
7284b2b2 435 } else // cluster has no partner; remove
436 associatedLay1[iC1] = 2;
437 } // end of loop over clusters in layer 1
438 }
439
440 // Step2: store tracklets; remove used clusters
441 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
de4c520e 442
7284b2b2 443 if (partners[iC2] == -1) continue;
7b116aa1 444
7284b2b2 445 if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (partners[iC2],iC2);
446
aba7fa71 447
448 if (fOverlapFlagClustersLay1[partners[iC2]] || fOverlapFlagClustersLay2[iC2]) continue;
449
fa9ed8e9 450 fTracklets[fNTracklets] = new Float_t[6];
451
7284b2b2 452 // use the theta from the clusters in the first layer
453 fTracklets[fNTracklets][0] = fClustersLay1[partners[iC2]][0];
454 // use the phi from the clusters in the first layer
455 fTracklets[fNTracklets][1] = fClustersLay1[partners[iC2]][1];
456 // store the difference between phi1 and phi2
457 fTracklets[fNTracklets][2] = fClustersLay1[partners[iC2]][1] - fClustersLay2[iC2][1];
458
459 // define dphi in the range [0,pi] with proper sign (track charge correlated)
460 if (fTracklets[fNTracklets][2] > TMath::Pi())
461 fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]-2.*TMath::Pi();
462 if (fTracklets[fNTracklets][2] < -TMath::Pi())
463 fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]+2.*TMath::Pi();
7b116aa1 464
7284b2b2 465 // store the difference between theta1 and theta2
466 fTracklets[fNTracklets][3] = fClustersLay1[partners[iC2]][0] - fClustersLay2[iC2][0];
467
468 if (fHistOn) {
469 fhClustersDPhiAcc->Fill(fTracklets[fNTracklets][2]);
470 fhClustersDThetaAcc->Fill(fTracklets[fNTracklets][3]);
471 fhDPhiVsDThetaAcc->Fill(fTracklets[fNTracklets][3],fTracklets[fNTracklets][2]);
ac903f1b 472 }
3ef75756 473
7284b2b2 474 // find label
475 // if equal label in both clusters found this label is assigned
476 // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
477 Int_t label1 = 0;
478 Int_t label2 = 0;
479 while (label2 < 3) {
480 if ((Int_t) fClustersLay1[partners[iC2]][3+label1] != -2 && (Int_t) fClustersLay1[partners[iC2]][3+label1] == (Int_t) fClustersLay2[iC2][3+label2])
481 break;
482 label1++;
483 if (label1 == 3) {
484 label1 = 0;
485 label2++;
486 }
487 }
488 if (label2 < 3) {
489 AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) fClustersLay1[partners[iC2]][3+label1], (Int_t) fClustersLay2[iC2][3+label2], fNTracklets));
490 fTracklets[fNTracklets][4] = fClustersLay1[partners[iC2]][3+label1];
491 fTracklets[fNTracklets][5] = fClustersLay2[iC2][3+label2];
492 } else {
493 AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", (Int_t) fClustersLay1[partners[iC2]][3], (Int_t) fClustersLay1[partners[iC2]][4], (Int_t) fClustersLay1[partners[iC2]][5], (Int_t) fClustersLay2[iC2][3], (Int_t) fClustersLay2[iC2][4], (Int_t) fClustersLay2[iC2][5], fNTracklets));
494 fTracklets[fNTracklets][4] = fClustersLay1[partners[iC2]][3];
495 fTracklets[fNTracklets][5] = fClustersLay2[iC2][3];
496 }
497
498 if (fHistOn) {
499 Float_t eta=fTracklets[fNTracklets][0];
500 eta= TMath::Tan(eta/2.);
501 eta=-TMath::Log(eta);
502 fhetaTracklets->Fill(eta);
503 fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
504 }
3ef75756 505
7284b2b2 506 AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
507 AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", partners[iC2], iC2));
508 fNTracklets++;
3ef75756 509
7284b2b2 510 associatedLay1[partners[iC2]] = 1;
511 }
512
513 // Delete the following else if you do not want to save Clusters!
514 // store the cluster
515 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
516 if (associatedLay1[iC1]==2||associatedLay1[iC1]==0) {
fa9ed8e9 517 fSClusters[fNSingleCluster] = new Float_t[2];
968e8539 518 fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
519 fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
de4c520e 520 AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
7284b2b2 521 fNSingleCluster, iC1));
968e8539 522 fNSingleCluster++;
3ef75756 523 }
7284b2b2 524 }
525
526 delete[] partners;
527 delete[] minDists;
528
529 for (Int_t i=0; i<fNClustersLay1; i++)
530 if (blacklist[i])
531 delete blacklist[i];
532 delete[] blacklist;
533
ac903f1b 534 AliDebug(1,Form("%d tracklets found", fNTracklets));
535}
536
537//____________________________________________________________________
538void
539AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
540 // This method
541 // - gets the clusters from the cluster tree
542 // - convert them into global coordinates
543 // - store them in the internal arrays
9b373e9a 544 // - count the number of cluster-fired chips
ac903f1b 545
9b373e9a 546 AliDebug(1,"Loading clusters and cluster-fired chips ...");
ac903f1b 547
548 fNClustersLay1 = 0;
549 fNClustersLay2 = 0;
9b373e9a 550 fNFiredChips[0] = 0;
551 fNFiredChips[1] = 0;
ac903f1b 552
b21c1af0 553 AliITSsegmentationSPD seg;
9b373e9a 554
b21c1af0 555 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
556 TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
557 if(!rpcont->IsSPDActive()){
558 AliWarning("No SPD rec points found, multiplicity not calculated");
559 return;
560 }
f606f16a 561 Float_t cluGlo[3]={0.,0.,0.};
ddced3c8 562
fa9ed8e9 563
564 // count clusters
b21c1af0 565 // loop over the SPD subdetectors
566 Int_t nSPDL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
567 for (Int_t iIts=0; iIts < nSPDL1; iIts++) {
568 itsClusters=rpcont->UncheckedGetClusters(iIts);
569 fNClustersLay1 += itsClusters->GetEntriesFast();
570 }
571 Int_t nSPDL2=AliITSgeomTGeo::GetModuleIndex(3,1,1);
572 for (Int_t iIts=nSPDL1; iIts < nSPDL2; iIts++) {
573 itsClusters=rpcont->UncheckedGetClusters(iIts);
574 fNClustersLay2 += itsClusters->GetEntriesFast();
fa9ed8e9 575 }
576
577 // create arrays
578 fClustersLay1 = new Float_t*[fNClustersLay1];
579 fDetectorIndexClustersLay1 = new Int_t[fNClustersLay1];
580 fOverlapFlagClustersLay1 = new Bool_t[fNClustersLay1];
581
582 fClustersLay2 = new Float_t*[fNClustersLay2];
583 fDetectorIndexClustersLay2 = new Int_t[fNClustersLay2];
584 fOverlapFlagClustersLay2 = new Bool_t[fNClustersLay2];
585
586 // no double association allowed
587 fTracklets = new Float_t*[TMath::Min(fNClustersLay1, fNClustersLay2)];
588 fSClusters = new Float_t*[fNClustersLay1];
589
590 for (Int_t i=0; i<fNClustersLay1; i++) {
591 fClustersLay1[i] = new Float_t[6];
592 fOverlapFlagClustersLay1[i] = kFALSE;
593 fSClusters[i] = 0;
594 }
595
596 for (Int_t i=0; i<fNClustersLay2; i++) {
597 fClustersLay2[i] = new Float_t[6];
598 fOverlapFlagClustersLay2[i] = kFALSE;
599 }
600
601 for (Int_t i=0; i<TMath::Min(fNClustersLay1, fNClustersLay2); i++)
602 fTracklets[i] = 0;
603
604 // fill clusters
ac903f1b 605 // loop over the its subdetectors
fa9ed8e9 606 fNClustersLay1 = 0; // reset to 0
607 fNClustersLay2 = 0;
b21c1af0 608 for (Int_t iIts=0; iIts < nSPDL2; iIts++) {
ac903f1b 609
b21c1af0 610 itsClusters=rpcont->UncheckedGetClusters(iIts);
ac903f1b 611
612 Int_t nClusters = itsClusters->GetEntriesFast();
9b373e9a 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;
ac903f1b 617
ac903f1b 618 // loop over clusters
619 while(nClusters--) {
de4c520e 620 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
ac903f1b 621
9b373e9a 622 layer = cluster->GetLayer();
623 if (layer>1) continue;
ac903f1b 624
f606f16a 625 cluster->GetGlobalXYZ(cluGlo);
626 Float_t x = cluGlo[0];
627 Float_t y = cluGlo[1];
628 Float_t z = cluGlo[2];
9b373e9a 629
630 // find the chip for the current cluster
631 Float_t locz = cluster->GetDetLocalZ();
b21c1af0 632 Int_t iChip = seg.GetChipFromLocal(0,locz);
9b373e9a 633 nClustersInChip[iChip]++;
ac903f1b 634
9b373e9a 635 if (layer==0) {
ac903f1b 636 fClustersLay1[fNClustersLay1][0] = x;
637 fClustersLay1[fNClustersLay1][1] = y;
638 fClustersLay1[fNClustersLay1][2] = z;
7b116aa1 639
640 fDetectorIndexClustersLay1[fNClustersLay1]=iIts;
641
de4c520e 642 for (Int_t i=0; i<3; i++)
643 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
ac903f1b 644 fNClustersLay1++;
645 }
9b373e9a 646 if (layer==1) {
ac903f1b 647 fClustersLay2[fNClustersLay2][0] = x;
648 fClustersLay2[fNClustersLay2][1] = y;
649 fClustersLay2[fNClustersLay2][2] = z;
7b116aa1 650
651 fDetectorIndexClustersLay2[fNClustersLay2]=iIts;
652
de4c520e 653 for (Int_t i=0; i<3; i++)
654 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
ac903f1b 655 fNClustersLay2++;
656 }
657
658 }// end of cluster loop
9b373e9a 659
660 // get number of fired chips in the current module
b21c1af0 661
9b373e9a 662 for(Int_t ifChip=0; ifChip<5; ifChip++) {
663 if(nClustersInChip[ifChip] >= 1) fNFiredChips[layer]++;
664 }
665
ac903f1b 666 } // end of its "subdetector" loop
9b373e9a 667
ac903f1b 668 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
9b373e9a 669 AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
670}
671//____________________________________________________________________
672void
673AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
674 // This method
675 // - gets the clusters from the cluster tree
676 // - counts the number of (cluster)fired chips
677
678 AliDebug(1,"Loading cluster-fired chips ...");
679
680 fNFiredChips[0] = 0;
681 fNFiredChips[1] = 0;
682
b21c1af0 683 AliITSsegmentationSPD seg;
684 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
685 TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
686 if(!rpcont->IsSPDActive()){
687 AliWarning("No SPD rec points found, multiplicity not calculated");
688 return;
689 }
9b373e9a 690
9b373e9a 691 // loop over the its subdetectors
b21c1af0 692 Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1);
693 for (Int_t iIts=0; iIts < nSPDmodules; iIts++) {
694 itsClusters=rpcont->UncheckedGetClusters(iIts);
9b373e9a 695 Int_t nClusters = itsClusters->GetEntriesFast();
696
697 // number of clusters in each chip of the current module
698 Int_t nClustersInChip[5] = {0,0,0,0,0};
699 Int_t layer = 0;
700
701 // loop over clusters
702 while(nClusters--) {
703 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
704
705 layer = cluster->GetLayer();
706 if (layer>1) continue;
707
708 // find the chip for the current cluster
709 Float_t locz = cluster->GetDetLocalZ();
b21c1af0 710 Int_t iChip = seg.GetChipFromLocal(0,locz);
9b373e9a 711 nClustersInChip[iChip]++;
712
713 }// end of cluster loop
714
715 // get number of fired chips in the current module
9b373e9a 716 for(Int_t ifChip=0; ifChip<5; ifChip++) {
717 if(nClustersInChip[ifChip] >= 1) fNFiredChips[layer]++;
718 }
719
720 } // end of its "subdetector" loop
721
b21c1af0 722
9b373e9a 723 AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
ac903f1b 724}
725//____________________________________________________________________
726void
727AliITSMultReconstructor::SaveHists() {
3ef75756 728 // This method save the histograms on the output file
729 // (only if fHistOn is TRUE).
ac903f1b 730
731 if (!fHistOn)
732 return;
733
ddced3c8 734 fhClustersDPhiAll->Write();
735 fhClustersDThetaAll->Write();
ac903f1b 736 fhDPhiVsDThetaAll->Write();
ddced3c8 737
738 fhClustersDPhiAcc->Write();
739 fhClustersDThetaAcc->Write();
ac903f1b 740 fhDPhiVsDThetaAcc->Write();
ddced3c8 741
742 fhetaTracklets->Write();
743 fhphiTracklets->Write();
744 fhetaClustersLay1->Write();
745 fhphiClustersLay1->Write();
ac903f1b 746}
7b116aa1 747
748//____________________________________________________________________
749void
750AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
751
752 Float_t distClSameMod=0.;
753 Float_t distClSameModMin=0.;
754 Int_t iClOverlap =0;
755 Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
756 Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
757
758 Float_t zproj1=0.;
759 Float_t zproj2=0.;
760 Float_t deZproj=0.;
761
762 // Loop on inner layer clusters
763 for (Int_t iiC1=0; iiC1<fNClustersLay1; iiC1++) {
764 if (!fOverlapFlagClustersLay1[iiC1]) {
765 // only for adjacent modules
766 if ((TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==4)||
767 (TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==76)) {
768 Float_t dePhi=TMath::Abs(fClustersLay1[iiC1][1]-fClustersLay1[iC1][1]);
769 if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
770
771 zproj1=meanRadiusLay1/TMath::Tan(fClustersLay1[iC1][0]);
772 zproj2=meanRadiusLay1/TMath::Tan(fClustersLay1[iiC1][0]);
773
774 deZproj=TMath::Abs(zproj1-zproj2);
775
776 distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
777 if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
778
779// if (distClSameMod<=1.) {
780// if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
781// distClSameModMin=distClSameMod;
782// iClOverlap=iiC1;
783// }
784// }
785
786
787 } // end adjacent modules
788 }
789 } // end Loop on inner layer clusters
790
791// if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
792
793 distClSameMod=0.;
794 distClSameModMin=0.;
795 iClOverlap =0;
796 // Loop on outer layer clusters
797 for (Int_t iiC2=0; iiC2<fNClustersLay2; iiC2++) {
798 if (!fOverlapFlagClustersLay2[iiC2]) {
799 // only for adjacent modules
800 if ((TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==4) ||
801 (TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==156)) {
802 Float_t dePhi=TMath::Abs(fClustersLay2[iiC2][1]-fClustersLay2[iC2WithBestDist][1]);
803 if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
804
805 zproj1=meanRadiusLay2/TMath::Tan(fClustersLay2[iC2WithBestDist][0]);
806 zproj2=meanRadiusLay2/TMath::Tan(fClustersLay2[iiC2][0]);
807
808 deZproj=TMath::Abs(zproj1-zproj2);
809 distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
810 if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE;
811
812// if (distClSameMod<=1.) {
813// if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
814// distClSameModMin=distClSameMod;
815// iClOverlap=iiC2;
816// }
817// }
818
819 } // end adjacent modules
820 }
821 } // end Loop on outer layer clusters
822
823// if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;
824
6b489238 825}