ROOT logo
#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;
}
 ITSAlignMille.C:1
 ITSAlignMille.C:2
 ITSAlignMille.C:3
 ITSAlignMille.C:4
 ITSAlignMille.C:5
 ITSAlignMille.C:6
 ITSAlignMille.C:7
 ITSAlignMille.C:8
 ITSAlignMille.C:9
 ITSAlignMille.C:10
 ITSAlignMille.C:11
 ITSAlignMille.C:12
 ITSAlignMille.C:13
 ITSAlignMille.C:14
 ITSAlignMille.C:15
 ITSAlignMille.C:16
 ITSAlignMille.C:17
 ITSAlignMille.C:18
 ITSAlignMille.C:19
 ITSAlignMille.C:20
 ITSAlignMille.C:21
 ITSAlignMille.C:22
 ITSAlignMille.C:23
 ITSAlignMille.C:24
 ITSAlignMille.C:25
 ITSAlignMille.C:26
 ITSAlignMille.C:27
 ITSAlignMille.C:28
 ITSAlignMille.C:29
 ITSAlignMille.C:30
 ITSAlignMille.C:31
 ITSAlignMille.C:32
 ITSAlignMille.C:33
 ITSAlignMille.C:34
 ITSAlignMille.C:35
 ITSAlignMille.C:36
 ITSAlignMille.C:37
 ITSAlignMille.C:38
 ITSAlignMille.C:39
 ITSAlignMille.C:40
 ITSAlignMille.C:41
 ITSAlignMille.C:42
 ITSAlignMille.C:43
 ITSAlignMille.C:44
 ITSAlignMille.C:45
 ITSAlignMille.C:46
 ITSAlignMille.C:47
 ITSAlignMille.C:48
 ITSAlignMille.C:49
 ITSAlignMille.C:50
 ITSAlignMille.C:51
 ITSAlignMille.C:52
 ITSAlignMille.C:53
 ITSAlignMille.C:54
 ITSAlignMille.C:55
 ITSAlignMille.C:56
 ITSAlignMille.C:57
 ITSAlignMille.C:58
 ITSAlignMille.C:59
 ITSAlignMille.C:60
 ITSAlignMille.C:61
 ITSAlignMille.C:62
 ITSAlignMille.C:63
 ITSAlignMille.C:64
 ITSAlignMille.C:65
 ITSAlignMille.C:66
 ITSAlignMille.C:67
 ITSAlignMille.C:68
 ITSAlignMille.C:69
 ITSAlignMille.C:70
 ITSAlignMille.C:71
 ITSAlignMille.C:72
 ITSAlignMille.C:73
 ITSAlignMille.C:74
 ITSAlignMille.C:75
 ITSAlignMille.C:76
 ITSAlignMille.C:77
 ITSAlignMille.C:78
 ITSAlignMille.C:79
 ITSAlignMille.C:80
 ITSAlignMille.C:81
 ITSAlignMille.C:82
 ITSAlignMille.C:83
 ITSAlignMille.C:84
 ITSAlignMille.C:85
 ITSAlignMille.C:86
 ITSAlignMille.C:87
 ITSAlignMille.C:88
 ITSAlignMille.C:89
 ITSAlignMille.C:90
 ITSAlignMille.C:91
 ITSAlignMille.C:92
 ITSAlignMille.C:93
 ITSAlignMille.C:94
 ITSAlignMille.C:95
 ITSAlignMille.C:96
 ITSAlignMille.C:97
 ITSAlignMille.C:98
 ITSAlignMille.C:99
 ITSAlignMille.C:100
 ITSAlignMille.C:101
 ITSAlignMille.C:102
 ITSAlignMille.C:103
 ITSAlignMille.C:104
 ITSAlignMille.C:105
 ITSAlignMille.C:106
 ITSAlignMille.C:107
 ITSAlignMille.C:108
 ITSAlignMille.C:109
 ITSAlignMille.C:110
 ITSAlignMille.C:111
 ITSAlignMille.C:112
 ITSAlignMille.C:113
 ITSAlignMille.C:114
 ITSAlignMille.C:115
 ITSAlignMille.C:116
 ITSAlignMille.C:117
 ITSAlignMille.C:118
 ITSAlignMille.C:119
 ITSAlignMille.C:120
 ITSAlignMille.C:121
 ITSAlignMille.C:122
 ITSAlignMille.C:123
 ITSAlignMille.C:124
 ITSAlignMille.C:125
 ITSAlignMille.C:126
 ITSAlignMille.C:127
 ITSAlignMille.C:128
 ITSAlignMille.C:129
 ITSAlignMille.C:130
 ITSAlignMille.C:131
 ITSAlignMille.C:132
 ITSAlignMille.C:133
 ITSAlignMille.C:134
 ITSAlignMille.C:135
 ITSAlignMille.C:136
 ITSAlignMille.C:137
 ITSAlignMille.C:138
 ITSAlignMille.C:139
 ITSAlignMille.C:140
 ITSAlignMille.C:141
 ITSAlignMille.C:142
 ITSAlignMille.C:143
 ITSAlignMille.C:144
 ITSAlignMille.C:145
 ITSAlignMille.C:146
 ITSAlignMille.C:147
 ITSAlignMille.C:148
 ITSAlignMille.C:149
 ITSAlignMille.C:150
 ITSAlignMille.C:151
 ITSAlignMille.C:152
 ITSAlignMille.C:153
 ITSAlignMille.C:154
 ITSAlignMille.C:155
 ITSAlignMille.C:156
 ITSAlignMille.C:157
 ITSAlignMille.C:158
 ITSAlignMille.C:159
 ITSAlignMille.C:160
 ITSAlignMille.C:161
 ITSAlignMille.C:162
 ITSAlignMille.C:163
 ITSAlignMille.C:164
 ITSAlignMille.C:165
 ITSAlignMille.C:166
 ITSAlignMille.C:167
 ITSAlignMille.C:168
 ITSAlignMille.C:169
 ITSAlignMille.C:170
 ITSAlignMille.C:171
 ITSAlignMille.C:172
 ITSAlignMille.C:173
 ITSAlignMille.C:174
 ITSAlignMille.C:175
 ITSAlignMille.C:176
 ITSAlignMille.C:177
 ITSAlignMille.C:178
 ITSAlignMille.C:179
 ITSAlignMille.C:180
 ITSAlignMille.C:181
 ITSAlignMille.C:182
 ITSAlignMille.C:183
 ITSAlignMille.C:184
 ITSAlignMille.C:185
 ITSAlignMille.C:186
 ITSAlignMille.C:187
 ITSAlignMille.C:188
 ITSAlignMille.C:189
 ITSAlignMille.C:190
 ITSAlignMille.C:191
 ITSAlignMille.C:192
 ITSAlignMille.C:193
 ITSAlignMille.C:194
 ITSAlignMille.C:195
 ITSAlignMille.C:196
 ITSAlignMille.C:197
 ITSAlignMille.C:198
 ITSAlignMille.C:199
 ITSAlignMille.C:200
 ITSAlignMille.C:201
 ITSAlignMille.C:202
 ITSAlignMille.C:203
 ITSAlignMille.C:204
 ITSAlignMille.C:205
 ITSAlignMille.C:206
 ITSAlignMille.C:207
 ITSAlignMille.C:208
 ITSAlignMille.C:209
 ITSAlignMille.C:210
 ITSAlignMille.C:211
 ITSAlignMille.C:212
 ITSAlignMille.C:213
 ITSAlignMille.C:214
 ITSAlignMille.C:215
 ITSAlignMille.C:216
 ITSAlignMille.C:217
 ITSAlignMille.C:218
 ITSAlignMille.C:219
 ITSAlignMille.C:220
 ITSAlignMille.C:221
 ITSAlignMille.C:222
 ITSAlignMille.C:223
 ITSAlignMille.C:224
 ITSAlignMille.C:225
 ITSAlignMille.C:226
 ITSAlignMille.C:227
 ITSAlignMille.C:228
 ITSAlignMille.C:229
 ITSAlignMille.C:230
 ITSAlignMille.C:231
 ITSAlignMille.C:232
 ITSAlignMille.C:233
 ITSAlignMille.C:234
 ITSAlignMille.C:235
 ITSAlignMille.C:236
 ITSAlignMille.C:237
 ITSAlignMille.C:238
 ITSAlignMille.C:239
 ITSAlignMille.C:240
 ITSAlignMille.C:241
 ITSAlignMille.C:242
 ITSAlignMille.C:243
 ITSAlignMille.C:244
 ITSAlignMille.C:245
 ITSAlignMille.C:246
 ITSAlignMille.C:247
 ITSAlignMille.C:248
 ITSAlignMille.C:249
 ITSAlignMille.C:250
 ITSAlignMille.C:251
 ITSAlignMille.C:252
 ITSAlignMille.C:253
 ITSAlignMille.C:254
 ITSAlignMille.C:255
 ITSAlignMille.C:256
 ITSAlignMille.C:257
 ITSAlignMille.C:258