+// for (Int_t iPart = 0; iPart < np; iPart++) {
+// p = gAlice->Particle(iPart);
+// printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
+// iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
+// }
+ return;
+}
+
+void AliMUONEventReconstructor::FillEvent()
+{
+// Create a new AliMUONRecoEvent, fill its track list, then add it as a
+// leaf in the Event branch of TreeRecoEvent tree
+ cout << "Enter FillEvent() ...\n";
+
+ if (!fRecoEvent) {
+ fRecoEvent = new AliMUONRecoEvent();
+ } else {
+ fRecoEvent->Clear();
+ }
+ //save current directory
+ TDirectory *current = gDirectory;
+ if (!fTreeFile) fTreeFile = new TFile("tree_reco.root", "RECREATE");
+ if (!fEventTree) fEventTree = new TTree("TreeRecoEvent", "MUON reconstructed events");
+ //AZif (fRecoEvent->MakeDumpTracks(fRecTracksPtr)) {
+ if (fRecoEvent->MakeDumpTracks(fMuons, fRecTracksPtr, this)) { //AZ
+ if (fPrintLevel > 1) fRecoEvent->EventInfo();
+ TBranch *branch = fEventTree->GetBranch("Event");
+ if (!branch) branch = fEventTree->Branch("Event", "AliMUONRecoEvent", &fRecoEvent, 64000);
+ branch->SetAutoDelete();
+ fTreeFile->cd();
+ fEventTree->Fill();
+ fTreeFile->Write();
+ }
+ // restore directory
+ current->cd();
+}
+
+//__________________________________________________________________________
+void AliMUONEventReconstructor::MakeTrackCandidatesK(void)
+{
+ // To make initial tracks for Kalman filter from the list of segments
+ Int_t istat, iseg;
+ AliMUONSegment *segment;
+ AliMUONTrackK *trackK;
+
+ if (fPrintLevel >= 1) cout << "enter MakeTrackCandidatesK" << endl;
+ // Reset the TClonesArray of reconstructed tracks
+ if (fRecTracksPtr) fRecTracksPtr->Delete();
+ // Delete in order that the Track destructors are called,
+ // hence the space for the TClonesArray of pointers to TrackHit's is freed
+ fNRecTracks = 0;
+
+ AliMUONTrackK a(this, fHitsForRecPtr); // bad idea ???
+ // Loop over stations(1...) 5 and 4
+ for (istat=4; istat>=3; istat--) {
+ // Loop over segments in the station
+ for (iseg=0; iseg<fNSegments[istat]; iseg++) {
+ // Transform segments to tracks and evaluate covariance matrix
+ segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]);
+ trackK = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrackK(segment);
+ fNRecTracks++;
+ } // for (iseg=0;...)
+ } // for (istat=4;...)
+ return;
+}
+
+//__________________________________________________________________________
+void AliMUONEventReconstructor::FollowTracksK(void)
+{
+ // Follow tracks using Kalman filter
+ Bool_t Ok;
+ Int_t icand, ichamBeg, ichamEnd, chamBits;
+ Double_t zDipole1, zDipole2;
+ AliMUONTrackK *trackK;
+ AliMUONHitForRec *hit;
+ AliMUONRawCluster *clus;
+ TClonesArray *rawclusters;
+ AliMUON *pMUON;
+ clus = 0; rawclusters = 0;
+
+ zDipole1 = GetSimpleBPosition() - GetSimpleBLength()/2;
+ zDipole2 = zDipole1 + GetSimpleBLength();
+
+ // Print hits
+ pMUON = (AliMUON*) gAlice->GetModule("MUON");
+ for (Int_t i1=0; i1<fNHitsForRec; i1++) {
+ hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]);
+ //if (hit->GetTHTrack() > 1 || hit->GetGeantSignal() == 0) continue;
+ /*
+ cout << " Hit #" << hit->GetChamberNumber() << " ";
+ cout << hit->GetBendingCoor() << " ";
+ cout << hit->GetNonBendingCoor() << " ";
+ cout << hit->GetZ() << " ";
+ cout << hit->GetGeantSignal() << " ";
+ cout << hit->GetTHTrack() << endl;
+ */
+ /*
+ printf(" Hit # %d %10.4f %10.4f %10.4f",
+ hit->GetChamberNumber(), hit->GetBendingCoor(),
+ hit->GetNonBendingCoor(), hit->GetZ());
+ if (fRecGeantHits) {
+ // from GEANT hits
+ printf(" %3d %3d \n", hit->GetGeantSignal(), hit->GetTHTrack());
+ } else {
+ // from raw clusters
+ rawclusters = pMUON->RawClustAddress(hit->GetChamberNumber());
+ clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
+ GetHitNumber());
+ printf("%3d", clus->fTracks[1]-1);
+ if (clus->fTracks[2] != 0) printf("%3d \n", clus->fTracks[2]-1);
+ else printf("\n");
+ }
+ */
+ }
+
+ icand = -1;
+ Int_t nSeeds = fNRecTracks; // starting number of seeds
+ // Loop over track candidates
+ while (icand < fNRecTracks-1) {
+ icand ++;
+ trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
+
+ // Discard candidate which will produce the double track
+ if (icand > 0) {
+ Ok = CheckCandidateK(icand,nSeeds);
+ if (!Ok) {
+ //trackK->SetRecover(-1); // mark candidate to be removed
+ //continue;
+ }
+ }
+
+ Ok = kTRUE;
+ if (trackK->GetRecover() == 0) hit = (AliMUONHitForRec*)
+ trackK->GetHitOnTrack()->Last(); // last hit
+ else hit = (AliMUONHitForRec*) (*trackK->GetHitOnTrack())[1]; // 2'nd hit
+ ichamBeg = hit->GetChamberNumber();
+ ichamEnd = 0;
+ // Check propagation direction
+ if (trackK->GetTrackDir() > 0) {
+ ichamEnd = 9; // forward propagation
+ Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
+ if (Ok) {
+ ichamBeg = ichamEnd;
+ ichamEnd = 6; // backward propagation
+ // Change weight matrix and zero fChi2 for backpropagation
+ trackK->StartBack();
+ Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
+ ichamBeg = ichamEnd;
+ ichamEnd = 0;
+ }
+ } else {
+ if (trackK->GetBPFlag()) {
+ // backpropagation
+ ichamEnd = 6; // backward propagation
+ // Change weight matrix and zero fChi2 for backpropagation
+ trackK->StartBack();
+ Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
+ ichamBeg = ichamEnd;
+ ichamEnd = 0;
+ }
+ }
+
+ if (Ok) {
+ trackK->SetTrackDir(-1);
+ trackK->SetBPFlag(kFALSE);
+ Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
+ }
+ if (Ok) trackK->SetTrackQuality(0); // compute "track quality"
+ else trackK->SetRecover(-1); // mark candidate to be removed
+
+ // Majority 3 of 4 in first 2 stations
+ chamBits = 0;
+ for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
+ hit = (AliMUONHitForRec*) (*trackK->GetHitOnTrack())[i];
+ chamBits |= BIT(hit->GetChamberNumber()-1);
+ }
+ //if (!((chamBits&3)==3 || (chamBits>>2&3)==3)) trackK->SetRecover(-1);
+ //mark candidate to be removed
+ } // while
+
+ for (Int_t i=0; i<fNRecTracks; i++) {
+ trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
+ if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);