HMPID/HMPIDFindAlign.C

   1 #include <TSystem.h>
   2 #include <TVector3.h>
   3 #include <TVector2.h>
   4 #include <TNtuple.h>
   5 #include <TFile.h>
   6 #include <TH2F.h>
   7 #include <TH1F.h>
   8 #include <TCanvas.h>
   9 #include <TGraphErrors.h>
  10 #include <TPaveText.h>
  11 #include <TText.h>
  12 #include <TF1.h>
  13 #include <AliHMPIDParam.h>
  14 #include <AliGeomManager.h>
  15 #include <AliCDBEntry.h>
  16 #include <TClonesArray.h>
  17
  18 TVector3 AlignAxis(Int_t ch,Int_t axis);
  19 Double_t FindChCenterD(Int_t ch);
  20 void FindAligners(Int_t ch, Double_t xShift,Double_t yShift, Double_t dist);
  21 TVector3 FindNewCenter(Int_t ch, Double_t xShLoc,Double_t yShLoc, Double_t zShLoc);
  22 void help();
  23
  24 //--------------------------------------------------
  25 void HMPIDFindAlign(Int_t cham=-1)
  26 {
  27   Printf("");
  28   Printf("******************************************");
  29   Printf("* You can type help() for instructions...*");
  30   Printf("******************************************");
  31   Printf("");
  32
  33   gSystem->Sleep(3000);
  34
  35   AliGeomManager::LoadGeometry("geometry.root");
  36
  37   TFile *file = TFile::Open("real_geometry.root");
  38   if(file) {
  39     AliCDBEntry *entry = (AliCDBEntry*)file->Get("AliCDBEntry");
  40     TClonesArray *array = (TClonesArray*)entry->GetObject();
  41     AliGeomManager::ApplyAlignObjsToGeom(*array);
  42 }
  43   AliHMPIDParam::Instance(); //just to print at the beginning...
  44
  45   Int_t ichMin = 0;
  46   Int_t ichMax = 6;
  47   if(cham>=0&&cham<=6) { ichMin = cham; ichMax = cham;}
  48
  49   for(Int_t ich=ichMin;ich<=ichMax;ich++) {
  50
  51     Printf("");
  52     Printf("*********SUMMARY FOR CHAMBER %i **********",ich);
  53     Printf ("Distance from IP of chamber %i = %f (cm)",ich,FindChCenterD(ich));
  54
  55     TVector3 xAl = AlignAxis(ich,0);
  56     TVector3 yAl = AlignAxis(ich,1);
  57
  58     Double_t p1 = 1/(xAl.Z()*xAl.Z());
  59     Double_t p2 = 1/(yAl.Z()*yAl.Z());
  60
  61     Double_t distMean = (p1*xAl.Y()+p2*yAl.Y())/(p1+p2);
  62
  63     FindAligners(ich,xAl.X(),yAl.X(),distMean);
  64   }
  65 }
  66 //--------------------------------------------------
  67 TVector3 AlignAxis(Int_t ch,Int_t axis)
  68 {
  69
  70 // Axis = 0  X Local Axis;
  71 // Axis = 1  Y Local Axis;
  72
  73
  74   TH1F* delta[6];
  75
  76   TFile* file = TFile::Open("histos.root","old");
  77   TNtuple *pNtuple;
  78   if(!(TNtuple*)file->Get("nTupla")) {
  79 // new version because Giacomo uses now a Task to strio the data...
  80   TList *hmpoutput = (TList*)(file->FindObjectAny("hmpoutput"));
  81   pNtuple = (TNtuple*)hmpoutput->At(0);
  82 } else {
  83   pNtuple = (TNtuple*)file->Get("nTupla");
  84 }
  85
  86   for(Int_t j=0;j<6;j++) {
  87     if(axis == 0) {
  88       delta[j] = new TH1F(Form("delta%i",j),Form("delta%i",j),100,-25.,25.);
  89     } else {
  90       delta[j] = new TH1F(Form("delta%i",j),Form("delta%i",j),100,-25.,25.);
  91     }
  92   }
  93   Double_t mean[6];
  94   Double_t sigma[6];
  95   Double_t siz[6];
  96   Double_t bin[6];
  97
  98   TCanvas *c = new TCanvas();
  99
 100   Int_t index = 0;
 101   for(Int_t j=0;j<6;j++) {
 102     Double_t xMin = j*20;
 103     Double_t xMax = xMin+20;
 104     if(axis == 0) {
 105       pNtuple->Draw(Form("Xtrk-Xmip>>delta%i",j),Form("ch==%i&&mipCharge>50&&Xmip>%f&&Xmip<%f",ch,xMin,xMax),"0");
 106     } else {
 107       pNtuple->Draw(Form("Ytrk-Ymip>>delta%i",j),Form("ch==%i&&mipCharge>50&&Ymip>%f&&Ymip<%f",ch,xMin,xMax),"0");
 108     }
 109     if(delta[j]->GetEntries()<50)  continue;
 110
 111     Double_t dIntFit = 5.;
 112     siz[index] = 0.5*(xMax-xMin);
 113     bin[index]= 0.5*(xMax+xMin);
 114     Double_t xMaxPos = delta[j]->GetBinCenter(delta[j]->FindFirstBinAbove(delta[j]->GetMaximum()-1));
 115     delta[j]->Fit("gaus","Q0","",xMaxPos-dIntFit,xMaxPos+dIntFit);
 116     mean[index]  = (delta[j]->GetFunction("gaus"))->GetParameter(1);
 117     sigma[index] = (delta[j]->GetFunction("gaus"))->GetParError(1);
 118     index++;
 119   }
 120   c->Clear();
 121   c->Divide(3,3);
 122
 123   for(Int_t j=0;j<6;j++) {
 124     c->cd(j+1);delta[j]->Draw();
 125   }
 126
 127   c->cd(8);
 128
 129   TGraphErrors *pGr = new TGraphErrors(index,bin,mean,siz,sigma);
 130
 131   pGr->SetTitle(Form("LRS for chamber %i",ch));
 132   if(axis == 0) {
 133     pGr->GetXaxis()->SetTitle("Xmip (cm)");
 134     pGr->GetYaxis()->SetTitle("deltaX (cm)");
 135   } else {
 136     pGr->GetXaxis()->SetTitle("Ymip (cm)");
 137     pGr->GetYaxis()->SetTitle("DeltaY (cm)");
 138   }
 139   pGr->SetMinimum(-5.);
 140   pGr->SetMaximum( 5.);
 141   //pGr->SetMinimum(TMath::MinElement(index,mean)-1.);
 142   //pGr->SetMaximum(TMath::MaxElement(index,mean)+1.);
 143   pGr->Draw("ALP");
 144   pGr->Fit("pol1","Q");
 145   Double_t  shift      = (pGr->GetFunction("pol1"))->GetParameter(0);
 146   Double_t  errshift   = (pGr->GetFunction("pol1"))->GetParError(0);
 147   Double_t  coefang    = (pGr->GetFunction("pol1"))->GetParameter(1);
 148   Double_t  errcoefang = (pGr->GetFunction("pol1"))->GetParError(1);
 149   Double_t  prob       = (pGr->GetFunction("pol1"))->GetProb();
 150   Double_t dist = FindChCenterD(ch);
 151
 152   Double_t rfShift    = -coefang/(1.+coefang)*dist;
 153   Double_t errrfShift = TMath::Abs(dist/((1.+coefang)*(1+coefang))*errcoefang);
 154   if (axis == 0) {
 155     Printf(" x shift = %5.3f +/- %5.3f dist shift = %5.3f +/- %5.3f prob chi2 %5.1f",shift,errshift,rfShift,errrfShift,prob*100);
 156   } else {
 157     Printf(" y shift = %5.3f +/- %5.3f dist shift = %5.3f +/- %5.3f prob chi2 %5.1f",shift,errshift,rfShift,errrfShift,prob*100);
 158   }
 159
 160   c->cd(7);
 161
 162   TPaveText *pPa = new TPaveText();
 163   pPa->DrawPave(0.2,0.2,0.8,0.8);
 164   TText *t = new TText();
 165   t->SetTextSize(0.08);
 166   t->DrawText(0.21,0.70,Form(" shift %5.3f +/- %5.3f",shift,errshift));
 167   t->DrawText(0.21,0.50,Form(" dist shift %5.3f +/- %5.3f",rfShift,errrfShift));
 168   t->DrawText(0.21,0.30,Form(" dist from IP %5.3f ",dist));
 169
 170   if (axis == 0) {
 171     c->SaveAs(Form("ShiftXChamber%i.gif",ch));
 172   } else {
 173     c->SaveAs(Form("ShiftYChamber%i.gif",ch));
 174   }
 175
 176   TVector3 align(shift,rfShift,errrfShift);
 177   return align;
 178 }
 179 //------------------------------------------------------------------
 180 Double_t FindChCenterD(Int_t ch)
 181 {
 182
 183   AliHMPIDParam *pParam = AliHMPIDParam::Instance();
 184   TVector3 xyz = pParam->Lors2Mars(ch,0.5*pParam->SizeAllX(),0.5*pParam->SizeAllY());
 185   return xyz.Mag();
 186 }
 187 //------------------------------------------------------------------
 188 void FindAligners(Int_t ch, Double_t xShift,Double_t yShift, Double_t dist)
 189 {
 190   TVector3 old = FindNewCenter(ch,0,0,0);
 191
 192   Printf("");
 193   Printf("Old X Y Z of the chamber n. %i",ch);
 194   Printf("XC = %f",old.X());
 195   Printf("YC = %f",old.Y());
 196   Printf("ZC = %f",old.Z());
 197   Printf("Dist. from IP = %f",old.Mag());
 198   Printf("");
 199
 200
 201   Printf("Local X shift = %f",xShift);
 202   Printf("Local Y shift = %f",yShift);
 203   Printf("RPhi    shift = %f",dist);
 204
 205   TVector3 shift = FindNewCenter(ch,xShift,yShift,dist);
 206
 207   Printf("");
 208   Printf("New X Y Z of the chamber n. %i",ch);
 209   Printf("");
 210   Printf("XC = %f YC = %f ZC = %f",shift.X(),shift.Y(),shift.Z());
 211   Printf("Dist. from IP = %f",shift.Mag());
 212   Printf("");
 213   Printf("(next line are the values to be inserted in the Align Object)");
 214   Printf("");
 215   Printf("XC shift = %5.2f YC shift = %5.2f ZC shift = %5.2f  for chamber %i",shift.X()-old.X(),shift.Y()-old.Y(),shift.Z()-old.Z(),ch);
 216   Printf("");
 217 }
 218 //------------------------------------------------------------------
 219 TVector3 FindNewCenter(Int_t ch, Double_t xShLoc,Double_t yShLoc, Double_t zShLoc)
 220 {
 221 //First shift in (LocX,LocY) and go in MRS
 222 //Then shift radially...
 223   AliHMPIDParam *pParam = AliHMPIDParam::Instance();
 224   TVector3 xyz = pParam->Lors2Mars(ch,0.5*pParam->SizeAllX()+xShLoc,0.5*pParam->SizeAllY()+yShLoc,-1);
 225   TVector3 xyzSh(1,1,1);
 226   xyzSh.SetTheta(xyz.Theta());
 227   xyzSh.SetPhi(xyz.Phi());
 228   xyzSh.SetMag(xyz.Mag()+zShLoc);
 229   return xyzSh;
 230 }
 231 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 232 void help()
 233 {
 234   Printf("Instructions to run this macros:");
 235   Printf("    1 - In the same dir where the macro runs it must exist the file geometry.root");
 236   Printf("    2 - assign symbolic links as it follows:");
 237   Printf("       a -  file from OCDB/HMPID/Align/Data linked to real_geometry.root");
 238   Printf("       b -  file with the NTupla of ESD Data (produced from AliHMPIDAnalysisTask.C) linked as histos.root");
 239   Printf("    ");
 240   Printf("   The macro have to be run in the following way:");
 241   Printf("    ");
 242   Printf("    > aliroot");
 243   Printf("    ");
 244   Printf("    AliRoot> .L HMPIDFindAlign.C");
 245   Printf("    AliRoot> HMPIDFindAlign(); summary.txt");
 246   Printf("    AliRoot> .q");
 247   Printf("    ");
 248   Printf("The results of the alignment procedures is in the text file summary.txt");
 249   Printf("In addition 14 Canvases (7 chamb. x X and Y) are produced to show the residual distrib. in Local Coord.");
 250   Printf("if some plots are funny, please try to run for a signle chamber modifying some parameters in the macro:");
 251   Printf("until you find a reasonable result.");
 252   Printf("    ");
 253   Printf("AliRoot> HMPIDFindAlign(chamber)");
 254   Printf("    ");
 255   Printf("Ciao!.");
 256 }