-#if !defined(__CINT__) || defined(__MAKECINT__)\r
-#include <Riostream.h>\r
-#include <TChain.h>\r
-#include <TClonesArray.h>\r
-#include <TFile.h>\r
-#include <TMath.h>\r
-#include <TStopwatch.h>\r
-#include "AliAlignObjParams.h"\r
-#include "AliTrackPointArray.h"\r
-#include "AliITSAlignMille.h"\r
-#endif\r
-//********************************************************************\r
-// Example to run ITS alignment using Millepede\r
-//\r
-// Read tracks from AliTrackPoints.N.root and fill Millepede.\r
-// Millepede configuration is taken from AliITSAlignMille.conf\r
-// Results are written as AliAlignObjParams in outfile.\r
-//\r
-// Origin: M. Lunardon\r
-//********************************************************************\r
-\r
-Bool_t SelectLayers(AliTrackPointArray *tpa, int *nreqpts_lay) {\r
- // selection on layers\r
- int npts=tpa->GetNPoints();\r
- int nlay[6]; \r
- for (int jj=0;jj<6;jj++) nlay[jj]=0;\r
- for (int ip=0; ip<npts; ip++) {\r
- int lay=-1;\r
- float r=TMath::Sqrt(tpa->GetX()[ip]*tpa->GetX()[ip] + tpa->GetY()[ip]*tpa->GetY()[ip]);\r
- if (r<5) lay=1;\r
- else if (r>5 && r<10) lay=2;\r
- else if (r>10 && r<18) lay=3;\r
- else if (r>18 && r<30) lay=4;\r
- else if (r>30 && r<40) lay=5;\r
- else if (r>40) lay=6;\r
- if (lay<0) continue;\r
- nlay[lay-1]++;\r
- }\r
- Bool_t isgood=1;\r
- for (int jj=0; jj<6; jj++) {\r
- if (nlay[jj]<nreqpts_lay[jj]) isgood=0;\r
- }\r
- return isgood;\r
-}\r
-//////////////////////////////////////\r
-\r
-int ITSAlignMille(int fromev=1, int toev=200, int maxentries=-1, char *outfile="ITSAlignMille.root") {\r
-\r
- // Read tracks from AliTrackPoints.N.root and fill Millepede.\r
- // Millepede results are written as AliAlignObjParams in outfile.\r
- \r
- int nreqpts=6;\r
- int nreqpts_lay[6]={0,0,0,0,0,0};\r
-\r
- TFile *fout=new TFile(outfile,"recreate");\r
- if (!fout->IsOpen()) {\r
- printf("\nCannot open output file!\n");\r
- return -1;\r
- }\r
-\r
- TChain *chain=new TChain("spTree"); \r
- char dir[100]="AliTrackPoints";\r
- char st[200];\r
- \r
- for (int iev=fromev; iev<=toev; iev++) {\r
- sprintf(st,"%s/AliTrackPoints.%d.root",dir,iev);\r
- chain->Add(st);\r
- }\r
- \r
- int nentries=chain->GetEntries();\r
- printf("There are %d entries in chain\n",nentries);\r
- \r
- if (maxentries>0 && maxentries<nentries) nentries=maxentries;\r
- \r
- AliTrackPointArray *tpa = 0;\r
- chain->SetBranchAddress("SP", &tpa);\r
-\r
- ////////////////////////////////////////////\r
- \r
- AliITSAlignMille *alig = new AliITSAlignMille("AliITSAlignMille.conf");\r
-\r
- Int_t nmod=alig->GetNModules();\r
-\r
- Double_t *parameters = new Double_t[nmod*6];\r
- Double_t *errors = new Double_t[nmod*6];\r
- Double_t *pulls = new Double_t[nmod*6];\r
-\r
- for(Int_t k=0;k<nmod*6;k++) {\r
- parameters[k]=0.;\r
- errors[k]=0.;\r
- pulls[k]=0.;\r
- }\r
-\r
- Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};\r
- alig->InitGlobalParameters(parameters);\r
- alig->Print();\r
-\r
- ////////////////////////////////////////////////////////////////////\r
-\r
- TStopwatch stimer;\r
- stimer.Start();\r
-\r
- int iTrackOk=0; // number of good passed track\r
- // loop on spTree entries\r
- // one entry = one track\r
- for (int i=0; i<nentries; i++) {\r
- chain->GetEvent(i);\r
- //////////////////////////////////////////////////////\r
- // track preselection \r
- int npts=tpa->GetNPoints();\r
- if (npts<nreqpts) continue;\r
- if (!SelectLayers(tpa,nreqpts_lay)) continue;\r
- //////////////////////////////////////////////////////\r
-\r
- if (!alig->ProcessTrack(tpa)) {\r
- if (!(iTrackOk%500)) \r
- printf("Calling LocalFit n. %d\n",iTrackOk);\r
- alig->LocalFit(iTrackOk++,trackParams,0);\r
- }\r
- }\r
- printf("\nMillepede fed with %d tracks\n\n",iTrackOk);\r
- \r
- alig->GlobalFit(parameters,errors,pulls);\r
- \r
- cout << "Done with GlobalFit " << endl;\r
- \r
- ////////////////////////////////////////////////////////////\r
- // output\r
- \r
-\r
- TClonesArray *array = new TClonesArray("AliAlignObjParams",4000);\r
- TClonesArray &alobj = *array;\r
-\r
- // first object: ITS\r
- new(alobj[0]) AliAlignObjParams("ITS", 0, 0., 0., 0., 0., 0., 0., kTRUE);\r
- \r
- UShort_t volid;\r
- const Char_t *symname;\r
- Double_t dx,dy,dz,dpsi,dtheta,dphi;\r
- Double_t corr[21];\r
-\r
- // all ITS modules\r
- for (int idx=0; idx<2198; idx++) {\r
- volid=alig->GetModuleVolumeID(idx);\r
- symname = AliGeomManager::SymName(volid);\r
- for (int jj=0;jj<21;jj++) corr[jj]=0.0;\r
-\r
- if (alig->CheckVolumeID(volid)) { // defined modules\r
- alig->SetCurrentModule(idx);\r
- int iidx=alig->GetCurrentModuleInternalIndex();\r
- dx = parameters[iidx*6+0]; corr[0] = errors[iidx*6+0]*errors[iidx*6+0];\r
- dy = parameters[iidx*6+1]; corr[2] = errors[iidx*6+1]*errors[iidx*6+1];\r
- dz = parameters[iidx*6+2]; corr[5] = errors[iidx*6+2]*errors[iidx*6+2];\r
- dpsi = parameters[iidx*6+3]; corr[9] = errors[iidx*6+3]*errors[iidx*6+3];\r
- dtheta= parameters[iidx*6+4]; corr[14]= errors[iidx*6+4]*errors[iidx*6+4];\r
- dphi = parameters[iidx*6+5]; corr[20]= errors[iidx*6+5]*errors[iidx*6+5];\r
- }\r
- else { // other modules\r
- dx=0.0; dy=0.0; dz=0.0; dpsi=0.0; dtheta=0.0; dphi=0.0;\r
- }\r
- \r
- new(alobj[idx+1]) AliAlignObjParams(symname, volid, dx, dy, dz, dpsi, dtheta, dphi, kFALSE);\r
- AliAlignObjParams* alo = (AliAlignObjParams*) array->UncheckedAt(idx+1);\r
- alo->SetCorrMatrix(corr);\r
- }\r
- \r
- fout->WriteObject(array,"ITSAlignObjs","kSingleKey");\r
- fout->Close();\r
-\r
- stimer.Stop();\r
- stimer.Print();\r
- \r
- return 0;\r
-}\r
+#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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff