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