| a9e2aefa |
1 | /************************************************************************** |
| 2 | * Copyright(c) 1998-1999, 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 | $Log$ |
| 9cbdf048 |
18 | Revision 1.3 2000/06/16 07:27:08 gosset |
| 19 | To remove problem in running RuleChecker, like in MUON-dev |
| 20 | |
| 79e1e601 |
21 | Revision 1.1.2.5 2000/06/16 07:00:26 gosset |
| 22 | To remove problem in running RuleChecker |
| 23 | |
| a9e2aefa |
24 | Revision 1.1.2.4 2000/06/12 08:00:07 morsch |
| 25 | Dummy streamer to solve CINT compilation problem (to be investigated !) |
| 26 | |
| 27 | Revision 1.1.2.3 2000/06/09 20:59:57 morsch |
| 28 | Make includes consistent with new file structure. |
| 29 | |
| 30 | Revision 1.1.2.2 2000/06/09 12:58:05 gosset |
| 31 | Removed comment beginnings in Log sections of .cxx files |
| 32 | Suppressed most violations of coding rules |
| 33 | |
| 34 | Revision 1.1.2.1 2000/06/07 14:44:53 gosset |
| 35 | Addition of files for track reconstruction in C++ |
| 36 | */ |
| 37 | |
| 38 | //__________________________________________________________________________ |
| 39 | // |
| 40 | // MUON event reconstructor in ALICE |
| 41 | // |
| 42 | // This class contains as data: |
| 43 | // * the parameters for the event reconstruction |
| 44 | // * a pointer to the array of hits to be reconstructed (the event) |
| 45 | // * a pointer to the array of segments made with these hits inside each station |
| 46 | // * a pointer to the array of reconstructed tracks |
| 47 | // |
| 48 | // It contains as methods, among others: |
| 49 | // * MakeEventToBeReconstructed to build the array of hits to be reconstructed |
| 50 | // * MakeSegments to build the segments |
| 51 | // * MakeTracks to build the tracks |
| 52 | //__________________________________________________________________________ |
| 53 | |
| 54 | #include <iostream.h> |
| 55 | |
| 56 | #include <TRandom.h> |
| 57 | #include <TFile.h> |
| 58 | |
| 59 | #include "AliCallf77.h" |
| 60 | #include "AliMUONEventReconstructor.h" |
| 61 | #include "AliMUON.h" |
| 62 | #include "AliMUONHitForRec.h" |
| 63 | #include "AliMUONSegment.h" |
| 64 | #include "AliMUONHit.h" |
| 65 | #include "AliMUONRawCluster.h" |
| 66 | #include "AliMUONTrack.h" |
| 67 | #include "AliMUONChamber.h" |
| 68 | #include "AliMUONTrackHit.h" |
| 69 | #include "AliRun.h" |
| 70 | |
| 71 | #ifndef WIN32 |
| 72 | # define initfield initfield_ |
| 73 | # define reco_gufld reco_gufld_ |
| 74 | #else |
| 75 | # define initfield INITFIELD |
| 76 | # define reco_gufld RECO_GUFLD |
| 77 | #endif |
| 78 | |
| 79 | extern "C" |
| 80 | { |
| 81 | void type_of_call initfield(); |
| 82 | void type_of_call reco_gufld(Double_t *Coor, Double_t *Field); |
| 83 | } |
| 84 | |
| 85 | //************* Defaults parameters for reconstruction |
| 9cbdf048 |
86 | static const Double_t kDefaultMinBendingMomentum = 3.0; |
| 87 | static const Double_t kDefaultMaxSigma2Distance = 16.0; |
| 88 | static const Double_t kDefaultBendingResolution = 0.01; |
| 89 | static const Double_t kDefaultNonBendingResolution = 0.144; |
| 90 | static const Double_t kDefaultChamberThicknessInX0 = 0.03; |
| a9e2aefa |
91 | // Simple magnetic field: |
| 92 | // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG |
| 93 | // Length and Position from reco_muon.F, with opposite sign: |
| 94 | // Length = ZMAGEND-ZCOIL |
| 95 | // Position = (ZMAGEND+ZCOIL)/2 |
| 96 | // to be ajusted differently from real magnetic field ???? |
| 9cbdf048 |
97 | static const Double_t kDefaultSimpleBValue = 7.0; |
| 98 | static const Double_t kDefaultSimpleBLength = 428.0; |
| 99 | static const Double_t kDefaultSimpleBPosition = 1019.0; |
| 100 | static const Int_t kDefaultRecGeantHits = 0; |
| 101 | static const Double_t kDefaultEfficiency = 0.95; |
| a9e2aefa |
102 | |
| 9cbdf048 |
103 | static const Int_t kDefaultPrintLevel = 0; |
| a9e2aefa |
104 | |
| 105 | ClassImp(AliMUONEventReconstructor) // Class implementation in ROOT context |
| 106 | |
| 107 | //__________________________________________________________________________ |
| 108 | AliMUONEventReconstructor::AliMUONEventReconstructor(void) |
| 109 | { |
| 110 | // Constructor for class AliMUONEventReconstructor |
| 111 | SetReconstructionParametersToDefaults(); |
| 112 | // Memory allocation for the TClonesArray of hits for reconstruction |
| 113 | // Is 10000 the right size ???? |
| 114 | fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000); |
| 115 | fNHitsForRec = 0; // really needed or GetEntriesFast sufficient ???? |
| 116 | // Memory allocation for the TClonesArray's of segments in stations |
| 117 | // Is 2000 the right size ???? |
| 9cbdf048 |
118 | for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) { |
| a9e2aefa |
119 | fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000); |
| 120 | fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ???? |
| 121 | } |
| 122 | // Memory allocation for the TClonesArray of reconstructed tracks |
| 123 | // Is 10 the right size ???? |
| 124 | fRecTracksPtr = new TClonesArray("AliMUONTrack", 10); |
| 125 | fNRecTracks = 0; // really needed or GetEntriesFast sufficient ???? |
| 126 | |
| 127 | // Initialize magnetic field |
| 128 | // using Fortran subroutine INITFIELD in "reco_muon.F". |
| 129 | // Should rather use AliMagF ???? and remove prototyping ... |
| 130 | initfield(); |
| 131 | // Impression de quelques valeurs |
| 132 | Double_t coor[3], field[3]; |
| 133 | coor[0] = 50.0; |
| 134 | coor[1] = 50.0; |
| 135 | coor[2] = 950.0; |
| 136 | reco_gufld(coor, field); |
| 137 | cout << "coor: " << coor[0] << ", " << coor[1] << ", " << coor[2] << endl; |
| 138 | cout << "field: " << field[0] << ", " << field[1] << ", " << field[2] << endl; |
| 139 | coor[2] = -950.0; |
| 140 | reco_gufld(coor, field); |
| 141 | cout << "coor: " << coor[0] << ", " << coor[1] << ", " << coor[2] << endl; |
| 142 | cout << "field: " << field[0] << ", " << field[1] << ", " << field[2] << endl; |
| 143 | coor[2] = -950.0; |
| 144 | |
| 145 | if (fPrintLevel >= 0) { |
| 146 | cout << "AliMUONEventReconstructor constructed with defaults" << endl; Dump();} |
| 147 | return; |
| 148 | } |
| 149 | |
| 9cbdf048 |
150 | AliMUONEventReconstructor::AliMUONEventReconstructor (const AliMUONEventReconstructor& Reconstructor) |
| 151 | { |
| 152 | // Dummy copy constructor |
| 153 | } |
| 154 | |
| 155 | AliMUONEventReconstructor & AliMUONEventReconstructor::operator=(const AliMUONEventReconstructor& Reconstructor) |
| 156 | { |
| 157 | // Dummy assignment operator |
| 158 | return *this; |
| 159 | } |
| 160 | |
| a9e2aefa |
161 | //__________________________________________________________________________ |
| 162 | AliMUONEventReconstructor::~AliMUONEventReconstructor(void) |
| 163 | { |
| 164 | // Destructor for class AliMUONEventReconstructor |
| 165 | delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ???? |
| 9cbdf048 |
166 | for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) |
| a9e2aefa |
167 | delete fSegmentsPtr[st]; // Correct destruction of everything ???? |
| 168 | return; |
| 169 | } |
| 170 | |
| 171 | //__________________________________________________________________________ |
| 172 | void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void) |
| 173 | { |
| 174 | // Set reconstruction parameters to default values |
| 175 | // Would be much more convenient with a structure (or class) ???? |
| 9cbdf048 |
176 | fMinBendingMomentum = kDefaultMinBendingMomentum; |
| 177 | fMaxSigma2Distance = kDefaultMaxSigma2Distance; |
| a9e2aefa |
178 | |
| 9cbdf048 |
179 | AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); |
| a9e2aefa |
180 | // ******** Parameters for making HitsForRec |
| 181 | // minimum radius, |
| 182 | // like in TRACKF_STAT: |
| 183 | // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3; |
| 184 | // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9 |
| 9cbdf048 |
185 | for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) { |
| 186 | if (ch < 4) fRMin[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) * |
| a9e2aefa |
187 | 2.0 * TMath::Pi() / 180.0; |
| 188 | else fRMin[ch] = 30.0; |
| 189 | } |
| 190 | // maximum radius |
| 191 | // like in TRACKF_STAT (10 degrees ????) |
| 192 | fRMax[0] = fRMax[1] = 91.5; |
| 193 | fRMax[2] = fRMax[3] = 122.5; |
| 194 | fRMax[4] = fRMax[5] = 158.3; |
| 195 | fRMax[6] = fRMax[7] = 260.0; |
| 196 | fRMax[8] = fRMax[9] = 260.0; |
| 197 | |
| 198 | // ******** Parameters for making segments |
| 199 | // should be parametrized ???? |
| 200 | // according to interval between chambers in a station ???? |
| 201 | // Maximum distance in non bending plane |
| 202 | // 5 * 0.22 just to remember the way it was made in TRACKF_STAT |
| 203 | // SIGCUT*DYMAX(IZ) |
| 9cbdf048 |
204 | for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) |
| a9e2aefa |
205 | fSegmentMaxDistNonBending[st] = 5. * 0.22; |
| 206 | // Maximum distance in bending plane |
| 207 | // values from TRACKF_STAT corresponding to (J psi 20cm) |
| 208 | fSegmentMaxDistBending[0] = 1.5; |
| 209 | fSegmentMaxDistBending[1] = 1.5; |
| 210 | fSegmentMaxDistBending[2] = 3.0; |
| 211 | fSegmentMaxDistBending[3] = 6.0; |
| 212 | fSegmentMaxDistBending[4] = 6.0; |
| 213 | |
| 9cbdf048 |
214 | fBendingResolution = kDefaultBendingResolution; |
| 215 | fNonBendingResolution = kDefaultNonBendingResolution; |
| 216 | fChamberThicknessInX0 = kDefaultChamberThicknessInX0; |
| 217 | fSimpleBValue = kDefaultSimpleBValue; |
| 218 | fSimpleBLength = kDefaultSimpleBLength; |
| 219 | fSimpleBPosition = kDefaultSimpleBPosition; |
| 220 | fRecGeantHits = kDefaultRecGeantHits; |
| 221 | fEfficiency = kDefaultEfficiency; |
| 222 | fPrintLevel = kDefaultPrintLevel; |
| a9e2aefa |
223 | return; |
| 224 | } |
| 225 | |
| 226 | //__________________________________________________________________________ |
| 227 | Double_t AliMUONEventReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum) |
| 228 | { |
| 229 | // Returns impact parameter at vertex in bending plane (cm), |
| 230 | // from the signed bending momentum "BendingMomentum" in bending plane (GeV/c), |
| 231 | // using simple values for dipole magnetic field. |
| 232 | // The sign is the sign of the charge. |
| 233 | return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition / |
| 234 | BendingMomentum); |
| 235 | } |
| 236 | |
| 237 | //__________________________________________________________________________ |
| 238 | Double_t AliMUONEventReconstructor::GetBendingMomentumFromImpactParam(Double_t ImpactParam) |
| 239 | { |
| 240 | // Returns signed bending momentum in bending plane (GeV/c), |
| 241 | // from the impact parameter "ImpactParam" at vertex in bending plane (cm), |
| 242 | // using simple values for dipole magnetic field. |
| 243 | // The sign is the sign of the charge. |
| 244 | return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition / |
| 245 | ImpactParam); |
| 246 | } |
| 247 | |
| 248 | //__________________________________________________________________________ |
| 249 | void AliMUONEventReconstructor::SetBkgGeantFile(Text_t *BkgGeantFileName) |
| 250 | { |
| 251 | // Set background file ... for GEANT hits |
| 252 | // Must be called after having loaded the firts signal event |
| 253 | if (fPrintLevel >= 0) { |
| 254 | cout << "Enter SetBkgGeantFile with BkgGeantFileName ``" |
| 255 | << BkgGeantFileName << "''" << endl;} |
| 256 | if (strlen(BkgGeantFileName)) { |
| 257 | // BkgGeantFileName not empty: try to open the file |
| 258 | if (fPrintLevel >= 2) {cout << "Before File(Bkg)" << endl; gDirectory->Dump();} |
| 259 | fBkgGeantFile = new TFile(BkgGeantFileName); |
| 260 | if (fPrintLevel >= 2) {cout << "After File(Bkg)" << endl; gDirectory->Dump();} |
| 261 | if (fBkgGeantFile-> IsOpen()) { |
| 262 | if (fPrintLevel >= 0) { |
| 263 | cout << "Background for GEANT hits in file: ``" << BkgGeantFileName |
| 264 | << "'' successfully opened" << endl;} |
| 265 | } |
| 266 | else { |
| 267 | cout << "Background for GEANT hits in file: " << BkgGeantFileName << endl; |
| 268 | cout << "NOT FOUND: EXIT" << endl; |
| 269 | exit(0); // right instruction for exit ???? |
| 270 | } |
| 271 | // Arrays for "particles" and "hits" |
| 272 | fBkgGeantParticles = new TClonesArray("TParticle", 200); |
| 273 | fBkgGeantHits = new TClonesArray("AliMUONHit", 2000); |
| 274 | // Event number to -1 for initialization |
| 275 | fBkgGeantEventNumber = -1; |
| 276 | // Back to the signal file: |
| 277 | // first signal event must have been loaded previously, |
| 278 | // otherwise, Segmentation violation at the next instruction |
| 279 | // How is it possible to do smething better ???? |
| 280 | ((gAlice->TreeK())->GetCurrentFile())->cd(); |
| 281 | if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();} |
| 282 | } |
| 283 | return; |
| 284 | } |
| 285 | |
| 286 | //__________________________________________________________________________ |
| 287 | void AliMUONEventReconstructor::NextBkgGeantEvent(void) |
| 288 | { |
| 289 | // Get next event in background file for GEANT hits |
| 290 | // Goes back to event number 0 when end of file is reached |
| 291 | char treeName[20]; |
| 292 | TBranch *branch; |
| 293 | if (fPrintLevel >= 0) { |
| 294 | cout << "Enter NextBkgGeantEvent" << endl;} |
| 295 | // Clean previous event |
| 296 | if(fBkgGeantTK) delete fBkgGeantTK; |
| 297 | fBkgGeantTK = NULL; |
| 298 | if(fBkgGeantParticles) fBkgGeantParticles->Clear(); |
| 299 | if(fBkgGeantTH) delete fBkgGeantTH; |
| 300 | fBkgGeantTH = NULL; |
| 301 | if(fBkgGeantHits) fBkgGeantHits->Clear(); |
| 302 | // Increment event number |
| 303 | fBkgGeantEventNumber++; |
| 304 | // Get access to Particles and Hits for event from background file |
| 305 | if (fPrintLevel >= 2) {cout << "Before cd(Bkg)" << endl; gDirectory->Dump();} |
| 306 | fBkgGeantFile->cd(); |
| 307 | if (fPrintLevel >= 2) {cout << "After cd(Bkg)" << endl; gDirectory->Dump();} |
| 308 | // Particles: TreeK for event and branch "Particles" |
| 309 | sprintf(treeName, "TreeK%d", fBkgGeantEventNumber); |
| 310 | fBkgGeantTK = (TTree*)gDirectory->Get(treeName); |
| 311 | if (!fBkgGeantTK) { |
| 312 | if (fPrintLevel >= 0) { |
| 313 | cout << "Cannot find Kine Tree for background event: " << |
| 314 | fBkgGeantEventNumber << endl; |
| 315 | cout << "Goes back to event 0" << endl; |
| 316 | } |
| 317 | fBkgGeantEventNumber = 0; |
| 318 | sprintf(treeName, "TreeK%d", fBkgGeantEventNumber); |
| 319 | fBkgGeantTK = (TTree*)gDirectory->Get(treeName); |
| 320 | if (!fBkgGeantTK) { |
| 321 | cout << "ERROR: cannot find Kine Tree for background event: " << |
| 322 | fBkgGeantEventNumber << endl; |
| 323 | exit(0); |
| 324 | } |
| 325 | } |
| 326 | if (fBkgGeantTK) |
| 327 | fBkgGeantTK->SetBranchAddress("Particles", &fBkgGeantParticles); |
| 328 | fBkgGeantTK->GetEvent(0); // why event 0 ???? necessary ???? |
| 329 | // Hits: TreeH for event and branch "MUON" |
| 330 | sprintf(treeName, "TreeH%d", fBkgGeantEventNumber); |
| 331 | fBkgGeantTH = (TTree*)gDirectory->Get(treeName); |
| 332 | if (!fBkgGeantTH) { |
| 333 | cout << "ERROR: cannot find Hits Tree for background event: " << |
| 334 | fBkgGeantEventNumber << endl; |
| 335 | exit(0); |
| 336 | } |
| 337 | if (fBkgGeantTH && fBkgGeantHits) { |
| 338 | branch = fBkgGeantTH->GetBranch("MUON"); |
| 339 | if (branch) branch->SetAddress(&fBkgGeantHits); |
| 340 | } |
| 341 | fBkgGeantTH->GetEntries(); // necessary ???? |
| 342 | // Back to the signal file |
| 343 | ((gAlice->TreeK())->GetCurrentFile())->cd(); |
| 344 | if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();} |
| 345 | return; |
| 346 | } |
| 347 | |
| 348 | //__________________________________________________________________________ |
| 349 | void AliMUONEventReconstructor::EventReconstruct(void) |
| 350 | { |
| 351 | // To reconstruct one event |
| 352 | if (fPrintLevel >= 1) cout << "enter EventReconstruct" << endl; |
| 353 | MakeEventToBeReconstructed(); |
| 354 | MakeSegments(); |
| 355 | MakeTracks(); |
| 356 | return; |
| 357 | } |
| 358 | |
| 359 | //__________________________________________________________________________ |
| 360 | void AliMUONEventReconstructor::ResetHitsForRec(void) |
| 361 | { |
| 362 | // To reset the array and the number of HitsForRec, |
| 363 | // and also the number of HitsForRec |
| 364 | // and the index of the first HitForRec per chamber |
| 365 | if (fHitsForRecPtr) fHitsForRecPtr->Clear(); |
| 366 | fNHitsForRec = 0; |
| 9cbdf048 |
367 | for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) |
| a9e2aefa |
368 | fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0; |
| 369 | return; |
| 370 | } |
| 371 | |
| 372 | //__________________________________________________________________________ |
| 373 | void AliMUONEventReconstructor::ResetSegments(void) |
| 374 | { |
| 375 | // To reset the TClonesArray of segments and the number of Segments |
| 376 | // for all stations |
| 9cbdf048 |
377 | for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) { |
| a9e2aefa |
378 | if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear(); |
| 379 | fNSegments[st] = 0; |
| 380 | } |
| 381 | return; |
| 382 | } |
| 383 | |
| 384 | //__________________________________________________________________________ |
| 385 | void AliMUONEventReconstructor::ResetTracks(void) |
| 386 | { |
| 387 | // To reset the TClonesArray of reconstructed tracks |
| 388 | if (fRecTracksPtr) fRecTracksPtr->Clear(); |
| 389 | fNRecTracks = 0; |
| 390 | return; |
| 391 | } |
| 392 | |
| 393 | //__________________________________________________________________________ |
| 394 | void AliMUONEventReconstructor::MakeEventToBeReconstructed(void) |
| 395 | { |
| 396 | // To make the list of hits to be reconstructed, |
| 397 | // either from the GEANT hits or from the raw clusters |
| 398 | // according to the parameter set for the reconstructor |
| 399 | if (fPrintLevel >= 1) cout << "enter MakeEventToBeReconstructed" << endl; |
| 400 | ResetHitsForRec(); |
| 401 | if (fRecGeantHits == 1) { |
| 402 | // Reconstruction from GEANT hits |
| 403 | // Back to the signal file |
| 404 | ((gAlice->TreeK())->GetCurrentFile())->cd(); |
| 405 | // Signal hits |
| 406 | // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? |
| 407 | // Security on MUON ???? |
| 408 | AddHitsForRecFromGEANT(gAlice->TreeH()); |
| 409 | // Background hits |
| 410 | AddHitsForRecFromBkgGEANT(fBkgGeantTH, fBkgGeantHits); |
| 411 | // Sort HitsForRec in increasing order with respect to chamber number |
| 412 | SortHitsForRecWithIncreasingChamber(); |
| 413 | } |
| 414 | else { |
| 415 | // Reconstruction from raw clusters |
| 416 | // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? |
| 417 | // Security on MUON ???? |
| 418 | // TreeR assumed to be be "prepared" in calling function |
| 419 | // by "MUON->GetTreeR(nev)" ???? |
| 9cbdf048 |
420 | TTree *treeR = gAlice->TreeR(); |
| 421 | AddHitsForRecFromRawClusters(treeR); |
| a9e2aefa |
422 | // No sorting: it is done automatically in the previous function |
| 423 | } |
| 424 | if (fPrintLevel >= 10) { |
| 425 | cout << "end of MakeEventToBeReconstructed" << endl; |
| 426 | cout << "NHitsForRec: " << fNHitsForRec << endl; |
| 9cbdf048 |
427 | for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) { |
| a9e2aefa |
428 | cout << "chamber(0...): " << ch |
| 429 | << " NHitsForRec: " << fNHitsForRecPerChamber[ch] |
| 430 | << " index(first HitForRec): " << fIndexOfFirstHitForRecPerChamber[ch] |
| 431 | << endl; |
| 432 | for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch]; |
| 433 | hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch]; |
| 434 | hit++) { |
| 435 | cout << "HitForRec index(0...): " << hit << endl; |
| 436 | ((*fHitsForRecPtr)[hit])->Dump(); |
| 437 | } |
| 438 | } |
| 439 | } |
| 440 | return; |
| 441 | } |
| 442 | |
| 443 | //__________________________________________________________________________ |
| 444 | void AliMUONEventReconstructor::AddHitsForRecFromGEANT(TTree *TH) |
| 445 | { |
| 446 | // To add to the list of hits for reconstruction |
| 447 | // the GEANT signal hits from a hit tree TH. |
| 448 | if (fPrintLevel >= 2) |
| 449 | cout << "enter AddHitsForRecFromGEANT with TH: " << TH << endl; |
| 450 | if (TH == NULL) return; |
| 9cbdf048 |
451 | AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? |
| a9e2aefa |
452 | // Security on MUON ???? |
| 453 | // See whether it could be the same for signal and background ???? |
| 454 | // Loop over tracks in tree |
| 455 | Int_t ntracks = (Int_t) TH->GetEntries(); |
| 456 | if (fPrintLevel >= 2) |
| 457 | cout << "ntracks: " << ntracks << endl; |
| 458 | for (Int_t track = 0; track < ntracks; track++) { |
| 459 | gAlice->ResetHits(); |
| 460 | TH->GetEvent(track); |
| 461 | // Loop over hits |
| 462 | Int_t hit = 0; |
| 9cbdf048 |
463 | for (AliMUONHit* mHit = (AliMUONHit*) pMUON->FirstHit(-1); |
| a9e2aefa |
464 | mHit; |
| 9cbdf048 |
465 | mHit = (AliMUONHit*) pMUON->NextHit(), hit++) { |
| a9e2aefa |
466 | NewHitForRecFromGEANT(mHit,track, hit, 1); |
| 467 | } // end of hit loop |
| 468 | } // end of track loop |
| 469 | return; |
| 470 | } |
| 471 | |
| 472 | //__________________________________________________________________________ |
| 473 | void AliMUONEventReconstructor::AddHitsForRecFromBkgGEANT(TTree *TH, TClonesArray *Hits) |
| 474 | { |
| 475 | // To add to the list of hits for reconstruction |
| 476 | // the GEANT background hits from a hit tree TH and a pointer Hits to a hit list. |
| 477 | // How to have only one function "AddHitsForRecFromGEANT" ???? |
| 478 | if (fPrintLevel >= 2) |
| 479 | cout << "enter AddHitsForRecFromBkgGEANT with TH: " << TH << endl; |
| 480 | if (TH == NULL) return; |
| 481 | // Loop over tracks in tree |
| 482 | Int_t ntracks = (Int_t) TH->GetEntries(); |
| 483 | if (fPrintLevel >= 2) |
| 484 | cout << "ntracks: " << ntracks << endl; |
| 485 | for (Int_t track = 0; track < ntracks; track++) { |
| 486 | if (Hits) Hits->Clear(); |
| 487 | TH->GetEvent(track); |
| 488 | // Loop over hits |
| 489 | for (Int_t hit = 0; hit < Hits->GetEntriesFast(); hit++) { |
| 490 | NewHitForRecFromGEANT((AliMUONHit*) (*Hits)[hit], track, hit, 0); |
| 491 | } // end of hit loop |
| 492 | } // end of track loop |
| 493 | return; |
| 494 | } |
| 495 | |
| 496 | //__________________________________________________________________________ |
| 497 | AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* Hit, Int_t TrackNumber, Int_t HitNumber, Int_t Signal) |
| 498 | { |
| 499 | // To make a new hit for reconstruction from a GEANT hit pointed to by "Hit", |
| 500 | // with hit number "HitNumber" in the track numbered "TrackNumber", |
| 501 | // either from signal ("Signal" = 1) or background ("Signal" = 0) event. |
| 502 | // Selects hits in tracking (not trigger) chambers. |
| 503 | // Takes into account the efficiency (fEfficiency) |
| 504 | // and the smearing from resolution (fBendingResolution and fNonBendingResolution). |
| 505 | // Adds a condition on the radius between RMin and RMax |
| 506 | // to better simulate the real chambers. |
| 507 | // Returns the pointer to the new hit for reconstruction, |
| 508 | // or NULL in case of inefficiency or non tracking chamber or bad radius. |
| 509 | // No condition on at most 20 cm from a muon from a resonance |
| 510 | // like in Fortran TRACKF_STAT. |
| 511 | AliMUONHitForRec* hitForRec; |
| 512 | Double_t bendCoor, nonBendCoor, radius; |
| 513 | Int_t chamber = Hit->fChamber - 1; // chamber(0...) |
| 514 | // only in tracking chambers (fChamber starts at 1) |
| 9cbdf048 |
515 | if (chamber >= kMaxMuonTrackingChambers) return NULL; |
| a9e2aefa |
516 | // only if hit is efficient (keep track for checking ????) |
| 517 | if (gRandom->Rndm() > fEfficiency) return NULL; |
| 518 | // only if radius between RMin and RMax |
| 519 | bendCoor = Hit->fY; |
| 520 | nonBendCoor = Hit->fX; |
| 521 | radius = TMath::Sqrt((bendCoor * bendCoor) + (nonBendCoor * nonBendCoor)); |
| 522 | if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL; |
| 523 | // new AliMUONHitForRec from GEANT hit and increment number of AliMUONHitForRec's |
| 524 | hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(Hit); |
| 525 | fNHitsForRec++; |
| 526 | // add smearing from resolution |
| 527 | hitForRec->SetBendingCoor(bendCoor + gRandom->Gaus(0., fBendingResolution)); |
| 528 | hitForRec->SetNonBendingCoor(nonBendCoor |
| 529 | + gRandom->Gaus(0., fNonBendingResolution)); |
| 530 | // more information into HitForRec |
| 531 | // resolution: angular effect to be added here ???? |
| 532 | hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution); |
| 533 | hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution); |
| 534 | // GEANT track info |
| 535 | hitForRec->SetHitNumber(HitNumber); |
| 536 | hitForRec->SetTHTrack(TrackNumber); |
| 537 | hitForRec->SetGeantSignal(Signal); |
| 538 | if (fPrintLevel >= 10) { |
| 539 | cout << "track: " << TrackNumber << " hit: " << HitNumber << endl; |
| 540 | Hit->Dump(); |
| 541 | cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl; |
| 542 | hitForRec->Dump();} |
| 543 | return hitForRec; |
| 544 | } |
| 545 | |
| 546 | //__________________________________________________________________________ |
| 547 | void AliMUONEventReconstructor::SortHitsForRecWithIncreasingChamber() |
| 548 | { |
| 549 | // Sort HitsForRec's in increasing order with respect to chamber number. |
| 550 | // Uses the function "Compare". |
| 551 | // Update the information for HitsForRec per chamber too. |
| 552 | Int_t ch, nhits, prevch; |
| 553 | fHitsForRecPtr->Sort(); |
| 9cbdf048 |
554 | for (ch = 0; ch < kMaxMuonTrackingChambers; ch++) { |
| a9e2aefa |
555 | fNHitsForRecPerChamber[ch] = 0; |
| 556 | fIndexOfFirstHitForRecPerChamber[ch] = 0; |
| 557 | } |
| 558 | prevch = 0; // previous chamber |
| 559 | nhits = 0; // number of hits in current chamber |
| 560 | // Loop over HitsForRec |
| 561 | for (Int_t hit = 0; hit < fNHitsForRec; hit++) { |
| 562 | // chamber number (0...) |
| 563 | ch = ((AliMUONHitForRec*) ((*fHitsForRecPtr)[hit]))->GetChamberNumber(); |
| 564 | // increment number of hits in current chamber |
| 565 | (fNHitsForRecPerChamber[ch])++; |
| 566 | // update index of first HitForRec in current chamber |
| 567 | // if chamber number different from previous one |
| 568 | if (ch != prevch) { |
| 569 | fIndexOfFirstHitForRecPerChamber[ch] = hit; |
| 570 | prevch = ch; |
| 571 | } |
| 572 | } |
| 573 | return; |
| 574 | } |
| 575 | |
| 576 | // //__________________________________________________________________________ |
| 577 | // void AliMUONEventReconstructor::AddHitsForRecFromCathodeCorrelations(TTree* TC) |
| 578 | // { |
| 579 | // // OLD VERSION WHEN ONE ONE WAS USING SO CALLED CATHODE CORRELATIONS |
| 580 | // // To add to the list of hits for reconstruction |
| 581 | // // the (cathode correlated) raw clusters |
| 582 | // // No condition added, like in Fortran TRACKF_STAT, |
| 583 | // // on the radius between RMin and RMax. |
| 584 | // AliMUONHitForRec *hitForRec; |
| 585 | // if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromCathodeCorrelations" << endl; |
| 586 | // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? |
| 587 | // // Security on MUON ???? |
| 588 | // // Loop over tracking chambers |
| 9cbdf048 |
589 | // for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) { |
| a9e2aefa |
590 | // // number of HitsForRec to 0 for the chamber |
| 591 | // fNHitsForRecPerChamber[ch] = 0; |
| 592 | // // index of first HitForRec for the chamber |
| 593 | // if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0; |
| 594 | // else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec; |
| 595 | // TClonesArray *reconst_hits = MUON->ReconstHitsAddress(ch); |
| 596 | // MUON->ResetReconstHits(); |
| 597 | // TC->GetEvent(); |
| 598 | // Int_t ncor = (Int_t)reconst_hits->GetEntries(); |
| 599 | // // Loop over (cathode correlated) raw clusters |
| 600 | // for (Int_t cor = 0; cor < ncor; cor++) { |
| 601 | // AliMUONReconstHit * mCor = |
| 602 | // (AliMUONReconstHit*) reconst_hits->UncheckedAt(cor); |
| 603 | // // new AliMUONHitForRec from (cathode correlated) raw cluster |
| 604 | // // and increment number of AliMUONHitForRec's (total and in chamber) |
| 605 | // hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(mCor); |
| 606 | // fNHitsForRec++; |
| 607 | // (fNHitsForRecPerChamber[ch])++; |
| 608 | // // more information into HitForRec |
| 609 | // hitForRec->SetChamberNumber(ch); |
| 610 | // hitForRec->SetHitNumber(cor); |
| 611 | // // Z coordinate of the chamber (cm) with sign opposite to GEANT sign |
| 612 | // // could (should) be more exact from chamber geometry ???? |
| 613 | // hitForRec->SetZ(-(&(MUON->Chamber(ch)))->Z()); |
| 614 | // if (fPrintLevel >= 10) { |
| 615 | // cout << "chamber (0...): " << ch << |
| 616 | // " cathcorrel (0...): " << cor << endl; |
| 617 | // mCor->Dump(); |
| 618 | // cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl; |
| 619 | // hitForRec->Dump();} |
| 620 | // } // end of cluster loop |
| 621 | // } // end of chamber loop |
| 622 | // return; |
| 623 | // } |
| 624 | |
| 625 | //__________________________________________________________________________ |
| 626 | void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR) |
| 627 | { |
| 628 | // To add to the list of hits for reconstruction all the raw clusters |
| 629 | // No condition added, like in Fortran TRACKF_STAT, |
| 630 | // on the radius between RMin and RMax. |
| 631 | AliMUONHitForRec *hitForRec; |
| 632 | AliMUONRawCluster *clus; |
| 633 | Int_t iclus, nclus; |
| 634 | TClonesArray *rawclusters; |
| 635 | if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromRawClusters" << endl; |
| 9cbdf048 |
636 | AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? |
| a9e2aefa |
637 | // Security on MUON ???? |
| 638 | // Loop over tracking chambers |
| 9cbdf048 |
639 | for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) { |
| a9e2aefa |
640 | // number of HitsForRec to 0 for the chamber |
| 641 | fNHitsForRecPerChamber[ch] = 0; |
| 642 | // index of first HitForRec for the chamber |
| 643 | if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0; |
| 644 | else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec; |
| 9cbdf048 |
645 | rawclusters = pMUON->RawClustAddress(ch); |
| 646 | pMUON->ResetRawClusters(); |
| a9e2aefa |
647 | TR->GetEvent((Int_t) (TR->GetEntries()) - 1); // to be checked ???? |
| 648 | nclus = (Int_t) (rawclusters->GetEntries()); |
| 649 | // Loop over (cathode correlated) raw clusters |
| 650 | for (iclus = 0; iclus < nclus; iclus++) { |
| 651 | clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus); |
| 652 | // new AliMUONHitForRec from raw cluster |
| 653 | // and increment number of AliMUONHitForRec's (total and in chamber) |
| 654 | hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus); |
| 655 | fNHitsForRec++; |
| 656 | (fNHitsForRecPerChamber[ch])++; |
| 657 | // more information into HitForRec |
| 658 | // resolution: info should be already in raw cluster and taken from it ???? |
| 659 | hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution); |
| 660 | hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution); |
| 661 | // original raw cluster |
| 662 | hitForRec->SetChamberNumber(ch); |
| 663 | hitForRec->SetHitNumber(iclus); |
| 664 | // Z coordinate of the chamber (cm) with sign opposite to GEANT sign |
| 665 | // could (should) be more exact from chamber geometry ???? |
| 9cbdf048 |
666 | hitForRec->SetZ(-(&(pMUON->Chamber(ch)))->Z()); |
| a9e2aefa |
667 | if (fPrintLevel >= 10) { |
| 668 | cout << "chamber (0...): " << ch << |
| 669 | " raw cluster (0...): " << iclus << endl; |
| 670 | clus->Dump(); |
| 671 | cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl; |
| 672 | hitForRec->Dump();} |
| 673 | } // end of cluster loop |
| 674 | } // end of chamber loop |
| 675 | return; |
| 676 | } |
| 677 | |
| 678 | //__________________________________________________________________________ |
| 679 | void AliMUONEventReconstructor::MakeSegments(void) |
| 680 | { |
| 681 | // To make the list of segments in all stations, |
| 682 | // from the list of hits to be reconstructed |
| 683 | if (fPrintLevel >= 1) cout << "enter MakeSegments" << endl; |
| 684 | ResetSegments(); |
| 685 | // Loop over stations |
| 9cbdf048 |
686 | for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) |
| a9e2aefa |
687 | MakeSegmentsPerStation(st); |
| 688 | if (fPrintLevel >= 10) { |
| 689 | cout << "end of MakeSegments" << endl; |
| 9cbdf048 |
690 | for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) { |
| a9e2aefa |
691 | cout << "station(0...): " << st |
| 692 | << " Segments: " << fNSegments[st] |
| 693 | << endl; |
| 694 | for (Int_t seg = 0; |
| 695 | seg < fNSegments[st]; |
| 696 | seg++) { |
| 697 | cout << "Segment index(0...): " << seg << endl; |
| 698 | ((*fSegmentsPtr[st])[seg])->Dump(); |
| 699 | } |
| 700 | } |
| 701 | } |
| 702 | return; |
| 703 | } |
| 704 | |
| 705 | //__________________________________________________________________________ |
| 706 | void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station) |
| 707 | { |
| 708 | // To make the list of segments in station number "Station" (0...) |
| 709 | // from the list of hits to be reconstructed. |
| 710 | // Updates "fNSegments"[Station]. |
| 711 | // Segments in stations 4 and 5 are sorted |
| 712 | // according to increasing absolute value of "impact parameter" |
| 713 | AliMUONHitForRec *hit1Ptr, *hit2Ptr; |
| 714 | AliMUONSegment *segment; |
| 715 | Bool_t last2st; |
| 716 | Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor, |
| 717 | impactParam, maxImpactParam; |
| 9cbdf048 |
718 | AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? |
| a9e2aefa |
719 | if (fPrintLevel >= 1) |
| 720 | cout << "enter MakeSegmentsPerStation (0...) " << Station << endl; |
| 721 | // first and second chambers (0...) in the station |
| 722 | Int_t ch1 = 2 * Station; |
| 723 | Int_t ch2 = ch1 + 1; |
| 724 | // variable true for stations downstream of the dipole: |
| 725 | // Station(0..4) equal to 3 or 4 |
| 726 | if ((Station == 3) || (Station == 4)) { |
| 727 | last2st = kTRUE; |
| 728 | // maximum impact parameter (cm) according to fMinBendingMomentum... |
| 729 | maxImpactParam = |
| 730 | TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum)); |
| 731 | } |
| 732 | else last2st = kFALSE; |
| 733 | // extrapolation factor from Z of first chamber to Z of second chamber |
| 734 | // dZ to be changed to take into account fine structure of chambers ???? |
| 735 | Double_t extrapFact = |
| 9cbdf048 |
736 | (&(pMUON->Chamber(ch2)))->Z() / (&(pMUON->Chamber(ch1)))->Z(); |
| a9e2aefa |
737 | // index for current segment |
| 738 | Int_t segmentIndex = 0; |
| 739 | // Loop over HitsForRec in the first chamber of the station |
| 740 | for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1]; |
| 741 | hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1]; |
| 742 | hit1++) { |
| 743 | // pointer to the HitForRec |
| 744 | hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]); |
| 745 | // extrapolation, |
| 746 | // on the straight line joining the HitForRec to the vertex (0,0,0), |
| 747 | // to the Z of the second chamber of the station |
| 748 | extBendCoor = extrapFact * hit1Ptr->GetBendingCoor(); |
| 749 | extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor(); |
| 750 | // Loop over HitsForRec in the second chamber of the station |
| 751 | for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2]; |
| 752 | hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2]; |
| 753 | hit2++) { |
| 754 | // pointer to the HitForRec |
| 755 | hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]); |
| 756 | // absolute values of distances, in bending and non bending planes, |
| 757 | // between the HitForRec in the second chamber |
| 758 | // and the previous extrapolation |
| 759 | distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor); |
| 760 | distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor); |
| 761 | if (last2st) { |
| 762 | // bending slope |
| 763 | bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) / |
| 764 | (hit1Ptr->GetZ() - hit2Ptr->GetZ()); |
| 765 | // absolute value of impact parameter |
| 766 | impactParam = |
| 767 | TMath::Abs(hit1Ptr->GetBendingCoor() - hit2Ptr->GetZ() * bendingSlope); |
| 768 | } |
| 769 | // check for distances not too large, |
| 770 | // and impact parameter not too big if stations downstream of the dipole. |
| 771 | // Conditions "distBend" and "impactParam" correlated for these stations ???? |
| 772 | if ((distBend < fSegmentMaxDistBending[Station]) && |
| 773 | (distNonBend < fSegmentMaxDistNonBending[Station]) && |
| 774 | (!last2st || (impactParam < maxImpactParam))) { |
| 775 | // make new segment |
| 776 | segment = new ((*fSegmentsPtr[Station])[segmentIndex]) |
| 777 | AliMUONSegment(hit1Ptr, hit2Ptr); |
| 778 | // update "link" to this segment from the hit in the first chamber |
| 779 | if (hit1Ptr->GetNSegments() == 0) |
| 780 | hit1Ptr->SetIndexOfFirstSegment(segmentIndex); |
| 781 | hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1); |
| 782 | if (fPrintLevel >= 10) { |
| 783 | cout << "segmentIndex(0...): " << segmentIndex |
| 784 | << " distBend: " << distBend |
| 785 | << " distNonBend: " << distNonBend |
| 786 | << endl; |
| 787 | segment->Dump(); |
| 788 | cout << "HitForRec in first chamber" << endl; |
| 789 | hit1Ptr->Dump(); |
| 790 | cout << "HitForRec in second chamber" << endl; |
| 791 | hit2Ptr->Dump(); |
| 792 | }; |
| 793 | // increment index for current segment |
| 794 | segmentIndex++; |
| 795 | } |
| 796 | } //for (Int_t hit2 |
| 797 | } // for (Int_t hit1... |
| 798 | fNSegments[Station] = segmentIndex; |
| 799 | // Sorting according to "impact parameter" if station(1..5) 4 or 5, |
| 800 | // i.e. Station(0..4) 3 or 4, using the function "Compare". |
| 801 | // After this sorting, it is impossible to use |
| 802 | // the "fNSegments" and "fIndexOfFirstSegment" |
| 803 | // of the HitForRec in the first chamber to explore all segments formed with it. |
| 804 | // Is this sorting really needed ???? |
| 805 | if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort(); |
| 806 | if (fPrintLevel >= 1) cout << "Station: " << Station << " NSegments: " |
| 807 | << fNSegments[Station] << endl; |
| 808 | return; |
| 809 | } |
| 810 | |
| 811 | //__________________________________________________________________________ |
| 812 | void AliMUONEventReconstructor::MakeTracks(void) |
| 813 | { |
| 814 | // To make the tracks, |
| 815 | // from the list of segments and points in all stations |
| 816 | if (fPrintLevel >= 1) cout << "enter MakeTracks" << endl; |
| 817 | ResetTracks(); |
| 818 | // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5 |
| 819 | MakeTrackCandidates(); |
| 820 | // Follow tracks in stations(1..) 3, 2 and 1 |
| 821 | FollowTracks(); |
| 822 | return; |
| 823 | } |
| 824 | |
| 825 | //__________________________________________________________________________ |
| 826 | Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment) |
| 827 | { |
| 828 | // To make track candidates with two segments in stations(1..) 4 and 5, |
| 829 | // the first segment being pointed to by "BegSegment". |
| 830 | // Returns the number of such track candidates. |
| 831 | Int_t endStation, iEndSegment, nbCan2Seg; |
| 832 | AliMUONSegment *endSegment, *extrapSegment; |
| 833 | AliMUONTrack *recTrack; |
| 834 | Double_t mcsFactor; |
| 835 | if (fPrintLevel >= 1) cout << "enter MakeTrackCandidatesWithTwoSegments" << endl; |
| 836 | // Station for the end segment |
| 837 | endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2; |
| 838 | // multiple scattering factor corresponding to one chamber |
| 839 | mcsFactor = 0.0136 / |
| 840 | GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact()); |
| 841 | mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor; |
| 842 | // linear extrapolation to end station |
| 843 | extrapSegment = |
| 844 | BegSegment->CreateSegmentFromLinearExtrapToStation(endStation, mcsFactor); |
| 845 | // number of candidates with 2 segments to 0 |
| 846 | nbCan2Seg = 0; |
| 847 | // Loop over segments in the end station |
| 848 | for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) { |
| 849 | // pointer to segment |
| 850 | endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]); |
| 851 | // test compatibility between current segment and "extrapSegment" |
| 852 | if ((endSegment-> |
| 853 | NormalizedChi2WithSegment(extrapSegment, |
| 854 | fMaxSigma2Distance)) <= 4.0) { |
| 855 | // both segments compatible: |
| 856 | // make track candidate from "begSegment" and "endSegment" |
| 857 | if (fPrintLevel >= 2) |
| 858 | cout << "TrackCandidate with Segment " << iEndSegment << |
| 859 | " in Station(0..) " << endStation << endl; |
| 860 | // flag for both segments in one track: |
| 861 | // to be done in track constructor ???? |
| 862 | BegSegment->SetInTrack(kTRUE); |
| 863 | endSegment->SetInTrack(kTRUE); |
| 864 | recTrack = new ((*fRecTracksPtr)[fNRecTracks]) |
| 865 | AliMUONTrack(BegSegment, endSegment, this); |
| 866 | fNRecTracks++; |
| 867 | if (fPrintLevel >= 10) recTrack->RecursiveDump(); |
| 868 | // increment number of track candidates with 2 segments |
| 869 | nbCan2Seg++; |
| 870 | } |
| 871 | } // for (iEndSegment = 0;... |
| 872 | delete extrapSegment; // should not delete HitForRec's it points to !!!! |
| 873 | return nbCan2Seg; |
| 874 | } |
| 875 | |
| 876 | //__________________________________________________________________________ |
| 877 | Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment) |
| 878 | { |
| 879 | // To make track candidates with one segment and one point |
| 880 | // in stations(1..) 4 and 5, |
| 881 | // the segment being pointed to by "BegSegment". |
| 882 | Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit; |
| 883 | AliMUONHitForRec *extrapHitForRec, *hit; |
| 884 | AliMUONTrack *recTrack; |
| 885 | Double_t mcsFactor; |
| 886 | if (fPrintLevel >= 1) |
| 887 | cout << "enter MakeTrackCandidatesWithOneSegmentAndOnePoint" << endl; |
| 888 | // station for the end point |
| 889 | endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2; |
| 890 | // multiple scattering factor corresponding to one chamber |
| 891 | mcsFactor = 0.0136 / |
| 892 | GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact()); |
| 893 | mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor; |
| 894 | // first and second chambers(0..) in the end station |
| 895 | ch1 = 2 * endStation; |
| 896 | ch2 = ch1 + 1; |
| 897 | // number of candidates to 0 |
| 898 | nbCan1Seg1Hit = 0; |
| 899 | // Loop over chambers of the end station |
| 900 | for (ch = ch2; ch >= ch1; ch--) { |
| 901 | // linear extrapolation to chamber |
| 902 | extrapHitForRec = |
| 903 | BegSegment->CreateHitForRecFromLinearExtrapToChamber(ch, mcsFactor); |
| 904 | // limits for the hit index in the loop |
| 905 | iHitMin = fIndexOfFirstHitForRecPerChamber[ch]; |
| 906 | iHitMax = iHitMin + fNHitsForRecPerChamber[ch]; |
| 907 | // Loop over HitForRec's in the chamber |
| 908 | for (iHit = iHitMin; iHit < iHitMax; iHit++) { |
| 909 | // pointer to HitForRec |
| 910 | hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]); |
| 911 | // test compatibility between current HitForRec and "extrapHitForRec" |
| 912 | if ((hit-> |
| 913 | NormalizedChi2WithHitForRec(extrapHitForRec, |
| 914 | fMaxSigma2Distance)) <= 2.0) { |
| 915 | // both HitForRec's compatible: |
| 916 | // make track candidate from begSegment and current HitForRec |
| 917 | if (fPrintLevel >= 2) |
| 918 | cout << "TrackCandidate with HitForRec " << iHit << |
| 919 | " in Chamber(0..) " << ch << endl; |
| 920 | // flag for beginning segments in one track: |
| 921 | // to be done in track constructor ???? |
| 922 | BegSegment->SetInTrack(kTRUE); |
| 923 | recTrack = new ((*fRecTracksPtr)[fNRecTracks]) |
| 924 | AliMUONTrack(BegSegment, hit, this); |
| 925 | // the right place to eliminate "double counting" ???? how ???? |
| 926 | fNRecTracks++; |
| 927 | if (fPrintLevel >= 10) recTrack->RecursiveDump(); |
| 928 | // increment number of track candidates |
| 929 | nbCan1Seg1Hit++; |
| 930 | } |
| 931 | } // for (iHit = iHitMin;... |
| 932 | delete extrapHitForRec; |
| 933 | } // for (ch = ch2;... |
| 934 | return nbCan1Seg1Hit; |
| 935 | } |
| 936 | |
| 937 | //__________________________________________________________________________ |
| 938 | void AliMUONEventReconstructor::MakeTrackCandidates(void) |
| 939 | { |
| 940 | // To make track candidates |
| 941 | // with at least 3 aligned points in stations(1..) 4 and 5 |
| 942 | // (two Segment's or one Segment and one HitForRec) |
| 943 | Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg; |
| 944 | AliMUONSegment *begSegment; |
| 945 | if (fPrintLevel >= 1) cout << "enter MakeTrackCandidates" << endl; |
| 946 | // Loop over stations(1..) 5 and 4 for the beginning segment |
| 947 | for (begStation = 4; begStation > 2; begStation--) { |
| 948 | // Loop over segments in the beginning station |
| 949 | for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) { |
| 950 | // pointer to segment |
| 951 | begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]); |
| 952 | if (fPrintLevel >= 2) |
| 953 | cout << "look for TrackCandidate's with Segment " << iBegSegment << |
| 954 | " in Station(0..) " << begStation << endl; |
| 955 | // Look for track candidates with two segments, |
| 956 | // "begSegment" and all compatible segments in other station. |
| 957 | // Only for beginning station(1..) 5 |
| 958 | // because candidates with 2 segments have to looked for only once. |
| 959 | if (begStation == 4) |
| 960 | nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment); |
| 961 | // Look for track candidates with one segments and one point, |
| 962 | // "begSegment" and all compatible HitForRec's in other station. |
| 963 | // Only if "begSegment" does not belong already to a track candidate. |
| 964 | // Is that a too strong condition ???? |
| 965 | if (!(begSegment->GetInTrack())) |
| 966 | nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment); |
| 967 | } // for (iBegSegment = 0;... |
| 968 | } // for (begStation = 4;... |
| 969 | return; |
| 970 | } |
| 971 | |
| 972 | //__________________________________________________________________________ |
| 973 | void AliMUONEventReconstructor::FollowTracks(void) |
| 974 | { |
| 975 | // Follow tracks in stations(1..) 3, 2 and 1 |
| 976 | AliMUONHitForRec *bestHit, *extrapHit, *hit; |
| 977 | AliMUONSegment *bestSegment, *extrapSegment, *segment; |
| 978 | AliMUONTrack *track; |
| 979 | AliMUONTrackParam *trackParam1, trackParam[2]; |
| 980 | Int_t ch, chBestHit, iHit, iSegment, station, trackIndex; |
| 981 | Double_t bestChi2, chi2, dZ1, dZ2, maxSigma2Distance, mcsFactor; |
| 9cbdf048 |
982 | AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? |
| a9e2aefa |
983 | maxSigma2Distance = 4. * fMaxSigma2Distance; // sigma2cut increased by 4 !!!! |
| 984 | if (fPrintLevel >= 2) |
| 985 | cout << "enter FollowTracks" << endl; |
| 986 | // Loop over track candidates |
| 987 | for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) { |
| 988 | if (fPrintLevel >= 2) |
| 989 | cout << "FollowTracks: track candidate(0..): " << trackIndex << endl; |
| 990 | // function for each track candidate ???? |
| 991 | track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]); |
| 992 | // Fit track candidate from vertex at X = Y = 0 |
| 993 | track->Fit(track->GetTrackParamAtVertex(), 3); |
| 994 | if (fPrintLevel >= 10) { |
| 995 | cout << "FollowTracks: track candidate(0..): " << trackIndex |
| 996 | << " after fit in stations(1..) 4 and 5" << endl; |
| 997 | track->RecursiveDump(); |
| 998 | } |
| 999 | // Loop over stations(1..) 3, 2 and 1 |
| 1000 | // something SPECIAL for stations 2 and 1 for majority coincidence ???? |
| 1001 | for (station = 2; station >= 0; station--) { |
| 1002 | // Track parameters at first track hit (smallest Z) |
| 1003 | trackParam1 = ((AliMUONTrackHit*) |
| 1004 | (track->GetTrackHitsPtr()->First()))->GetTrackParam(); |
| 1005 | // extrapolation to station |
| 1006 | trackParam1->ExtrapToStation(station, trackParam); |
| 79e1e601 |
1007 | extrapSegment = new AliMUONSegment(); // empty segment |
| a9e2aefa |
1008 | // multiple scattering factor corresponding to one chamber |
| 1009 | // and momentum in bending plane (not total) |
| 1010 | mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum(); |
| 1011 | mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor; |
| 1012 | // Z difference from previous station |
| 9cbdf048 |
1013 | dZ1 = (&(pMUON->Chamber(2 * station)))->Z() - |
| 1014 | (&(pMUON->Chamber(2 * station + 2)))->Z(); |
| a9e2aefa |
1015 | // Z difference between the two previous stations |
| 9cbdf048 |
1016 | dZ2 = (&(pMUON->Chamber(2 * station + 2)))->Z() - |
| 1017 | (&(pMUON->Chamber(2 * station + 4)))->Z(); |
| a9e2aefa |
1018 | extrapSegment->SetBendingCoorReso2(fBendingResolution); |
| 1019 | extrapSegment->SetNonBendingCoorReso2(fNonBendingResolution); |
| 1020 | extrapSegment->UpdateFromStationTrackParam(trackParam, mcsFactor, dZ1, dZ2); |
| 1021 | bestChi2 = 5.0; |
| 1022 | bestSegment = NULL; |
| 1023 | if (fPrintLevel >= 10) { |
| 1024 | cout << "FollowTracks: track candidate(0..): " << trackIndex |
| 1025 | << " Look for segment in station(0..): " << station << endl; |
| 1026 | } |
| 1027 | // Loop over segments in station |
| 1028 | for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) { |
| 1029 | // Look for best compatible Segment in station |
| 1030 | // should consider all possibilities ???? |
| 1031 | // multiple scattering ???? |
| 1032 | // separation in 2 functions: Segment and HitForRec ???? |
| 1033 | segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]); |
| 1034 | chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance); |
| 1035 | if (chi2 < bestChi2) { |
| 1036 | // update best Chi2 and Segment if better found |
| 1037 | bestSegment = segment; |
| 1038 | bestChi2 = chi2; |
| 1039 | } |
| 1040 | } |
| 1041 | if (bestSegment) { |
| 1042 | // best segment found: add it to track candidate |
| 1043 | track->AddSegment(bestSegment); |
| 1044 | // set track parameters at these two TrakHit's |
| 1045 | track->SetTrackParamAtHit(track->GetNTrackHits() - 2, &(trackParam[0])); |
| 1046 | track->SetTrackParamAtHit(track->GetNTrackHits() - 1, &(trackParam[1])); |
| 1047 | if (fPrintLevel >= 10) { |
| 1048 | cout << "FollowTracks: track candidate(0..): " << trackIndex |
| 1049 | << " Added segment in station(0..): " << station << endl; |
| 1050 | track->RecursiveDump(); |
| 1051 | } |
| 1052 | } |
| 1053 | else { |
| 1054 | // No best segment found: |
| 1055 | // Look for best compatible HitForRec in station: |
| 1056 | // should consider all possibilities ???? |
| 1057 | // multiple scattering ???? do about like for extrapSegment !!!! |
| 79e1e601 |
1058 | extrapHit = new AliMUONHitForRec(); // empty hit |
| a9e2aefa |
1059 | bestChi2 = 3.0; |
| 1060 | bestHit = NULL; |
| 1061 | if (fPrintLevel >= 10) { |
| 1062 | cout << "FollowTracks: track candidate(0..): " << trackIndex |
| 1063 | << " Segment not found, look for hit in station(0..): " << station |
| 1064 | << endl; |
| 1065 | } |
| 1066 | // Loop over chambers of the station |
| 1067 | for (ch = 0; ch < 2; ch++) { |
| 1068 | // coordinates of extrapolated hit |
| 1069 | extrapHit->SetBendingCoor((&(trackParam[ch]))->GetBendingCoor()); |
| 1070 | extrapHit->SetNonBendingCoor((&(trackParam[ch]))->GetNonBendingCoor()); |
| 1071 | // resolutions from "extrapSegment" |
| 1072 | extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2()); |
| 1073 | extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2()); |
| 1074 | // Loop over hits in the chamber |
| 1075 | for (iHit = fIndexOfFirstHitForRecPerChamber[ch]; |
| 1076 | iHit < fIndexOfFirstHitForRecPerChamber[ch] + |
| 1077 | fNHitsForRecPerChamber[ch]; |
| 1078 | iHit++) { |
| 1079 | hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]); |
| 1080 | // condition for hit not already in segment ???? |
| 1081 | chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance); |
| 1082 | if (chi2 < bestChi2) { |
| 1083 | // update best Chi2 and HitForRec if better found |
| 1084 | bestHit = hit; |
| 1085 | bestChi2 = chi2; |
| 1086 | chBestHit = ch; |
| 1087 | } |
| 1088 | } |
| 1089 | } |
| 1090 | if (bestHit) { |
| 1091 | // best hit found: add it to track candidate |
| 1092 | track->AddHitForRec(bestHit); |
| 1093 | // set track parameters at these two TrakHit's |
| 1094 | track->SetTrackParamAtHit(track->GetNTrackHits() - 1, |
| 1095 | &(trackParam[chBestHit])); |
| 1096 | if (fPrintLevel >= 10) { |
| 1097 | cout << "FollowTracks: track candidate(0..): " << trackIndex |
| 1098 | << " Added hit in station(0..): " << station << endl; |
| 1099 | track->RecursiveDump(); |
| 1100 | } |
| 1101 | } |
| 1102 | else { |
| 1103 | fRecTracksPtr->RemoveAt(trackIndex); // Remove candidate |
| 1104 | break; // stop the search for this candidate: |
| 1105 | // exit from the loop over station |
| 1106 | } |
| 1107 | } |
| 1108 | // Sort track hits according to increasing Z |
| 1109 | track->GetTrackHitsPtr()->Sort(); |
| 1110 | // Update track parameters at first track hit (smallest Z) |
| 1111 | trackParam1 = ((AliMUONTrackHit*) |
| 1112 | (track->GetTrackHitsPtr()->First()))->GetTrackParam(); |
| 1113 | // Track fit from first track hit varying X and Y |
| 1114 | track->Fit(trackParam1, 5); |
| 1115 | if (fPrintLevel >= 10) { |
| 1116 | cout << "FollowTracks: track candidate(0..): " << trackIndex |
| 1117 | << " after fit from station(0..): " << station << " to 4" << endl; |
| 1118 | track->RecursiveDump(); |
| 1119 | } |
| 1120 | delete extrapSegment; |
| 1121 | } // for (station = 2;... |
| 1122 | } // for (trackIndex = 0;... |
| 1123 | return; |
| 1124 | } |
| 1125 | |
| 1126 | void AliMUONEventReconstructor::Streamer(TBuffer &R__b) |
| 1127 | { |
| 1128 | ; |
| 1129 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame