]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AlignITSmacro3.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New version compatible with the current HEAD
[u/mrichter/AliRoot.git] / ITS / AlignITSmacro3.C
CommitLineData
e8189707 1//////////////////////////////////////////////////////////////////////////
2// Alice ITS first detector alignment program. //
3// //
4// version: 0.0.0 Draft. //
5// Date: April 18 1999 //
6// By: Bjorn S. Nilsen //
7// //
8//////////////////////////////////////////////////////////////////////////
9
10#include <fstream.h>
11
12// Data structure to hold averaged clusts.
13struct ClustAl_sl{
14 Int_t lay,lad,det;
15 Float_t xg,yg,zg,xl,yl,zl;
16};
17struct ClustAl_tl{
18 Int_t track,nclust; // track number and number of data points.
19 ClustAl_sl *clust; // data points to fit.
20 Float_t a,b,c,d,a0,b0,c0,d0,qual; // fit parameters and fit quality.
21 Float_t px,py,pz,p,pt;
22 // x=a+bz and y=c+dz;
23 // x=a0+b0*y and z=c0+d0*y in coordinate system of clust[0].lay,lad,det
24};
25
26//
27void AlignITSmacro3(const char *Rfilename="galice_ITS_B0.root",
28 const char *sfile="Align_ITS_B0",
29 Float_t sigma1=0.0,Float_t sigma2=0.0,Float_t sigma3=0.0,
30 Int_t evNumber=0) {
31/////////////////////////////////////////////////////////////////////////
32// This macro is a small example of a ROOT macro illustrating how to
33// read the output of GALICE and fill some histograms.
34//
35// Root > .L AlignITSmacro2.C // this loads the macro in memory
36// Root > AlignITSmacro2(); // by default process first event
37// Root > AlignITSmacro2("galice.root"); // process file galice.root
38// Root > AlignITSmacro2("galice.root",3); // process file galice.root
39// the first 4 events
40// or
41// Root > .x AlignITSmacro2.C; //by default process first event
42// Root > .x AlignITSmacro2.C("galice.root"); // process file galice.root
43// Root > .x AlignITSmacro2.C("galice.root",3); // process file galice.root
44// the first 4 events
45/////////////////////////////////////////////////////////////////////////
46//
47 gROOT->Reset(); // Reset root to it's default state
48// Dynamically link some shared libs
49 if (gClassTable->GetID("AliRun") < 0) {
50 gROOT->LoadMacro("loadlibs.C");
51 loadlibs();
52 } // end if gClassTable...
53
54 // Connect the Root Galice file containing Geometry, Kine and Clusts
55 TFile *Rfile = (TFile*)gROOT->GetListOfFiles()->FindObject(Rfilename);
56 if(!Rfile) Rfile = new TFile(Rfilename);
57 printf("reading from file %s\n",Rfilename);
58
59 // Get AliRun object from file or create it if not on file
60 if(!gAlice) {
61 gAlice = (AliRun*)Rfile->Get("gAlice");
62 if( gAlice) printf("AliRun object found on file\n");
63 if(!gAlice) gAlice = new AliRun("gAlice","Alice test program");
64 } /* end if gAlice */
65
66 Float_t v0[3] = {0.0,0.0,0.0};
67 Float_t tran[3] = {0.0,0.0,0.0},rot[3] = {0.0,0.0,0.0};
68 Float_t trans[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
69 2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
70 2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // cm
71 Float_t rots[15] ={0.0E-0,1.0E-4,4.0E-4,7.0E-4,1.0E-3,
72 2.0E-3,4.0E-3,6.0E-3,8.0E-3,1.0E-2,
73 2.0E-2,3.0E-2,5.0E-2,7.5E-2,1.0E-1}; // rad
74 Int_t nparticles,Ntrkp,ntrk,Itimes,Isigmas;
75 AliITS *ITS = 0;
76 TTree *TH = 0;
77 AliITSgeom *gm,gm2;
78 char Hfilename[80];
79// char Gfilename[80];
80 TFile *Hfile = 0;
81 ClustAl_tl *trk = 0;
82 Float_t Rdta[6];
83 Float_t Fdta[24],Fdta0[24];
84 Int_t Ndta[12];
85 FILE *fp1,*fp2;
86
87 for(Int_t evnt=0;evnt<=evNumber;evnt++){
88 // define some variables for later use.
89 nparticles = gAlice->GetEvent(evnt);
90 printf("nparticles %d\n",nparticles);
91 if (nparticles <= 0) continue; /* get next event */
92
93 // Get pointers to Alice detectors and Clusts containers
94 ITS = (AliITS*)gAlice->GetDetector("ITS");
95 if(!ITS) return; /* error no ITS data exit */
96 TH = gAlice->TreeH();
97 Ntrkp = TH->GetEntries();
98 gm = ITS->GetITSgeom();
99
100 // Array (stucture) of clusts for the first and second layer
101 // this should be replaced with either clusters or digits
102 // when they are proporly defined.
103 trk = new ClustAl_tl[Ntrkp];
104
105 printf("Ntrkp=%d\n",Ntrkp);
106
107 HitsToClustAl(trk,ntrk,Ntrkp,TH,ITS,1.0);
108 printf("Filled data structures ntrk=%d fitting lines next\n",ntrk);
109//
110 for(Itimes=1;Itimes<2;Itimes++){
111 if(Itimes==0){
112 fp1 = fopen("RvariationsR.csv","w");
113 fp2 = fopen("RvariationsN.csv","w");
114 fprintf(fp1,"Rr,Rs,Rrx1,Rerx1,Rsrx1,Rserx1,Rrz1,Rerz1,Rsrz1,Rserz1,"
115 "Rrx2,Rerx2,Rsrx2,Rserx2,Rrz2,Rerz2,Rsrz2,Rserz2,"
116 "Rrx3,Rerx3,Rsrx3,Rserx3,Rrz3,Rerz3,Rsrz3,Rserz3,"
117 "Rrx4,Rerx4,Rsrx4,Rserx4,Rrz4,Rerz4,Rsrz4,Rserz4,"
118 "Rrx5,Rerx5,Rsrx5,Rserx5,Rrz5,Rerz5,Rsrz5,Rserz5,"
119 "Rrx6,Rerx6,Rsrx6,Rserx6,Rrz6,Rerz6,Rsrz6,Rserz6");
120 fprintf(fp2,"Rr,Rs,RN1,RN2,RN3,RN4,RN5,RN6,RN7,RN8,RN9,RN10,"
121 "RN50,Rcut");
122 }else if(Itimes==1){
123 fp1 = fopen("RPhivariationsR.csv","w");
124 fp2 = fopen("RPhivariationsN.csv","w");
125 fprintf(fp1,"Fr,Fs,Frx1,Ferx1,Fsrx1,Fserx1,Frz1,Ferz1,Fsrz1,Fserz1,"
126 "Frx2,Ferx2,Fsrx2,Fserx2,Frz2,Ferz2,Fsrz2,Fserz2,"
127 "Frx3,Ferx3,Fsrx3,Fserx3,Frz3,Ferz3,Fsrz3,Fserz3,"
128 "Frx4,Ferx4,Fsrx4,Fserx4,Frz4,Ferz4,Fsrz4,Fserz4,"
129 "Frx5,Ferx5,Fsrx5,Fserx5,Frz5,Ferz5,Fsrz5,Fserz5,"
130 "Frx6,Ferx6,Fsrx6,Fserx6,Frz6,Ferz6,Fsrz6,Fserz6");
131 fprintf(fp2,"Fr,Fs,FN1,FN2,FN3,FN4,FN5,FN6,FN7,FN8,FN9,FN10,"
132 "FN50,Fcut");
133 }else if(Itimes==2){
134 fp1 = fopen("ZvariationsR.csv","w");
135 fp2 = fopen("ZvariationsN.csv","w");
136 fprintf(fp1,"Zr,Zs,Zrx1,Zerx1,Zsrx1,Zserx1,Zrz1,Zerz1,Zsrz1,Zserz1,"
137 "Zrx2,Zerx2,Zsrx2,Zserx2,Zrz2,Zerz2,Zsrz2,Zserz2,"
138 "Zrx3,Zerx3,Zsrx3,Zserx3,Zrz3,Zerz3,Zsrz3,Zserz3,"
139 "Zrx4,Zerx4,Zsrx4,Zserx4,Zrz4,Zerz4,Zsrz4,Zserz4,"
140 "Zrx5,Zerx5,Zsrx5,Zserx5,Zrz5,Zerz5,Zsrz5,Zserz5,"
141 "Zrx6,Zerx6,Zsrx6,Zserx6,Zrz6,Zerz6,Zsrz6,Zserz6");
142 fprintf(fp2,"Zr,Zs,ZN1,ZN2,ZN3,ZN4,ZN5,ZN6,ZN7,ZN8,ZN9,ZN10,"
143 "ZN50,Zcut");
144 }else if(Itimes==3){
145 fp1 = fopen("AvariationsR.csv","w");
146 fp2 = fopen("AvariationsN.csv","w");
147 fprintf(fp1,"Ar,As,Arx1,Aerx1,Asrx1,Aserx1,Arz1,Aerz1,Asrz1,Aserz1,"
148 "Arx2,Aerx2,Asrx2,Aserx2,Arz2,Aerz2,Asrz2,Aserz2,"
149 "Arx3,Aerx3,Asrx3,Aserx3,Arz3,Aerz3,Asrz3,Aserz3,"
150 "Arx4,Aerx4,Asrx4,Aserx4,Arz4,Aerz4,Asrz4,Aserz4,"
151 "Arx5,Aerx5,Asrx5,Aserx5,Arz5,Aerz5,Asrz5,Aserz5,"
152 "Arx6,Aerx6,Asrx6,Aserx6,Arz6,Aerz6,Asrz6,Aserz6");
153 fprintf(fp2,"Ar,As,AN1,AN2,AN3,AN4,AN5,AN6,AN7,AN8,AN9,AN10,"
154 "AN50,Acut");
155 }else if(Itimes==4){
156 fp1 = fopen("BvariationsR.csv","w");
157 fp2 = fopen("BvariationsN.csv","w");
158 fprintf(fp1,"Br,Bs,Brx1,Berx1,Bsrx1,Bserx1,Brz1,Berz1,Bsrz1,Bserz1,"
159 "Brx2,Berx2,Bsrx2,Bserx2,Brz2,Berz2,Bsrz2,Bserz2,"
160 "Brx3,Berx3,Bsrx3,Bserx3,Brz3,Berz3,Bsrz3,Bserz3,"
161 "Brx4,Berx4,Bsrx4,Bserx4,Brz4,Berz4,Bsrz4,Bserz4,"
162 "Brx5,Berx5,Bsrx5,Bserx5,Brz5,Berz5,Bsrz5,Bserz5,"
163 "Brx6,Berx6,Bsrx6,Bserx6,Brz6,Berz6,Bsrz6,Bserz6");
164 fprintf(fp2,"Br,Bs,BN1,BN2,BN3,BN4,BN5,BN6,BN7,BN8,BN9,BN10,"
165 "BN50,Bcut");
166 }else if(Itimes==5){
167 fp1 = fopen("CvariationsR.csv","w");
168 fp2 = fopen("CvariationsN.csv","w");
169 fprintf(fp1,"Cr,Cs,Crx1,Cerx1,Csrx1,Cserx1,Crz1,Cerz1,Csrz1,Cserz1,"
170 "Crx2,Cerx2,Csrx2,Cserx2,Crz2,Cerz2,Csrz2,Cserz2,"
171 "Crx3,Cerx3,Csrx3,Cserx3,Crz3,Cerz3,Csrz3,Cserz3,"
172 "Crx4,Cerx4,Csrx4,Cserx4,Crz4,Cerz4,Csrz4,Cserz4,"
173 "Crx5,Cerx5,Csrx5,Cserx5,Crz5,Cerz5,Csrz5,Cserz5,"
174 "Crx6,Cerx6,Csrx6,Cserx6,Crz6,Cerz6,Csrz6,Cserz6");
175 fprintf(fp2,"Cr,Cs,CN1,CN2,CN3,CN4,CN5,CN6,CN7,CN8,CN9,CN10,"
176 "CN50,Ccut");
177 } // end if
178 fprintf(fp1,"\n");
179 fprintf(fp2,"\n");
180 for(Isigmas=0;Isigmas<15;Isigmas++){
181//
182// tran[0] = sigma1;
183// tran[1] = sigma2;
184// tran[2] = sigma3;
185 if(Itimes==0){ tran[0] = trans[Isigmas];
186 }else tran[0] = 0.0;
187 if(Itimes==1){ tran[1] = trans[Isigmas];
188 }else tran[1] = 0.0;
189 if(Itimes==2){ tran[2] = trans[Isigmas];
190 }else tran[2] = 0.0;
191 if(Itimes==3){ rot[0] = rots[Isigmas];
192 }else rot[0] = 0.0;
193 if(Itimes==4){ rot[1] = rots[Isigmas];
194 }else rot[1] = 0.0;
195 if(Itimes==5){ rot[2] = rots[Isigmas];
196 }else rot[2] = 0.0;
197 printf("tran= %e %e %e (cm), rot=%e %e %e (rad)\n",
198 tran[0],tran[1],tran[2],rot[0],rot[1],rot[2]);
199//
200 gm2 = *gm;
201 gm2.RandomCylindericalChange(tran,rot);
202
203 FillGlobalPositions(trk,ntrk,&gm2);
204 // setup to save all created histograms.
205 sprintf(Hfilename,"%s_%04.0fr%04.0fp%04.0fz%04.0fx%04.0fy%04.0fz.root",
206 sfile,
207 10000.*tran[0],10000.*tran[1],10000.*tran[2],
208 10000.* rot[0],10000.* rot[1],10000.* rot[2]);
209 Hfile = (TFile*)gROOT->GetListOfFiles()->FindObject(Hfilename);
210 if(Hfile) Hfile->Close();
211 Hfile = new TFile(Hfilename,"RECREATE","Histograms from AlignITS.C");
212 printf("Histograms saved to file %s\n",Hfilename);
213 //
214 PlotGeomChanges(gm,&gm2,Hfile,Rdta);
215 //
216 // fit all tracks and do a track quality hist.
217 FitAllTracks(trk,ntrk,v0,&gm2,sfile,Hfile,Fdta,Ndta);
218//
219 if(Isigmas==0){
220 for(Int_t fp=0;fp<24;fp++) Fdta0[fp] = Fdta[fp];
221 } // end if Itimes==0&&Isigmas==0
222 if(Itimes==0) fprintf(fp1,"%e,%e,",tran[0],Rdta[0]);
223 else if(Itimes==1) fprintf(fp1,"%e,%e,",tran[1],Rdta[1]);
224 else if(Itimes==2) fprintf(fp1,"%e,%e,",tran[2],Rdta[2]);
225 else if(Itimes==3) fprintf(fp1,"%e,%e,",ran[0],Rdta[3]);
226 else if(Itimes==4) fprintf(fp1,"%e,%e,",ran[1],Rdta[4]);
227 else if(Itimes==5) fprintf(fp1,"%e,%e,",ran[2],Rdta[5]);
228 if(Itimes==0) fprintf(fp2,"%e,%e,",tran[0],Rdta[0]);
229 else if(Itimes==1) fprintf(fp2,"%e,%e,",tran[1],Rdta[1]);
230 else if(Itimes==2) fprintf(fp2,"%e,%e,",tran[2],Rdta[2]);
231 else if(Itimes==3) fprintf(fp2,"%e,%e,",ran[0],Rdta[3]);
232 else if(Itimes==4) fprintf(fp2,"%e,%e,",ran[1],Rdta[4]);
233 else if(Itimes==5) fprintf(fp2,"%e,%e,",ran[2],Rdta[5]);
234 fprintf(fp1,"%e,%e,%e,%e,",Fdta[0],Fdta[1],Fdta[0]/Fdta0[0],Fdta[0]/Fdta0[0]);
235 fprintf(fp1,"%e,%e,%e,%e,",Fdta[2],Fdta[3],Fdta[2]/Fdta0[2],Fdta[3]/Fdta0[2]);
236 fprintf(fp1,"%e,%e,%e,%e,",Fdta[4],Fdta[5],Fdta[4]/Fdta0[4],Fdta[5]/Fdta0[4]);
237 fprintf(fp1,"%e,%e,%e,%e,",Fdta[6],Fdta[7],Fdta[6]/Fdta0[6],Fdta[7]/Fdta0[6]);
238 fprintf(fp1,"%e,%e,%e,%e,",Fdta[8],Fdta[9],Fdta[8]/Fdta0[8],Fdta[9]/Fdta0[8]);
239 fprintf(fp1,"%e,%e,%e,%e,",Fdta[10],Fdta[11],Fdta[10]/Fdta0[10],Fdta[11]/Fdta0[10]);
240 fprintf(fp1,"%e,%e,%e,%e,",Fdta[12],Fdta[12],Fdta[12]/Fdta0[12],Fdta[13]/Fdta0[12]);
241 fprintf(fp1,"%e,%e,%e,%e,",Fdta[14],Fdta[13],Fdta[14]/Fdta0[14],Fdta[15]/Fdta0[14]);
242 fprintf(fp1,"%e,%e,%e,%e,",Fdta[16],Fdta[15],Fdta[16]/Fdta0[16],Fdta[17]/Fdta0[16]);
243 fprintf(fp1,"%e,%e,%e,%e,",Fdta[18],Fdta[17],Fdta[18]/Fdta0[18],Fdta[19]/Fdta0[18]);
244 fprintf(fp1,"%e,%e,%e,%e,",Fdta[20],Fdta[19],Fdta[20]/Fdta0[20],Fdta[21]/Fdta0[20]);
245 fprintf(fp1,"%e,%e,%e,%e\n",Fdta[22],Fdta[21],Fdta[22]/Fdta0[22],Fdta[23]/Fdta0[22]);
246 fprintf(fp2,"%d,%d,%d,%d,",Ndta[0],Ndta[1],Ndta[2],Ndta[3]);
247 fprintf(fp2,"%d,%d,%d,%d,",Ndta[4],Ndta[5],Ndta[6],Ndta[7]);
248 fprintf(fp2,"%d,%d,%d,%d\n",Ndta[8],Ndta[9],Ndta[10],Ndta[11]);
249// FitVertexAll(trk,ntrk,sfile,Hfile);
250 // Find all 2 track vertecies and hist. them.
251 Hfile->Close();
252 //
253 } // end for Isigmas
254 fclose(fp1);
255 fclose(fp2);
256 } // end for Itimes
257//
258 printf("Event %d done\n",evnt);
259//
260 deleteClustAl(trk,ntrk); // subrotine to delet memory allocated
261 // inside HitsToclustAl.
262 delete[] trk; // now delet memory allocated above.
263 } // end for evnt
264 Rfile->Close();
265 return;
266}