EVE/alice-macros/muon_tracks.C

   1 void muon_tracks()
   2 {
   3
   4   AliRunLoader* rl =  Alieve::Event::AssertRunLoader();
   5
   6   rl->LoadTracks("MUON");
   7
   8   TTree* tt = rl->GetTreeT("MUON", false);
   9
  10   Alieve::MUONDigitsInfo* di = new Alieve::MUONDigitsInfo();
  11   di->SetTTree(tt);
  12
  13   rl->LoadgAlice();
  14   gAlice = rl->GetAliRun();
  15   AliMUON *pMUON = (AliMUON*)gAlice->GetModule("MUON");
  16   const AliMUONGeometryTransformer* kGeomTransformer = pMUON->GetGeometryTransformer();
  17
  18   gStyle->SetPalette(1, 0);
  19
  20   char name[128];
  21   char title[128];
  22
  23   /* TRACKS */
  24
  25   /* Get the z positions of the trigger chambers */
  26
  27   Float_t xg, yg, zg[4];
  28
  29   // MT11 = ST6CH1, detElemId = 1100 ... 1117
  30   // MT12 = ST6CH2, detElemId = 1200 ... 1217
  31   // MT21 = ST7CH1, detElemId = 1300 ... 1317
  32   // MT22 = ST7CH2, detElemId = 1400 ... 1417
  33
  34   kGeomTransformer->Local2Global(1100, 0, 0, 0, xg, yg, zg[0]);
  35   kGeomTransformer->Local2Global(1200, 0, 0, 0, xg, yg, zg[1]);
  36   kGeomTransformer->Local2Global(1300, 0, 0, 0, xg, yg, zg[2]);
  37   kGeomTransformer->Local2Global(1400, 0, 0, 0, xg, yg, zg[3]);
  38
  39   /* from tracking chambers */
  40
  41   TMatrixD smatrix(2,2);
  42   TMatrixD sums(2,1);
  43   TMatrixD res(2,1);
  44
  45   TClonesArray *tracks = 0;
  46   tt->SetBranchAddress("MUONTrack",&tracks);
  47   tt->GetEntry(0);  // load event 0
  48
  49   Int_t ntracks = tracks->GetEntriesFast();
  50   printf("Found %d tracks. \n",ntracks);
  51
  52   Reve::TrackList* cont = new Reve::TrackList("M-Tracks");
  53   cont->SetMainColor(Color_t(6));
  54
  55   gReve->AddRenderElement(cont);
  56
  57   Float_t xRec, xRec0;
  58   Float_t yRec, yRec0;
  59   Float_t zRec, zRec0;
  60
  61   AliMUONTrack *mt;
  62   Reve::RecTrack  rt;
  63   for (Int_t n = 0; n < ntracks; n++) {
  64
  65     mt = (AliMUONTrack*) tracks->At(n);
  66
  67     printf("Match trigger %d \n",mt->GetMatchTrigger());
  68
  69     rt.label = n;
  70
  71     Int_t count = 0;
  72     Reve::Track* track = new Reve::Track(&rt, cont->GetRnrStyle());
  73
  74     if (mt->GetMatchTrigger()) {
  75       track->SetName(Form("MUONTrack %2d (MT)", rt.label));
  76       track->SetLineColor(7);
  77     } else {
  78       track->SetName(Form("MUONTrack %2d     ", rt.label));
  79       track->SetLineColor(6);
  80     }
  81
  82     AliMUONTrackParam *trackParam = mt->GetTrackParamAtVertex();
  83     xRec0  = trackParam->GetNonBendingCoor();
  84     yRec0  = trackParam->GetBendingCoor();
  85     zRec0  = trackParam->GetZ();
  86
  87     track->SetPoint(count,zRec0,yRec0,xRec0);
  88     count++;
  89
  90     Float_t xr[20], yr[20], zr[20];
  91     for (Int_t i = 0; i < 10; i++) xr[i]=yr[i]=zr[i]=0.0;
  92
  93     Int_t nTrackHits = mt->GetNTrackHits();
  94     printf("Nhits = %d \n",nTrackHits);
  95     for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
  96       trackParamAtHit = mt->GetTrackParamAtHit();
  97       trackParam = (AliMUONTrackParam*) trackParamAtHit->At(iHit);
  98       xRec  = trackParam->GetNonBendingCoor();
  99       yRec  = trackParam->GetBendingCoor();
 100       zRec  = trackParam->GetZ();
 101
 102       //printf("Hit %d x %f y %f z %f \n",iHit,xRec,yRec,zRec);
 103
 104       xr[iHit] = xRec;
 105       yr[iHit] = yRec;
 106       zr[iHit] = zRec;
 107
 108       track->SetPoint(count,zRec,yRec,xRec);
 109       count++;
 110
 111     }
 112
 113     Float_t xrc[20], yrc[20], zrc[20];
 114     Int_t nrc = 0;
 115     if (mt->GetMatchTrigger()) {
 116
 117       for (Int_t i = 0; i < nTrackHits; i++) {
 118         if (TMath::Abs(zr[i]) > 1000.0) {
 119           //printf("Hit %d x %f y %f z %f \n",iHit,xr[i],yr[i],zr[i]);
 120           xrc[nrc] = xr[i];
 121           yrc[nrc] = yr[i];
 122           zrc[nrc] = zr[i];
 123           nrc++;
 124         }
 125       }
 126
 127       Double_t xv, yv;
 128       Float_t ax, bx, ay, by;
 129
 130       // fit x-z
 131       smatrix.Zero();
 132       sums.Zero();
 133       for (Int_t i = 0; i < nrc; i++) {
 134         xv = (Double_t)zrc[i];
 135         yv = (Double_t)xrc[i];
 136         //printf("x-z: xv %f yv %f \n",xv,yv);
 137         smatrix(0,0) += 1.0;
 138         smatrix(1,1) += xv*xv;
 139         smatrix(0,1) += xv;
 140         smatrix(1,0) += xv;
 141         sums(0,0)    += yv;
 142         sums(1,0)    += xv*yv;
 143       }
 144       res = smatrix.Invert() * sums;
 145       ax = res(0,0);
 146       bx = res(1,0);
 147
 148       // fit y-z
 149       smatrix.Zero();
 150       sums.Zero();
 151       for (Int_t i = 0; i < nrc; i++) {
 152         xv = (Double_t)zrc[i];
 153         yv = (Double_t)yrc[i];
 154         //printf("y-z: xv %f yv %f \n",xv,yv);
 155         smatrix(0,0) += 1.0;
 156         smatrix(1,1) += xv*xv;
 157         smatrix(0,1) += xv;
 158         smatrix(1,0) += xv;
 159         sums(0,0)    += yv;
 160         sums(1,0)    += xv*yv;
 161       }
 162       res = smatrix.Invert() * sums;
 163       ay = res(0,0);
 164       by = res(1,0);
 165
 166       Float_t xtc, ytc, ztc;
 167       for (Int_t itc = 0; itc < 4; itc++) {
 168
 169         ztc = zg[itc];
 170         ytc = ay+by*zg[itc];
 171         xtc = ax+bx*zg[itc];
 172
 173         //printf("tc: x %f y %f z %f \n",xtc,ytc,ztc);
 174
 175         track->SetPoint(count,ztc,ytc,xtc);
 176         count++;
 177
 178       }
 179
 180     }
 181
 182     gReve->AddRenderElement(cont, track);
 183
 184   }
 185
 186   gReve->Redraw3D();
 187
 188 }