--- /dev/null
+#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
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff