]>
Commit | Line | Data |
---|---|---|
d3603b4d | 1 | #if !defined(__CINT__) || defined(__MAKECINT__) |
2 | #include <Riostream.h> | |
3 | #include <TChain.h> | |
4 | #include <TClonesArray.h> | |
5 | #include <TFile.h> | |
6 | #include <TMath.h> | |
7 | #include <TStopwatch.h> | |
8 | #include "AliAlignObjParams.h" | |
9 | #include "AliTrackPointArray.h" | |
10 | #include "AliITSAlignMille.h" | |
11 | #endif | |
12 | ||
13 | //******************************************************************** | |
14 | // Example to run ITS alignment using Millepede | |
15 | // | |
16 | // Read tracks from AliTrackPoints.N.root and fill Millepede. | |
17 | // Millepede configuration is taken from AliITSAlignMille.conf | |
18 | // Results are written as AliAlignObjParams in outfile. | |
19 | // | |
20 | // Origin: M. Lunardon | |
21 | //******************************************************************** | |
22 | ||
23 | ||
24 | int ITSAlignMille(int fromev=1, int toev=1, int fromentry=0, int nentries=-1, char *outfile="ITSAlignMille.root", char *confile="AliITSAlignMille.conf",int nreqpts=3) { | |
25 | ||
26 | //AliLog::SetGlobalLogLevel(6); | |
27 | ||
28 | TFile *fout=new TFile(outfile,"recreate"); | |
29 | if (!fout->IsOpen()) { | |
30 | printf("\nCannot open output file!\n"); | |
31 | return -1; | |
32 | } | |
33 | ||
34 | TChain *chain=new TChain("spTree"); | |
35 | char dir[100]="AliTrackPoints"; | |
36 | char st[200]; | |
37 | ||
38 | for (int iev=fromev; iev<=toev; iev++) { | |
39 | sprintf(st,"%s/AliTrackPoints.%d.root",dir,iev); | |
40 | chain->Add(st); | |
41 | } | |
42 | ||
43 | int toentry=chain->GetEntries(); | |
44 | printf("There are %d entries in chain\n",toentry); | |
45 | ||
46 | if (nentries>0 && nentries<(toentry-fromentry)) toentry=fromentry+nentries; | |
47 | ||
48 | AliTrackPointArray *tpa = 0; | |
49 | chain->SetBranchAddress("SP", &tpa); | |
50 | ||
51 | //////////////////////////////////////////// | |
52 | ||
53 | AliITSAlignMille *alig = new AliITSAlignMille(confile); | |
54 | if (!alig->IsConfigured()) return -3; | |
55 | ||
56 | Int_t nmod=alig->GetNModules(); | |
57 | alig->SetMinNPtsPerTrack(nreqpts); | |
58 | ||
59 | // require 4 points in SPD (one per layer, up and down) | |
60 | if (nreqpts>3) { | |
61 | alig->SetRequiredPoint("LAYER",1,1,1); | |
62 | alig->SetRequiredPoint("LAYER",1,-1,1); | |
63 | alig->SetRequiredPoint("LAYER",2,1,1); | |
64 | alig->SetRequiredPoint("LAYER",2,-1,1); | |
65 | } | |
66 | ||
67 | // correction for SSD bug (1) : inverisione sensor 18 e 19 | |
68 | // alig->SetBug(1); | |
69 | ||
70 | Double_t *parameters = new Double_t[nmod*6]; | |
71 | Double_t *errors = new Double_t[nmod*6]; | |
72 | Double_t *pulls = new Double_t[nmod*6]; | |
73 | ||
74 | for(Int_t k=0;k<nmod*6;k++) { | |
75 | parameters[k]=0.; | |
76 | errors[k]=0.; | |
77 | pulls[k]=0.; | |
78 | } | |
79 | ||
80 | // need array with size fNLocal*2 | |
81 | Double_t trackParams[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; | |
82 | alig->InitGlobalParameters(parameters); | |
83 | alig->Print(); | |
84 | ||
85 | Double_t sigmatra=alig->GetParSigTranslations(); | |
86 | Double_t sigmarot=alig->GetParSigRotations(); | |
87 | //////////////////////////////////////////////////////////////////// | |
88 | ||
89 | TStopwatch stimer; | |
90 | stimer.Start(); | |
91 | ||
92 | ///////////////////// | |
93 | ||
94 | int iTrackOk=0; // number of good passed track | |
95 | // loop on spTree entries | |
96 | // one entry = one track | |
97 | for (int i=fromentry; i<toentry; i++) { | |
98 | ||
99 | chain->GetEvent(i); | |
100 | if (!alig->ProcessTrack(tpa)) { | |
101 | if (!(iTrackOk%500)) | |
102 | printf("Calling LocalFit n. %d\n",iTrackOk); | |
103 | alig->LocalFit(iTrackOk++,trackParams,0); | |
104 | } | |
105 | } | |
106 | printf("\nMillepede fed with %d tracks\n",iTrackOk); | |
107 | printf("Total number of rejected points because of bad EqLoc : %d\n\n",alig->GetTotBadLocEqPoints()); | |
108 | ||
109 | stimer.Print(); | |
110 | stimer.Continue(); | |
111 | ||
112 | // make global fit | |
113 | alig->GlobalFit(parameters,errors,pulls); | |
114 | ||
115 | cout << "Done with GlobalFit " << endl; | |
116 | ||
117 | //////////////////////////////////////////////////////////// | |
118 | // output | |
119 | ||
120 | printf("\nProcessed points statistics:\n"); | |
121 | Int_t maxstat=0; | |
122 | ||
123 | printf("\nOutput values and point statistics:\n\n"); | |
124 | for (int i=0; i<nmod; i++) { | |
125 | Int_t statm=alig->GetProcessedPoints()[i]; | |
126 | if (statm>maxstat) maxstat=statm; | |
127 | printf("index: %-4d stat: %-7d pars: %f %f %f %f %f %f\n",alig->GetModuleIndexArray()[i], statm, | |
128 | parameters[i*6+0], | |
129 | parameters[i*6+1], | |
130 | parameters[i*6+2], | |
131 | parameters[i*6+3], | |
132 | parameters[i*6+4], | |
133 | parameters[i*6+5]); | |
134 | ||
135 | } | |
136 | FILE *pf=fopen("ITSAlignMille.out","w"); | |
137 | if (pf) { | |
138 | fprintf(pf,"# param dx dy dz dpsi dtheta dphi\n"); | |
139 | for (int i=0; i<nmod; i++) { | |
140 | fprintf(pf," %-5d %-10f %-10f %-10f %-10f %-10f %-10f\n",i, | |
141 | parameters[i*6+0], | |
142 | parameters[i*6+1], | |
143 | parameters[i*6+2], | |
144 | parameters[i*6+3], | |
145 | parameters[i*6+4], | |
146 | parameters[i*6+5]); | |
147 | ||
148 | } | |
149 | fclose(pf); | |
150 | } | |
151 | ||
152 | printf("Max statistics = %d\n",maxstat); | |
153 | if (maxstat<1) { | |
154 | printf("No points for alignment! quitting now...\n"); | |
155 | return -1; | |
156 | } | |
157 | ||
158 | TClonesArray *array = new TClonesArray("AliAlignObjParams",4000); | |
159 | TClonesArray &alobj = *array; | |
160 | ||
161 | // first object: ITS | |
162 | new(alobj[0]) AliAlignObjParams("ITS", 0, 0., 0., 0., 0., 0., 0., kTRUE); | |
163 | ||
164 | UShort_t volid; | |
165 | Char_t *symname; | |
166 | Double_t dx,dy,dz,dpsi,dtheta,dphi; | |
167 | Double_t corr[21]; | |
168 | Double_t deltalocal[6]; | |
169 | Double_t t[3],r[3]; | |
170 | ||
171 | AliAlignObjParams *tmpa=NULL; | |
172 | ||
173 | // quality factor: 0 = NOT ALIGNED or OFF | |
174 | // N = ALIGNED with statistics = N | |
175 | UInt_t QF; | |
176 | ||
177 | // all ITS sensitive modules | |
178 | for (int idx=0; idx<2198; idx++) { | |
179 | volid=AliITSAlignMilleModule::GetVolumeIDFromIndex(idx); | |
180 | symname = AliGeomManager::SymName(volid); | |
181 | // default null misalignment | |
182 | for (int jj=0;jj<21;jj++) corr[jj]=0.0; | |
183 | for (int jj=0;jj<3;jj++) {t[jj]=0.0;r[jj]=0.0;} | |
184 | QF=0; | |
185 | ||
186 | if (alig->IsContained(volid)>=0) { // defined modules or inside a supermodule | |
187 | alig->SetCurrentModule(idx); // set the supermodule that contains idx | |
188 | int iidx=alig->GetCurrentModuleInternalIndex(); | |
189 | ||
190 | // check if all 6 parameters have been evaluated by millepede | |
191 | Bool_t isoff=0; | |
192 | Double_t parsig=0; | |
193 | for (int kk=0; kk<6; kk++) { | |
194 | parsig=sigmatra; | |
195 | if (kk>2) parsig=sigmarot; | |
196 | if (pulls[iidx*6+kk]==0.0 && errors[iidx*6+kk]==parsig) isoff=1; | |
197 | } | |
198 | ||
199 | // check if module was fixed | |
200 | Bool_t isfixed=1; | |
201 | for (int kk=0; kk<6; kk++) { | |
202 | if (parameters[iidx*6+kk]!=0.0 || pulls[iidx*6+kk]!=0.0 || errors[iidx*6+kk]!=0.0) isfixed=0; | |
203 | } | |
204 | ||
205 | if (!isoff && !isfixed) { // OK, has been evaluated | |
206 | deltalocal[0] = parameters[iidx*6+0]; | |
207 | deltalocal[1] = parameters[iidx*6+1]; | |
208 | deltalocal[2] = parameters[iidx*6+2]; | |
209 | deltalocal[3] = parameters[iidx*6+3]; | |
210 | deltalocal[4] = parameters[iidx*6+4]; | |
211 | deltalocal[5] = parameters[iidx*6+5]; | |
212 | tmpa = alig->GetCurrentModule()->GetSensitiveVolumeMisalignment(volid,deltalocal); | |
213 | if (!tmpa) { | |
214 | printf("error transforming millepede parameters for module %d\n",idx); | |
215 | continue; | |
216 | } | |
217 | tmpa->GetPars(t,r); | |
218 | // at the moment sigma of supermodule is given to sensitive modules. to be fixed... | |
219 | corr[0] = errors[iidx*6+0]*errors[iidx*6+0]; | |
220 | corr[2] = errors[iidx*6+1]*errors[iidx*6+1]; | |
221 | corr[5] = errors[iidx*6+2]*errors[iidx*6+2]; | |
222 | corr[9] = errors[iidx*6+3]*errors[iidx*6+3]; | |
223 | corr[14]= errors[iidx*6+4]*errors[iidx*6+4]; | |
224 | corr[20]= errors[iidx*6+5]*errors[iidx*6+5]; | |
225 | Int_t statm=alig->GetProcessedPoints()[iidx]; | |
226 | QF=statm; | |
227 | } | |
228 | else { | |
229 | if (isoff) { | |
230 | printf("Module %d is OFF!\n",idx); | |
231 | QF=0; | |
232 | } | |
233 | if (isfixed) { | |
234 | printf("Module %d was FIXED!\n",idx); | |
235 | if (alig->GetPreAlignmentQualityFactor(idx)>0) | |
236 | QF=alig->GetPreAlignmentQualityFactor(idx); | |
237 | else | |
238 | QF=0; | |
239 | } | |
240 | } | |
241 | ||
242 | } | |
243 | ||
244 | new(alobj[idx+1]) AliAlignObjParams(symname,volid,t[0],t[1],t[2],r[0],r[1],r[2],kFALSE); | |
245 | AliAlignObjParams* alo = (AliAlignObjParams*) array->UncheckedAt(idx+1); | |
246 | alo->SetCorrMatrix(corr); | |
247 | alo->SetUniqueID(QF); | |
248 | } | |
249 | ||
250 | fout->WriteObject(array,"ITSAlignObjs","kSingleKey"); | |
251 | fout->Close(); | |
252 | ||
253 | stimer.Stop(); | |
254 | stimer.Print(); | |
255 | ||
256 | return 0; | |
257 | } |