MUON/AliMUONRecoEvent.cxx

   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 /* $Id$ */
  17
  18 //Authors: Mihaela Gheata, Andrei Gheata 09/10/00
  19 ////////////////////////////////////////////////////////////////////
  20 //                                                                //
  21 // AliMUONRecoEvent, AliMUONRecoTrack (and AliMUONRecoDisplay)    //
  22 //                                                                //
  23 // Theses classes are used to store and retrieve                  //
  24 // MUON reconstructed events.                                     //
  25 // The corresponding tree is constructed and filled               //
  26 // during the FillEvent() method of AliMUONEventReconstructor,    //
  27 // when all reconstruction information is available.              //
  28 //                                                                //
  29 ////////////////////////////////////////////////////////////////////
  30 ////////////////////////////////////////////////////////////////////
  31 //                                                                //
  32 // AliMUONRecoEvent                                               //
  33 //                                                                //
  34 // This class handles an array of reconstructed tracks.           //
  35 // It provides :                                                  //
  36 //      - filling the tracks array according to the information   //
  37 //        stored in AliMUONEventReconstructor class ;             //
  38 //      - printing event and track informations : event numer,    //
  39 //        number of tracks, hits positions, reconstr. mometum.    //
  40 //                                                                //
  41 ////////////////////////////////////////////////////////////////////
  42 ////////////////////////////////////////////////////////////////////
  43 //                                                                //
  44 // AliMUONRecoTrack                                               //
  45 //                                                                //
  46 // This class represents a reconstructed muon track               //
  47 // in the tree of reconstructed events.                           //
  48 //                                                                //
  49 ////////////////////////////////////////////////////////////////////
  50
  51 #include <Riostream.h>
  52 #include <AliRun.h>
  53 #include <TClonesArray.h>
  54 #include <TClass.h>
  55
  56 #include <TFile.h>
  57 #include <TMatrixD.h>
  58 #include <TParticle.h>
  59
  60 #include "AliMUONRecoEvent.h"
  61 #include "AliMUONEventReconstructor.h"
  62 #include "AliMUONTrack.h"
  63 #include "AliMUONTrackK.h"
  64 #include "AliMUONTrackParam.h"
  65 #include "AliMUONHitForRec.h"
  66 #include "AliMUONTrackHit.h"
  67 #include "AliHeader.h"
  68
  69 ClassImp(AliMUONRecoTrack)
  70 ClassImp(AliMUONRecoEvent)
  71
  72 //-------------------------------------------------------------------
  73 AliMUONRecoEvent::AliMUONRecoEvent(Int_t eventNo)
  74 {
  75 // Reconstructed event constructor
  76    fTracks      = new TClonesArray("AliMUONRecoTrack",200);
  77    fNevr        = eventNo;
  78    fNtracks = 0;
  79 }
  80
  81 //-------------------------------------------------------------------
  82 AliMUONRecoEvent::~AliMUONRecoEvent()
  83 {
  84 // Destructor of AliMUONRecoEvent
  85    fTracks->Delete();
  86    delete fTracks;
  87    fTracks = 0;
  88 }
  89
  90 //-------------------------------------------------------------------
  91 AliMUONRecoTrack* AliMUONRecoEvent::AddEmptyTrack()
  92 {
  93 // Add a empty AliMUONRecoTrackObject to the track list
  94    TClonesArray &dumptracks = *fTracks;
  95    return (new(dumptracks[fNtracks++])AliMUONRecoTrack(kTRUE));
  96 }
  97
  98 //-------------------------------------------------------------------
  99 void AliMUONRecoEvent::Clear(Option_t *option)
 100 {
 101 // Clears all track pointers from the list
 102 //   fTracks->Clear(option);
 103    fTracks->Delete();
 104    fNtracks=0;
 105 }
 106
 107 //-------------------------------------------------------------------
 108 void AliMUONRecoEvent::EventInfo()
 109 {
 110 // Prints reconstructed event information
 111    cout << "*********************Reco Dumper**************************" << endl;
 112    cout << "Event number : " << fNevr << endl;
 113    cout << "   Number of tracks : " << fNtracks << endl;
 114    AliMUONRecoTrack *currentTrack =0;
 115    Int_t trackIndex = 0;
 116    for(trackIndex=0; trackIndex<fNtracks; trackIndex++) {
 117       currentTrack = (AliMUONRecoTrack*)fTracks->UncheckedAt(trackIndex);
 118       cout << "Track : " << trackIndex << endl;
 119       cout << "   Sign : " << currentTrack->GetSign() << endl;
 120       cout << "   Vertex position    : " << currentTrack->GetVertexPos() << endl;
 121       Double_t momreco[3];
 122       for (Int_t axis=0; axis<3; axis++) {
 123          momreco[axis] = currentTrack->GetMomReconstr(axis);
 124       }
 125       cout << "   Reconstructed mom. : " << "Px=" << momreco[0] << "Py=" << momreco[1] << "Pz=" << momreco[2] << endl;
 126       cout << "   Chi squared        : " << currentTrack->GetChi2r() << endl;
 127       cout << "   Hits positions     : " << endl;
 128       Double_t xhit, yhit, zhit;
 129       for (Int_t chamber=0; chamber<10; chamber++) {
 130          xhit = currentTrack->GetPosX(chamber);
 131          yhit = currentTrack->GetPosY(chamber);
 132          zhit = currentTrack->GetPosZ(chamber);
 133 //         cout <<"      chamber" << chamber << " X=" << xhit << " Y=" << yhit << " Z=" << zhit << endl;
 134       }
 135    }
 136    cout << "**********************************************************" << endl;
 137 }
 138
 139 //-------------------------------------------------------------------
 140 Bool_t AliMUONRecoEvent::MakeDumpTracks(Int_t muons, TClonesArray *tracksPtr,
 141   AliMUONEventReconstructor *EventReco)
 142 {
 143 // This method takes the pointer of the list of reconstructed tracks from
 144 // AliMUONEventReconstructor and fill the reconstructed AliMUONRecoEvent
 145 // fields.
 146
 147         cout << "Enter MakeDumpTracks..." << endl;
 148    Int_t nTracks = tracksPtr->GetEntriesFast();
 149    cout << "nTracks = "<< nTracks << endl;
 150    if (nTracks == 0) {
 151       cout << "AliMUONRecoEvent::MakeDumpTracks: Number of tracks is zero !" << endl;
 152       //AZ return kFALSE;
 153    }
 154    cout << tracksPtr << endl;
 155    if (!tracksPtr) {
 156       cout << "AliMUONRecoEvent::MakeDumpTracks() : You should call SetRecoTracksPtr() first..." << endl;
 157       return kFALSE;
 158    }
 159         // Get event number
 160    Int_t noEvent = gAlice->GetHeader()->GetEvent();
 161    cout << "noEvent = "<< nTracks << endl;
 162    tracksPtr->Compress();  // simple loop
 163    AliMUONRecoTrack *currentTrack;
 164    Int_t trackIndex, nTrackHits = 0;
 165    Double_t z, pYZ, bendingSlope, nonBendingSlope;
 166    Double_t pX, pY, pZ;                 // reconstructed momentum components
 167    Int_t isign, flag=0;         // charge sign, flag of reconstructed track
 168    Double_t alpha, beta;
 169    TObjArray *hitsOnTrack = 0;
 170    AliMUONTrackHit *trackHit = 0;
 171    AliMUONTrack *track = 0;
 172    AliMUONTrackK *trackK = 0;
 173    TMatrixD *trackParamK; //AZnon
 174    AliMUONTrackParam *trackParam = 0;
 175    // Fill event number and number of tracks
 176    fNevr = noEvent;
 177    fMuons = muons; //AZ - number of muons within acceptance
 178    // Loop over reconstructed tracks
 179    for (trackIndex=0; trackIndex<nTracks; trackIndex++) {
 180       cout << " trackIndex = " << trackIndex << endl;
 181       currentTrack = AddEmptyTrack();
 182       cout << " currentTrack = " << currentTrack << endl;
 183
 184       if (EventReco->GetTrackMethod() == 2) { // Kalman
 185
 186         trackK = (AliMUONTrackK*) ((*tracksPtr)[trackIndex]);
 187         nTrackHits = trackK->GetNTrackHits();
 188         trackParamK = trackK->GetTrackParameters();
 189         isign = Int_t(TMath::Sign(1., (*trackParamK)(4,0)));
 190         z = trackK->GetZ();
 191         alpha = (*trackParamK)(2,0);
 192         beta = (*trackParamK)(3,0);
 193         pYZ = TMath::Cos(beta)/TMath::Abs((*trackParamK)(4,0));
 194         pZ = pYZ/TMath::Cos(alpha);
 195         pX = TMath::Sin(beta)/TMath::Abs((*trackParamK)(4,0));
 196         pY = pYZ*TMath::Sin(alpha);
 197
 198         currentTrack->SetVertexPos(z);
 199         currentTrack->SetMomReconstr(pX,pY,pZ);
 200         currentTrack->SetSign(isign);
 201         currentTrack->SetChi2r(trackK->GetTrackQuality());
 202
 203         // Check hits on the track
 204         hitsOnTrack = trackK->GetHitOnTrack();
 205         Float_t signal = 0;
 206         Float_t tht = 0;
 207         for (int ihit = 0; ihit < nTrackHits; ihit++) {
 208           signal += ((AliMUONHitForRec*)((*hitsOnTrack)[ihit]))->GetGeantSignal();
 209           tht += TMath::Min (1,((AliMUONHitForRec*)((*hitsOnTrack)[ihit]))->GetTHTrack());
 210         }
 211         signal /= nTrackHits;
 212         tht /= nTrackHits;
 213         flag = 0;
 214         if (TMath::Nint(signal) > 0) { // signal muon
 215           for (int ihit = 0; ihit < nTrackHits ; ihit++) {
 216             if (((AliMUONHitForRec*)((*hitsOnTrack)[ihit]))->GetTHTrack() != TMath::Nint(tht)) flag++;
 217           }
 218         } else flag = -9; // background track
 219         //cout << TMath::Nint(signal) << " " << TMath::Nint(tht) << " " << recTrackNt->fFlag << endl;
 220         currentTrack->SetFlag(flag);
 221       } else { // default tracking
 222
 223         track = (AliMUONTrack*) ((*tracksPtr)[trackIndex]);
 224         nTrackHits = track->GetNTrackHits();
 225         // track parameters at Vertex
 226         trackParam = track->GetTrackParamAtVertex();
 227         bendingSlope = trackParam->GetBendingSlope();
 228         nonBendingSlope = trackParam->GetNonBendingSlope();
 229
 230         z = trackParam->GetZ();
 231         pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
 232         pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
 233         pX = pZ * nonBendingSlope;
 234         pY = pZ * bendingSlope;
 235
 236         if (trackParam->GetInverseBendingMomentum()<0) isign=-1; else isign=1;
 237         currentTrack->SetVertexPos(z);
 238         currentTrack->SetMomReconstr(pX,pY,pZ);
 239         currentTrack->SetSign(isign);
 240         //         currentTrack->SetChi2r(trackParam->GetChi2());
 241         currentTrack->SetChi2r(0);
 242
 243         // Check hits on the track
 244         hitsOnTrack = track->GetTrackHitsPtr();
 245         Float_t signal = 0;
 246         Float_t tht = 0;
 247         AliMUONHitForRec *hitForRec = 0;
 248         for (int ihit = 0; ihit < nTrackHits; ihit++) {
 249           hitForRec = ((AliMUONTrackHit*)(*hitsOnTrack)[ihit])->GetHitForRecPtr();
 250           signal += hitForRec->GetGeantSignal();
 251           tht += TMath::Min (1,hitForRec->GetTHTrack());
 252         }
 253         signal /= nTrackHits;
 254         tht /= nTrackHits;
 255         flag = 0;
 256         if (TMath::Nint(signal) > 0) { // signal muon
 257           for (int ihit = 0; ihit < nTrackHits ; ihit++) {
 258             hitForRec = ((AliMUONTrackHit*)(*hitsOnTrack)[ihit])->GetHitForRecPtr();
 259             if (hitForRec->GetTHTrack() != TMath::Nint(tht)) flag++;
 260           }
 261         } else flag = -9; // background track
 262         //cout << TMath::Nint(signal) << " " << TMath::Nint(tht) << " " << recTrackNt->fFlag << endl;
 263         currentTrack->SetFlag(flag);
 264       }
 265
 266       Double_t xhit,yhit,zhit;
 267       // Loop over track hits
 268       for (Int_t trackHitIndex = 0; trackHitIndex < nTrackHits; trackHitIndex++) {
 269         if (EventReco->GetTrackMethod() == 2) { // Kalman
 270           xhit = ((AliMUONHitForRec*)((*hitsOnTrack)[trackHitIndex]))->GetNonBendingCoor();
 271           yhit = ((AliMUONHitForRec*)((*hitsOnTrack)[trackHitIndex]))->GetBendingCoor();
 272           zhit = ((AliMUONHitForRec*)((*hitsOnTrack)[trackHitIndex]))->GetZ();
 273         } else {
 274           trackHit = (AliMUONTrackHit*) (*(track->GetTrackHitsPtr()))[trackHitIndex];
 275           xhit = trackHit->GetHitForRecPtr()->GetNonBendingCoor();
 276           yhit = trackHit->GetHitForRecPtr()->GetBendingCoor();
 277           zhit = trackHit->GetHitForRecPtr()->GetZ();
 278         }
 279         if (trackHitIndex >= 0 && trackHitIndex < 10) {
 280           currentTrack->SetHitPosition(trackHitIndex,xhit,yhit,zhit);
 281         } else { cout << "track " << trackIndex << " hit out of range" << endl;}
 282       }
 283
 284    }
 285    cout << "Leave MakeDumpTracks..." << endl;
 286    return kTRUE;
 287 }
 288
 289 //-------------------------------------------------------------------
 290 void AliMUONRecoEvent::Streamer(TBuffer &R__b)
 291 {
 292 // Streams an object of class AliMUONRecoEvent
 293    if (R__b.IsReading()) {
 294       fTracks->Clear();
 295       AliMUONRecoEvent::Class()->ReadBuffer(R__b, this);
 296    } else {
 297       cout << "...writing event to file...\n";
 298       AliMUONRecoEvent::Class()->WriteBuffer(R__b, this);
 299    }
 300 }
 301
 302 //-------------------------------------------------------------------
 303 AliMUONRecoTrack::AliMUONRecoTrack(Bool_t active)
 304 {
 305 //Constructor of AliMUONRecoTrack
 306    fSign  = 0;
 307    fZvr   = 0.0;
 308    fChi2r = 0.0;
 309    if (active) {
 310         for (Int_t axis=0; axis<3; axis++) {
 311         fPr[axis] = 0.0;
 312         }
 313         for (Int_t chamber=0; chamber<10; chamber++) {
 314         fPosX[chamber] = 0.0;
 315         fPosY[chamber] = 0.0;
 316         fPosZ[chamber] = 0.0;
 317         }
 318    }
 319 }
 320
 321 //-------------------------------------------------------------------
 322 const Double_t AliMUONRecoTrack::Phi()
 323 {
 324 // Return trach phi angle
 325         return TMath::ATan2(fPr[2], fPr[1]);
 326 }
 327
 328 //-------------------------------------------------------------------
 329 const Double_t AliMUONRecoTrack::Theta()
 330 {
 331 // Return trach theta angle
 332    return TMath::ACos(fPr[2] / P());
 333 }
 334
 335 //-------------------------------------------------------------------
 336 void AliMUONRecoTrack::SetMomReconstr(Double_t px, Double_t py, Double_t pz)
 337 {
 338 // Set the track reconstructed momentum
 339    fPr[0] = px;
 340    fPr[1] = py;
 341    fPr[2] = pz;
 342 }
 343
 344 //-------------------------------------------------------------------
 345 void AliMUONRecoTrack::SetHitPosition(Int_t chamber, Double_t x, Double_t y, Double_t z)
 346 {
 347 // Set hit coordinates in chamber[0..9]
 348    fPosX[chamber] = x;
 349    fPosY[chamber] = y;
 350    fPosZ[chamber] = z;
 351 }
 352 //-------------------------------------------------------------------
 353 void AliMUONRecoTrack::TrackInfo()
 354 {
 355 // Prints momentum info for this track
 356    cout << "Px=" << GetMomReconstr(0) << " Py=" << GetMomReconstr(1) <<
 357           " Pz=" << GetMomReconstr(2) << " P=" << P() << endl;
 358 }