EVE/alice-macros/muon_trackRefs.C

   1 // $Id$
   2
   3 /**************************************************************************
   4  * Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. *
   5  * See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for          *
   6  * full copyright notice.                                                 *
   7  **************************************************************************/
   8
   9 /// \ingroup evemacros
  10 /// \file muon_trackRefs.C
  11 /// \brief Macro to visualise trackRef in MUON spectrometer
  12 /// (both tracker and trigger).
  13 ///
  14 /// Use muon_trackRefs(Bool_t showSimClusters) in order to run it
  15 ///
  16 /// Needs that alieve_init() is already called
  17 ///
  18 /// \author P. Pillot, L. Aphecetche; Subatech
  19
  20 #if !defined(__CINT__) || defined(__MAKECINT__)
  21
  22 #include "AliMUONClusterStoreV2.h"
  23 #include "AliMUONRawClusterV2.h"
  24 #include "AliMUONVCluster.h"
  25 #include "AliMUONConstants.h"
  26 #include "AliMUONRecoParam.h"
  27 #include "AliMUONCDB.h"
  28
  29 #include "AliRunLoader.h"
  30 #include "AliStack.h"
  31 #include "AliTrackReference.h"
  32
  33 #include "EveBase/AliEveMagField.h"
  34 #include "EveBase/AliEveTrack.h"
  35 #include "EveBase/AliEveEventManager.h"
  36
  37 #include <TEveManager.h>
  38 #include <TEveUtil.h>
  39 #include <TEveTrack.h>
  40 #include <TEvePointSet.h>
  41 #include <TEveVSDStructs.h>
  42 #include <TEveTrackPropagator.h>
  43
  44 #include <TClonesArray.h>
  45 #include <TTree.h>
  46 #include <TParticle.h>
  47 #include <TMath.h>
  48 #include <TROOT.h>
  49
  50 #endif
  51
  52 //______________________________________________________________________________
  53 void muon_trackRef_propagator_setup(TEveTrackPropagator* trkProp, Bool_t showVertex)
  54 {
  55   // set magnetic field
  56   trkProp->SetMagFieldObj(new AliEveMagField);
  57   trkProp->SetStepper(TEveTrackPropagator::kRungeKutta);
  58
  59   // set propagation range
  60   trkProp->SetMaxR(500.);
  61   trkProp->SetMaxZ(-AliMUONConstants::DefaultChamberZ(13)+3.5);
  62
  63   // go through pathmarks
  64   trkProp->SetFitDaughters(kFALSE);
  65   trkProp->SetFitReferences(kTRUE);
  66   trkProp->SetFitDecay(kTRUE);
  67   trkProp->SetFitCluster2Ds(kFALSE);
  68
  69   // Render first vertex if required
  70   if (showVertex)
  71   {
  72     trkProp->SetRnrFV(kTRUE);
  73     trkProp->RefFVAtt().SetMarkerSize(0.5);
  74     trkProp->RefFVAtt().SetMarkerColor(kRed);
  75   }
  76 }
  77
  78 //______________________________________________________________________________
  79 Bool_t isReconstructible(Bool_t* chHit)
  80 {
  81   // load recoParam from OCDB
  82   static AliMUONRecoParam* gRecoParam = 0x0;
  83   static UInt_t gRequestedStationMask = 0;
  84   static Bool_t gRequest2ChInSameSt45 = kFALSE;
  85   if (!gRecoParam)
  86   {
  87     gRecoParam = AliMUONCDB::LoadRecoParam();
  88     if (!gRecoParam) exit(-1);
  89     // compute the mask of requested stations
  90     gRequestedStationMask = 0;
  91     for (Int_t i = 0; i < 5; i++) if (gRecoParam->RequestStation(i)) gRequestedStationMask |= ( 1 << i );
  92     // get whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
  93     gRequest2ChInSameSt45 = !gRecoParam->MakeMoreTrackCandidates();
  94   }
  95
  96   // check which chambers are hit
  97   UInt_t presentStationMask(0);
  98   Int_t nChHitInSt4 = 0, nChHitInSt5 = 0;
  99   for (Int_t ich=0; ich<10; ich++)
 100   {
 101     if (!chHit[ich]) continue;
 102     Int_t ist = ich/2;
 103     presentStationMask |= ( 1 << ist );
 104     if (ist == 3) nChHitInSt4++;
 105     if (ist == 4) nChHitInSt5++;
 106   }
 107
 108   // at least one cluster per requested station
 109   if ((gRequestedStationMask & presentStationMask) != gRequestedStationMask) return kFALSE;
 110
 111   // 2 chambers hit in the same station (4 or 5)
 112   if (gRequest2ChInSameSt45) return (nChHitInSt4 == 2 || nChHitInSt5 == 2);
 113   // or 2 chambers hit in station 4 & 5 together
 114   else return (nChHitInSt4+nChHitInSt5 >= 2);
 115 }
 116
 117 //______________________________________________________________________________
 118 void add_muon_trackRefs(AliStack* stack, TTree* treeTR, TEveTrackList* reco, TEveTrackList* other,
 119                         TEvePointSet* RecoClusters, TEvePointSet* OtherClusters)
 120 {
 121   TClonesArray* trackRefs = 0;
 122   treeTR->SetBranchAddress("TrackReferences", &trackRefs);
 123   Int_t nSimTracks = stack->GetNtrack();
 124   AliMUONClusterStoreV2 clusters;
 125   TClonesArray clustersTrigger("AliMUONRawClusterV2", 10);
 126
 127   // loop over simulated track
 128   for (Int_t itr = 0; itr < nSimTracks; itr++)
 129   {
 130     treeTR->GetEntry(stack->TreeKEntry(itr));
 131     Int_t nTrackRefs = trackRefs->GetEntriesFast();
 132     if (nTrackRefs <= 0) continue;
 133
 134     TEveTrack* track = 0x0;
 135     Bool_t chHit[10];
 136     for (Int_t ich=0; ich<10; ich++) chHit[ich] = kFALSE;
 137
 138     // loop over simulated track hits
 139     for (Int_t itrR = 0; itrR < nTrackRefs; ++itrR)
 140     {
 141       AliTrackReference* atr = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(itrR));
 142
 143       // skip trackRefs not in MUON
 144       if (atr->DetectorId() != AliTrackReference::kMUON) continue;
 145
 146       // record chamber hit
 147       Int_t detElemId = atr->UserId();
 148       Int_t chamberId = detElemId / 100 - 1;
 149       if (chamberId < 0) continue;
 150       if (chamberId < 10) chHit[chamberId] = kTRUE;
 151
 152       // produce eve track if not already done
 153       if (!track)
 154       {
 155         TParticle* p = stack->Particle(itr);
 156         track = new AliEveTrack(p, itr, 0x0);
 157         track->SetName(Form("%s [%d]", p->GetName(), itr));
 158         track->SetStdTitle();
 159         track->SetSourceObject(p);
 160
 161         clusters.Clear();
 162         clustersTrigger.Clear("C");
 163       }
 164
 165       // add path mark
 166       track->AddPathMark(TEvePathMark(TEvePathMark::kReference,
 167                                       TEveVector(atr->X(),  atr->Y(),  atr->Z()),
 168                                       TEveVector(atr->Px(), atr->Py(), atr->Pz()),
 169                                       atr->GetTime()));
 170
 171       // produce clusters if required
 172       if (RecoClusters || OtherClusters)
 173       {
 174         // from tracker
 175         if (chamberId < 10)
 176         {
 177           // produce a new cluster on that DE or update existing one
 178           AliMUONVCluster* cl = 0x0;
 179           Int_t clNum = -1;
 180           do cl = clusters.FindObject(AliMUONVCluster::BuildUniqueID(chamberId, detElemId, ++clNum));
 181           while (cl && cl->GetNDigits() == 2);
 182           if (cl) {
 183             cl->SetXYZ((cl->GetX() + atr->X()) / 2., (cl->GetY() + atr->Y()) / 2., (cl->GetZ() + atr->Z()) / 2.);
 184           }
 185           else
 186           {
 187             cl = clusters.Add(chamberId, detElemId, clNum);
 188             cl->SetXYZ(atr->X(), atr->Y(), atr->Z());
 189           }
 190           cl->AddDigitId((UInt_t)itrR);
 191         }
 192         else // from trigger
 193         {
 194           AliMUONVCluster* cl = new(clustersTrigger[clustersTrigger.GetLast()+1]) AliMUONRawClusterV2();
 195           cl->SetXYZ(atr->X(), atr->Y(), atr->Z());
 196         }
 197       }
 198     }
 199
 200     if (track)
 201     {
 202       track->SortPathMarksByTime();
 203       // stop track propagation at last path mark
 204       track->RefPathMarks().back().fType = TEvePathMark::kDecay;
 205
 206       // add the track and trackRefs to proper lists
 207       if (isReconstructible(chHit)) {
 208         track->SetPropagator(reco->GetPropagator());
 209         track->SetAttLineAttMarker(reco);
 210         reco->AddElement(track);
 211
 212         // trackRefs
 213         if (RecoClusters)
 214         {
 215           // tracker
 216           TIter next(clusters.CreateIterator());
 217           gROOT->ProcessLine(Form("add_muon_clusters((TIter*)%p, (TEvePointSet*)%p);",&next, RecoClusters));
 218           // trigger
 219           TIter next2(clustersTrigger.MakeIterator());
 220           gROOT->ProcessLine(Form("add_muon_clusters((TIter*)%p, (TEvePointSet*)%p);",&next2, RecoClusters));
 221         }
 222       }
 223       else {
 224         track->SetPropagator(other->GetPropagator());
 225         track->SetAttLineAttMarker(other);
 226         other->AddElement(track);
 227
 228         // trackRefs
 229         if (OtherClusters)
 230         {
 231           // tracker
 232           TIter next(clusters.CreateIterator());
 233           gROOT->ProcessLine(Form("add_muon_clusters((TIter*)%p, (TEvePointSet*)%p);",&next, OtherClusters));
 234           // trigger
 235           TIter next2(clustersTrigger.MakeIterator());
 236           gROOT->ProcessLine(Form("add_muon_clusters((TIter*)%p, (TEvePointSet*)%p);",&next2, OtherClusters));
 237         }
 238       }
 239     }
 240   }
 241   delete trackRefs;
 242 }
 243
 244 //______________________________________________________________________________
 245 void muon_trackRefs(Bool_t showSimClusters)
 246 {
 247   // load kinematics and trackRefs
 248   AliRunLoader* rl =  AliEveEventManager::AssertRunLoader();
 249   rl->LoadKinematics();
 250   AliStack* stack = rl->Stack();
 251   if (!stack) return;
 252   rl->LoadTrackRefs();
 253   TTree* treeTR = rl->TreeTR();
 254   if (!treeTR) return;
 255
 256   // track containers
 257   TEveElementList* trackCont = new TEveElementList("Sim MUON Tracks");
 258
 259   TEveTrackList* reco = new TEveTrackList("reconstructible");
 260   reco->SetRnrPoints(kFALSE);
 261   reco->SetRnrLine(kTRUE);
 262   reco->SetLineColor(kRed);
 263   reco->SetLineStyle(2);
 264   muon_trackRef_propagator_setup(reco->GetPropagator(), kTRUE);
 265   trackCont->AddElement(reco);
 266
 267   TEveTrackList* other = new TEveTrackList("others");
 268   other->SetRnrPoints(kFALSE);
 269   other->SetRnrLine(kTRUE);
 270   other->SetLineColor(kWhite);
 271   other->SetLineStyle(3);
 272   muon_trackRef_propagator_setup(other->GetPropagator(), kFALSE);
 273   trackCont->AddElement(other);
 274
 275   // cluster container
 276   TEveElementList* clusterCont = 0x0;
 277   TEvePointSet* RecoClusters = 0x0;
 278   TEvePointSet* OtherClusters = 0x0;
 279   if (showSimClusters)
 280   {
 281     clusterCont = new TEveElementList("Sim MUON Clusters");
 282
 283     RecoClusters = new TEvePointSet(1000);
 284     RecoClusters->SetName("Reconstructibles");
 285     RecoClusters->SetPickable(kFALSE);
 286     RecoClusters->SetMarkerStyle(2);
 287     RecoClusters->SetMarkerColor(kRed);
 288     RecoClusters->SetMarkerSize(0.5);
 289     clusterCont->AddElement(RecoClusters);
 290
 291     OtherClusters = new TEvePointSet(10000);
 292     OtherClusters->SetName("Others");
 293     OtherClusters->SetPickable(kFALSE);
 294     OtherClusters->SetMarkerStyle(2);
 295     OtherClusters->SetMarkerColor(kWhite);
 296     OtherClusters->SetMarkerSize(0.5);
 297     clusterCont->AddElement(OtherClusters);
 298
 299     TEveUtil::LoadMacro("muon_clusters.C+");
 300   }
 301
 302   // add tracks to the proper list and propagate them. Add also clusters if required.
 303   add_muon_trackRefs(stack, treeTR, reco, other, RecoClusters, OtherClusters);
 304   reco->MakeTracks();
 305   other->MakeTracks();
 306   trackCont->SetTitle(Form("N=%d", reco->NumChildren()+other->NumChildren()));
 307   reco->SetTitle(Form("N=%d",reco->NumChildren()));
 308   other->SetTitle(Form("N=%d",other->NumChildren()));
 309   if (showSimClusters)
 310   {
 311     clusterCont->SetTitle(Form("N=%d",RecoClusters->GetLastPoint()+OtherClusters->GetLastPoint()+2));
 312     RecoClusters->SetTitle(Form("N=%d",RecoClusters->GetLastPoint()+1));
 313     OtherClusters->SetTitle(Form("N=%d",OtherClusters->GetLastPoint()+1));
 314   }
 315
 316   // add graphic containers
 317   gEve->DisableRedraw();
 318   gEve->AddElement(trackCont);
 319   if (clusterCont) gEve->AddElement(clusterCont);
 320   gEve->EnableRedraw();
 321   gEve->Redraw3D();
 322 }