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