Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
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 : //_________________________________________________________________________
17 : //
18 : // Implementation of the ITS-SPD trackleter class
19 : //
20 : // It retrieves clusters in the pixels (theta and phi) and finds tracklets.
21 : // These can be used to extract charged particle multiplicity from the ITS.
22 : //
23 : // A tracklet consists of two ITS clusters, one in the first pixel layer and
24 : // one in the second. The clusters are associated if the differences in
25 : // Phi (azimuth) and Theta (polar angle) are within fiducial windows.
26 : // In case of multiple candidates the candidate with minimum
27 : // distance is selected.
28 : //
29 : // Two methods return the number of tracklets and the number of unassociated
30 : // clusters (i.e. not used in any tracklet) in the first SPD layer
31 : // (GetNTracklets and GetNSingleClusters)
32 : //
33 : // The cuts on phi and theta depend on the interacting system (p-p or Pb-Pb)
34 : // and can be set via AliITSRecoParam class
35 : // (SetPhiWindow and SetThetaWindow)
36 : //
37 : // Origin: Tiziano Virgili
38 : //
39 : // Current support and development:
40 : // Domenico Elia, Maria Nicassio (INFN Bari)
41 : // Domenico.Elia@ba.infn.it, Maria.Nicassio@ba.infn.it
42 : //
43 : // Most recent updates:
44 : // - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE)
45 : // - phi definition changed to ALICE convention (0,2*TMath::pi())
46 : // - cluster coordinates taken with GetGlobalXYZ()
47 : // - fGeometry removed
48 : // - number of fired chips on the two layers
49 : // - option to cut duplicates in the overlaps
50 : // - options and fiducial cuts via AliITSRecoParam
51 : // - move from DeltaZeta to DeltaTheta cut
52 : // - update to the new algorithm by Mariella and Jan Fiete
53 : // - store also DeltaTheta in the ESD
54 : // - less new and delete calls when creating the needed arrays
55 : //
56 : // - RS: to decrease the number of new/deletes the clusters data are stored
57 : // not in float[6] attached to float**, but in 1-D array.
58 : // - RS: Clusters are sorted in Z in roder to have the same numbering as in the ITS reco
59 : // - RS: Clusters used by ESDtrack are flagged, this information is passed to AliMulitiplicity object
60 : // when storing the tracklets and single cluster info
61 : // - MN: first MC label of single clusters stored
62 : //_________________________________________________________________________
63 :
64 : #include <TClonesArray.h>
65 : #include <TH1F.h>
66 : #include <TH2F.h>
67 : #include <TTree.h>
68 : #include <TBits.h>
69 : #include <TArrayI.h>
70 : #include <string.h>
71 :
72 : #include "AliITSMultReconstructor.h"
73 : #include "AliITSReconstructor.h"
74 : #include "AliITSRecPoint.h"
75 : #include "AliITSRecPointContainer.h"
76 : #include "AliITSgeom.h"
77 : #include "AliITSgeomTGeo.h"
78 : #include "AliITSDetTypeRec.h"
79 : #include "AliESDEvent.h"
80 : #include "AliESDVertex.h"
81 : #include "AliESDtrack.h"
82 : #include "AliMultiplicity.h"
83 : #include "AliLog.h"
84 : #include "TGeoGlobalMagField.h"
85 : #include "AliMagF.h"
86 : #include "AliESDv0.h"
87 : #include "AliV0.h"
88 : #include "AliKFParticle.h"
89 : #include "AliKFVertex.h"
90 : #include "AliRefArray.h"
91 :
92 : //____________________________________________________________________
93 118 : ClassImp(AliITSMultReconstructor)
94 :
95 :
96 : //____________________________________________________________________
97 8 : AliITSMultReconstructor::AliITSMultReconstructor():
98 8 : fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fTreeRPMix(0),
99 8 : fTracklets(0),
100 8 : fSClusters(0),
101 8 : fNTracklets(0),
102 8 : fNSingleCluster(0),
103 8 : fNSingleClusterSPD2(0),
104 8 : fDPhiWindow(0),
105 8 : fDThetaWindow(0),
106 8 : fPhiShift(0),
107 8 : fRemoveClustersFromOverlaps(0),
108 8 : fPhiOverlapCut(0),
109 8 : fZetaOverlapCut(0),
110 8 : fPhiRotationAngle(0),
111 8 : fScaleDTBySin2T(0),
112 8 : fNStdDev(1.0),
113 8 : fNStdDevSq(1.0),
114 : //
115 8 : fCutPxDrSPDin(0.1),
116 8 : fCutPxDrSPDout(0.15),
117 8 : fCutPxDz(0.2),
118 8 : fCutDCArz(0.5),
119 8 : fCutMinElectronProbTPC(0.5),
120 8 : fCutMinElectronProbESD(0.1),
121 8 : fCutMinP(0.05),
122 8 : fCutMinRGamma(2.),
123 8 : fCutMinRK0(1.),
124 8 : fCutMinPointAngle(0.98),
125 8 : fCutMaxDCADauther(0.5),
126 8 : fCutMassGamma(0.03),
127 8 : fCutMassGammaNSigma(5.),
128 8 : fCutMassK0(0.03),
129 8 : fCutMassK0NSigma(5.),
130 8 : fCutChi2cGamma(2.),
131 8 : fCutChi2cK0(2.),
132 8 : fCutGammaSFromDecay(-10.),
133 8 : fCutK0SFromDecay(-10.),
134 8 : fCutMaxDCA(1.),
135 : //
136 8 : fHistOn(0),
137 8 : fhClustersDPhiAcc(0),
138 8 : fhClustersDThetaAcc(0),
139 8 : fhClustersDPhiAll(0),
140 8 : fhClustersDThetaAll(0),
141 8 : fhDPhiVsDThetaAll(0),
142 8 : fhDPhiVsDThetaAcc(0),
143 8 : fhetaTracklets(0),
144 8 : fhphiTracklets(0),
145 8 : fhetaClustersLay1(0),
146 8 : fhphiClustersLay1(0),
147 : //
148 8 : fDPhiShift(0),
149 8 : fDPhiWindow2(0),
150 8 : fDThetaWindow2(0),
151 8 : fPartners(0),
152 8 : fAssociatedLay1(0),
153 8 : fMinDists(0),
154 8 : fBlackList(0),
155 : //
156 8 : fCreateClustersCopy(0),
157 8 : fClustersLoaded(0),
158 8 : fRecoDone(0),
159 8 : fBuildRefs(kTRUE),
160 8 : fStoreSPD2SingleCl(kFALSE),
161 8 : fSPDSeg()
162 40 : {
163 : // default c-tor
164 48 : for (int i=0;i<2;i++) {
165 16 : fNFiredChips[i] = 0;
166 16 : fClArr[i] = 0;
167 96 : for (int j=0;j<2;j++) fUsedClusLay[i][j] = 0;
168 16 : fDetectorIndexClustersLay[i] = 0;
169 16 : fClusterCopyIndex[i] = 0;
170 16 : fOverlapFlagClustersLay[i] = 0;
171 16 : fNClustersLay[i] = 0;
172 16 : fClustersLay[i] = 0;
173 : }
174 : // Method to reconstruct the charged particles multiplicity with the
175 : // SPD (tracklets).
176 :
177 8 : SetHistOn();
178 :
179 16 : if (AliITSReconstructor::GetRecoParam()) {
180 16 : SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow());
181 16 : SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow());
182 16 : SetPhiShift(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiShift());
183 16 : SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps());
184 16 : SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut());
185 16 : SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut());
186 16 : SetPhiRotationAngle(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiRotationAngle());
187 16 : SetNStdDev(AliITSReconstructor::GetRecoParam()->GetTrackleterNStdDevCut());
188 16 : SetScaleDThetaBySin2T(AliITSReconstructor::GetRecoParam()->GetTrackleterScaleDThetaBySin2T());
189 16 : SetBuildRefs(AliITSReconstructor::GetRecoParam()->GetTrackleterBuildCl2TrkRefs());
190 16 : SetStoreSPD2SingleCl(AliITSReconstructor::GetRecoParam()->GetTrackleterStoreSPD2SingleCl());
191 : //
192 16 : SetCutPxDrSPDin(AliITSReconstructor::GetRecoParam()->GetMultCutPxDrSPDin());
193 16 : SetCutPxDrSPDout(AliITSReconstructor::GetRecoParam()->GetMultCutPxDrSPDout());
194 16 : SetCutPxDz(AliITSReconstructor::GetRecoParam()->GetMultCutPxDz());
195 16 : SetCutDCArz(AliITSReconstructor::GetRecoParam()->GetMultCutDCArz());
196 16 : SetCutMinElectronProbTPC(AliITSReconstructor::GetRecoParam()->GetMultCutMinElectronProbTPC());
197 16 : SetCutMinElectronProbESD(AliITSReconstructor::GetRecoParam()->GetMultCutMinElectronProbESD());
198 16 : SetCutMinP(AliITSReconstructor::GetRecoParam()->GetMultCutMinP());
199 16 : SetCutMinRGamma(AliITSReconstructor::GetRecoParam()->GetMultCutMinRGamma());
200 16 : SetCutMinRK0(AliITSReconstructor::GetRecoParam()->GetMultCutMinRK0());
201 16 : SetCutMinPointAngle(AliITSReconstructor::GetRecoParam()->GetMultCutMinPointAngle());
202 16 : SetCutMaxDCADauther(AliITSReconstructor::GetRecoParam()->GetMultCutMaxDCADauther());
203 16 : SetCutMassGamma(AliITSReconstructor::GetRecoParam()->GetMultCutMassGamma());
204 16 : SetCutMassGammaNSigma(AliITSReconstructor::GetRecoParam()->GetMultCutMassGammaNSigma());
205 16 : SetCutMassK0(AliITSReconstructor::GetRecoParam()->GetMultCutMassK0());
206 16 : SetCutMassK0NSigma(AliITSReconstructor::GetRecoParam()->GetMultCutMassK0NSigma());
207 16 : SetCutChi2cGamma(AliITSReconstructor::GetRecoParam()->GetMultCutChi2cGamma());
208 16 : SetCutChi2cK0(AliITSReconstructor::GetRecoParam()->GetMultCutChi2cK0());
209 16 : SetCutGammaSFromDecay(AliITSReconstructor::GetRecoParam()->GetMultCutGammaSFromDecay());
210 16 : SetCutK0SFromDecay(AliITSReconstructor::GetRecoParam()->GetMultCutK0SFromDecay());
211 16 : SetCutMaxDCA(AliITSReconstructor::GetRecoParam()->GetMultCutMaxDCA());
212 : //
213 8 : } else {
214 0 : SetPhiWindow();
215 0 : SetThetaWindow();
216 0 : SetPhiShift();
217 0 : SetRemoveClustersFromOverlaps();
218 0 : SetPhiOverlapCut();
219 0 : SetZetaOverlapCut();
220 0 : SetPhiRotationAngle();
221 :
222 : //
223 0 : SetCutPxDrSPDin();
224 0 : SetCutPxDrSPDout();
225 0 : SetCutPxDz();
226 0 : SetCutDCArz();
227 0 : SetCutMinElectronProbTPC();
228 0 : SetCutMinElectronProbESD();
229 0 : SetCutMinP();
230 0 : SetCutMinRGamma();
231 0 : SetCutMinRK0();
232 0 : SetCutMinPointAngle();
233 0 : SetCutMaxDCADauther();
234 0 : SetCutMassGamma();
235 0 : SetCutMassGammaNSigma();
236 0 : SetCutMassK0();
237 0 : SetCutMassK0NSigma();
238 0 : SetCutChi2cGamma();
239 0 : SetCutChi2cK0();
240 0 : SetCutGammaSFromDecay();
241 0 : SetCutK0SFromDecay();
242 0 : SetCutMaxDCA();
243 : }
244 : //
245 8 : fTracklets = 0;
246 8 : fSClusters = 0;
247 : //
248 : // definition of histograms
249 16 : Bool_t oldStatus = TH1::AddDirectoryStatus();
250 8 : TH1::AddDirectory(kFALSE);
251 :
252 24 : fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,-0.1,0.1);
253 24 : fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
254 :
255 24 : fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
256 :
257 24 : fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
258 24 : fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,0.0,0.5);
259 :
260 24 : fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,0.,0.5,100,0.,0.5);
261 :
262 24 : fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
263 24 : fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
264 24 : fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
265 24 : fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
266 48 : for (int i=2;i--;) fStoreRefs[i][0] = fStoreRefs[i][1] = kFALSE;
267 8 : TH1::AddDirectory(oldStatus);
268 16 : }
269 : /*
270 : //______________________________________________________________________
271 : AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) :
272 : AliTrackleter(mr),
273 : fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fTreeRPMix(0),
274 : fTracklets(0),
275 : fSClusters(0),
276 : fNTracklets(0),
277 : fNSingleCluster(0),
278 : fNSingleClusterSPD2(0),
279 : fDPhiWindow(0),
280 : fDThetaWindow(0),
281 : fPhiShift(0),
282 : fRemoveClustersFromOverlaps(0),
283 : fPhiOverlapCut(0),
284 : fZetaOverlapCut(0),
285 : fPhiRotationAngle(0),
286 : fScaleDTBySin2T(0),
287 : fNStdDev(1.0),
288 : fNStdDevSq(1.0),
289 : //
290 : fCutPxDrSPDin(0.1),
291 : fCutPxDrSPDout(0.15),
292 : fCutPxDz(0.2),
293 : fCutDCArz(0.5),
294 : fCutMinElectronProbTPC(0.5),
295 : fCutMinElectronProbESD(0.1),
296 : fCutMinP(0.05),
297 : fCutMinRGamma(2.),
298 : fCutMinRK0(1.),
299 : fCutMinPointAngle(0.98),
300 : fCutMaxDCADauther(0.5),
301 : fCutMassGamma(0.03),
302 : fCutMassGammaNSigma(5.),
303 : fCutMassK0(0.03),
304 : fCutMassK0NSigma(5.),
305 : fCutChi2cGamma(2.),
306 : fCutChi2cK0(2.),
307 : fCutGammaSFromDecay(-10.),
308 : fCutK0SFromDecay(-10.),
309 : fCutMaxDCA(1.),
310 : //
311 : fHistOn(0),
312 : fhClustersDPhiAcc(0),
313 : fhClustersDThetaAcc(0),
314 : fhClustersDPhiAll(0),
315 : fhClustersDThetaAll(0),
316 : fhDPhiVsDThetaAll(0),
317 : fhDPhiVsDThetaAcc(0),
318 : fhetaTracklets(0),
319 : fhphiTracklets(0),
320 : fhetaClustersLay1(0),
321 : fhphiClustersLay1(0),
322 : fDPhiShift(0),
323 : fDPhiWindow2(0),
324 : fDThetaWindow2(0),
325 : fPartners(0),
326 : fAssociatedLay1(0),
327 : fMinDists(0),
328 : fBlackList(0),
329 : //
330 : fCreateClustersCopy(0),
331 : fClustersLoaded(0),
332 : fRecoDone(0),
333 : fBuildRefs(kTRUE),
334 : fStoreSPD2SingleCl(kFALSE),
335 : fSPDSeg()
336 : {
337 : // Copy constructor :!!! RS ATTENTION: old c-tor reassigned the pointers instead of creating a new copy -> would crash on delete
338 : AliError("May not use");
339 : }
340 :
341 : //______________________________________________________________________
342 : AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
343 : // Assignment operator
344 : if (this != &mr) {
345 : this->~AliITSMultReconstructor();
346 : new(this) AliITSMultReconstructor(mr);
347 : }
348 : return *this;
349 : }
350 : */
351 :
352 : //______________________________________________________________________
353 48 : AliITSMultReconstructor::~AliITSMultReconstructor(){
354 : // Destructor
355 :
356 : // delete histograms
357 16 : delete fhClustersDPhiAcc;
358 16 : delete fhClustersDThetaAcc;
359 16 : delete fhClustersDPhiAll;
360 16 : delete fhClustersDThetaAll;
361 16 : delete fhDPhiVsDThetaAll;
362 16 : delete fhDPhiVsDThetaAcc;
363 16 : delete fhetaTracklets;
364 16 : delete fhphiTracklets;
365 16 : delete fhetaClustersLay1;
366 16 : delete fhphiClustersLay1;
367 : //
368 : // delete arrays
369 208 : for(Int_t i=0; i<fNTracklets; i++) delete [] fTracklets[i];
370 :
371 176 : for(Int_t i=0; i<fNSingleCluster; i++) delete [] fSClusters[i];
372 :
373 : //
374 48 : for (int i=0;i<2;i++) {
375 32 : delete[] fClustersLay[i];
376 32 : delete[] fDetectorIndexClustersLay[i];
377 16 : delete[] fClusterCopyIndex[i];
378 32 : delete[] fOverlapFlagClustersLay[i];
379 16 : delete fClArr[i];
380 160 : for (int j=0;j<2;j++) delete fUsedClusLay[i][j];
381 : }
382 16 : delete [] fTracklets;
383 16 : delete [] fSClusters;
384 : //
385 24 : delete[] fPartners; fPartners = 0;
386 24 : delete[] fMinDists; fMinDists = 0;
387 24 : delete fBlackList; fBlackList = 0;
388 : //
389 24 : }
390 :
391 : //____________________________________________________________________
392 : void AliITSMultReconstructor::Reconstruct(AliESDEvent* esd, TTree* treeRP)
393 : {
394 16 : if (!treeRP) { AliError(" Invalid ITS cluster tree !\n"); return; }
395 8 : if (!esd) {AliError("ESDEvent is not available, use old reconstructor"); return;}
396 : // reset counters
397 16 : if (fMult) delete fMult; fMult = 0;
398 8 : fNClustersLay[0] = 0;
399 8 : fNClustersLay[1] = 0;
400 8 : fNTracklets = 0;
401 8 : fNSingleCluster = 0;
402 8 : fNSingleClusterSPD2 = 0;
403 : //
404 8 : fESDEvent = esd;
405 8 : fTreeRP = treeRP;
406 : //
407 : // >>>> RS: this part is equivalent to former AliITSVertexer::FindMultiplicity
408 : //
409 : // see if there is a SPD vertex
410 : Bool_t isVtxOK=kTRUE, isCosmics=kFALSE;
411 8 : AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
412 16 : if (!vtx || vtx->GetNContributors()<1) isVtxOK = kFALSE;
413 16 : if (vtx && strstr(vtx->GetTitle(),"cosmics")) {
414 : isVtxOK = kFALSE;
415 : isCosmics = kTRUE;
416 0 : }
417 : //
418 8 : if (!isVtxOK) {
419 0 : if (!isCosmics) {
420 0 : AliDebug(1,"Tracklets multiplicity not determined because the primary vertex was not found");
421 0 : AliDebug(1,"Just counting the number of cluster-fired chips on the SPD layers");
422 : }
423 : vtx = 0;
424 0 : }
425 8 : if(vtx){
426 8 : float vtxf[3] = {static_cast<float>(vtx->GetX()),static_cast<float>(vtx->GetY()),static_cast<float>(vtx->GetZ())};
427 8 : FindTracklets(vtxf);
428 8 : }
429 : else {
430 0 : FindTracklets(0);
431 : }
432 : //
433 8 : CreateMultiplicityObject();
434 16 : }
435 :
436 : //____________________________________________________________________
437 : void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
438 : //
439 : // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder and in analysis mode
440 :
441 0 : if (fMult) delete fMult; fMult = 0;
442 0 : fNClustersLay[0] = 0;
443 0 : fNClustersLay[1] = 0;
444 0 : fNTracklets = 0;
445 0 : fNSingleCluster = 0;
446 0 : fNSingleClusterSPD2 = 0;
447 : //
448 0 : if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
449 : //
450 0 : fESDEvent = 0;
451 0 : SetTreeRP(clusterTree);
452 : //
453 0 : FindTracklets(vtx);
454 : //
455 0 : }
456 :
457 :
458 : //____________________________________________________________________
459 : void AliITSMultReconstructor::ReconstructMix(TTree* clusterTree, TTree* clusterTreeMix, const Float_t* vtx, Float_t*)
460 : {
461 : //
462 : // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder and in analysis mode
463 :
464 0 : if (fMult) delete fMult; fMult = 0;
465 0 : fNClustersLay[0] = 0;
466 0 : fNClustersLay[1] = 0;
467 0 : fNTracklets = 0;
468 0 : fNSingleCluster = 0;
469 0 : fNSingleClusterSPD2 = 0;
470 : //
471 0 : if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
472 0 : if (!clusterTreeMix) { AliError(" Invalid ITS cluster tree 2nd event !\n"); return; }
473 : //
474 0 : fESDEvent = 0;
475 0 : SetTreeRP(clusterTree);
476 0 : SetTreeRPMix(clusterTreeMix);
477 : //
478 0 : FindTracklets(vtx);
479 : //
480 0 : }
481 :
482 :
483 : //____________________________________________________________________
484 : void AliITSMultReconstructor::FindTracklets(const Float_t *vtx)
485 : {
486 : // - calls LoadClusterArrays that finds the position of the clusters
487 : // (in global coord)
488 :
489 : // - convert the cluster coordinates to theta, phi (seen from the
490 : // interaction vertex). Clusters in the inner layer can be now
491 : // rotated for combinatorial studies
492 : // - makes an array of tracklets
493 : //
494 : // After this method has been called, the clusters of the two layers
495 : // and the tracklets can be retrieved by calling the Get'er methods.
496 :
497 :
498 : // Find tracklets converging to vertex
499 : //
500 16 : LoadClusterArrays(fTreeRP,fTreeRPMix);
501 : // flag clusters used by ESD tracks
502 16 : if (fESDEvent) ProcessESDTracks();
503 8 : fRecoDone = kTRUE;
504 :
505 8 : if (!vtx) return;
506 :
507 8 : InitAux();
508 :
509 : // find the tracklets
510 24 : AliDebug(1,"Looking for tracklets... ");
511 :
512 8 : ClusterPos2Angles(vtx); // convert cluster position to angles wrt vtx
513 : //
514 : // Step1: find all tracklets allowing double assocation:
515 : int found = 1;
516 48 : while (found > 0) {
517 : found = 0;
518 384 : for (Int_t iC1=0; iC1<fNClustersLay[0]; iC1++) found += AssociateClusterOfL1(iC1);
519 : }
520 : //
521 : // Step2: store tracklets; remove used clusters
522 132 : for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) StoreTrackletForL2Cluster(iC2);
523 : //
524 : // store unused single clusters of L1 (optionally for L2 too)
525 8 : StoreL1Singles();
526 : //
527 24 : AliDebug(1,Form("%d tracklets found", fNTracklets));
528 16 : }
529 :
530 : //____________________________________________________________________
531 : void AliITSMultReconstructor::CreateMultiplicityObject()
532 : {
533 : // create AliMultiplicity object and store it in the ESD event
534 : //
535 16 : TBits fastOrFiredMap,firedChipMap;
536 8 : if (fDetTypeRec) {
537 24 : fastOrFiredMap = fDetTypeRec->GetFastOrFiredMap();
538 24 : firedChipMap = fDetTypeRec->GetFiredChipMap(fTreeRP);
539 8 : }
540 : //
541 24 : fMult = new AliMultiplicity(fNTracklets,fNSingleCluster,fNFiredChips[0],fNFiredChips[1],fastOrFiredMap);
542 8 : fMult->SetMultTrackRefs( fBuildRefs );
543 8 : fMult->SetSPD2SinglesStored(fStoreSPD2SingleCl);
544 8 : fMult->SetNumberOfSingleClustersSPD2(fNSingleClusterSPD2);
545 : // store some details of reco:
546 8 : fMult->SetScaleDThetaBySin2T(fScaleDTBySin2T);
547 8 : fMult->SetDPhiWindow2(fDPhiWindow2);
548 8 : fMult->SetDThetaWindow2(fDThetaWindow2);
549 8 : fMult->SetDPhiShift(fDPhiShift);
550 8 : fMult->SetNStdDev(fNStdDev);
551 : //
552 8 : fMult->SetFiredChipMap(firedChipMap);
553 8 : AliITSRecPointContainer* rcont = AliITSRecPointContainer::Instance();
554 16 : fMult->SetITSClusters(0,rcont->GetNClustersInLayer(1,fTreeRP));
555 176 : for(Int_t kk=2;kk<=6;kk++) fMult->SetITSClusters(kk-1,rcont->GetNClustersInLayerFast(kk));
556 : //
557 8 : UInt_t shared[100];
558 8 : AliRefArray *refs[2][2] = {{0,0},{0,0}};
559 8 : if (fBuildRefs) {
560 48 : for (int il=2;il--;)
561 64 : for (int it=2;it--;) // tracklet_clusters->track references to stor
562 80 : if (fStoreRefs[il][it]) refs[il][it] = new AliRefArray(fNTracklets,0);
563 8 : }
564 : //
565 64 : for (int i=fNTracklets;i--;) {
566 48 : float* tlInfo = fTracklets[i];
567 48 : fMult->SetTrackletData(i,tlInfo);
568 : //
569 48 : if (!fBuildRefs) continue; // do we need references?
570 288 : for (int itp=0;itp<2;itp++) {
571 576 : for (int ilr=0;ilr<2;ilr++) {
572 192 : if (!fStoreRefs[ilr][itp]) continue; // nothing to store
573 96 : int clID = int(tlInfo[ilr ? kClID2:kClID1]);
574 96 : int nref = fUsedClusLay[ilr][itp]->GetReferences(clID,shared,100);
575 96 : if (!nref) continue;
576 192 : else if (nref==1) refs[ilr][itp]->AddReference(i,shared[0]);
577 0 : else refs[ilr][itp]->AddReferences(i,shared,nref);
578 96 : }
579 : }
580 48 : }
581 16 : if (fBuildRefs) fMult->AttachTracklet2TrackRefs(refs[0][0],refs[0][1],refs[1][0],refs[1][1]);
582 : //
583 8 : AliRefArray *refsc[2] = {0,0};
584 88 : if (fBuildRefs) for (int it=2;it--;) if (fStoreRefs[0][it]) refsc[it] = new AliRefArray(fNClustersLay[0]);
585 56 : for (int i=fNSingleCluster;i--;) {
586 40 : float* clInfo = fSClusters[i];
587 40 : fMult->SetSingleClusterData(i,clInfo);
588 : //
589 40 : if (!fBuildRefs) continue; // do we need references?
590 40 : int ilr = i>=(fNSingleCluster-fNSingleClusterSPD2) ? 1:0;
591 40 : int clID = int(clInfo[kSCID]);
592 240 : for (int itp=0;itp<2;itp++) {
593 80 : if (!fStoreRefs[ilr][itp]) continue;
594 40 : int nref = fUsedClusLay[ilr][itp]->GetReferences(clID,shared,100);
595 52 : if (!nref) continue;
596 56 : else if (nref==1) refsc[itp]->AddReference(i,shared[0]);
597 0 : else refsc[itp]->AddReferences(i,shared,nref);
598 28 : }
599 40 : }
600 : //
601 16 : if (fBuildRefs) fMult->AttachCluster2TrackRefs(refsc[0],refsc[1]);
602 8 : fMult->CompactBits();
603 : //
604 8 : }
605 :
606 :
607 : //____________________________________________________________________
608 : void AliITSMultReconstructor::LoadClusterArrays(TTree* tree, TTree* treeMix)
609 : {
610 : // load cluster info and prepare tracklets arrays
611 : //
612 16 : if (AreClustersLoaded()) {AliInfo("Clusters are already loaded"); return;}
613 8 : LoadClusterArrays(tree,0);
614 8 : LoadClusterArrays(treeMix ? treeMix:tree,1);
615 8 : int nmaxT = TMath::Min(fNClustersLay[0], fNClustersLay[1]);
616 8 : if (fTracklets) delete[] fTracklets;
617 8 : fTracklets = new Float_t*[nmaxT];
618 8 : memset(fTracklets,0,nmaxT*sizeof(Float_t*));
619 : //
620 8 : if (fSClusters) delete[] fSClusters;
621 16 : int nSlots = GetStoreSPD2SingleCl() ? fNClustersLay[0]+fNClustersLay[1] : fNClustersLay[0];
622 8 : fSClusters = new Float_t*[nSlots];
623 8 : memset(fSClusters,0,nSlots*sizeof(Float_t*));
624 : //
625 24 : AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay[0],fNClustersLay[1]));
626 24 : AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
627 8 : SetClustersLoaded();
628 16 : }
629 :
630 : //____________________________________________________________________
631 : void AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree, int il)
632 : {
633 : // This method
634 : // - gets the clusters from the cluster tree for layer il
635 : // - convert them into global coordinates
636 : // - store them in the internal arrays
637 : // - count the number of cluster-fired chips
638 : //
639 : // RS: This method was strongly modified wrt original. In order to have the same numbering
640 : // of clusters as in the ITS reco I had to introduce sorting in Z
641 : // Also note that now the clusters data are stored not in float[6] attached to float**, but in 1-D array
642 64 : AliDebug(1,Form("Loading clusters and cluster-fired chips for layer %d",il));
643 : //
644 16 : fNClustersLay[il] = 0;
645 16 : fNFiredChips[il] = 0;
646 96 : for (int i=2;i--;) fStoreRefs[il][i] = kFALSE;
647 : //
648 : AliITSRecPointContainer* rpcont = 0;
649 22 : static TClonesArray statITSrec("AliITSRecPoint");
650 22 : static TObjArray clArr(100);
651 : TBranch* branch = 0;
652 16 : TClonesArray* itsClusters = 0;
653 : //
654 16 : if (!fCreateClustersCopy) {
655 16 : rpcont=AliITSRecPointContainer::Instance();
656 16 : itsClusters = rpcont->FetchClusters(0,itsClusterTree);
657 16 : if(!rpcont->IsSPDActive()){
658 0 : AliWarning("No SPD rec points found, multiplicity not calculated");
659 0 : return;
660 : }
661 : }
662 : else {
663 0 : itsClusters = &statITSrec;
664 0 : branch = itsClusterTree->GetBranch("ITSRecPoints");
665 0 : branch->SetAddress(&itsClusters);
666 0 : if (!fClArr[il]) fClArr[il] = new TClonesArray("AliITSRecPoint",100);
667 0 : delete[] fClusterCopyIndex[il];
668 : }
669 : //
670 : // count clusters
671 : // loop over the SPD subdetectors
672 : int nclLayer = 0;
673 16 : int detMin = TMath::Max(0,AliITSgeomTGeo::GetModuleIndex(il+1,1,1));
674 16 : int detMax = AliITSgeomTGeo::GetModuleIndex(il+2,1,1);
675 3872 : for (int idt=detMin;idt<detMax;idt++) {
676 3840 : if (!fCreateClustersCopy) itsClusters = rpcont->UncheckedGetClusters(idt);
677 0 : else branch->GetEvent(idt);
678 1920 : int nClusters = itsClusters->GetEntriesFast();
679 3710 : if (!nClusters) continue;
680 130 : Int_t nClustersInChip[5] = {0,0,0,0,0};
681 406 : while(nClusters--) {
682 146 : AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
683 146 : if (!cluster) continue;
684 146 : if (fCreateClustersCopy) cluster = new ((*fClArr[il])[nclLayer]) AliITSRecPoint(*cluster);
685 146 : clArr.AddAtAndExpand(cluster,nclLayer++);
686 146 : Int_t chipNo = fSPDSeg.GetChipFromLocal(0,cluster->GetDetLocalZ());
687 292 : if(chipNo>=0)nClustersInChip[ chipNo ]++;
688 146 : }
689 1698 : for(Int_t ifChip=5;ifChip--;) if (nClustersInChip[ifChip]) fNFiredChips[il]++;
690 130 : }
691 : // sort the clusters in Z (to have the same numbering as in ITS reco
692 16 : Float_t *z = new Float_t[nclLayer];
693 16 : Int_t *index = new Int_t[nclLayer];
694 324 : for (int ic=0;ic<nclLayer;ic++) z[ic] = ((AliITSRecPoint*)clArr[ic])->GetZ();
695 16 : TMath::Sort(nclLayer,z,index,kFALSE);
696 16 : Float_t* clustersLay = new Float_t[nclLayer*kClNPar];
697 16 : Int_t* detectorIndexClustersLay = new Int_t[nclLayer];
698 16 : Bool_t* overlapFlagClustersLay = new Bool_t[nclLayer];
699 16 : if (fCreateClustersCopy) fClusterCopyIndex[il] = new Int_t[nclLayer];
700 : //
701 324 : for (int ic=0;ic<nclLayer;ic++) {
702 146 : AliITSRecPoint* cluster = (AliITSRecPoint*)clArr[index[ic]];
703 146 : float* clPar = &clustersLay[ic*kClNPar];
704 : //
705 146 : cluster->GetGlobalXYZ( clPar );
706 146 : detectorIndexClustersLay[ic] = cluster->GetDetectorIndex();
707 146 : overlapFlagClustersLay[ic] = kFALSE;
708 1168 : for (Int_t i=3;i--;) clPar[kClMC0+i] = cluster->GetLabel(i);
709 146 : if (fCreateClustersCopy) fClusterCopyIndex[il][ic] = index[ic];
710 : }
711 16 : clArr.Clear();
712 32 : delete[] z;
713 32 : delete[] index;
714 : //
715 16 : if (fOverlapFlagClustersLay[il]) delete[] fOverlapFlagClustersLay[il];
716 16 : fOverlapFlagClustersLay[il] = overlapFlagClustersLay;
717 : //
718 16 : if (fDetectorIndexClustersLay[il]) delete[] fDetectorIndexClustersLay[il];
719 16 : fDetectorIndexClustersLay[il] = detectorIndexClustersLay;
720 : //
721 16 : if (fBuildRefs) {
722 96 : for (int it=0;it<2;it++) {
723 32 : if (fUsedClusLay[il][it]) delete fUsedClusLay[il][it];
724 64 : fUsedClusLay[il][it] = new AliRefArray(nclLayer);
725 : }
726 16 : }
727 : //
728 16 : if (fClustersLay[il]) delete[] fClustersLay[il];
729 16 : fClustersLay[il] = clustersLay;
730 16 : fNClustersLay[il] = nclLayer;
731 : //
732 32 : }
733 :
734 : //____________________________________________________________________
735 : void AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
736 : // This method
737 : // - gets the clusters from the cluster tree
738 : // - counts the number of (cluster)fired chips
739 :
740 0 : AliDebug(1,"Loading cluster-fired chips ...");
741 :
742 0 : fNFiredChips[0] = 0;
743 0 : fNFiredChips[1] = 0;
744 :
745 0 : AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
746 : TClonesArray* itsClusters=NULL;
747 0 : rpcont->FetchClusters(0,itsClusterTree);
748 0 : if(!rpcont->IsSPDActive()){
749 0 : AliWarning("No SPD rec points found, multiplicity not calculated");
750 0 : return;
751 : }
752 :
753 : // loop over the its subdetectors
754 0 : Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1);
755 0 : for (Int_t iIts=0; iIts < nSPDmodules; iIts++) {
756 0 : itsClusters=rpcont->UncheckedGetClusters(iIts);
757 0 : Int_t nClusters = itsClusters->GetEntriesFast();
758 :
759 : // number of clusters in each chip of the current module
760 0 : Int_t nClustersInChip[5] = {0,0,0,0,0};
761 0 : Int_t layer = 0;
762 0 : Int_t ladder=0;
763 0 : Int_t det=0;
764 0 : AliITSgeomTGeo::GetModuleId(iIts,layer,ladder,det);
765 0 : --layer; // layer is from 1 to 6 in AliITSgeomTGeo, but from 0 to 5 here
766 0 : if(layer<0 || layer >1)continue;
767 :
768 : // loop over clusters
769 0 : while(nClusters--) {
770 0 : AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
771 :
772 : // find the chip for the current cluster
773 0 : Float_t locz = cluster->GetDetLocalZ();
774 0 : Int_t iChip = fSPDSeg.GetChipFromLocal(0,locz);
775 0 : if (iChip>=0) nClustersInChip[iChip]++;
776 :
777 : }// end of cluster loop
778 :
779 : // get number of fired chips in the current module
780 0 : for(Int_t ifChip=0; ifChip<5; ifChip++) {
781 0 : if(nClustersInChip[ifChip] >= 1) fNFiredChips[layer]++;
782 : }
783 :
784 0 : } // end of its "subdetector" loop
785 :
786 :
787 0 : AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
788 0 : }
789 : //____________________________________________________________________
790 : void
791 : AliITSMultReconstructor::SaveHists() {
792 : // This method save the histograms on the output file
793 : // (only if fHistOn is TRUE).
794 :
795 0 : if (!fHistOn)
796 : return;
797 :
798 0 : fhClustersDPhiAll->Write();
799 0 : fhClustersDThetaAll->Write();
800 0 : fhDPhiVsDThetaAll->Write();
801 :
802 0 : fhClustersDPhiAcc->Write();
803 0 : fhClustersDThetaAcc->Write();
804 0 : fhDPhiVsDThetaAcc->Write();
805 :
806 0 : fhetaTracklets->Write();
807 0 : fhphiTracklets->Write();
808 0 : fhetaClustersLay1->Write();
809 0 : fhphiClustersLay1->Write();
810 0 : }
811 :
812 : //____________________________________________________________________
813 : void AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist)
814 : {
815 : // Flags clusters in the overlapping regions
816 : Float_t distClSameMod=0.;
817 : Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
818 : Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
819 :
820 : Float_t zproj1=0.;
821 : Float_t zproj2=0.;
822 : Float_t deZproj=0.;
823 0 : Float_t* clPar1 = GetClusterLayer1(iC1);
824 0 : Float_t* clPar2B = GetClusterLayer2(iC2WithBestDist);
825 : // Loop on inner layer clusters
826 0 : for (Int_t iiC1=0; iiC1<fNClustersLay[0]; iiC1++) {
827 0 : if (!fOverlapFlagClustersLay[0][iiC1]) {
828 : // only for adjacent modules
829 0 : if ((TMath::Abs(fDetectorIndexClustersLay[0][iC1]-fDetectorIndexClustersLay[0][iiC1])==4)||
830 0 : (TMath::Abs(fDetectorIndexClustersLay[0][iC1]-fDetectorIndexClustersLay[0][iiC1])==76)) {
831 0 : Float_t *clPar11 = GetClusterLayer1(iiC1);
832 0 : Float_t dePhi=TMath::Abs(clPar11[kClPh]-clPar1[kClPh]);
833 0 : if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
834 :
835 0 : zproj1=meanRadiusLay1/TMath::Tan(clPar1[kClTh]);
836 0 : zproj2=meanRadiusLay1/TMath::Tan(clPar11[kClTh]);
837 :
838 0 : deZproj=TMath::Abs(zproj1-zproj2);
839 :
840 0 : distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
841 0 : if (distClSameMod<=1.) fOverlapFlagClustersLay[0][iiC1]=kTRUE;
842 :
843 0 : } // end adjacent modules
844 : }
845 : } // end Loop on inner layer clusters
846 :
847 :
848 : distClSameMod=0.;
849 : // Loop on outer layer clusters
850 0 : for (Int_t iiC2=0; iiC2<fNClustersLay[1]; iiC2++) {
851 0 : if (!fOverlapFlagClustersLay[1][iiC2]) {
852 : // only for adjacent modules
853 0 : Float_t *clPar2 = GetClusterLayer2(iiC2);
854 0 : if ((TMath::Abs(fDetectorIndexClustersLay[1][iC2WithBestDist]-fDetectorIndexClustersLay[1][iiC2])==4) ||
855 0 : (TMath::Abs(fDetectorIndexClustersLay[1][iC2WithBestDist]-fDetectorIndexClustersLay[1][iiC2])==156)) {
856 0 : Float_t dePhi=TMath::Abs(clPar2[kClPh]-clPar2B[kClPh]);
857 0 : if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
858 :
859 0 : zproj1=meanRadiusLay2/TMath::Tan(clPar2B[kClTh]);
860 0 : zproj2=meanRadiusLay2/TMath::Tan(clPar2[kClTh]);
861 :
862 0 : deZproj=TMath::Abs(zproj1-zproj2);
863 0 : distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
864 0 : if (distClSameMod<=1.) fOverlapFlagClustersLay[1][iiC2]=kTRUE;
865 :
866 0 : } // end adjacent modules
867 0 : }
868 : } // end Loop on outer layer clusters
869 :
870 0 : }
871 :
872 : //____________________________________________________________________
873 : void AliITSMultReconstructor::InitAux()
874 : {
875 : // init arrays/parameters for tracklet reconstruction
876 :
877 : // dPhi shift is field dependent, get average magnetic field
878 : Float_t bz = 0;
879 : AliMagF* field = 0;
880 40 : if (TGeoGlobalMagField::Instance()) field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
881 8 : if (!field) {
882 0 : AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor.");
883 0 : }
884 8 : else bz = TMath::Abs(field->SolenoidField());
885 8 : fDPhiShift = fPhiShift / 5 * bz;
886 24 : AliDebug(1, Form("Using phi shift of %f", fDPhiShift));
887 : //
888 16 : if (fPartners) delete[] fPartners; fPartners = new Int_t[fNClustersLay[1]];
889 16 : if (fMinDists) delete[] fMinDists; fMinDists = new Float_t[fNClustersLay[1]];
890 16 : if (fAssociatedLay1) delete[] fAssociatedLay1; fAssociatedLay1 = new Int_t[fNClustersLay[0]];
891 : //
892 24 : if (fBlackList) delete fBlackList; fBlackList = new AliRefArray(fNClustersLay[0]);
893 : //
894 : // Printf("Vertex in find tracklets...%f %f %f",vtx[0],vtx[1],vtx[2]);
895 132 : for (Int_t i=0; i<fNClustersLay[1]; i++) {
896 58 : fPartners[i] = -1;
897 58 : fMinDists[i] = 2*fNStdDev;
898 : }
899 8 : memset(fAssociatedLay1,0,fNClustersLay[0]*sizeof(Int_t));
900 : //
901 8 : }
902 :
903 : //____________________________________________________________________
904 : void AliITSMultReconstructor::ClusterPos2Angles(const Float_t *vtx)
905 : {
906 : // convert cluster coordinates to angles wrt vertex
907 56 : for (int ilr=0;ilr<2;ilr++) {
908 324 : for (Int_t iC=0; iC<fNClustersLay[ilr]; iC++) {
909 146 : float* clPar = GetClusterOfLayer(ilr,iC);
910 146 : CalcThetaPhi(clPar[kClTh]-vtx[0],clPar[kClPh]-vtx[1],clPar[kClZ]-vtx[2],clPar[kClTh],clPar[kClPh]);
911 146 : if (ilr==0) {
912 88 : clPar[kClPh] = clPar[kClPh] + fPhiRotationAngle; // rotation of inner layer for comb studies
913 88 : if (fHistOn) {
914 0 : Float_t eta = clPar[kClTh];
915 0 : eta= TMath::Tan(eta/2.);
916 0 : eta=-TMath::Log(eta);
917 0 : fhetaClustersLay1->Fill(eta);
918 0 : fhphiClustersLay1->Fill(clPar[kClPh]);
919 0 : }
920 : }
921 : }
922 : }
923 : //
924 8 : }
925 :
926 : //____________________________________________________________________
927 : Int_t AliITSMultReconstructor::AssociateClusterOfL1(Int_t iC1)
928 : {
929 : // search association of cluster iC1 of L1 with all clusters of L2
930 440 : if (fAssociatedLay1[iC1] != 0) return 0;
931 : Int_t iC2WithBestDist = -1; // reset
932 88 : Double_t minDist = 2*fNStdDev; // reset
933 88 : float* clPar1 = GetClusterLayer1(iC1);
934 1452 : for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) {
935 : //
936 638 : if (fBlackList->IsReferred(iC1,iC2)) continue;
937 638 : float* clPar2 = GetClusterLayer2(iC2);
938 : //
939 : // find the difference in angles
940 638 : Double_t dTheta = TMath::Abs(clPar2[kClTh] - clPar1[kClTh]);
941 638 : Double_t dPhi = TMath::Abs(clPar2[kClPh] - clPar1[kClPh]);
942 : // Printf("detheta %f dephi %f", dTheta,dPhi);
943 : //
944 712 : if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi; // take into account boundary condition
945 : //
946 638 : if (fHistOn) {
947 0 : fhClustersDPhiAll->Fill(dPhi);
948 0 : fhClustersDThetaAll->Fill(dTheta);
949 0 : fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
950 0 : }
951 638 : Float_t d = CalcDist(dPhi,clPar2[kClTh] - clPar1[kClTh],clPar1[kClTh]); // make "elliptical" cut in Phi and Theta!
952 : // look for the minimum distance: the minimum is in iC2WithBestDist
953 734 : if (d<fNStdDev && d<minDist) { minDist=d; iC2WithBestDist = iC2; }
954 638 : }
955 : //
956 88 : if (minDist<fNStdDev) { // This means that a cluster in layer 2 was found that matches with iC1
957 : //
958 48 : if (fMinDists[iC2WithBestDist] > minDist) {
959 48 : Int_t oldPartner = fPartners[iC2WithBestDist];
960 48 : fPartners[iC2WithBestDist] = iC1;
961 48 : fMinDists[iC2WithBestDist] = minDist;
962 : //
963 48 : fAssociatedLay1[iC1] = 1; // mark as assigned
964 : //
965 48 : if (oldPartner != -1) {
966 : // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its fBlackList
967 0 : fBlackList->AddReference(oldPartner,iC2WithBestDist);
968 0 : fAssociatedLay1[oldPartner] = 0; // mark as free
969 0 : }
970 48 : } else {
971 : // try again to find a cluster without considering iC2WithBestDist
972 0 : fBlackList->AddReference(iC1,iC2WithBestDist);
973 : }
974 : //
975 : }
976 40 : else fAssociatedLay1[iC1] = 2;// cluster has no partner; remove
977 : //
978 : return 1;
979 176 : }
980 :
981 : //____________________________________________________________________
982 : Int_t AliITSMultReconstructor::StoreTrackletForL2Cluster(Int_t iC2)
983 : {
984 : // build tracklet for cluster iC2 of layer 2
985 126 : if (fPartners[iC2] == -1) return 0;
986 48 : if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (fPartners[iC2],iC2);
987 : // Printf("saving tracklets");
988 96 : if (fOverlapFlagClustersLay[0][fPartners[iC2]] || fOverlapFlagClustersLay[1][iC2]) return 0;
989 48 : float* clPar2 = GetClusterLayer2(iC2);
990 48 : float* clPar1 = GetClusterLayer1(fPartners[iC2]);
991 : //
992 48 : Float_t* tracklet = fTracklets[fNTracklets] = new Float_t[kTrNPar]; // RS Add also the cluster id's
993 : //
994 48 : tracklet[kTrTheta] = clPar1[kClTh]; // use the theta from the clusters in the first layer
995 48 : tracklet[kTrPhi] = clPar1[kClPh]; // use the phi from the clusters in the first layer
996 48 : tracklet[kTrDPhi] = clPar1[kClPh] - clPar2[kClPh]; // store the difference between phi1 and phi2
997 : //
998 : // define dphi in the range [0,pi] with proper sign (track charge correlated)
999 48 : if (tracklet[kTrDPhi] > TMath::Pi()) tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi();
1000 48 : if (tracklet[kTrDPhi] < -TMath::Pi()) tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi();
1001 : //
1002 48 : tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh]; // store the theta1-theta2
1003 : //
1004 48 : if (fHistOn) {
1005 0 : fhClustersDPhiAcc->Fill(tracklet[kTrDPhi]);
1006 0 : fhClustersDThetaAcc->Fill(tracklet[kTrDTheta]);
1007 0 : fhDPhiVsDThetaAcc->Fill(tracklet[kTrDTheta],tracklet[kTrDPhi]);
1008 0 : }
1009 : //
1010 : // find label
1011 : // if equal label in both clusters found this label is assigned
1012 : // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
1013 : Int_t label1=0,label2=0;
1014 313 : while (label2 < 3) {
1015 266 : if ( int(clPar1[kClMC0+label1])!=-2 && int(clPar1[kClMC0+label1])==int(clPar2[kClMC0+label2])) break;
1016 289 : if (++label1 == 3) { label1 = 0; label2++; }
1017 : }
1018 48 : if (label2 < 3) {
1019 120 : AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n",
1020 : (Int_t) clPar1[kClMC0+label1], (Int_t) clPar1[kClMC0+label2], fNTracklets));
1021 24 : tracklet[kTrLab1] = tracklet[kTrLab2] = clPar1[kClMC0+label1];
1022 24 : } else {
1023 72 : AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n",
1024 : (Int_t) clPar1[kClMC0], (Int_t) clPar1[kClMC1], (Int_t) clPar1[kClMC2],
1025 : (Int_t) clPar2[kClMC0], (Int_t) clPar2[kClMC1], (Int_t) clPar2[kClMC2], fNTracklets));
1026 24 : tracklet[kTrLab1] = clPar1[kClMC0];
1027 24 : tracklet[kTrLab2] = clPar2[kClMC0];
1028 : }
1029 : //
1030 48 : if (fHistOn) {
1031 0 : Float_t eta = tracklet[kTrTheta];
1032 0 : eta= TMath::Tan(eta/2.);
1033 0 : eta=-TMath::Log(eta);
1034 0 : fhetaTracklets->Fill(eta);
1035 0 : fhphiTracklets->Fill(tracklet[kTrPhi]);
1036 0 : }
1037 : //
1038 48 : tracklet[kClID1] = fPartners[iC2];
1039 48 : tracklet[kClID2] = iC2;
1040 : //
1041 : // Printf("Adding tracklet candidate");
1042 144 : AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
1043 144 : AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", fPartners[iC2], iC2));
1044 48 : fNTracklets++;
1045 48 : fAssociatedLay1[fPartners[iC2]] = 1;
1046 : //
1047 : return 1;
1048 58 : }
1049 :
1050 : //____________________________________________________________________
1051 : void AliITSMultReconstructor::StoreL1Singles()
1052 : {
1053 : // Printf("saving single clusters...");
1054 200 : for (Int_t iC1=0; iC1<fNClustersLay[0]; iC1++) {
1055 88 : float* clPar1 = GetClusterLayer1(iC1);
1056 136 : if (fAssociatedLay1[iC1]==2||fAssociatedLay1[iC1]==0) {
1057 40 : fSClusters[fNSingleCluster] = new Float_t[kClNPar];
1058 40 : fSClusters[fNSingleCluster][kSCTh] = clPar1[kClTh];
1059 40 : fSClusters[fNSingleCluster][kSCPh] = clPar1[kClPh];
1060 40 : fSClusters[fNSingleCluster][kSCLab] = clPar1[kClMC0];
1061 40 : fSClusters[fNSingleCluster][kSCID] = iC1;
1062 120 : AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
1063 : fNSingleCluster, iC1));
1064 40 : fNSingleCluster++;
1065 40 : }
1066 : }
1067 : //
1068 8 : if (GetStoreSPD2SingleCl()) {
1069 0 : for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) {
1070 0 : if (fPartners[iC2]<0 || (fOverlapFlagClustersLay[0][fPartners[iC2]] || fOverlapFlagClustersLay[1][iC2])) {
1071 0 : float* clPar2 = GetClusterLayer2(iC2);
1072 0 : fSClusters[fNSingleCluster] = new Float_t[kClNPar];
1073 0 : fSClusters[fNSingleCluster][kSCTh] = clPar2[kClTh];
1074 0 : fSClusters[fNSingleCluster][kSCPh] = clPar2[kClPh];
1075 0 : fSClusters[fNSingleCluster][kSCLab] = clPar2[kClMC0];
1076 0 : fSClusters[fNSingleCluster][kSCID] = iC2;
1077 0 : AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 2)",
1078 : fNSingleCluster, iC2));
1079 0 : fNSingleCluster++;
1080 0 : fNSingleClusterSPD2++;
1081 0 : }
1082 : }
1083 0 : }
1084 : //
1085 8 : }
1086 :
1087 : //____________________________________________________________________
1088 : void AliITSMultReconstructor::ProcessESDTracks()
1089 : {
1090 : // Flag the clusters used by ESD tracks
1091 : // Flag primary tracks to be used for multiplicity counting
1092 : //
1093 24 : if (!fESDEvent || !fBuildRefs) return;
1094 8 : AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexTracks();
1095 16 : if (!vtx || vtx->GetNContributors()<1) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
1096 16 : if (!vtx || vtx->GetNContributors()<1) {
1097 0 : AliDebug(1,"No primary vertex: cannot flag primary tracks");
1098 0 : return;
1099 : }
1100 8 : Int_t ntracks = fESDEvent->GetNumberOfTracks();
1101 300 : for(Int_t itr=0; itr<ntracks; itr++) {
1102 142 : AliESDtrack* track = fESDEvent->GetTrack(itr);
1103 186 : if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx
1104 98 : FlagTrackClusters(itr);
1105 98 : FlagIfSecondary(track,vtx);
1106 98 : }
1107 8 : FlagV0s(vtx);
1108 : //
1109 16 : }
1110 :
1111 : //____________________________________________________________________
1112 : void AliITSMultReconstructor::FlagTrackClusters(Int_t id)
1113 : {
1114 : // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation
1115 : //
1116 196 : const AliESDtrack* track = fESDEvent->GetTrack(id);
1117 98 : Int_t idx[12];
1118 98 : if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit
1119 98 : Int_t itsType = track->IsOn(AliESDtrack::kITSpureSA) ? 1:0;
1120 :
1121 784 : for (int i=6/*AliESDfriendTrack::kMaxITScluster*/;i--;) { // ignore extras: note: i>=6 is for extra clusters
1122 588 : if (idx[i]<0) continue;
1123 514 : int layID= (idx[i] & 0xf0000000) >> 28;
1124 894 : if (layID>1) continue; // SPD only
1125 134 : int clID = (idx[i] & 0x0fffffff);
1126 134 : fUsedClusLay[layID][itsType]->AddReference(clID,id);
1127 134 : fStoreRefs[layID][itsType] = kTRUE;
1128 134 : }
1129 : //
1130 196 : }
1131 :
1132 : //____________________________________________________________________
1133 : void AliITSMultReconstructor::FlagIfSecondary(AliESDtrack* track, const AliVertex* vtx)
1134 : {
1135 : // RS: check if the track is primary and set the flag
1136 414 : double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? fCutPxDrSPDin:fCutPxDrSPDout;
1137 98 : float xz[2];
1138 98 : track->GetDZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), fESDEvent->GetMagneticField(), xz);
1139 278 : if (TMath::Abs(xz[0]*track->P())>cut || TMath::Abs(xz[1]*track->P())>fCutPxDz ||
1140 180 : TMath::Abs(xz[0])>fCutDCArz || TMath::Abs(xz[1])>fCutDCArz)
1141 8 : track->SetStatus(AliESDtrack::kMultSec);
1142 90 : else track->ResetStatus(AliESDtrack::kMultSec);
1143 98 : }
1144 :
1145 : //____________________________________________________________________
1146 : void AliITSMultReconstructor::FlagV0s(const AliESDVertex *vtx)
1147 : {
1148 : // flag tracks belonging to v0s
1149 : //
1150 : const double kK0Mass = 0.4976;
1151 : //
1152 16 : AliV0 pvertex;
1153 8 : AliKFVertex vertexKF;
1154 56 : AliKFParticle epKF0,epKF1,pipmKF0,piKF0,piKF1,gammaKF,k0KF;
1155 8 : Double_t mass,massErr,chi2c;
1156 : enum {kKFIni=BIT(14)};
1157 : //
1158 8 : double recVtx[3];
1159 8 : float recVtxF[3];
1160 8 : vtx->GetXYZ(recVtx);
1161 64 : for (int i=3;i--;) recVtxF[i] = recVtx[i];
1162 : //
1163 8 : int ntracks = fESDEvent->GetNumberOfTracks();
1164 8 : if (ntracks<2) return;
1165 : //
1166 16 : vertexKF.X() = recVtx[0];
1167 16 : vertexKF.Y() = recVtx[1];
1168 16 : vertexKF.Z() = recVtx[2];
1169 16 : vertexKF.Covariance(0,0) = vtx->GetXRes()*vtx->GetXRes();
1170 16 : vertexKF.Covariance(1,2) = vtx->GetYRes()*vtx->GetYRes();
1171 16 : vertexKF.Covariance(2,2) = vtx->GetZRes()*vtx->GetZRes();
1172 : //
1173 : AliESDtrack *trc0,*trc1;
1174 300 : for (int it0=0;it0<ntracks;it0++) {
1175 142 : trc0 = fESDEvent->GetTrack(it0);
1176 284 : if (trc0->IsOn(AliESDtrack::kMultInV0)) continue;
1177 284 : if (!trc0->IsOn(AliESDtrack::kITSin)) continue;
1178 196 : Bool_t isSAP = trc0->IsPureITSStandalone();
1179 196 : Int_t q0 = trc0->Charge();
1180 196 : Bool_t testGamma = CanBeElectron(trc0);
1181 98 : epKF0.ResetBit(kKFIni);
1182 98 : piKF0.ResetBit(kKFIni);
1183 : double bestChi2=1e16;
1184 : int bestID = -1;
1185 : //
1186 2236 : for (int it1=it0+1;it1<ntracks;it1++) {
1187 1020 : trc1 = fESDEvent->GetTrack(it1);
1188 2040 : if (trc1->IsOn(AliESDtrack::kMultInV0)) continue;
1189 2040 : if (!trc1->IsOn(AliESDtrack::kITSin)) continue;
1190 1140 : if (trc1->IsPureITSStandalone() != isSAP) continue; // pair separately ITS_SA_Pure tracks and TPC/ITS+ITS_SA
1191 1140 : if ( (q0+trc1->Charge())!=0 ) continue; // don't pair like signs
1192 : //
1193 306 : pvertex.SetParamN(q0<0 ? *trc0:*trc1);
1194 306 : pvertex.SetParamP(q0>0 ? *trc0:*trc1);
1195 306 : pvertex.Update(recVtxF);
1196 612 : if (pvertex.P()<fCutMinP) continue;
1197 306 : if (pvertex.GetV0CosineOfPointingAngle()<fCutMinPointAngle) continue;
1198 54 : if (pvertex.GetDcaV0Daughters()>fCutMaxDCADauther) continue;
1199 76 : double d = pvertex.GetD(recVtx[0],recVtx[1],recVtx[2]);
1200 38 : if (d>fCutMaxDCA) continue;
1201 114 : double dx=recVtx[0]-pvertex.Xv(), dy=recVtx[1]-pvertex.Yv();
1202 38 : double rv = TMath::Sqrt(dx*dx+dy*dy);
1203 : //
1204 : // check gamma conversion hypothesis ----------------------------------------------------------->>>
1205 : Bool_t gammaOK = kFALSE;
1206 114 : while (testGamma && CanBeElectron(trc1)) {
1207 38 : if (rv<fCutMinRGamma) break;
1208 8 : if (!epKF0.TestBit(kKFIni)) {
1209 8 : new(&epKF0) AliKFParticle(*trc0,q0>0 ? kPositron:kElectron);
1210 4 : epKF0.SetBit(kKFIni);
1211 4 : }
1212 16 : new(&epKF1) AliKFParticle(*trc1,q0<0 ? kPositron:kElectron);
1213 8 : gammaKF.Initialize();
1214 8 : gammaKF += epKF0;
1215 8 : gammaKF += epKF1;
1216 8 : gammaKF.SetProductionVertex(vertexKF);
1217 8 : gammaKF.GetMass(mass,massErr);
1218 16 : if (mass>fCutMassGamma || (massErr>0&&(mass>massErr*fCutMassGammaNSigma))) break;
1219 8 : if (gammaKF.GetS()<fCutGammaSFromDecay) break;
1220 4 : gammaKF.SetMassConstraint(0.,0.001);
1221 24 : chi2c = (gammaKF.GetNDF()!=0) ? gammaKF.GetChi2()/gammaKF.GetNDF() : 1000;
1222 4 : if (chi2c>fCutChi2cGamma) break;
1223 : gammaOK = kTRUE;
1224 0 : if (chi2c>bestChi2) break;
1225 : bestChi2 = chi2c;
1226 : bestID = it1;
1227 0 : break;
1228 : }
1229 38 : if (gammaOK) continue;
1230 : // check gamma conversion hypothesis -----------------------------------------------------------<<<
1231 : // check K0 conversion hypothesis ----------------------------------------------------------->>>
1232 : while (1) {
1233 38 : if (rv<fCutMinRK0) break;
1234 10 : if (!piKF0.TestBit(kKFIni)) {
1235 12 : new(&piKF0) AliKFParticle(*trc0,q0>0 ? kPiPlus:kPiMinus);
1236 6 : piKF0.SetBit(kKFIni);
1237 6 : }
1238 20 : new(&piKF1) AliKFParticle(*trc1,q0<0 ? kPiPlus:kPiMinus);
1239 10 : k0KF.Initialize();
1240 10 : k0KF += piKF0;
1241 10 : k0KF += piKF1;
1242 10 : k0KF.SetProductionVertex(vertexKF);
1243 10 : k0KF.GetMass(mass,massErr);
1244 10 : mass -= kK0Mass;
1245 10 : if (TMath::Abs(mass)>fCutMassK0 || (massErr>0&&(std::abs(mass)>massErr*fCutMassK0NSigma))) break;
1246 0 : if (k0KF.GetS()<fCutK0SFromDecay) break;
1247 0 : k0KF.SetMassConstraint(kK0Mass,0.001);
1248 0 : chi2c = (k0KF.GetNDF()!=0) ? k0KF.GetChi2()/k0KF.GetNDF() : 1000;
1249 0 : if (chi2c>fCutChi2cK0) break;
1250 0 : if (chi2c>bestChi2) break;
1251 : bestChi2 = chi2c;
1252 : bestID = it1;
1253 0 : break;
1254 : }
1255 : // check K0 conversion hypothesis -----------------------------------------------------------<<<
1256 38 : }
1257 : //
1258 98 : if (bestID>=0) {
1259 0 : trc0->SetStatus(AliESDtrack::kMultInV0);
1260 0 : fESDEvent->GetTrack(bestID)->SetStatus(AliESDtrack::kMultInV0);
1261 : }
1262 98 : }
1263 : //
1264 16 : }
1265 :
1266 : //____________________________________________________________________
1267 : Bool_t AliITSMultReconstructor::CanBeElectron(const AliESDtrack* trc) const
1268 : {
1269 : // check if the track can be electron
1270 272 : Double_t pid[AliPID::kSPECIES];
1271 272 : if (!trc->IsOn(AliESDtrack::kESDpid)) return kTRUE;
1272 0 : trc->GetESDpid(pid);
1273 0 : return (trc->IsOn(AliESDtrack::kTPCpid)) ?
1274 0 : pid[AliPID::kElectron]>fCutMinElectronProbTPC :
1275 0 : pid[AliPID::kElectron]>fCutMinElectronProbESD;
1276 : //
1277 136 : }
1278 :
1279 : //____________________________________________________________________
1280 : AliITSRecPoint* AliITSMultReconstructor::GetRecPoint(Int_t lr, Int_t n) const
1281 : {
1282 : // return a cluster of lr corresponding to orderer cluster index n
1283 0 : if (fClArr[lr] && fClusterCopyIndex[lr] && n<fNClustersLay[lr])
1284 0 : return (AliITSRecPoint*) fClArr[lr]->At(fClusterCopyIndex[lr][n]);
1285 : else {
1286 0 : AliError("To access the clusters SetCreateClustersCopy should have been called");
1287 0 : return 0;
1288 : }
1289 0 : }
|
