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 // Implementation of the ITS-SPD trackleter class
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.
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.
31 // Two methods return the number of tracklets and the number of unassociated
32 // clusters (i.e. not used in any tracklet) in the first SPD layer
33 // (GetNTracklets and GetNSingleClusters)
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)
39 // Origin: Tiziano Virgili
41 // Current support and development:
42 // Domenico Elia, Maria Nicassio (INFN Bari)
43 // Domenico.Elia@ba.infn.it, Maria.Nicassio@ba.infn.it
45 // Most recent updates:
46 // - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE)
47 // - phi definition changed to ALICE convention (0,2*TMath::pi())
48 // - cluster coordinates taken with GetGlobalXYZ()
49 // - fGeometry removed
50 // - number of fired chips on the two layers
51 // - option to cut duplicates in the overlaps
52 // - options and fiducial cuts via AliITSRecoParam
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
57 //_________________________________________________________________________
59 #include <TClonesArray.h>
65 #include "AliITSMultReconstructor.h"
66 #include "AliITSReconstructor.h"
67 #include "AliITSsegmentationSPD.h"
68 #include "AliITSRecPoint.h"
69 #include "AliITSRecPointContainer.h"
70 #include "AliITSgeom.h"
71 #include "AliITSgeomTGeo.h"
73 #include "TGeoGlobalMagField.h"
76 //____________________________________________________________________
77 ClassImp(AliITSMultReconstructor)
80 //____________________________________________________________________
81 AliITSMultReconstructor::AliITSMultReconstructor():
85 fDetectorIndexClustersLay1(0),
86 fDetectorIndexClustersLay2(0),
87 fOverlapFlagClustersLay1(0),
88 fOverlapFlagClustersLay2(0),
98 fRemoveClustersFromOverlaps(0),
102 fhClustersDPhiAcc(0),
103 fhClustersDThetaAcc(0),
104 fhClustersDPhiAll(0),
105 fhClustersDThetaAll(0),
106 fhDPhiVsDThetaAll(0),
107 fhDPhiVsDThetaAcc(0),
110 fhetaClustersLay1(0),
111 fhphiClustersLay1(0){
116 // Method to reconstruct the charged particles multiplicity with the
122 if(AliITSReconstructor::GetRecoParam()) {
123 SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow());
124 SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow());
125 SetPhiShift(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiShift());
126 SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps());
127 SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut());
128 SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut());
133 SetRemoveClustersFromOverlaps();
141 fDetectorIndexClustersLay1 = 0;
142 fDetectorIndexClustersLay2 = 0;
143 fOverlapFlagClustersLay1 = 0;
144 fOverlapFlagClustersLay2 = 0;
148 // definition of histograms
149 Bool_t oldStatus = TH1::AddDirectoryStatus();
150 TH1::AddDirectory(kFALSE);
152 fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,-0.1,0.1);
153 fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
155 fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
157 fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
158 fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,0.0,0.5);
160 fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,0.,0.5,100,0.,0.5);
162 fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
163 fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
164 fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
165 fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
167 TH1::AddDirectory(oldStatus);
170 //______________________________________________________________________
171 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
172 fClustersLay1(mr.fClustersLay1),
173 fClustersLay2(mr.fClustersLay2),
174 fDetectorIndexClustersLay1(mr.fDetectorIndexClustersLay1),
175 fDetectorIndexClustersLay2(mr.fDetectorIndexClustersLay2),
176 fOverlapFlagClustersLay1(mr.fOverlapFlagClustersLay1),
177 fOverlapFlagClustersLay2(mr.fOverlapFlagClustersLay2),
178 fTracklets(mr.fTracklets),
179 fSClusters(mr.fSClusters),
180 fNClustersLay1(mr.fNClustersLay1),
181 fNClustersLay2(mr.fNClustersLay2),
182 fNTracklets(mr.fNTracklets),
183 fNSingleCluster(mr.fNSingleCluster),
184 fPhiWindow(mr.fPhiWindow),
185 fThetaWindow(mr.fThetaWindow),
186 fPhiShift(mr.fPhiShift),
187 fRemoveClustersFromOverlaps(mr.fRemoveClustersFromOverlaps),
188 fPhiOverlapCut(mr.fPhiOverlapCut),
189 fZetaOverlapCut(mr.fZetaOverlapCut),
191 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
192 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
193 fhClustersDPhiAll(mr.fhClustersDPhiAll),
194 fhClustersDThetaAll(mr.fhClustersDThetaAll),
195 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
196 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
197 fhetaTracklets(mr.fhetaTracklets),
198 fhphiTracklets(mr.fhphiTracklets),
199 fhetaClustersLay1(mr.fhetaClustersLay1),
200 fhphiClustersLay1(mr.fhphiClustersLay1) {
205 //______________________________________________________________________
206 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
207 // Assignment operator
208 this->~AliITSMultReconstructor();
209 new(this) AliITSMultReconstructor(mr);
213 //______________________________________________________________________
214 AliITSMultReconstructor::~AliITSMultReconstructor(){
218 delete fhClustersDPhiAcc;
219 delete fhClustersDThetaAcc;
220 delete fhClustersDPhiAll;
221 delete fhClustersDThetaAll;
222 delete fhDPhiVsDThetaAll;
223 delete fhDPhiVsDThetaAcc;
224 delete fhetaTracklets;
225 delete fhphiTracklets;
226 delete fhetaClustersLay1;
227 delete fhphiClustersLay1;
230 for(Int_t i=0; i<fNClustersLay1; i++)
231 delete [] fClustersLay1[i];
233 for(Int_t i=0; i<fNClustersLay2; i++)
234 delete [] fClustersLay2[i];
236 for(Int_t i=0; i<fNTracklets; i++)
237 delete [] fTracklets[i];
239 for(Int_t i=0; i<fNSingleCluster; i++)
240 delete [] fSClusters[i];
242 delete [] fClustersLay1;
243 delete [] fClustersLay2;
244 delete [] fDetectorIndexClustersLay1;
245 delete [] fDetectorIndexClustersLay2;
246 delete [] fOverlapFlagClustersLay1;
247 delete [] fOverlapFlagClustersLay2;
248 delete [] fTracklets;
249 delete [] fSClusters;
252 //____________________________________________________________________
253 void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
255 // - calls LoadClusterArray that finds the position of the clusters
257 // - convert the cluster coordinates to theta, phi (seen from the
258 // interaction vertex).
259 // - makes an array of tracklets
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.
270 // loading the clusters
271 LoadClusterArrays(clusterTree);
273 const Double_t pi = TMath::Pi();
275 // dPhi shift is field dependent
276 // get average magnetic field
279 if (TGeoGlobalMagField::Instance())
280 field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
283 AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor.")
286 bz = TMath::Abs(field->SolenoidField());
288 const Double_t dPhiShift = fPhiShift / 5 * bz;
289 AliDebug(1, Form("Using phi shift of %f", dPhiShift));
291 const Double_t dPhiWindow2 = fPhiWindow * fPhiWindow;
292 const Double_t dThetaWindow2 = fThetaWindow * fThetaWindow;
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];
299 for (Int_t i=0; i<fNClustersLay2; i++) {
303 for (Int_t i=0; i<fNClustersLay1; i++)
304 associatedLay1[i] = 0;
305 for (Int_t i=0; i<fNClustersLay1; i++)
308 // find the tracklets
309 AliDebug(1,"Looking for tracklets... ");
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];
318 Float_t r = TMath::Sqrt(x*x + y*y + z*z);
320 fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
321 fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
324 Float_t eta=fClustersLay1[iC1][0];
325 eta= TMath::Tan(eta/2.);
326 eta=-TMath::Log(eta);
327 fhetaClustersLay1->Fill(eta);
328 fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
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];
338 Float_t r = TMath::Sqrt(x*x + y*y + z*z);
340 fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
341 fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
344 //###########################################################
349 // Step1: find all tracklets allowing double assocation
351 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
353 // already used or in the overlap ?
354 if (associatedLay1[iC1] != 0 || fOverlapFlagClustersLay1[iC1]) continue;
358 // reset of variables for multiple candidates
359 Int_t iC2WithBestDist = -1; // reset
360 Double_t minDist = 2; // reset
363 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
366 if (fOverlapFlagClustersLay2[iC2]) continue;
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) {
376 if (blacklisted) continue;
379 // find the difference in angles
380 Double_t dTheta = TMath::Abs(fClustersLay2[iC2][0] - fClustersLay1[iC1][0]);
381 Double_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
382 // take into account boundary condition
383 if (dPhi>pi) dPhi=2.*pi-dPhi;
386 fhClustersDPhiAll->Fill(dPhi);
387 fhClustersDThetaAll->Fill(dTheta);
388 fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
393 // make "elliptical" cut in Phi and Theta!
394 Float_t d = dPhi*dPhi/dPhiWindow2 + dTheta*dTheta/dThetaWindow2;
396 // look for the minimum distance: the minimum is in iC2WithBestDist
397 if (d<1 && d<minDist) {
399 iC2WithBestDist = iC2;
401 } // end of loop over clusters in layer 2
403 if (minDist<1) { // This means that a cluster in layer 2 was found that matches with iC1
405 if (minDists[iC2WithBestDist] > minDist) {
406 Int_t oldPartner = partners[iC2WithBestDist];
407 partners[iC2WithBestDist] = iC1;
408 minDists[iC2WithBestDist] = minDist;
411 associatedLay1[iC1] = 1;
413 if (oldPartner != -1) {
414 // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its blacklist
415 if (blacklist[oldPartner] == 0) {
416 blacklist[oldPartner] = new TArrayI(1);
417 } else blacklist[oldPartner]->Set(blacklist[oldPartner]->GetSize()+1);
419 blacklist[oldPartner]->AddAt(iC2WithBestDist, blacklist[oldPartner]->GetSize()-1);
422 associatedLay1[oldPartner] = 0;
425 // try again to find a cluster without considering iC2WithBestDist
426 if (blacklist[iC1] == 0) {
427 blacklist[iC1] = new TArrayI(1);
430 blacklist[iC1]->Set(blacklist[iC1]->GetSize()+1);
432 blacklist[iC1]->AddAt(iC2WithBestDist, blacklist[iC1]->GetSize()-1);
435 } else // cluster has no partner; remove
436 associatedLay1[iC1] = 2;
437 } // end of loop over clusters in layer 1
440 // Step2: store tracklets; remove used clusters
441 for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
443 if (partners[iC2] == -1) continue;
445 if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (partners[iC2],iC2);
448 if (fOverlapFlagClustersLay1[partners[iC2]] || fOverlapFlagClustersLay2[iC2]) continue;
450 fTracklets[fNTracklets] = new Float_t[6];
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];
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();
465 // store the difference between theta1 and theta2
466 fTracklets[fNTracklets][3] = fClustersLay1[partners[iC2]][0] - fClustersLay2[iC2][0];
469 fhClustersDPhiAcc->Fill(fTracklets[fNTracklets][2]);
470 fhClustersDThetaAcc->Fill(fTracklets[fNTracklets][3]);
471 fhDPhiVsDThetaAcc->Fill(fTracklets[fNTracklets][3],fTracklets[fNTracklets][2]);
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
480 if ((Int_t) fClustersLay1[partners[iC2]][3+label1] != -2 && (Int_t) fClustersLay1[partners[iC2]][3+label1] == (Int_t) fClustersLay2[iC2][3+label2])
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];
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];
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]);
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));
510 associatedLay1[partners[iC2]] = 1;
513 // Delete the following else if you do not want to save Clusters!
515 for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
516 if (associatedLay1[iC1]==2||associatedLay1[iC1]==0) {
517 fSClusters[fNSingleCluster] = new Float_t[2];
518 fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
519 fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
520 AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
521 fNSingleCluster, iC1));
529 for (Int_t i=0; i<fNClustersLay1; i++)
534 AliDebug(1,Form("%d tracklets found", fNTracklets));
537 //____________________________________________________________________
539 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
541 // - gets the clusters from the cluster tree
542 // - convert them into global coordinates
543 // - store them in the internal arrays
544 // - count the number of cluster-fired chips
546 AliDebug(1,"Loading clusters and cluster-fired chips ...");
553 AliITSsegmentationSPD seg;
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");
561 Float_t cluGlo[3]={0.,0.,0.};
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();
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();
578 fClustersLay1 = new Float_t*[fNClustersLay1];
579 fDetectorIndexClustersLay1 = new Int_t[fNClustersLay1];
580 fOverlapFlagClustersLay1 = new Bool_t[fNClustersLay1];
582 fClustersLay2 = new Float_t*[fNClustersLay2];
583 fDetectorIndexClustersLay2 = new Int_t[fNClustersLay2];
584 fOverlapFlagClustersLay2 = new Bool_t[fNClustersLay2];
586 // no double association allowed
587 fTracklets = new Float_t*[TMath::Min(fNClustersLay1, fNClustersLay2)];
588 fSClusters = new Float_t*[fNClustersLay1];
590 for (Int_t i=0; i<fNClustersLay1; i++) {
591 fClustersLay1[i] = new Float_t[6];
592 fOverlapFlagClustersLay1[i] = kFALSE;
596 for (Int_t i=0; i<fNClustersLay2; i++) {
597 fClustersLay2[i] = new Float_t[6];
598 fOverlapFlagClustersLay2[i] = kFALSE;
601 for (Int_t i=0; i<TMath::Min(fNClustersLay1, fNClustersLay2); i++)
605 // loop over the its subdetectors
606 fNClustersLay1 = 0; // reset to 0
608 for (Int_t iIts=0; iIts < nSPDL2; iIts++) {
610 itsClusters=rpcont->UncheckedGetClusters(iIts);
612 Int_t nClusters = itsClusters->GetEntriesFast();
614 // number of clusters in each chip of the current module
615 Int_t nClustersInChip[5] = {0,0,0,0,0};
618 // loop over clusters
620 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
622 layer = cluster->GetLayer();
623 if (layer>1) continue;
625 cluster->GetGlobalXYZ(cluGlo);
626 Float_t x = cluGlo[0];
627 Float_t y = cluGlo[1];
628 Float_t z = cluGlo[2];
630 // find the chip for the current cluster
631 Float_t locz = cluster->GetDetLocalZ();
632 Int_t iChip = seg.GetChipFromLocal(0,locz);
633 nClustersInChip[iChip]++;
636 fClustersLay1[fNClustersLay1][0] = x;
637 fClustersLay1[fNClustersLay1][1] = y;
638 fClustersLay1[fNClustersLay1][2] = z;
640 fDetectorIndexClustersLay1[fNClustersLay1]=iIts;
642 for (Int_t i=0; i<3; i++)
643 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
647 fClustersLay2[fNClustersLay2][0] = x;
648 fClustersLay2[fNClustersLay2][1] = y;
649 fClustersLay2[fNClustersLay2][2] = z;
651 fDetectorIndexClustersLay2[fNClustersLay2]=iIts;
653 for (Int_t i=0; i<3; i++)
654 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
658 }// end of cluster loop
660 // get number of fired chips in the current module
662 for(Int_t ifChip=0; ifChip<5; ifChip++) {
663 if(nClustersInChip[ifChip] >= 1) fNFiredChips[layer]++;
666 } // end of its "subdetector" loop
668 AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
669 AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
671 //____________________________________________________________________
673 AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
675 // - gets the clusters from the cluster tree
676 // - counts the number of (cluster)fired chips
678 AliDebug(1,"Loading cluster-fired chips ...");
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");
691 // loop over the its subdetectors
692 Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1);
693 for (Int_t iIts=0; iIts < nSPDmodules; iIts++) {
694 itsClusters=rpcont->UncheckedGetClusters(iIts);
695 Int_t nClusters = itsClusters->GetEntriesFast();
697 // number of clusters in each chip of the current module
698 Int_t nClustersInChip[5] = {0,0,0,0,0};
701 // loop over clusters
703 AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
705 layer = cluster->GetLayer();
706 if (layer>1) continue;
708 // find the chip for the current cluster
709 Float_t locz = cluster->GetDetLocalZ();
710 Int_t iChip = seg.GetChipFromLocal(0,locz);
711 nClustersInChip[iChip]++;
713 }// end of cluster loop
715 // get number of fired chips in the current module
716 for(Int_t ifChip=0; ifChip<5; ifChip++) {
717 if(nClustersInChip[ifChip] >= 1) fNFiredChips[layer]++;
720 } // end of its "subdetector" loop
723 AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
725 //____________________________________________________________________
727 AliITSMultReconstructor::SaveHists() {
728 // This method save the histograms on the output file
729 // (only if fHistOn is TRUE).
734 fhClustersDPhiAll->Write();
735 fhClustersDThetaAll->Write();
736 fhDPhiVsDThetaAll->Write();
738 fhClustersDPhiAcc->Write();
739 fhClustersDThetaAcc->Write();
740 fhDPhiVsDThetaAcc->Write();
742 fhetaTracklets->Write();
743 fhphiTracklets->Write();
744 fhetaClustersLay1->Write();
745 fhphiClustersLay1->Write();
748 //____________________________________________________________________
750 AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
752 Float_t distClSameMod=0.;
753 Float_t distClSameModMin=0.;
755 Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
756 Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
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;
771 zproj1=meanRadiusLay1/TMath::Tan(fClustersLay1[iC1][0]);
772 zproj2=meanRadiusLay1/TMath::Tan(fClustersLay1[iiC1][0]);
774 deZproj=TMath::Abs(zproj1-zproj2);
776 distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
777 if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
779 // if (distClSameMod<=1.) {
780 // if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
781 // distClSameModMin=distClSameMod;
787 } // end adjacent modules
789 } // end Loop on inner layer clusters
791 // if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
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;
805 zproj1=meanRadiusLay2/TMath::Tan(fClustersLay2[iC2WithBestDist][0]);
806 zproj2=meanRadiusLay2/TMath::Tan(fClustersLay2[iiC2][0]);
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;
812 // if (distClSameMod<=1.) {
813 // if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
814 // distClSameModMin=distClSameMod;
819 } // end adjacent modules
821 } // end Loop on outer layer clusters
823 // if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;