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