1 #if !defined(__CINT__) || defined(__MAKECINT__)
\r
2 #include <Riostream.h>
\r
4 #include <TClonesArray.h>
\r
7 #include <TStopwatch.h>
\r
8 #include "AliAlignObjParams.h"
\r
9 #include "AliTrackPointArray.h"
\r
10 #include "AliITSAlignMille.h"
\r
12 //********************************************************************
\r
13 // Example to run ITS alignment using Millepede
\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
19 // Origin: M. Lunardon
\r
20 //********************************************************************
\r
22 Bool_t SelectLayers(AliTrackPointArray *tpa, int *nreqpts_lay) {
\r
23 // selection on layers
\r
24 int npts=tpa->GetNPoints();
\r
26 for (int jj=0;jj<6;jj++) nlay[jj]=0;
\r
27 for (int ip=0; ip<npts; ip++) {
\r
29 float r=TMath::Sqrt(tpa->GetX()[ip]*tpa->GetX()[ip] + tpa->GetY()[ip]*tpa->GetY()[ip]);
\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
40 for (int jj=0; jj<6; jj++) {
\r
41 if (nlay[jj]<nreqpts_lay[jj]) isgood=0;
\r
45 //////////////////////////////////////
\r
47 int ITSAlignMille(int fromev=1, int toev=200, int maxentries=-1, char *outfile="ITSAlignMille.root") {
\r
49 // Read tracks from AliTrackPoints.N.root and fill Millepede.
\r
50 // Millepede results are written as AliAlignObjParams in outfile.
\r
53 int nreqpts_lay[6]={0,0,0,0,0,0};
\r
55 TFile *fout=new TFile(outfile,"recreate");
\r
56 if (!fout->IsOpen()) {
\r
57 printf("\nCannot open output file!\n");
\r
61 TChain *chain=new TChain("spTree");
\r
62 char dir[100]="AliTrackPoints";
\r
65 for (int iev=fromev; iev<=toev; iev++) {
\r
66 sprintf(st,"%s/AliTrackPoints.%d.root",dir,iev);
\r
70 int nentries=chain->GetEntries();
\r
71 printf("There are %d entries in chain\n",nentries);
\r
73 if (maxentries>0 && maxentries<nentries) nentries=maxentries;
\r
75 AliTrackPointArray *tpa = 0;
\r
76 chain->SetBranchAddress("SP", &tpa);
\r
78 ////////////////////////////////////////////
\r
80 AliITSAlignMille *alig = new AliITSAlignMille("AliITSAlignMille.conf");
\r
82 Int_t nmod=alig->GetNModules();
\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
88 for(Int_t k=0;k<nmod*6;k++) {
\r
94 Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
\r
95 alig->InitGlobalParameters(parameters);
\r
98 ////////////////////////////////////////////////////////////////////
\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
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
121 printf("\nMillepede fed with %d tracks\n\n",iTrackOk);
\r
123 alig->GlobalFit(parameters,errors,pulls);
\r
125 cout << "Done with GlobalFit " << endl;
\r
127 ////////////////////////////////////////////////////////////
\r
131 TClonesArray *array = new TClonesArray("AliAlignObjParams",4000);
\r
132 TClonesArray &alobj = *array;
\r
134 // first object: ITS
\r
135 new(alobj[0]) AliAlignObjParams("ITS", 0, 0., 0., 0., 0., 0., 0., kTRUE);
\r
138 const Char_t *symname;
\r
139 Double_t dx,dy,dz,dpsi,dtheta,dphi;
\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
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
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
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
167 fout->WriteObject(array,"ITSAlignObjs","kSingleKey");
\r