EVE/alice-macros/MUON_digits.C

   1 #include "TGLViewer.h"
   2
   3 namespace Alieve {
   4 class MUONData;
   5 class Event;
   6 }
   7
   8 Alieve::MUONData* g_muon_data = 0;
   9 Alieve::Event*    g_muon_last_event = 0;
  10
  11 Int_t g_currentEvent = -1;
  12 Bool_t g_fromRaw = kFALSE;
  13
  14 void MUON_digits(Bool_t fromRaw = kFALSE, Bool_t showTracks = kTRUE)
  15 {
  16
  17   TTree* dt = 0;
  18
  19   if (Alieve::gEvent == 0) {
  20     printf("No alieve event: use alieve_init(...) \n");
  21     return;
  22   }
  23
  24   if (g_currentEvent == Alieve::gEvent->GetEventId()) {
  25     if (g_fromRaw == fromRaw) {
  26       printf("Same event... \n");
  27       return;
  28     } else {
  29       if (g_fromRaw) {
  30         printf("Same event with digits.\n");
  31         Alieve::gEvent->GotoEvent(g_currentEvent);
  32       } else {
  33         printf("Same event with raw.\n");
  34         Alieve::gEvent->GotoEvent(g_currentEvent);
  35       }
  36     }
  37   }
  38
  39   g_fromRaw = fromRaw;
  40
  41   AliRunLoader* rl =  Alieve::Event::AssertRunLoader();
  42   TString fileName = rl->GetFileName();
  43   Int_t length = fileName.Length();
  44   fileName.Resize(length-11);
  45   fileName.Append("raw.root");
  46
  47   g_muon_data = new Alieve::MUONData;
  48
  49   if (!fromRaw) {
  50     rl->LoadDigits("MUON");
  51     dt = rl->GetTreeD("MUON", false);
  52     if (dt == 0) {
  53       cout << "No digits produced!" << endl;
  54     } else {
  55       cout << "Display aliroot digits!" << endl;
  56       g_muon_data->LoadDigits(dt);
  57     }
  58   } else {
  59     if (gSystem->AccessPathName(fileName.Data(),kFileExists)) {
  60       cout << "No raw data produced!" << endl;
  61     } else {
  62       cout << "Display raw digits!" << endl;
  63       g_muon_data->LoadRaw(fileName.Data());
  64     }
  65   }
  66
  67   g_muon_last_event = Alieve::gEvent;
  68
  69   g_currentEvent = g_muon_last_event->GetEventId();
  70
  71   gStyle->SetPalette(1, 0);
  72
  73   gReve->DisableRedraw();
  74
  75   Reve::RenderElementList* l = new Reve::RenderElementList("MUONChambers");
  76   l->SetTitle("MUON chambers");
  77   l->SetMainColor(Color_t(2));
  78   gReve->AddRenderElement(l);
  79
  80   for (Int_t ic = 0; ic < 14; ic++) {
  81
  82     Alieve::MUONChamber* mucha = new Alieve::MUONChamber();
  83
  84     mucha->SetFrameColor(2);
  85     mucha->SetChamberID(ic);
  86     mucha->SetDataSource(g_muon_data);
  87
  88     gReve->AddRenderElement(l,mucha);
  89
  90   }
  91
  92   if (showTracks) MUON_tracks();
  93
  94   gReve->EnableRedraw();
  95
  96   TGLViewer* view = dynamic_cast<TGLViewer*>(gReve->GetGLCanvas()->GetViewer3D());
  97   view->ResetCamerasAfterNextUpdate();
  98   gReve->GetGLCanvas()->Modified();
  99   gReve->GetGLCanvas()->Update();
 100
 101 }
 102
 103 //_____________________________________________________________________________
 104 void MUON_tracks() {
 105
 106   AliRunLoader* rl =  Alieve::Event::AssertRunLoader();
 107   rl->LoadTracks("MUON");
 108   TTree* tt = rl->GetTreeT("MUON", false);
 109
 110   TClonesArray *tracks = 0;
 111   tt->SetBranchAddress("MUONTrack",&tracks);
 112   tt->GetEntry(0);
 113
 114   Int_t ntracks = tracks->GetEntriesFast();
 115   //printf("Found %d tracks. \n",ntracks);
 116
 117   Reve::TrackList* lt = new Reve::TrackList("M-Tracks");
 118   lt->SetMainColor(Color_t(6));
 119
 120   gReve->AddRenderElement(lt);
 121
 122   TMatrixD smatrix(2,2);
 123   TMatrixD sums(2,1);
 124   TMatrixD res(2,1);
 125
 126   Float_t xRec, xRec0;
 127   Float_t yRec, yRec0;
 128   Float_t zRec, zRec0;
 129
 130   Float_t zg[4] = { -1603.5, -1620.5, -1703.5, -1720.5 };
 131
 132   AliMUONTrack *mt;
 133   Reve::RecTrack  rt;
 134   Int_t count;
 135   for (Int_t n = 0; n < ntracks; n++) {
 136
 137     count = 0;
 138
 139     mt = (AliMUONTrack*) tracks->At(n);
 140
 141     //printf("Match trigger %d \n",mt->GetMatchTrigger());
 142
 143     rt.label = n;
 144
 145     Reve::Track* track = new Reve::Track(&rt, lt->GetRnrStyle());
 146
 147     //PH The line below is replaced waiting for a fix in Root
 148     //PH which permits to use variable siza arguments in CINT
 149     //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
 150     char form[1000];
 151     if (mt->GetMatchTrigger()) {
 152       //PH      track->SetName(Form("MUONTrack %2d (MT)", rt.label));
 153       sprintf(form,"MUONTrack %2d (MT)", rt.label);
 154       track->SetName(form);
 155       track->SetLineColor(7);
 156     } else {
 157       //PH      track->SetName(Form("MUONTrack %2d     ", rt.label));
 158       sprintf(form,"MUONTrack %2d     ", rt.label);
 159       track->SetName(form);
 160       track->SetLineColor(6);
 161     }
 162
 163     AliMUONTrackParam *trackParam = mt->GetTrackParamAtVertex();
 164     xRec0  = trackParam->GetNonBendingCoor();
 165     yRec0  = trackParam->GetBendingCoor();
 166     zRec0  = trackParam->GetZ();
 167
 168     track->SetPoint(count,xRec0,yRec0,zRec0);
 169     count++;
 170
 171     Float_t xr[20], yr[20], zr[20];
 172     for (Int_t i = 0; i < 10; i++) xr[i]=yr[i]=zr[i]=0.0;
 173
 174     Int_t nTrackHits = mt->GetNTrackHits();
 175     //printf("Nhits = %d \n",nTrackHits);
 176     TClonesArray* trackParamAtHit;
 177     for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
 178       trackParamAtHit = mt->GetTrackParamAtHit();
 179       trackParam = (AliMUONTrackParam*) trackParamAtHit->At(iHit);
 180       xRec  = trackParam->GetNonBendingCoor();
 181       yRec  = trackParam->GetBendingCoor();
 182       zRec  = trackParam->GetZ();
 183
 184       //printf("Hit %d x %f y %f z %f \n",iHit,xRec,yRec,zRec);
 185
 186       xr[iHit] = xRec;
 187       yr[iHit] = yRec;
 188       zr[iHit] = zRec;
 189
 190       track->SetPoint(count,xRec,yRec,zRec);
 191       count++;
 192
 193     }
 194
 195     Float_t xrc[20], yrc[20], zrc[20];
 196     Int_t nrc = 0;
 197     if (mt->GetMatchTrigger() && 1) {
 198
 199       for (Int_t i = 0; i < nTrackHits; i++) {
 200         if (TMath::Abs(zr[i]) > 1000.0) {
 201           //printf("Hit %d x %f y %f z %f \n",iHit,xr[i],yr[i],zr[i]);
 202           xrc[nrc] = xr[i];
 203           yrc[nrc] = yr[i];
 204           zrc[nrc] = zr[i];
 205           nrc++;
 206         }
 207       }
 208
 209       Double_t xv, yv;
 210       Float_t ax, bx, ay, by;
 211
 212       if (nrc < 2) continue;
 213
 214       // fit x-z
 215       smatrix.Zero();
 216       sums.Zero();
 217       for (Int_t i = 0; i < nrc; i++) {
 218         xv = (Double_t)zrc[i];
 219         yv = (Double_t)xrc[i];
 220         //printf("x-z: xv %f yv %f \n",xv,yv);
 221         smatrix(0,0) += 1.0;
 222         smatrix(1,1) += xv*xv;
 223         smatrix(0,1) += xv;
 224         smatrix(1,0) += xv;
 225         sums(0,0)    += yv;
 226         sums(1,0)    += xv*yv;
 227       }
 228       res = smatrix.Invert() * sums;
 229       ax = res(0,0);
 230       bx = res(1,0);
 231
 232       // fit y-z
 233       smatrix.Zero();
 234       sums.Zero();
 235       for (Int_t i = 0; i < nrc; i++) {
 236         xv = (Double_t)zrc[i];
 237         yv = (Double_t)yrc[i];
 238         //printf("y-z: xv %f yv %f \n",xv,yv);
 239         smatrix(0,0) += 1.0;
 240         smatrix(1,1) += xv*xv;
 241         smatrix(0,1) += xv;
 242         smatrix(1,0) += xv;
 243         sums(0,0)    += yv;
 244         sums(1,0)    += xv*yv;
 245       }
 246       res = smatrix.Invert() * sums;
 247       ay = res(0,0);
 248       by = res(1,0);
 249
 250       Float_t xtc, ytc, ztc;
 251       for (Int_t ii = 0; ii < 4; ii++) {
 252
 253         ztc = zg[ii];
 254         ytc = ay+by*zg[ii];
 255         xtc = ax+bx*zg[ii];
 256
 257         //printf("tc: x %f y %f z %f \n",xtc,ytc,ztc);
 258
 259         track->SetPoint(count,xtc,ytc,ztc);
 260         count++;
 261
 262       }
 263
 264     }  // end match trigger
 265
 266     gReve->AddRenderElement(lt, track);
 267
 268   }
 269
 270 }
 271