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