ITS/ITSAlignMille.C

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