Alignment.C 13 KB
Newer Older
1
2
3
4
5
6
7
#include "Alignment.h"

#include <TVirtualFitter.h>

Alignment::Alignment(bool debugging)
: Algorithm("Alignment"){
  debug = debugging;
CLICdp user's avatar
CLICdp user committed
8
9
  m_numberOfTracksForAlignment = 20000;
  nIterations = 3;
10
11
12
}

// Global container declarations
13
Tracks globalTracks;
14
15
16
17
18
string detectorToAlign;
Parameters* globalParameters;
int detNum;

void Alignment::initialise(Parameters* par){
19
  // Pick up the global parameters
20
21
22
  parameters = par;
}

23
// During run, just pick up tracks and save them till the end
Daniel Hynds's avatar
Daniel Hynds committed
24
StatusCode Alignment::run(Clipboard* clipboard){
25
26
 
  // Get the tracks
27
  Tracks* tracks = (Tracks*)clipboard->get("tracks");
28
  if(tracks == NULL){
Daniel Hynds's avatar
Daniel Hynds committed
29
    return Success;
30
31
32
33
  }
  
  // Make a local copy and store it 
  for(int iTrack=0; iTrack<tracks->size(); iTrack++){
34
35
    Track* track = (*tracks)[iTrack];
    Track* alignmentTrack = new Track(track);
36
37
38
    m_alignmenttracks.push_back(alignmentTrack);
  }
  
39
  // If we have enough tracks for the alignment, tell the event loop to finish
Daniel Hynds's avatar
Daniel Hynds committed
40
  if(m_alignmenttracks.size() >= m_numberOfTracksForAlignment) return Failure;
41
42
  
  // Otherwise keep going
Daniel Hynds's avatar
Daniel Hynds committed
43
  return Success;
44
45
46

}

47
48
49
50
51

// ========================================
//  Minimisation functions for Minuit
// ========================================

52
// METHOD 0
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
// This method will move the detector in question, refit all of the tracks, and try to minimise the
// track chi2. If there were no clusters from this detector on any tracks then it would do nothing!
void MinimiseTrackChi2(Int_t &npar, Double_t *grad, Double_t &result, Double_t *par, Int_t flag) {
  
  // Pick up new alignment conditions
  globalParameters->detector[detectorToAlign]->displacementX(par[detNum*6 + 0]);
  globalParameters->detector[detectorToAlign]->displacementY(par[detNum*6 + 1]);
  globalParameters->detector[detectorToAlign]->displacementZ(par[detNum*6 + 2]);
  globalParameters->detector[detectorToAlign]->rotationX(par[detNum*6 + 3]);
  globalParameters->detector[detectorToAlign]->rotationY(par[detNum*6 + 4]);
  globalParameters->detector[detectorToAlign]->rotationZ(par[detNum*6 + 5]);
  
  // Apply new alignment conditions
  globalParameters->detector[detectorToAlign]->update();
  
  // The chi2 value to be returned
  result = 0.;
  
  // Loop over all tracks
  for(int iTrack=0; iTrack<globalTracks.size(); iTrack++){
    // Get the track
74
    Track* track = globalTracks[iTrack];
75
    // Get all clusters on the track
76
    Clusters trackClusters = track->clusters();
77
78
    // Find the cluster that needs to have its position recalculated
    for(int iTrackCluster=0; iTrackCluster<trackClusters.size(); iTrackCluster++){
79
      Cluster* trackCluster = trackClusters[iTrackCluster];
80
      string detectorID = trackCluster->detectorID();
81
      if(detectorID != detectorToAlign) continue;
82
83
84
85
86
      // Recalculate the global position from the local
      PositionVector3D<Cartesian3D<double> > positionLocal(trackCluster->localX(),trackCluster->localY(),trackCluster->localZ());
      PositionVector3D<Cartesian3D<double> > positionGlobal = *(globalParameters->detector[detectorID]->m_localToGlobal) * positionLocal;
      trackCluster->setClusterCentre(positionGlobal.X(), positionGlobal.Y(),positionGlobal.Z());
    }
87

88
    // Refit the track
89
    track->fit();
90
    
91
92
93
94
95
96
    // Add the new chi2
    result += track->chi2();
  }
  
}

97
// METHOD 1
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
// This method will move the detector in question and try to minimise the (unbiased) residuals. It uses
// the associated cluster container on the track (no refitting of the track)
void MinimiseResiduals(Int_t &npar, Double_t *grad, Double_t &result, Double_t *par, Int_t flag) {
  
  // Pick up new alignment conditions
  globalParameters->detector[globalParameters->detectorToAlign]->displacementX(par[0]);
  globalParameters->detector[globalParameters->detectorToAlign]->displacementY(par[1]);
  globalParameters->detector[globalParameters->detectorToAlign]->displacementZ(par[2]);
  globalParameters->detector[globalParameters->detectorToAlign]->rotationX(par[3]);
  globalParameters->detector[globalParameters->detectorToAlign]->rotationY(par[4]);
  globalParameters->detector[globalParameters->detectorToAlign]->rotationZ(par[5]);
  
  // Apply new alignment conditions
  globalParameters->detector[globalParameters->detectorToAlign]->update();
  
  // The chi2 value to be returned
  result = 0.;
  
  // Loop over all tracks
  for(int iTrack=0; iTrack<globalTracks.size(); iTrack++){
    // Get the track
119
    Track* track = globalTracks[iTrack];
120
    // Get all clusters on the track
121
    Clusters associatedClusters = track->associatedClusters();
122
123
124

    // Find the cluster that needs to have its position recalculated
    for(int iAssociatedCluster=0; iAssociatedCluster<associatedClusters.size(); iAssociatedCluster++){
125
      Cluster* associatedCluster = associatedClusters[iAssociatedCluster];
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
      string detectorID = associatedCluster->detectorID();
      if(detectorID != globalParameters->detectorToAlign) continue;
      // Recalculate the global position from the local
      PositionVector3D<Cartesian3D<double> > positionLocal(associatedCluster->localX(),associatedCluster->localY(),associatedCluster->localZ());
      PositionVector3D<Cartesian3D<double> > positionGlobal = *(globalParameters->detector[globalParameters->detectorToAlign]->m_localToGlobal) * positionLocal;
      // Get the track intercept with the detector
      ROOT::Math::XYZPoint intercept = track->intercept(positionGlobal.Z());
      // Calculate the residuals
      double residualX = intercept.X() - positionGlobal.X();
      double residualY = intercept.Y() - positionGlobal.Y();
      double error = associatedCluster->error();
      // Add the new residual2
      result += ( (residualX*residualX + residualY*residualY) / (error*error)) ;
    }
  }
}

Daniel Hynds's avatar
Daniel Hynds committed
143
144
145
146
147

// ==================================================================
//  The finalise function - effectively the brains of the alignment!
// ==================================================================

148
149
void Alignment::finalise(){
  
150
  // If not enough tracks were produced, do nothing
151
  //if(m_alignmenttracks.size() < m_numberOfTracksForAlignment) return;
152

153
154
155
156
  // Make the fitting object
  TVirtualFitter* residualFitter = TVirtualFitter::Fitter(0,50);
  
  // Tell it what to minimise
157
158
  if(parameters->alignmentMethod == 0)residualFitter->SetFCN(MinimiseTrackChi2);
  if(parameters->alignmentMethod == 1)residualFitter->SetFCN(MinimiseResiduals);
159
160
161
162
163
164
165
166
167
168
169
170
171
172
  
  // Set the global parameters
  globalTracks = m_alignmenttracks;
  globalParameters = parameters;

  // Set the printout arguments of the fitter
  Double_t arglist[10];
  arglist[0] = 3;
  residualFitter->ExecuteCommand("SET PRINT",arglist,1);
  
  // Set some fitter parameters
  arglist[0] = 1000; // number of function calls
  arglist[1] = 0.001; // tolerance
  
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
  // This has been inserted in a temporary way. If the alignment method is 1 then it will align the single detector and then
  // return. This should be made into separate functions.
  if(parameters->alignmentMethod == 1){
    
    // Get the name of the detector to align
    string detectorID = parameters->detectorToAlign;
    
    // Add the parameters to the fitter (z displacement not allowed to move!)
    residualFitter->SetParameter(0,(detectorID+"_displacementX").c_str(),parameters->detector[detectorID]->displacementX(),0.01,-50,50);
    residualFitter->SetParameter(1,(detectorID+"_displacementY").c_str(),parameters->detector[detectorID]->displacementY(),0.01,-50,50);
    residualFitter->SetParameter(2,(detectorID+"_displacementZ").c_str(),parameters->detector[detectorID]->displacementZ(),0,-10,500);
    residualFitter->SetParameter(3,(detectorID+"_rotationX").c_str(),parameters->detector[detectorID]->rotationX(),0.001,-6.30,6.30);
    residualFitter->SetParameter(4,(detectorID+"_rotationY").c_str(),parameters->detector[detectorID]->rotationY(),0.001,-6.30,6.30);
    residualFitter->SetParameter(5,(detectorID+"_rotationZ").c_str(),parameters->detector[detectorID]->rotationZ(),0.001,-6.30,6.30);
    
Daniel Hynds's avatar
Daniel Hynds committed
188
189
190
191
    for(int iteration=0;iteration<nIterations;iteration++){
    
      // Fit this plane (minimising global track chi2)
    	residualFitter->ExecuteCommand("MIGRAD",arglist,2);
192
   
Daniel Hynds's avatar
Daniel Hynds committed
193
194
195
196
197
198
199
      // Set the alignment parameters of this plane to be the optimised values from the alignment
      parameters->detector[detectorID]->displacementX(residualFitter->GetParameter(0));
      parameters->detector[detectorID]->displacementY(residualFitter->GetParameter(1));
      parameters->detector[detectorID]->displacementZ(residualFitter->GetParameter(2));
      parameters->detector[detectorID]->rotationX(residualFitter->GetParameter(3));
      parameters->detector[detectorID]->rotationY(residualFitter->GetParameter(4));
      parameters->detector[detectorID]->rotationZ(residualFitter->GetParameter(5));
200

Daniel Hynds's avatar
Daniel Hynds committed
201
202
    }
    
203
204
205
206
207
208
    // Write the output alignment file
    parameters->writeConditions();

    return;
  }
  
209
210
211
212
213
  // Loop over all planes. For each plane, set the plane alignment parameters which will be varied, and
  // then minimise the track chi2 (sum of biased residuals). This means that tracks are refitted with
  // each minimisation step.
  
  int det = 0;
Daniel Hynds's avatar
Daniel Hynds committed
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
  for(int iteration=0;iteration<nIterations;iteration++){

    det = 0;
    for(int ndet = 0; ndet<parameters->nDetectors; ndet++){
      string detectorID = parameters->detectors[ndet];
      // Do not align the reference plane
      if(detectorID == parameters->reference) continue;
      if(parameters->excludedFromTracking.count(detectorID) != 0) continue;
      // Say that this is the detector we align
      detectorToAlign = detectorID;
      detNum=det;
      // Add the parameters to the fitter (z displacement not allowed to move!)
      residualFitter->SetParameter(det*6+0,(detectorID+"_displacementX").c_str(),parameters->detector[detectorID]->displacementX(),0.01,-50,50);
      residualFitter->SetParameter(det*6+1,(detectorID+"_displacementY").c_str(),parameters->detector[detectorID]->displacementY(),0.01,-50,50);
      residualFitter->SetParameter(det*6+2,(detectorID+"_displacementZ").c_str(),parameters->detector[detectorID]->displacementZ(),0,-10,500);
      residualFitter->SetParameter(det*6+3,(detectorID+"_rotationX").c_str(),parameters->detector[detectorID]->rotationX(),0.001,-6.30,6.30);
      residualFitter->SetParameter(det*6+4,(detectorID+"_rotationY").c_str(),parameters->detector[detectorID]->rotationY(),0.001,-6.30,6.30);
      residualFitter->SetParameter(det*6+5,(detectorID+"_rotationZ").c_str(),parameters->detector[detectorID]->rotationZ(),0.001,-6.30,6.30);
      
      // Fit this plane (minimising global track chi2)
      residualFitter->ExecuteCommand("MIGRAD",arglist,2);
      
      // Now that this device is fitted, set parameter errors to 0 so that they are not fitted again
      residualFitter->SetParameter(det*6+0,(detectorID+"_displacementX").c_str(),residualFitter->GetParameter(det*6+0),0,-50,50);
      residualFitter->SetParameter(det*6+1,(detectorID+"_displacementY").c_str(),residualFitter->GetParameter(det*6+1),0,-50,50);
      residualFitter->SetParameter(det*6+2,(detectorID+"_displacementZ").c_str(),residualFitter->GetParameter(det*6+2),0,-10,500);
      residualFitter->SetParameter(det*6+3,(detectorID+"_rotationX").c_str(),residualFitter->GetParameter(det*6+3),0,-6.30,6.30);
      residualFitter->SetParameter(det*6+4,(detectorID+"_rotationY").c_str(),residualFitter->GetParameter(det*6+4),0,-6.30,6.30);
      residualFitter->SetParameter(det*6+5,(detectorID+"_rotationZ").c_str(),residualFitter->GetParameter(det*6+5),0,-6.30,6.30);
      
      // Set the alignment parameters of this plane to be the optimised values from the alignment
      parameters->detector[detectorID]->displacementX(residualFitter->GetParameter(det*6+0));
      parameters->detector[detectorID]->displacementY(residualFitter->GetParameter(det*6+1));
      parameters->detector[detectorID]->displacementZ(residualFitter->GetParameter(det*6+2));
      parameters->detector[detectorID]->rotationX(residualFitter->GetParameter(det*6+3));
      parameters->detector[detectorID]->rotationY(residualFitter->GetParameter(det*6+4));
      parameters->detector[detectorID]->rotationZ(residualFitter->GetParameter(det*6+5));
      parameters->detector[detectorID]->update();
      det++;
    }
254
  
Daniel Hynds's avatar
Daniel Hynds committed
255
  }
256
  det = 0;
257

258
259
260
261
262
  // Now list the new alignment parameters
  for(int ndet = 0; ndet<parameters->nDetectors; ndet++){
    string detectorID = parameters->detectors[ndet];
    // Do not align the reference plane
    if(detectorID == parameters->reference) continue;
263
264
    if(parameters->excludedFromTracking.count(detectorID) != 0) continue;

265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
    // Get the alignment parameters
    double displacementX = residualFitter->GetParameter(det*6+0);
    double displacementY = residualFitter->GetParameter(det*6+1);
    double displacementZ = residualFitter->GetParameter(det*6+2);
    double rotationX = residualFitter->GetParameter(det*6+3);
    double rotationY = residualFitter->GetParameter(det*6+4);
    double rotationZ = residualFitter->GetParameter(det*6+5);

    tcout<<" Detector "<<detectorID<<" new alignment parameters: T("<<displacementX<<","<<displacementY<<","<<displacementZ<<") R("<<rotationX<<","<<rotationY<<","<<rotationZ<<")"<<endl;
    
    det++;
  }
  
  // Write the output alignment file
  parameters->writeConditions();
}