| | 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 | ClassImp(AliITSMultReconstructor) |
| | 94 | |
| | 95 | |
| | 96 | //____________________________________________________________________ |
| | 97 | AliITSMultReconstructor::AliITSMultReconstructor(): |
| | 98 | fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fTreeRPMix(0), |
| | 99 | fTracklets(0), |
| | 100 | fSClusters(0), |
| | 101 | fNTracklets(0), |
| | 102 | fNSingleCluster(0), |
| | 103 | fNSingleClusterSPD2(0), |
| | 104 | fDPhiWindow(0), |
| | 105 | fDThetaWindow(0), |
| | 106 | fPhiShift(0), |
| | 107 | fRemoveClustersFromOverlaps(0), |
| | 108 | fPhiOverlapCut(0), |
| | 109 | fZetaOverlapCut(0), |
| | 110 | fPhiRotationAngle(0), |
| | 111 | fScaleDTBySin2T(0), |
| | 112 | fNStdDev(1.0), |
| | 113 | fNStdDevSq(1.0), |
| | 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 | // |
| | 136 | fHistOn(0), |
| | 137 | fhClustersDPhiAcc(0), |
| | 138 | fhClustersDThetaAcc(0), |
| | 139 | fhClustersDPhiAll(0), |
| | 140 | fhClustersDThetaAll(0), |
| | 141 | fhDPhiVsDThetaAll(0), |
| | 142 | fhDPhiVsDThetaAcc(0), |
| | 143 | fhetaTracklets(0), |
| | 144 | fhphiTracklets(0), |
| | 145 | fhetaClustersLay1(0), |
| | 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), |
| | 159 | fBuildRefs(kTRUE), |
| | 160 | fStoreSPD2SingleCl(kFALSE), |
| | 161 | fSPDSeg() |
| | 162 | { |
| | 163 | // default c-tor |
| | 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; |
| | 169 | fClusterCopyIndex[i] = 0; |
| | 170 | fOverlapFlagClustersLay[i] = 0; |
| | 171 | fNClustersLay[i] = 0; |
| | 172 | fClustersLay[i] = 0; |
| | 173 | } |
| | 174 | // Method to reconstruct the charged particles multiplicity with the |
| | 175 | // SPD (tracklets). |
| | 176 | |
| | 177 | SetHistOn(); |
| | 178 | |
| | 179 | if (AliITSReconstructor::GetRecoParam()) { |
| | 180 | SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow()); |
| | 181 | SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow()); |
| | 182 | SetPhiShift(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiShift()); |
| | 183 | SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps()); |
| | 184 | SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut()); |
| | 185 | SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut()); |
| | 186 | SetPhiRotationAngle(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiRotationAngle()); |
| | 187 | SetNStdDev(AliITSReconstructor::GetRecoParam()->GetTrackleterNStdDevCut()); |
| | 188 | SetScaleDThetaBySin2T(AliITSReconstructor::GetRecoParam()->GetTrackleterScaleDThetaBySin2T()); |
| | 189 | SetBuildRefs(AliITSReconstructor::GetRecoParam()->GetTrackleterBuildCl2TrkRefs()); |
| | 190 | SetStoreSPD2SingleCl(AliITSReconstructor::GetRecoParam()->GetTrackleterStoreSPD2SingleCl()); |
| | 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 | // |
| | 213 | } else { |
| | 214 | SetPhiWindow(); |
| | 215 | SetThetaWindow(); |
| | 216 | SetPhiShift(); |
| | 217 | SetRemoveClustersFromOverlaps(); |
| | 218 | SetPhiOverlapCut(); |
| | 219 | SetZetaOverlapCut(); |
| | 220 | SetPhiRotationAngle(); |
| | 221 | |
| | 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(); |
| | 243 | } |
| | 244 | // |
| | 245 | fTracklets = 0; |
| | 246 | fSClusters = 0; |
| | 247 | // |
| | 248 | // definition of histograms |
| | 249 | Bool_t oldStatus = TH1::AddDirectoryStatus(); |
| | 250 | TH1::AddDirectory(kFALSE); |
| | 251 | |
| | 252 | fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,-0.1,0.1); |
| | 253 | fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1); |
| | 254 | |
| | 255 | fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1); |
| | 256 | |
| | 257 | fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5); |
| | 258 | fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,0.0,0.5); |
| | 259 | |
| | 260 | fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,0.,0.5,100,0.,0.5); |
| | 261 | |
| | 262 | fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.); |
| | 263 | fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi()); |
| | 264 | fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.); |
| | 265 | fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi()); |
| | 266 | for (int i=2;i--;) fStoreRefs[i][0] = fStoreRefs[i][1] = kFALSE; |
| | 267 | TH1::AddDirectory(oldStatus); |
| | 268 | } |
| | 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 | AliITSMultReconstructor::~AliITSMultReconstructor(){ |
| | 353 | // Destructor |
| | 354 | |
| | 355 | // delete histograms |
| | 356 | delete fhClustersDPhiAcc; |
| | 357 | delete fhClustersDThetaAcc; |
| | 358 | delete fhClustersDPhiAll; |
| | 359 | delete fhClustersDThetaAll; |
| | 360 | delete fhDPhiVsDThetaAll; |
| | 361 | delete fhDPhiVsDThetaAcc; |
| | 362 | delete fhetaTracklets; |
| | 363 | delete fhphiTracklets; |
| | 364 | delete fhetaClustersLay1; |
| | 365 | delete fhphiClustersLay1; |
| | 366 | // |
| | 367 | // delete arrays |
| | 368 | for(Int_t i=0; i<fNTracklets; i++) delete [] fTracklets[i]; |
| | 369 | |
| | 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]; |
| | 376 | delete[] fClusterCopyIndex[i]; |
| | 377 | delete[] fOverlapFlagClustersLay[i]; |
| | 378 | delete fClArr[i]; |
| | 379 | for (int j=0;j<2;j++) delete fUsedClusLay[i][j]; |
| | 380 | } |
| | 381 | delete [] fTracklets; |
| | 382 | delete [] fSClusters; |
| | 383 | // |
| | 384 | delete[] fPartners; fPartners = 0; |
| | 385 | delete[] fMinDists; fMinDists = 0; |
| | 386 | delete fBlackList; fBlackList = 0; |
| | 387 | // |
| | 388 | } |
| | 389 | |
| | 390 | //____________________________________________________________________ |
| | 391 | void AliITSMultReconstructor::Reconstruct(AliESDEvent* esd, TTree* treeRP) |
| | 392 | { |
| | 393 | if (!treeRP) { AliError(" Invalid ITS cluster tree !\n"); return; } |
| | 394 | if (!esd) {AliError("ESDEvent is not available, use old reconstructor"); return;} |
| | 395 | // reset counters |
| | 396 | if (fMult) delete fMult; fMult = 0; |
| | 397 | fNClustersLay[0] = 0; |
| | 398 | fNClustersLay[1] = 0; |
| | 399 | fNTracklets = 0; |
| | 400 | fNSingleCluster = 0; |
| | 401 | fNSingleClusterSPD2 = 0; |
| | 402 | // |
| | 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(); |
| | 411 | if (!vtx || vtx->GetNContributors()<1) isVtxOK = kFALSE; |
| | 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 | } |
| | 424 | if(vtx){ |
| | 425 | float vtxf[3] = {vtx->GetX(),vtx->GetY(),vtx->GetZ()}; |
| | 426 | FindTracklets(vtxf); |
| | 427 | } |
| | 428 | else { |
| | 429 | FindTracklets(0); |
| | 430 | } |
| | 431 | // |
| | 432 | CreateMultiplicityObject(); |
| | 433 | } |
| | 434 | |
| | 435 | //____________________________________________________________________ |
| | 436 | void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) { |
| | 437 | // |
| | 438 | // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder and in analysis mode |
| | 439 | |
| | 440 | if (fMult) delete fMult; fMult = 0; |
| | 441 | fNClustersLay[0] = 0; |
| | 442 | fNClustersLay[1] = 0; |
| | 443 | fNTracklets = 0; |
| | 444 | fNSingleCluster = 0; |
| | 445 | fNSingleClusterSPD2 = 0; |
| | 446 | // |
| | 447 | if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; } |
| | 448 | // |
| | 449 | fESDEvent = 0; |
| | 450 | SetTreeRP(clusterTree); |
| | 451 | // |
| | 452 | FindTracklets(vtx); |
| | 453 | // |
| | 454 | } |
| | 455 | |
| | 456 | |
| | 457 | //____________________________________________________________________ |
| | 458 | void AliITSMultReconstructor::ReconstructMix(TTree* clusterTree, TTree* clusterTreeMix, const Float_t* vtx, Float_t*) |
| | 459 | { |
| | 460 | // |
| | 461 | // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder and in analysis mode |
| | 462 | |
| | 463 | if (fMult) delete fMult; fMult = 0; |
| | 464 | fNClustersLay[0] = 0; |
| | 465 | fNClustersLay[1] = 0; |
| | 466 | fNTracklets = 0; |
| | 467 | fNSingleCluster = 0; |
| | 468 | fNSingleClusterSPD2 = 0; |
| | 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 | { |
| | 485 | // - calls LoadClusterArrays that finds the position of the clusters |
| | 486 | // (in global coord) |
| | 487 | |
| | 488 | // - convert the cluster coordinates to theta, phi (seen from the |
| | 489 | // interaction vertex). Clusters in the inner layer can be now |
| | 490 | // rotated for combinatorial studies |
| | 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 | |
| | 497 | // Find tracklets converging to vertex |
| | 498 | // |
| | 499 | LoadClusterArrays(fTreeRP,fTreeRPMix); |
| | 500 | // flag clusters used by ESD tracks |
| | 501 | if (fESDEvent) ProcessESDTracks(); |
| | 502 | fRecoDone = kTRUE; |
| | 503 | |
| | 504 | if (!vtx) return; |
| | 505 | |
| | 506 | InitAux(); |
| | 507 | |
| | 508 | // find the tracklets |
| | 509 | AliDebug(1,"Looking for tracklets... "); |
| | 510 | |
| | 511 | ClusterPos2Angles(vtx); // convert cluster position to angles wrt vtx |
| | 512 | // |
| | 513 | // Step1: find all tracklets allowing double assocation: |
| | 514 | int found = 1; |
| | 515 | while (found > 0) { |
| | 516 | found = 0; |
| | 517 | for (Int_t iC1=0; iC1<fNClustersLay[0]; iC1++) found += AssociateClusterOfL1(iC1); |
| | 518 | } |
| | 519 | // |
| | 520 | // Step2: store tracklets; remove used clusters |
| | 521 | for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) StoreTrackletForL2Cluster(iC2); |
| | 522 | // |
| | 523 | // store unused single clusters of L1 (optionally for L2 too) |
| | 524 | StoreL1Singles(); |
| | 525 | // |
| | 526 | AliDebug(1,Form("%d tracklets found", fNTracklets)); |
| | 527 | } |
| | 528 | |
| | 529 | //____________________________________________________________________ |
| | 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); |
| | 541 | fMult->SetMultTrackRefs( fBuildRefs ); |
| | 542 | fMult->SetSPD2SinglesStored(fStoreSPD2SingleCl); |
| | 543 | fMult->SetNumberOfSingleClustersSPD2(fNSingleClusterSPD2); |
| | 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 | // |
| | 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 | // |
| | 556 | UInt_t shared[100]; |
| | 557 | AliRefArray *refs[2][2] = {{0,0},{0,0}}; |
| | 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 | } |
| | 563 | // |
| | 564 | for (int i=fNTracklets;i--;) { |
| | 565 | float* tlInfo = fTracklets[i]; |
| | 566 | fMult->SetTrackletData(i,tlInfo); |
| | 567 | // |
| | 568 | if (!fBuildRefs) continue; // do we need references? |
| | 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 | } |
| | 579 | } |
| | 580 | if (fBuildRefs) fMult->AttachTracklet2TrackRefs(refs[0][0],refs[0][1],refs[1][0],refs[1][1]); |
| | 581 | // |
| | 582 | AliRefArray *refsc[2] = {0,0}; |
| | 583 | if (fBuildRefs) for (int it=2;it--;) if (fStoreRefs[0][it]) refsc[it] = new AliRefArray(fNClustersLay[0]); |
| | 584 | for (int i=fNSingleCluster;i--;) { |
| | 585 | float* clInfo = fSClusters[i]; |
| | 586 | fMult->SetSingleClusterData(i,clInfo); |
| | 587 | // |
| | 588 | if (!fBuildRefs) continue; // do we need references? |
| | 589 | int ilr = i>=(fNSingleCluster-fNSingleClusterSPD2) ? 1:0; |
| | 590 | int clID = int(clInfo[kSCID]); |
| | 591 | for (int itp=0;itp<2;itp++) { |
| | 592 | if (!fStoreRefs[ilr][itp]) continue; |
| | 593 | int nref = fUsedClusLay[ilr][itp]->GetReferences(clID,shared,100); |
| | 594 | if (!nref) continue; |
| | 595 | else if (nref==1) refsc[itp]->AddReference(i,shared[0]); |
| | 596 | else refsc[itp]->AddReferences(i,shared,nref); |
| | 597 | } |
| | 598 | } |
| | 599 | // |
| | 600 | if (fBuildRefs) fMult->AttachCluster2TrackRefs(refsc[0],refsc[1]); |
| | 601 | fMult->CompactBits(); |
| | 602 | // |
| | 603 | } |
| | 604 | |
| | 605 | |
| | 606 | //____________________________________________________________________ |
| | 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; |
| | 620 | int nSlots = GetStoreSPD2SingleCl() ? fNClustersLay[0]+fNClustersLay[1] : fNClustersLay[0]; |
| | 621 | fSClusters = new Float_t*[nSlots]; |
| | 622 | memset(fSClusters,0,nSlots*sizeof(Float_t*)); |
| | 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) |
| | 631 | { |
| | 632 | // This method |
| | 633 | // - gets the clusters from the cluster tree for layer il |
| | 634 | // - convert them into global coordinates |
| | 635 | // - store them in the internal arrays |
| | 636 | // - count the number of cluster-fired chips |
| | 637 | // |
| | 638 | // RS: This method was strongly modified wrt original. In order to have the same numbering |
| | 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 |
| | 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); |
| | 666 | delete[] fClusterCopyIndex[il]; |
| | 667 | } |
| | 668 | // |
| | 669 | // count clusters |
| | 670 | // loop over the SPD subdetectors |
| | 671 | int nclLayer = 0; |
| | 672 | int detMin = TMath::Max(0,AliITSgeomTGeo::GetModuleIndex(il+1,1,1)); |
| | 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; |
| | 683 | if (fCreateClustersCopy) cluster = new ((*fClArr[il])[nclLayer]) AliITSRecPoint(*cluster); |
| | 684 | clArr.AddAtAndExpand(cluster,nclLayer++); |
| | 685 | Int_t chipNo = fSPDSeg.GetChipFromLocal(0,cluster->GetDetLocalZ()); |
| | 686 | if(chipNo>=0)nClustersInChip[ chipNo ]++; |
| | 687 | } |
| | 688 | for(Int_t ifChip=5;ifChip--;) if (nClustersInChip[ifChip]) fNFiredChips[il]++; |
| | 689 | } |
| | 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]; |
| | 698 | if (fCreateClustersCopy) fClusterCopyIndex[il] = new Int_t[nclLayer]; |
| | 699 | // |
| | 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); |
| | 708 | if (fCreateClustersCopy) fClusterCopyIndex[il][ic] = index[ic]; |
| | 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 | // |
| | 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 | } |
| | 725 | } |
| | 726 | // |
| | 727 | if (fClustersLay[il]) delete[] fClustersLay[il]; |
| | 728 | fClustersLay[il] = clustersLay; |
| | 729 | fNClustersLay[il] = nclLayer; |
| | 730 | // |
| | 731 | } |
| | 732 | |
| | 733 | //____________________________________________________________________ |
| | 734 | void AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) { |
| | 735 | // This method |
| | 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 | |
| | 744 | AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance(); |
| | 745 | TClonesArray* itsClusters=NULL; |
| | 746 | rpcont->FetchClusters(0,itsClusterTree); |
| | 747 | if(!rpcont->IsSPDActive()){ |
| | 748 | AliWarning("No SPD rec points found, multiplicity not calculated"); |
| | 749 | return; |
| | 750 | } |
| | 751 | |
| | 752 | // loop over the its subdetectors |
| | 753 | Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1); |
| | 754 | for (Int_t iIts=0; iIts < nSPDmodules; iIts++) { |
| | 755 | itsClusters=rpcont->UncheckedGetClusters(iIts); |
| | 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; |
| | 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; |
| | 766 | |
| | 767 | // loop over clusters |
| | 768 | while(nClusters--) { |
| | 769 | AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters); |
| | 770 | |
| | 771 | // find the chip for the current cluster |
| | 772 | Float_t locz = cluster->GetDetLocalZ(); |
| | 773 | Int_t iChip = fSPDSeg.GetChipFromLocal(0,locz); |
| | 774 | if (iChip>=0) nClustersInChip[iChip]++; |
| | 775 | |
| | 776 | }// end of cluster loop |
| | 777 | |
| | 778 | // get number of fired chips in the current module |
| | 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 | |
| | 785 | |
| | 786 | AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1])); |
| | 787 | } |
| | 788 | //____________________________________________________________________ |
| | 789 | void |
| | 790 | AliITSMultReconstructor::SaveHists() { |
| | 791 | // This method save the histograms on the output file |
| | 792 | // (only if fHistOn is TRUE). |
| | 793 | |
| | 794 | if (!fHistOn) |
| | 795 | return; |
| | 796 | |
| | 797 | fhClustersDPhiAll->Write(); |
| | 798 | fhClustersDThetaAll->Write(); |
| | 799 | fhDPhiVsDThetaAll->Write(); |
| | 800 | |
| | 801 | fhClustersDPhiAcc->Write(); |
| | 802 | fhClustersDThetaAcc->Write(); |
| | 803 | fhDPhiVsDThetaAcc->Write(); |
| | 804 | |
| | 805 | fhetaTracklets->Write(); |
| | 806 | fhphiTracklets->Write(); |
| | 807 | fhetaClustersLay1->Write(); |
| | 808 | fhphiClustersLay1->Write(); |
| | 809 | } |
| | 810 | |
| | 811 | //____________________________________________________________________ |
| | 812 | void AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) |
| | 813 | { |
| | 814 | // Flags clusters in the overlapping regions |
| | 815 | Float_t distClSameMod=0.; |
| | 816 | Float_t distClSameModMin=0.; |
| | 817 | Int_t iClOverlap =0; |
| | 818 | Float_t meanRadiusLay1 = 3.99335; // average radius inner layer |
| | 819 | Float_t meanRadiusLay2 = 7.37935; // average radius outer layer; |
| | 820 | |
| | 821 | Float_t zproj1=0.; |
| | 822 | Float_t zproj2=0.; |
| | 823 | Float_t deZproj=0.; |
| | 824 | Float_t* clPar1 = GetClusterLayer1(iC1); |
| | 825 | Float_t* clPar2B = GetClusterLayer2(iC2WithBestDist); |
| | 826 | // Loop on inner layer clusters |
| | 827 | for (Int_t iiC1=0; iiC1<fNClustersLay[0]; iiC1++) { |
| | 828 | if (!fOverlapFlagClustersLay[0][iiC1]) { |
| | 829 | // only for adjacent modules |
| | 830 | if ((TMath::Abs(fDetectorIndexClustersLay[0][iC1]-fDetectorIndexClustersLay[0][iiC1])==4)|| |
| | 831 | (TMath::Abs(fDetectorIndexClustersLay[0][iC1]-fDetectorIndexClustersLay[0][iiC1])==76)) { |
| | 832 | Float_t *clPar11 = GetClusterLayer1(iiC1); |
| | 833 | Float_t dePhi=TMath::Abs(clPar11[kClPh]-clPar1[kClPh]); |
| | 834 | if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi; |
| | 835 | |
| | 836 | zproj1=meanRadiusLay1/TMath::Tan(clPar1[kClTh]); |
| | 837 | zproj2=meanRadiusLay1/TMath::Tan(clPar11[kClTh]); |
| | 838 | |
| | 839 | deZproj=TMath::Abs(zproj1-zproj2); |
| | 840 | |
| | 841 | distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2)); |
| | 842 | if (distClSameMod<=1.) fOverlapFlagClustersLay[0][iiC1]=kTRUE; |
| | 843 | |
| | 844 | // if (distClSameMod<=1.) { |
| | 845 | // if (distClSameModMin==0. || distClSameMod<distClSameModMin) { |
| | 846 | // distClSameModMin=distClSameMod; |
| | 847 | // iClOverlap=iiC1; |
| | 848 | // } |
| | 849 | // } |
| | 850 | |
| | 851 | |
| | 852 | } // end adjacent modules |
| | 853 | } |
| | 854 | } // end Loop on inner layer clusters |
| | 855 | |
| | 856 | // if (distClSameModMin!=0.) fOverlapFlagClustersLay[0][iClOverlap]=kTRUE; |
| | 857 | |
| | 858 | distClSameMod=0.; |
| | 859 | distClSameModMin=0.; |
| | 860 | iClOverlap =0; |
| | 861 | // Loop on outer layer clusters |
| | 862 | for (Int_t iiC2=0; iiC2<fNClustersLay[1]; iiC2++) { |
| | 863 | if (!fOverlapFlagClustersLay[1][iiC2]) { |
| | 864 | // only for adjacent modules |
| | 865 | Float_t *clPar2 = GetClusterLayer2(iiC2); |
| | 866 | if ((TMath::Abs(fDetectorIndexClustersLay[1][iC2WithBestDist]-fDetectorIndexClustersLay[1][iiC2])==4) || |
| | 867 | (TMath::Abs(fDetectorIndexClustersLay[1][iC2WithBestDist]-fDetectorIndexClustersLay[1][iiC2])==156)) { |
| | 868 | Float_t dePhi=TMath::Abs(clPar2[kClPh]-clPar2B[kClPh]); |
| | 869 | if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi; |
| | 870 | |
| | 871 | zproj1=meanRadiusLay2/TMath::Tan(clPar2B[kClTh]); |
| | 872 | zproj2=meanRadiusLay2/TMath::Tan(clPar2[kClTh]); |
| | 873 | |
| | 874 | deZproj=TMath::Abs(zproj1-zproj2); |
| | 875 | distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2)); |
| | 876 | if (distClSameMod<=1.) fOverlapFlagClustersLay[1][iiC2]=kTRUE; |
| | 877 | |
| | 878 | // if (distClSameMod<=1.) { |
| | 879 | // if (distClSameModMin==0. || distClSameMod<distClSameModMin) { |
| | 880 | // distClSameModMin=distClSameMod; |
| | 881 | // iClOverlap=iiC2; |
| | 882 | // } |
| | 883 | // } |
| | 884 | |
| | 885 | } // end adjacent modules |
| | 886 | } |
| | 887 | } // end Loop on outer layer clusters |
| | 888 | |
| | 889 | // if (distClSameModMin!=0.) fOverlapFlagClustersLay[1][iClOverlap]=kTRUE; |
| | 890 | |
| | 891 | } |
| | 892 | |
| | 893 | //____________________________________________________________________ |
| | 894 | void AliITSMultReconstructor::InitAux() |
| | 895 | { |
| | 896 | // init arrays/parameters for tracklet reconstruction |
| | 897 | |
| | 898 | // dPhi shift is field dependent, get average magnetic field |
| | 899 | Float_t bz = 0; |
| | 900 | AliMagF* field = 0; |
| | 901 | if (TGeoGlobalMagField::Instance()) field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField()); |
| | 902 | if (!field) { |
| | 903 | AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor."); |
| | 904 | } |
| | 905 | else bz = TMath::Abs(field->SolenoidField()); |
| | 906 | fDPhiShift = fPhiShift / 5 * bz; |
| | 907 | AliDebug(1, Form("Using phi shift of %f", fDPhiShift)); |
| | 908 | // |
| | 909 | if (fPartners) delete[] fPartners; fPartners = new Int_t[fNClustersLay[1]]; |
| | 910 | if (fMinDists) delete[] fMinDists; fMinDists = new Float_t[fNClustersLay[1]]; |
| | 911 | if (fAssociatedLay1) delete[] fAssociatedLay1; fAssociatedLay1 = new Int_t[fNClustersLay[0]]; |
| | 912 | // |
| | 913 | if (fBlackList) delete fBlackList; fBlackList = new AliRefArray(fNClustersLay[0]); |
| | 914 | // |
| | 915 | // Printf("Vertex in find tracklets...%f %f %f",vtx[0],vtx[1],vtx[2]); |
| | 916 | for (Int_t i=0; i<fNClustersLay[1]; i++) { |
| | 917 | fPartners[i] = -1; |
| | 918 | fMinDists[i] = 2*fNStdDev; |
| | 919 | } |
| | 920 | memset(fAssociatedLay1,0,fNClustersLay[0]*sizeof(Int_t)); |
| | 921 | // |
| | 922 | } |
| | 923 | |
| | 924 | //____________________________________________________________________ |
| | 925 | void AliITSMultReconstructor::ClusterPos2Angles(const Float_t *vtx) |
| | 926 | { |
| | 927 | // convert cluster coordinates to angles wrt vertex |
| | 928 | for (int ilr=0;ilr<2;ilr++) { |
| | 929 | for (Int_t iC=0; iC<fNClustersLay[ilr]; iC++) { |
| | 930 | float* clPar = GetClusterOfLayer(ilr,iC); |
| | 931 | CalcThetaPhi(clPar[kClTh]-vtx[0],clPar[kClPh]-vtx[1],clPar[kClZ]-vtx[2],clPar[kClTh],clPar[kClPh]); |
| | 932 | if (ilr==0) { |
| | 933 | clPar[kClPh] = clPar[kClPh] + fPhiRotationAngle; // rotation of inner layer for comb studies |
| | 934 | if (fHistOn) { |
| | 935 | Float_t eta = clPar[kClTh]; |
| | 936 | eta= TMath::Tan(eta/2.); |
| | 937 | eta=-TMath::Log(eta); |
| | 938 | fhetaClustersLay1->Fill(eta); |
| | 939 | fhphiClustersLay1->Fill(clPar[kClPh]); |
| | 940 | } |
| | 941 | } |
| | 942 | } |
| | 943 | } |
| | 944 | // |
| | 945 | } |
| | 946 | |
| | 947 | //____________________________________________________________________ |
| | 948 | Int_t AliITSMultReconstructor::AssociateClusterOfL1(Int_t iC1) |
| | 949 | { |
| | 950 | // search association of cluster iC1 of L1 with all clusters of L2 |
| | 951 | if (fAssociatedLay1[iC1] != 0) return 0; |
| | 952 | Int_t iC2WithBestDist = -1; // reset |
| | 953 | Double_t minDist = 2*fNStdDev; // reset |
| | 954 | float* clPar1 = GetClusterLayer1(iC1); |
| | 955 | for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) { |
| | 956 | // |
| | 957 | if (fBlackList->IsReferred(iC1,iC2)) continue; |
| | 958 | float* clPar2 = GetClusterLayer2(iC2); |
| | 959 | // |
| | 960 | // find the difference in angles |
| | 961 | Double_t dTheta = TMath::Abs(clPar2[kClTh] - clPar1[kClTh]); |
| | 962 | Double_t dPhi = TMath::Abs(clPar2[kClPh] - clPar1[kClPh]); |
| | 963 | // Printf("detheta %f dephi %f", dTheta,dPhi); |
| | 964 | // |
| | 965 | if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi; // take into account boundary condition |
| | 966 | // |
| | 967 | if (fHistOn) { |
| | 968 | fhClustersDPhiAll->Fill(dPhi); |
| | 969 | fhClustersDThetaAll->Fill(dTheta); |
| | 970 | fhDPhiVsDThetaAll->Fill(dTheta, dPhi); |
| | 971 | } |
| | 972 | Float_t d = CalcDist(dPhi,dTheta,clPar1[kClTh]); // make "elliptical" cut in Phi and Theta! |
| | 973 | // look for the minimum distance: the minimum is in iC2WithBestDist |
| | 974 | if (d<fNStdDev && d<minDist) { minDist=d; iC2WithBestDist = iC2; } |
| | 975 | } |
| | 976 | // |
| | 977 | if (minDist<fNStdDev) { // This means that a cluster in layer 2 was found that matches with iC1 |
| | 978 | // |
| | 979 | if (fMinDists[iC2WithBestDist] > minDist) { |
| | 980 | Int_t oldPartner = fPartners[iC2WithBestDist]; |
| | 981 | fPartners[iC2WithBestDist] = iC1; |
| | 982 | fMinDists[iC2WithBestDist] = minDist; |
| | 983 | // |
| | 984 | fAssociatedLay1[iC1] = 1; // mark as assigned |
| | 985 | // |
| | 986 | if (oldPartner != -1) { |
| | 987 | // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its fBlackList |
| | 988 | fBlackList->AddReference(oldPartner,iC2WithBestDist); |
| | 989 | fAssociatedLay1[oldPartner] = 0; // mark as free |
| | 990 | } |
| | 991 | } else { |
| | 992 | // try again to find a cluster without considering iC2WithBestDist |
| | 993 | fBlackList->AddReference(iC1,iC2WithBestDist); |
| | 994 | } |
| | 995 | // |
| | 996 | } |
| | 997 | else fAssociatedLay1[iC1] = 2;// cluster has no partner; remove |
| | 998 | // |
| | 999 | return 1; |
| | 1000 | } |
| | 1001 | |
| | 1002 | //____________________________________________________________________ |
| | 1003 | Int_t AliITSMultReconstructor::StoreTrackletForL2Cluster(Int_t iC2) |
| | 1004 | { |
| | 1005 | // build tracklet for cluster iC2 of layer 2 |
| | 1006 | if (fPartners[iC2] == -1) return 0; |
| | 1007 | if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (fPartners[iC2],iC2); |
| | 1008 | // Printf("saving tracklets"); |
| | 1009 | if (fOverlapFlagClustersLay[0][fPartners[iC2]] || fOverlapFlagClustersLay[1][iC2]) return 0; |
| | 1010 | float* clPar2 = GetClusterLayer2(iC2); |
| | 1011 | float* clPar1 = GetClusterLayer1(fPartners[iC2]); |
| | 1012 | // |
| | 1013 | Float_t* tracklet = fTracklets[fNTracklets] = new Float_t[kTrNPar]; // RS Add also the cluster id's |
| | 1014 | // |
| | 1015 | tracklet[kTrTheta] = clPar1[kClTh]; // use the theta from the clusters in the first layer |
| | 1016 | tracklet[kTrPhi] = clPar1[kClPh]; // use the phi from the clusters in the first layer |
| | 1017 | tracklet[kTrDPhi] = clPar1[kClPh] - clPar2[kClPh]; // store the difference between phi1 and phi2 |
| | 1018 | // |
| | 1019 | // define dphi in the range [0,pi] with proper sign (track charge correlated) |
| | 1020 | if (tracklet[kTrDPhi] > TMath::Pi()) tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi(); |
| | 1021 | if (tracklet[kTrDPhi] < -TMath::Pi()) tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi(); |
| | 1022 | // |
| | 1023 | tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh]; // store the theta1-theta2 |
| | 1024 | // |
| | 1025 | if (fHistOn) { |
| | 1026 | fhClustersDPhiAcc->Fill(tracklet[kTrDPhi]); |
| | 1027 | fhClustersDThetaAcc->Fill(tracklet[kTrDTheta]); |
| | 1028 | fhDPhiVsDThetaAcc->Fill(tracklet[kTrDTheta],tracklet[kTrDPhi]); |
| | 1029 | } |
| | 1030 | // |
| | 1031 | // find label |
| | 1032 | // if equal label in both clusters found this label is assigned |
| | 1033 | // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned |
| | 1034 | Int_t label1=0,label2=0; |
| | 1035 | while (label2 < 3) { |
| | 1036 | if ( int(clPar1[kClMC0+label1])!=-2 && int(clPar1[kClMC0+label1])==int(clPar2[kClMC0+label2])) break; |
| | 1037 | if (++label1 == 3) { label1 = 0; label2++; } |
| | 1038 | } |
| | 1039 | if (label2 < 3) { |
| | 1040 | AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", |
| | 1041 | (Int_t) clPar1[kClMC0+label1], (Int_t) clPar1[kClMC0+label2], fNTracklets)); |
| | 1042 | tracklet[kTrLab1] = tracklet[kTrLab2] = clPar1[kClMC0+label1]; |
| | 1043 | } else { |
| | 1044 | AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", |
| | 1045 | (Int_t) clPar1[kClMC0], (Int_t) clPar1[kClMC1], (Int_t) clPar1[kClMC2], |
| | 1046 | (Int_t) clPar2[kClMC0], (Int_t) clPar2[kClMC1], (Int_t) clPar2[kClMC2], fNTracklets)); |
| | 1047 | tracklet[kTrLab1] = clPar1[kClMC0]; |
| | 1048 | tracklet[kTrLab2] = clPar2[kClMC0]; |
| | 1049 | } |
| | 1050 | // |
| | 1051 | if (fHistOn) { |
| | 1052 | Float_t eta = tracklet[kTrTheta]; |
| | 1053 | eta= TMath::Tan(eta/2.); |
| | 1054 | eta=-TMath::Log(eta); |
| | 1055 | fhetaTracklets->Fill(eta); |
| | 1056 | fhphiTracklets->Fill(tracklet[kTrPhi]); |
| | 1057 | } |
| | 1058 | // |
| | 1059 | tracklet[kClID1] = fPartners[iC2]; |
| | 1060 | tracklet[kClID2] = iC2; |
| | 1061 | // |
| | 1062 | // Printf("Adding tracklet candidate"); |
| | 1063 | AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets)); |
| | 1064 | AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", fPartners[iC2], iC2)); |
| | 1065 | fNTracklets++; |
| | 1066 | fAssociatedLay1[fPartners[iC2]] = 1; |
| | 1067 | // |
| | 1068 | return 1; |
| | 1069 | } |
| | 1070 | |
| | 1071 | //____________________________________________________________________ |
| | 1072 | void AliITSMultReconstructor::StoreL1Singles() |
| | 1073 | { |
| | 1074 | // Printf("saving single clusters..."); |
| | 1075 | for (Int_t iC1=0; iC1<fNClustersLay[0]; iC1++) { |
| | 1076 | float* clPar1 = GetClusterLayer1(iC1); |
| | 1077 | if (fAssociatedLay1[iC1]==2||fAssociatedLay1[iC1]==0) { |
| | 1078 | fSClusters[fNSingleCluster] = new Float_t[kClNPar]; |
| | 1079 | fSClusters[fNSingleCluster][kSCTh] = clPar1[kClTh]; |
| | 1080 | fSClusters[fNSingleCluster][kSCPh] = clPar1[kClPh]; |
| | 1081 | fSClusters[fNSingleCluster][kSCLab] = clPar1[kClMC0]; |
| | 1082 | fSClusters[fNSingleCluster][kSCID] = iC1; |
| | 1083 | AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)", |
| | 1084 | fNSingleCluster, iC1)); |
| | 1085 | fNSingleCluster++; |
| | 1086 | } |
| | 1087 | } |
| | 1088 | // |
| | 1089 | if (GetStoreSPD2SingleCl()) { |
| | 1090 | for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) { |
| | 1091 | if (fPartners[iC2]<0 || (fOverlapFlagClustersLay[0][fPartners[iC2]] || fOverlapFlagClustersLay[1][iC2])) { |
| | 1092 | float* clPar2 = GetClusterLayer2(iC2); |
| | 1093 | fSClusters[fNSingleCluster] = new Float_t[kClNPar]; |
| | 1094 | fSClusters[fNSingleCluster][kSCTh] = clPar2[kClTh]; |
| | 1095 | fSClusters[fNSingleCluster][kSCPh] = clPar2[kClPh]; |
| | 1096 | fSClusters[fNSingleCluster][kSCLab] = clPar2[kClMC0]; |
| | 1097 | fSClusters[fNSingleCluster][kSCID] = iC2; |
| | 1098 | AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 2)", |
| | 1099 | fNSingleCluster, iC2)); |
| | 1100 | fNSingleCluster++; |
| | 1101 | fNSingleClusterSPD2++; |
| | 1102 | } |
| | 1103 | } |
| | 1104 | } |
| | 1105 | // |
| | 1106 | } |
| | 1107 | |
| | 1108 | //____________________________________________________________________ |
| | 1109 | void AliITSMultReconstructor::ProcessESDTracks() |
| | 1110 | { |
| | 1111 | // Flag the clusters used by ESD tracks |
| | 1112 | // Flag primary tracks to be used for multiplicity counting |
| | 1113 | // |
| | 1114 | if (!fESDEvent || !fBuildRefs) return; |
| | 1115 | AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexTracks(); |
| | 1116 | if (!vtx || vtx->GetNContributors()<1) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD(); |
| | 1117 | if (!vtx || vtx->GetNContributors()<1) { |
| | 1118 | AliDebug(1,"No primary vertex: cannot flag primary tracks"); |
| | 1119 | return; |
| | 1120 | } |
| | 1121 | Int_t ntracks = fESDEvent->GetNumberOfTracks(); |
| | 1122 | for(Int_t itr=0; itr<ntracks; itr++) { |
| | 1123 | AliESDtrack* track = fESDEvent->GetTrack(itr); |
| | 1124 | if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx |
| | 1125 | FlagTrackClusters(itr); |
| | 1126 | FlagIfSecondary(track,vtx); |
| | 1127 | } |
| | 1128 | FlagV0s(vtx); |
| | 1129 | // |
| | 1130 | } |
| | 1131 | |
| | 1132 | //____________________________________________________________________ |
| | 1133 | void AliITSMultReconstructor::FlagTrackClusters(Int_t id) |
| | 1134 | { |
| | 1135 | // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation |
| | 1136 | // |
| | 1137 | const AliESDtrack* track = fESDEvent->GetTrack(id); |
| | 1138 | Int_t idx[12]; |
| | 1139 | if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit |
| | 1140 | Int_t itsType = track->IsOn(AliESDtrack::kITSpureSA) ? 1:0; |
| | 1141 | |
| | 1142 | for (int i=6/*AliESDfriendTrack::kMaxITScluster*/;i--;) { // ignore extras: note: i>=6 is for extra clusters |
| | 1143 | if (idx[i]<0) continue; |
| | 1144 | int layID= (idx[i] & 0xf0000000) >> 28; |
| | 1145 | if (layID>1) continue; // SPD only |
| | 1146 | int clID = (idx[i] & 0x0fffffff); |
| | 1147 | fUsedClusLay[layID][itsType]->AddReference(clID,id); |
| | 1148 | fStoreRefs[layID][itsType] = kTRUE; |
| | 1149 | } |
| | 1150 | // |
| | 1151 | } |
| | 1152 | |
| | 1153 | //____________________________________________________________________ |
| | 1154 | void AliITSMultReconstructor::FlagIfSecondary(AliESDtrack* track, const AliVertex* vtx) |
| | 1155 | { |
| | 1156 | // RS: check if the track is primary and set the flag |
| | 1157 | double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? fCutPxDrSPDin:fCutPxDrSPDout; |
| | 1158 | float xz[2]; |
| | 1159 | track->GetDZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), fESDEvent->GetMagneticField(), xz); |
| | 1160 | if (TMath::Abs(xz[0]*track->P())>cut || TMath::Abs(xz[1]*track->P())>fCutPxDz || |
| | 1161 | TMath::Abs(xz[0])>fCutDCArz || TMath::Abs(xz[1])>fCutDCArz) |
| | 1162 | track->SetStatus(AliESDtrack::kMultSec); |
| | 1163 | else track->ResetStatus(AliESDtrack::kMultSec); |
| | 1164 | } |
| | 1165 | |
| | 1166 | //____________________________________________________________________ |
| | 1167 | void AliITSMultReconstructor::FlagV0s(const AliESDVertex *vtx) |
| | 1168 | { |
| | 1169 | // flag tracks belonging to v0s |
| | 1170 | // |
| | 1171 | const double kK0Mass = 0.4976; |
| | 1172 | // |
| | 1173 | AliV0 pvertex; |
| | 1174 | AliKFVertex vertexKF; |
| | 1175 | AliKFParticle epKF0,epKF1,pipmKF0,piKF0,piKF1,gammaKF,k0KF; |
| | 1176 | Double_t mass,massErr,chi2c; |
| | 1177 | enum {kKFIni=BIT(14)}; |
| | 1178 | // |
| | 1179 | double recVtx[3]; |
| | 1180 | float recVtxF[3]; |
| | 1181 | vtx->GetXYZ(recVtx); |
| | 1182 | for (int i=3;i--;) recVtxF[i] = recVtx[i]; |
| | 1183 | // |
| | 1184 | int ntracks = fESDEvent->GetNumberOfTracks(); |
| | 1185 | if (ntracks<2) return; |
| | 1186 | // |
| | 1187 | vertexKF.X() = recVtx[0]; |
| | 1188 | vertexKF.Y() = recVtx[1]; |
| | 1189 | vertexKF.Z() = recVtx[2]; |
| | 1190 | vertexKF.Covariance(0,0) = vtx->GetXRes()*vtx->GetXRes(); |
| | 1191 | vertexKF.Covariance(1,2) = vtx->GetYRes()*vtx->GetYRes(); |
| | 1192 | vertexKF.Covariance(2,2) = vtx->GetZRes()*vtx->GetZRes(); |
| | 1193 | // |
| | 1194 | AliESDtrack *trc0,*trc1; |
| | 1195 | for (int it0=0;it0<ntracks;it0++) { |
| | 1196 | trc0 = fESDEvent->GetTrack(it0); |
| | 1197 | if (trc0->IsOn(AliESDtrack::kMultInV0)) continue; |
| | 1198 | if (!trc0->IsOn(AliESDtrack::kITSin)) continue; |
| | 1199 | Bool_t isSAP = trc0->IsPureITSStandalone(); |
| | 1200 | Int_t q0 = trc0->Charge(); |
| | 1201 | Bool_t testGamma = CanBeElectron(trc0); |
| | 1202 | epKF0.ResetBit(kKFIni); |
| | 1203 | piKF0.ResetBit(kKFIni); |
| | 1204 | double bestChi2=1e16; |
| | 1205 | int bestID = -1; |
| | 1206 | // |
| | 1207 | for (int it1=it0+1;it1<ntracks;it1++) { |
| | 1208 | trc1 = fESDEvent->GetTrack(it1); |
| | 1209 | if (trc1->IsOn(AliESDtrack::kMultInV0)) continue; |
| | 1210 | if (!trc1->IsOn(AliESDtrack::kITSin)) continue; |
| | 1211 | if (trc1->IsPureITSStandalone() != isSAP) continue; // pair separately ITS_SA_Pure tracks and TPC/ITS+ITS_SA |
| | 1212 | if ( (q0+trc1->Charge())!=0 ) continue; // don't pair like signs |
| | 1213 | // |
| | 1214 | pvertex.SetParamN(q0<0 ? *trc0:*trc1); |
| | 1215 | pvertex.SetParamP(q0>0 ? *trc0:*trc1); |
| | 1216 | pvertex.Update(recVtxF); |
| | 1217 | if (pvertex.P()<fCutMinP) continue; |
| | 1218 | if (pvertex.GetV0CosineOfPointingAngle()<fCutMinPointAngle) continue; |
| | 1219 | if (pvertex.GetDcaV0Daughters()>fCutMaxDCADauther) continue; |
| | 1220 | double d = pvertex.GetD(recVtx[0],recVtx[1],recVtx[2]); |
| | 1221 | if (d>fCutMaxDCA) continue; |
| | 1222 | double dx=recVtx[0]-pvertex.Xv(), dy=recVtx[1]-pvertex.Yv(); |
| | 1223 | double rv = TMath::Sqrt(dx*dx+dy*dy); |
| | 1224 | // |
| | 1225 | // check gamma conversion hypothesis ----------------------------------------------------------->>> |
| | 1226 | Bool_t gammaOK = kFALSE; |
| | 1227 | while (testGamma && CanBeElectron(trc1)) { |
| | 1228 | if (rv<fCutMinRGamma) break; |
| | 1229 | if (!epKF0.TestBit(kKFIni)) { |
| | 1230 | new(&epKF0) AliKFParticle(*trc0,q0>0 ? kPositron:kElectron); |
| | 1231 | epKF0.SetBit(kKFIni); |
| | 1232 | } |
| | 1233 | new(&epKF1) AliKFParticle(*trc1,q0<0 ? kPositron:kElectron); |
| | 1234 | gammaKF.Initialize(); |
| | 1235 | gammaKF += epKF0; |
| | 1236 | gammaKF += epKF1; |
| | 1237 | gammaKF.SetProductionVertex(vertexKF); |
| | 1238 | gammaKF.GetMass(mass,massErr); |
| | 1239 | if (mass>fCutMassGamma || (massErr>0&&(mass>massErr*fCutMassGammaNSigma))) break; |
| | 1240 | if (gammaKF.GetS()<fCutGammaSFromDecay) break; |
| | 1241 | gammaKF.SetMassConstraint(0.,0.001); |
| | 1242 | chi2c = (gammaKF.GetNDF()!=0) ? gammaKF.GetChi2()/gammaKF.GetNDF() : 1000; |
| | 1243 | if (chi2c>fCutChi2cGamma) break; |
| | 1244 | gammaOK = kTRUE; |
| | 1245 | if (chi2c>bestChi2) break; |
| | 1246 | bestChi2 = chi2c; |
| | 1247 | bestID = it1; |
| | 1248 | break; |
| | 1249 | } |
| | 1250 | if (gammaOK) continue; |
| | 1251 | // check gamma conversion hypothesis -----------------------------------------------------------<<< |
| | 1252 | // check K0 conversion hypothesis ----------------------------------------------------------->>> |
| | 1253 | while (1) { |
| | 1254 | if (rv<fCutMinRK0) break; |
| | 1255 | if (!piKF0.TestBit(kKFIni)) { |
| | 1256 | new(&piKF0) AliKFParticle(*trc0,q0>0 ? kPiPlus:kPiMinus); |
| | 1257 | piKF0.SetBit(kKFIni); |
| | 1258 | } |
| | 1259 | new(&piKF1) AliKFParticle(*trc1,q0<0 ? kPiPlus:kPiMinus); |
| | 1260 | k0KF.Initialize(); |
| | 1261 | k0KF += piKF0; |
| | 1262 | k0KF += piKF1; |
| | 1263 | k0KF.SetProductionVertex(vertexKF); |
| | 1264 | k0KF.GetMass(mass,massErr); |
| | 1265 | mass -= kK0Mass; |
| | 1266 | if (TMath::Abs(mass)>fCutMassK0 || (massErr>0&&(abs(mass)>massErr*fCutMassK0NSigma))) break; |
| | 1267 | if (k0KF.GetS()<fCutK0SFromDecay) break; |
| | 1268 | k0KF.SetMassConstraint(kK0Mass,0.001); |
| | 1269 | chi2c = (k0KF.GetNDF()!=0) ? k0KF.GetChi2()/k0KF.GetNDF() : 1000; |
| | 1270 | if (chi2c>fCutChi2cK0) break; |
| | 1271 | if (chi2c>bestChi2) break; |
| | 1272 | bestChi2 = chi2c; |
| | 1273 | bestID = it1; |
| | 1274 | break; |
| | 1275 | } |
| | 1276 | // check K0 conversion hypothesis -----------------------------------------------------------<<< |
| | 1277 | } |
| | 1278 | // |
| | 1279 | if (bestID>=0) { |
| | 1280 | trc0->SetStatus(AliESDtrack::kMultInV0); |
| | 1281 | fESDEvent->GetTrack(bestID)->SetStatus(AliESDtrack::kMultInV0); |
| | 1282 | } |
| | 1283 | } |
| | 1284 | // |
| | 1285 | } |
| | 1286 | |
| | 1287 | //____________________________________________________________________ |
| | 1288 | Bool_t AliITSMultReconstructor::CanBeElectron(const AliESDtrack* trc) const |
| | 1289 | { |
| | 1290 | // check if the track can be electron |
| | 1291 | Double_t pid[AliPID::kSPECIES]; |
| | 1292 | if (!trc->IsOn(AliESDtrack::kESDpid)) return kTRUE; |
| | 1293 | trc->GetESDpid(pid); |
| | 1294 | return (trc->IsOn(AliESDtrack::kTPCpid)) ? |
| | 1295 | pid[AliPID::kElectron]>fCutMinElectronProbTPC : |
| | 1296 | pid[AliPID::kElectron]>fCutMinElectronProbESD; |
| | 1297 | // |
| | 1298 | } |
| | 1299 | |
| | 1300 | //____________________________________________________________________ |
| | 1301 | AliITSRecPoint* AliITSMultReconstructor::GetRecPoint(Int_t lr, Int_t n) const |
| | 1302 | { |
| | 1303 | // return a cluster of lr corresponding to orderer cluster index n |
| | 1304 | if (fClArr[lr] && fClusterCopyIndex[lr] && n<fNClustersLay[lr]) |
| | 1305 | return (AliITSRecPoint*) fClArr[lr]->At(fClusterCopyIndex[lr][n]); |
| | 1306 | else { |
| | 1307 | AliError("To access the clusters SetCreateClustersCopy should have been called"); |
| | 1308 | return 0; |
| | 1309 | } |
| | 1310 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental