Updated macro
[u/mrichter/AliRoot.git] / ITS / ITSAlignMille.C
CommitLineData
95958115 1#if !defined(__CINT__) || defined(__MAKECINT__)\r
2#include <Riostream.h>\r
3#include <TChain.h>\r
4#include <TClonesArray.h>\r
5#include <TFile.h>\r
6#include <TMath.h>\r
7#include <TStopwatch.h>\r
8#include "AliAlignObjParams.h"\r
9#include "AliTrackPointArray.h"\r
10#include "AliITSAlignMille.h"\r
11#endif\r
12//********************************************************************\r
13// Example to run ITS alignment using Millepede\r
14//\r
15// Read tracks from AliTrackPoints.N.root and fill Millepede.\r
16// Millepede configuration is taken from AliITSAlignMille.conf\r
17// Results are written as AliAlignObjParams in outfile.\r
18//\r
19// Origin: M. Lunardon\r
20//********************************************************************\r
21\r
22Bool_t SelectLayers(AliTrackPointArray *tpa, int *nreqpts_lay) {\r
23 // selection on layers\r
24 int npts=tpa->GetNPoints();\r
25 int nlay[6]; \r
26 for (int jj=0;jj<6;jj++) nlay[jj]=0;\r
27 for (int ip=0; ip<npts; ip++) {\r
28 int lay=-1;\r
29 float r=TMath::Sqrt(tpa->GetX()[ip]*tpa->GetX()[ip] + tpa->GetY()[ip]*tpa->GetY()[ip]);\r
30 if (r<5) lay=1;\r
31 else if (r>5 && r<10) lay=2;\r
32 else if (r>10 && r<18) lay=3;\r
33 else if (r>18 && r<30) lay=4;\r
34 else if (r>30 && r<40) lay=5;\r
35 else if (r>40) lay=6;\r
36 if (lay<0) continue;\r
37 nlay[lay-1]++;\r
38 }\r
39 Bool_t isgood=1;\r
40 for (int jj=0; jj<6; jj++) {\r
41 if (nlay[jj]<nreqpts_lay[jj]) isgood=0;\r
42 }\r
43 return isgood;\r
44}\r
45//////////////////////////////////////\r
46\r
47int ITSAlignMille(int fromev=1, int toev=200, int maxentries=-1, char *outfile="ITSAlignMille.root") {\r
48\r
49 // Read tracks from AliTrackPoints.N.root and fill Millepede.\r
50 // Millepede results are written as AliAlignObjParams in outfile.\r
51 \r
52 int nreqpts=6;\r
53 int nreqpts_lay[6]={0,0,0,0,0,0};\r
54\r
55 TFile *fout=new TFile(outfile,"recreate");\r
56 if (!fout->IsOpen()) {\r
57 printf("\nCannot open output file!\n");\r
58 return -1;\r
59 }\r
60\r
61 TChain *chain=new TChain("spTree"); \r
62 char dir[100]="AliTrackPoints";\r
63 char st[200];\r
64 \r
65 for (int iev=fromev; iev<=toev; iev++) {\r
66 sprintf(st,"%s/AliTrackPoints.%d.root",dir,iev);\r
67 chain->Add(st);\r
68 }\r
69 \r
70 int nentries=chain->GetEntries();\r
71 printf("There are %d entries in chain\n",nentries);\r
72 \r
73 if (maxentries>0 && maxentries<nentries) nentries=maxentries;\r
74 \r
75 AliTrackPointArray *tpa = 0;\r
76 chain->SetBranchAddress("SP", &tpa);\r
77\r
78 ////////////////////////////////////////////\r
79 \r
80 AliITSAlignMille *alig = new AliITSAlignMille("AliITSAlignMille.conf");\r
81\r
82 Int_t nmod=alig->GetNModules();\r
83\r
84 Double_t *parameters = new Double_t[nmod*6];\r
85 Double_t *errors = new Double_t[nmod*6];\r
86 Double_t *pulls = new Double_t[nmod*6];\r
87\r
88 for(Int_t k=0;k<nmod*6;k++) {\r
89 parameters[k]=0.;\r
90 errors[k]=0.;\r
91 pulls[k]=0.;\r
92 }\r
93\r
94 Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};\r
95 alig->InitGlobalParameters(parameters);\r
96 alig->Print();\r
97\r
98 ////////////////////////////////////////////////////////////////////\r
99\r
100 TStopwatch stimer;\r
101 stimer.Start();\r
102\r
103 int iTrackOk=0; // number of good passed track\r
104 // loop on spTree entries\r
105 // one entry = one track\r
106 for (int i=0; i<nentries; i++) {\r
107 chain->GetEvent(i);\r
108 //////////////////////////////////////////////////////\r
109 // track preselection \r
110 int npts=tpa->GetNPoints();\r
111 if (npts<nreqpts) continue;\r
112 if (!SelectLayers(tpa,nreqpts_lay)) continue;\r
113 //////////////////////////////////////////////////////\r
114\r
115 if (!alig->ProcessTrack(tpa)) {\r
116 if (!(iTrackOk%500)) \r
117 printf("Calling LocalFit n. %d\n",iTrackOk);\r
118 alig->LocalFit(iTrackOk++,trackParams,0);\r
119 }\r
120 }\r
121 printf("\nMillepede fed with %d tracks\n\n",iTrackOk);\r
122 \r
123 alig->GlobalFit(parameters,errors,pulls);\r
124 \r
125 cout << "Done with GlobalFit " << endl;\r
126 \r
127 ////////////////////////////////////////////////////////////\r
128 // output\r
129 \r
130\r
131 TClonesArray *array = new TClonesArray("AliAlignObjParams",4000);\r
132 TClonesArray &alobj = *array;\r
133\r
134 // first object: ITS\r
135 new(alobj[0]) AliAlignObjParams("ITS", 0, 0., 0., 0., 0., 0., 0., kTRUE);\r
136 \r
137 UShort_t volid;\r
138 const Char_t *symname;\r
139 Double_t dx,dy,dz,dpsi,dtheta,dphi;\r
140 Double_t corr[21];\r
141\r
142 // all ITS modules\r
143 for (int idx=0; idx<2198; idx++) {\r
144 volid=alig->GetModuleVolumeID(idx);\r
145 symname = AliGeomManager::SymName(volid);\r
146 for (int jj=0;jj<21;jj++) corr[jj]=0.0;\r
147\r
148 if (alig->CheckVolumeID(volid)) { // defined modules\r
149 alig->SetCurrentModule(idx);\r
150 int iidx=alig->GetCurrentModuleInternalIndex();\r
151 dx = parameters[iidx*6+0]; corr[0] = errors[iidx*6+0]*errors[iidx*6+0];\r
152 dy = parameters[iidx*6+1]; corr[2] = errors[iidx*6+1]*errors[iidx*6+1];\r
153 dz = parameters[iidx*6+2]; corr[5] = errors[iidx*6+2]*errors[iidx*6+2];\r
154 dpsi = parameters[iidx*6+3]; corr[9] = errors[iidx*6+3]*errors[iidx*6+3];\r
155 dtheta= parameters[iidx*6+4]; corr[14]= errors[iidx*6+4]*errors[iidx*6+4];\r
156 dphi = parameters[iidx*6+5]; corr[20]= errors[iidx*6+5]*errors[iidx*6+5];\r
157 }\r
158 else { // other modules\r
159 dx=0.0; dy=0.0; dz=0.0; dpsi=0.0; dtheta=0.0; dphi=0.0;\r
160 }\r
161 \r
162 new(alobj[idx+1]) AliAlignObjParams(symname, volid, dx, dy, dz, dpsi, dtheta, dphi, kFALSE);\r
163 AliAlignObjParams* alo = (AliAlignObjParams*) array->UncheckedAt(idx+1);\r
164 alo->SetCorrMatrix(corr);\r
165 }\r
166 \r
167 fout->WriteObject(array,"ITSAlignObjs","kSingleKey");\r
168 fout->Close();\r
169\r
170 stimer.Stop();\r
171 stimer.Print();\r
172 \r
173 return 0;\r
174}\r