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