Eff C++ warning removal (Marian)

[u/mrichter/AliRoot.git] / ITS / ITSAlignMille.C

#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