ITS/ITSAlignMille.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)\r
   2 #include <Riostream.h>\r
   3 #include <TChain.h>\r
   4 #include <TClonesArray.h>\r
   5 #include <TFile.h>\r
   6 #include <TMath.h>\r
   7 #include <TStopwatch.h>\r
   8 #include "AliAlignObjParams.h"\r
   9 #include "AliTrackPointArray.h"\r
  10 #include "AliITSAlignMille.h"\r
  11 #endif\r
  12 //********************************************************************\r
  13 //  Example to run ITS alignment using Millepede\r
  14 //\r
  15 //  Read tracks from AliTrackPoints.N.root and fill Millepede.\r
  16 //  Millepede configuration is taken from AliITSAlignMille.conf\r
  17 //  Results are written as AliAlignObjParams in outfile.\r
  18 //\r
  19 //      Origin: M. Lunardon\r
  20 //********************************************************************\r
  21 \r
  22 Bool_t SelectLayers(AliTrackPointArray *tpa, int *nreqpts_lay) {\r
  23   // selection on layers\r
  24   int npts=tpa->GetNPoints();\r
  25   int nlay[6]; \r
  26   for (int jj=0;jj<6;jj++) nlay[jj]=0;\r
  27   for (int ip=0; ip<npts; ip++) {\r
  28     int lay=-1;\r
  29     float r=TMath::Sqrt(tpa->GetX()[ip]*tpa->GetX()[ip] + tpa->GetY()[ip]*tpa->GetY()[ip]);\r
  30     if (r<5) lay=1;\r
  31     else if (r>5 && r<10) lay=2;\r
  32     else if (r>10 && r<18) lay=3;\r
  33     else if (r>18 && r<30) lay=4;\r
  34     else if (r>30 && r<40) lay=5;\r
  35     else if (r>40) lay=6;\r
  36     if (lay<0) continue;\r
  37     nlay[lay-1]++;\r
  38   }\r
  39   Bool_t isgood=1;\r
  40   for (int jj=0; jj<6; jj++) {\r
  41     if (nlay[jj]<nreqpts_lay[jj]) isgood=0;\r
  42   }\r
  43   return isgood;\r
  44 }\r
  45 //////////////////////////////////////\r
  46 \r
  47 int ITSAlignMille(int fromev=1, int toev=200, int maxentries=-1, char *outfile="ITSAlignMille.root") {\r
  48 \r
  49   // Read tracks from AliTrackPoints.N.root and fill Millepede.\r
  50   // Millepede results are written as AliAlignObjParams in outfile.\r
  51   \r
  52   int nreqpts=6;\r
  53   int nreqpts_lay[6]={0,0,0,0,0,0};\r
  54 \r
  55   TFile *fout=new TFile(outfile,"recreate");\r
  56   if (!fout->IsOpen()) {\r
  57     printf("\nCannot open output file!\n");\r
  58     return -1;\r
  59   }\r
  60 \r
  61   TChain *chain=new TChain("spTree");  \r
  62   char dir[100]="AliTrackPoints";\r
  63   char st[200];\r
  64   \r
  65   for (int iev=fromev; iev<=toev; iev++) {\r
  66     sprintf(st,"%s/AliTrackPoints.%d.root",dir,iev);\r
  67     chain->Add(st);\r
  68   }\r
  69   \r
  70   int nentries=chain->GetEntries();\r
  71   printf("There are %d entries in chain\n",nentries);\r
  72   \r
  73   if (maxentries>0 && maxentries<nentries) nentries=maxentries;\r
  74   \r
  75   AliTrackPointArray *tpa = 0;\r
  76   chain->SetBranchAddress("SP", &tpa);\r
  77 \r
  78   ////////////////////////////////////////////\r
  79   \r
  80   AliITSAlignMille *alig = new AliITSAlignMille("AliITSAlignMille.conf");\r
  81 \r
  82   Int_t nmod=alig->GetNModules();\r
  83 \r
  84   Double_t *parameters = new Double_t[nmod*6];\r
  85   Double_t *errors = new Double_t[nmod*6];\r
  86   Double_t *pulls = new Double_t[nmod*6];\r
  87 \r
  88   for(Int_t k=0;k<nmod*6;k++) {\r
  89     parameters[k]=0.;\r
  90     errors[k]=0.;\r
  91     pulls[k]=0.;\r
  92   }\r
  93 \r
  94   Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};\r
  95   alig->InitGlobalParameters(parameters);\r
  96   alig->Print();\r
  97 \r
  98   ////////////////////////////////////////////////////////////////////\r
  99 \r
 100   TStopwatch stimer;\r
 101   stimer.Start();\r
 102 \r
 103   int iTrackOk=0; // number of good passed track\r
 104   // loop on spTree entries\r
 105   // one entry = one track\r
 106   for (int i=0; i<nentries; i++) {\r
 107     chain->GetEvent(i);\r
 108     //////////////////////////////////////////////////////\r
 109     // track preselection    \r
 110     int npts=tpa->GetNPoints();\r
 111     if (npts<nreqpts) continue;\r
 112     if (!SelectLayers(tpa,nreqpts_lay)) continue;\r
 113     //////////////////////////////////////////////////////\r
 114 \r
 115     if (!alig->ProcessTrack(tpa)) {\r
 116       if (!(iTrackOk%500)) \r
 117         printf("Calling LocalFit n. %d\n",iTrackOk);\r
 118       alig->LocalFit(iTrackOk++,trackParams,0);\r
 119     }\r
 120   }\r
 121   printf("\nMillepede fed with %d tracks\n\n",iTrackOk);\r
 122   \r
 123   alig->GlobalFit(parameters,errors,pulls);\r
 124   \r
 125   cout << "Done with GlobalFit " << endl;\r
 126   \r
 127   ////////////////////////////////////////////////////////////\r
 128   // output\r
 129   \r
 130 \r
 131   TClonesArray *array = new TClonesArray("AliAlignObjParams",4000);\r
 132   TClonesArray &alobj = *array;\r
 133 \r
 134   // first object: ITS\r
 135   new(alobj[0]) AliAlignObjParams("ITS", 0, 0., 0., 0., 0., 0., 0., kTRUE);\r
 136   \r
 137   UShort_t volid;\r
 138   const Char_t *symname;\r
 139   Double_t dx,dy,dz,dpsi,dtheta,dphi;\r
 140   Double_t corr[21];\r
 141 \r
 142   // all ITS modules\r
 143   for (int idx=0; idx<2198; idx++) {\r
 144     volid=alig->GetModuleVolumeID(idx);\r
 145     symname = AliGeomManager::SymName(volid);\r
 146     for (int jj=0;jj<21;jj++) corr[jj]=0.0;\r
 147 \r
 148     if (alig->CheckVolumeID(volid)) { // defined modules\r
 149       alig->SetCurrentModule(idx);\r
 150       int iidx=alig->GetCurrentModuleInternalIndex();\r
 151       dx    = parameters[iidx*6+0];  corr[0] = errors[iidx*6+0]*errors[iidx*6+0];\r
 152       dy    = parameters[iidx*6+1];  corr[2] = errors[iidx*6+1]*errors[iidx*6+1];\r
 153       dz    = parameters[iidx*6+2];  corr[5] = errors[iidx*6+2]*errors[iidx*6+2];\r
 154       dpsi  = parameters[iidx*6+3];  corr[9] = errors[iidx*6+3]*errors[iidx*6+3];\r
 155       dtheta= parameters[iidx*6+4];  corr[14]= errors[iidx*6+4]*errors[iidx*6+4];\r
 156       dphi  = parameters[iidx*6+5];  corr[20]= errors[iidx*6+5]*errors[iidx*6+5];\r
 157     }\r
 158     else { // other modules\r
 159       dx=0.0; dy=0.0; dz=0.0; dpsi=0.0; dtheta=0.0; dphi=0.0;\r
 160     }\r
 161     \r
 162     new(alobj[idx+1]) AliAlignObjParams(symname, volid, dx, dy, dz, dpsi, dtheta, dphi, kFALSE);\r
 163     AliAlignObjParams* alo = (AliAlignObjParams*) array->UncheckedAt(idx+1);\r
 164     alo->SetCorrMatrix(corr);\r
 165   }\r
 166   \r
 167   fout->WriteObject(array,"ITSAlignObjs","kSingleKey");\r
 168   fout->Close();\r
 169 \r
 170   stimer.Stop();\r
 171   stimer.Print();\r
 172    \r
 173   return 0;\r
 174 }\r