[u/mrichter/AliRoot.git] / ITS / ITSAlignMille.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TFile.h>
#include <TMath.h>
#include <TStopwatch.h>
#include "AliAlignObjParams.h"
#include "AliTrackPointArray.h"
#include "AliITSAlignMille.h"
#endif

//********************************************************************
//  Example to run ITS alignment using Millepede
//
//  Read tracks from AliTrackPoints.N.root and fill Millepede.
//  Millepede configuration is taken from AliITSAlignMille.conf
//  Results are written as AliAlignObjParams in outfile.
//
//      Origin: M. Lunardon
//********************************************************************


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) {

  //AliLog::SetGlobalLogLevel(6);
  
  TFile *fout=new TFile(outfile,"recreate");
  if (!fout->IsOpen()) {
    printf("\nCannot open output file!\n");
    return -1;
  }

  TChain *chain=new TChain("spTree");  
  char dir[100]="AliTrackPoints";
  char st[200];
  
  for (int iev=fromev; iev<=toev; iev++) {
    sprintf(st,"%s/AliTrackPoints.%d.root",dir,iev);
    chain->Add(st);
  }
  
  int toentry=chain->GetEntries();
  printf("There are %d entries in chain\n",toentry);
  
  if (nentries>0 && nentries<(toentry-fromentry)) toentry=fromentry+nentries;
  
  AliTrackPointArray *tpa = 0;
  chain->SetBranchAddress("SP", &tpa);

  ////////////////////////////////////////////
  
  AliITSAlignMille *alig = new AliITSAlignMille(confile);
  if (!alig->IsConfigured()) return -3;

  Int_t nmod=alig->GetNModules();
  alig->SetMinNPtsPerTrack(nreqpts);

  // require 4 points in SPD (one per layer, up and down)
  if (nreqpts>3) {
    alig->SetRequiredPoint("LAYER",1,1,1);
    alig->SetRequiredPoint("LAYER",1,-1,1);
    alig->SetRequiredPoint("LAYER",2,1,1);
    alig->SetRequiredPoint("LAYER",2,-1,1);
  }

  // correction for SSD bug (1) : inverisione sensor 18 e 19
  // alig->SetBug(1);

  Double_t *parameters = new Double_t[nmod*6];
  Double_t *errors = new Double_t[nmod*6];
  Double_t *pulls = new Double_t[nmod*6];

  for(Int_t k=0;k<nmod*6;k++) {
    parameters[k]=0.;
    errors[k]=0.;
    pulls[k]=0.;
  }

  // need array with size fNLocal*2
  Double_t trackParams[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  alig->InitGlobalParameters(parameters);
  alig->Print();

  Double_t sigmatra=alig->GetParSigTranslations();
  Double_t sigmarot=alig->GetParSigRotations();
  ////////////////////////////////////////////////////////////////////

  TStopwatch stimer;
  stimer.Start();

  /////////////////////

  int iTrackOk=0; // number of good passed track
  // loop on spTree entries
  // one entry = one track
  for (int i=fromentry; i<toentry; i++) {

    chain->GetEvent(i);
    if (!alig->ProcessTrack(tpa)) {
      if (!(iTrackOk%500)) 
	printf("Calling LocalFit n. %d\n",iTrackOk);
      alig->LocalFit(iTrackOk++,trackParams,0);
    }
  }
  printf("\nMillepede fed with %d tracks\n",iTrackOk);
  printf("Total number of rejected points because of bad EqLoc : %d\n\n",alig->GetTotBadLocEqPoints());
  
  stimer.Print();
  stimer.Continue(); 

  // make global fit
  alig->GlobalFit(parameters,errors,pulls);
  
  cout << "Done with GlobalFit " << endl;
  
  ////////////////////////////////////////////////////////////
  // output

  printf("\nProcessed points statistics:\n");
  Int_t maxstat=0;

  printf("\nOutput values and point statistics:\n\n");
  for (int i=0; i<nmod; i++) {
    Int_t statm=alig->GetProcessedPoints()[i];
    if (statm>maxstat) maxstat=statm;
    printf("index: %-4d   stat: %-7d   pars: %f   %f   %f   %f   %f   %f\n",alig->GetModuleIndexArray()[i], statm,
	   parameters[i*6+0],
	   parameters[i*6+1],
	   parameters[i*6+2],
	   parameters[i*6+3],
	   parameters[i*6+4],
	   parameters[i*6+5]);
    
  }
  FILE *pf=fopen("ITSAlignMille.out","w");
  if (pf) {
    fprintf(pf,"# param     dx         dy         dz         dpsi       dtheta     dphi\n");
    for (int i=0; i<nmod; i++) {
      fprintf(pf,"   %-5d %-10f %-10f %-10f %-10f %-10f %-10f\n",i,
	      parameters[i*6+0],
	      parameters[i*6+1],
	      parameters[i*6+2],
	      parameters[i*6+3],
	      parameters[i*6+4],
	      parameters[i*6+5]);
      
    }
    fclose(pf);
  }

  printf("Max statistics = %d\n",maxstat);
  if (maxstat<1) {
    printf("No points for alignment! quitting now...\n");
    return -1;
  }  

  TClonesArray *array = new TClonesArray("AliAlignObjParams",4000);
  TClonesArray &alobj = *array;

  // first object: ITS
  new(alobj[0]) AliAlignObjParams("ITS", 0, 0., 0., 0., 0., 0., 0., kTRUE);
  
  UShort_t volid;
  Char_t *symname;
  Double_t dx,dy,dz,dpsi,dtheta,dphi;
  Double_t corr[21];
  Double_t deltalocal[6];
  Double_t t[3],r[3];

  AliAlignObjParams *tmpa=NULL;

  // quality factor: 0 = NOT ALIGNED or OFF
  //                 N = ALIGNED with statistics = N  
  UInt_t QF; 

  // all ITS sensitive modules
  for (int idx=0; idx<2198; idx++) {
    volid=AliITSAlignMilleModule::GetVolumeIDFromIndex(idx);
    symname = AliGeomManager::SymName(volid);
    // default null misalignment
    for (int jj=0;jj<21;jj++) corr[jj]=0.0;
    for (int jj=0;jj<3;jj++) {t[jj]=0.0;r[jj]=0.0;}
    QF=0;

    if (alig->IsContained(volid)>=0) { // defined modules or inside a supermodule
      alig->SetCurrentModule(idx); // set the supermodule that contains idx
      int iidx=alig->GetCurrentModuleInternalIndex();
      
      // check if all 6 parameters have been evaluated by millepede
      Bool_t isoff=0;
      Double_t parsig=0;
      for (int kk=0; kk<6; kk++) {
	parsig=sigmatra;
	if (kk>2) parsig=sigmarot;
	if (pulls[iidx*6+kk]==0.0 && errors[iidx*6+kk]==parsig) isoff=1;
      }

      // check if module was fixed
      Bool_t isfixed=1;
      for (int kk=0; kk<6; kk++) {
	if (parameters[iidx*6+kk]!=0.0 || pulls[iidx*6+kk]!=0.0 || errors[iidx*6+kk]!=0.0) isfixed=0;
      }

      if (!isoff && !isfixed) { // OK, has been evaluated
	deltalocal[0] = parameters[iidx*6+0];  
	deltalocal[1] = parameters[iidx*6+1]; 
	deltalocal[2] = parameters[iidx*6+2]; 
	deltalocal[3] = parameters[iidx*6+3]; 
	deltalocal[4] = parameters[iidx*6+4];
	deltalocal[5] = parameters[iidx*6+5]; 
	tmpa = alig->GetCurrentModule()->GetSensitiveVolumeMisalignment(volid,deltalocal);
	if (!tmpa) {
	  printf("error transforming millepede parameters for module %d\n",idx);
	  continue;
	}
	tmpa->GetPars(t,r);
	// at the moment sigma of supermodule is given to sensitive modules. to be fixed...
	corr[0] = errors[iidx*6+0]*errors[iidx*6+0];
	corr[2] = errors[iidx*6+1]*errors[iidx*6+1];
	corr[5] = errors[iidx*6+2]*errors[iidx*6+2];
	corr[9] = errors[iidx*6+3]*errors[iidx*6+3];
	corr[14]= errors[iidx*6+4]*errors[iidx*6+4];
	corr[20]= errors[iidx*6+5]*errors[iidx*6+5];
	Int_t statm=alig->GetProcessedPoints()[iidx];
	QF=statm;
      }
      else {
	if (isoff) {
	  printf("Module %d is OFF!\n",idx);
	  QF=0;
	}
	if (isfixed) {
	  printf("Module %d was FIXED!\n",idx);
	  if (alig->GetPreAlignmentQualityFactor(idx)>0) 
	    QF=alig->GetPreAlignmentQualityFactor(idx);
	  else 
	    QF=0;
	}
      }

    }

    new(alobj[idx+1]) AliAlignObjParams(symname,volid,t[0],t[1],t[2],r[0],r[1],r[2],kFALSE);
    AliAlignObjParams* alo = (AliAlignObjParams*) array->UncheckedAt(idx+1);
    alo->SetCorrMatrix(corr);
    alo->SetUniqueID(QF);
  }
  
  fout->WriteObject(array,"ITSAlignObjs","kSingleKey");
  fout->Close();

  stimer.Stop();
  stimer.Print();

  return 0;
}
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[nmod6];
71	Double_t errors = new Double_t[nmod6];
72	Double_t pulls = new Double_t[nmod6];
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[iidx6+kk]==0.0 && errors[iidx6+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[iidx6+kk]!=0.0 \|\| pulls[iidx6+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[iidx6+0]errors[iidx*6+0];
220	corr[2] = errors[iidx6+1]errors[iidx*6+1];
221	corr[5] = errors[iidx6+2]errors[iidx*6+2];
222	corr[9] = errors[iidx6+3]errors[iidx*6+3];
223	corr[14]= errors[iidx6+4]errors[iidx*6+4];
224	corr[20]= errors[iidx6+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	}