Corrected UInt_t <-> Int_t conversion
[u/mrichter/AliRoot.git] / TPC / AliTPCExBEffectiveSector.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 ////////////////////////////////////////////////////////////////////////////
17 //                                                                        //
18 // AliTPCExBEffectiveSector class                                                   //
19 // Correct for the rest of ExB effect which are not covered yet by physical models
20 //
21 // Motivation:
22 //   ExB correction: 
23 //      dr    =  c0* integral(Er/Ez) + c1* integral(Erphi/Ez)
24 //      drphi = -c1* integral(Er/Ez) + c0* integral(Erphi/Ez)
25 //   Where:      
26 //   wt = Bz*(k*vdrift/E)           ~ 0.3 at B=0.5 T 
27 //   c0 = 1/(1+T2*T2*wt*wt) 
28 //   c1 = T1*wt/(1+T1*T1*wt*wt)
29 //   
30 //  
31 //  3 correction maps 0 implemented as histogram used
32 //  R-Phi correction map obtained minimizing residuals betwee the track
33 //        and space points (AliTPCcalibAlign class). Track is defined using
34 //        the points from the refernce plain at the middle of the TPC
35 //        and vertex
36 //        Corrected primar tracks straight pointing to the primary vertex
37 //
38 //  R distortion - obtained using the cluster residuals in the setup with
39 //                 plus and minus field 
40 //                 Only high momenta tracks used for this calibration (1 GeV threshold)
41 //     drphi_plus-drphi_minus=-2*c1 integral(Er/Ez)
42 //               - Erphi/Ez cancels   
43 ////////////////////////////////////////////////////////////////////////////
44 #include "AliMagF.h"
45 #include "TGeoGlobalMagField.h"
46 #include "AliTPCcalibDB.h"
47 #include "AliTPCParam.h"
48 #include "AliLog.h"
49
50 #include "TMath.h"
51 #include "AliTPCROC.h"
52 #include "TFile.h"
53 #include "TAxis.h"
54 #include "TTree.h"
55 #include "TTreeStream.h"
56 #include "THnSparse.h"
57 #include "TProfile.h"
58 #include "TH2F.h"
59 #include "TH3F.h"
60 #include "TROOT.h"
61 #include "AliTPCExBEffectiveSector.h"
62 ClassImp(AliTPCExBEffectiveSector)
63
64 AliTPCExBEffectiveSector::AliTPCExBEffectiveSector()
65   : AliTPCCorrection("ExB_effectiveSector","ExB effective sector"),
66     fC0(1.),fC1(0.), 
67     fCorrectionR(0),        // radial correction
68     fCorrectionRPhi(0),     // r-phi correction
69     fCorrectionZ(0)        // z correction
70 {
71   //
72   // default constructor
73   //
74 }
75
76 AliTPCExBEffectiveSector::~AliTPCExBEffectiveSector() {
77   //
78   // default destructor
79   //
80   delete fCorrectionR;        // radial correction
81   delete fCorrectionRPhi;     // r-phi correction
82   delete fCorrectionZ;        // z correction
83 }
84
85
86
87 void AliTPCExBEffectiveSector::Init() {
88   //
89   // Initialization funtion
90   //
91   
92   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
93   if (!magF) AliError("Magneticd field - not initialized");
94   Double_t bzField = magF->SolenoidField()/10.; //field in T
95   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
96   if (!param) AliError("Parameters - not initialized");
97   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
98   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
99   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
100   // Correction Terms for effective omegaTau; obtained by a laser calibration run
101   SetOmegaTauT1T2(wt,fT1,fT2);
102 }
103
104 void AliTPCExBEffectiveSector::Update(const TTimeStamp &/*timeStamp*/) {
105   //
106   // Update function 
107   //
108   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
109   if (!magF) AliError("Magneticd field - not initialized");
110   Double_t bzField = magF->SolenoidField()/10.; //field in T
111   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
112   if (!param) AliError("Parameters - not initialized");
113   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
114   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
115   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
116   // Correction Terms for effective omegaTau; obtained by a laser calibration run
117   SetOmegaTauT1T2(wt,fT1,fT2);
118 }
119
120
121
122 void AliTPCExBEffectiveSector::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
123   //
124   // Calculates the correction using the lookup table (histogram) of distortion
125   // The histogram is created as poscl - postrack
126   //   
127   dx[0]=0;
128   dx[1]=0;
129   dx[2]=0;
130   if (!fCorrectionRPhi) return;
131   Double_t phi      = TMath::ATan2(x[1],x[0]);
132   Double_t r        = TMath::Sqrt(x[1]*x[1]+x[0]*x[0]);
133   Double_t sector   = 9.*phi/TMath::Pi();
134   if (sector<0) sector+=18.;
135   Double_t        kZ=x[2]/r;
136   //
137   if (kZ>1.2)        kZ= 1.2;
138   if (kZ<-1.2)       kZ= -1.2;
139   if (roc%36<18)  kZ= TMath::Abs(kZ);
140   if (roc%36>=18) kZ=-TMath::Abs(kZ);
141   if (TMath::Abs(kZ)<0.15){
142     kZ = (roc%36<18) ? 0.15:-0.15;
143   }  
144   //
145   Double_t dlR=0;
146   Double_t dlRPhi=0;
147   Double_t dlZ=0;
148   Double_t rr=TMath::Max(r,fCorrectionRPhi->GetYaxis()->GetXmin()+0.01);
149   rr=TMath::Min(rr,fCorrectionRPhi->GetYaxis()->GetXmax()-0.01);
150   Double_t kZZ=TMath::Max(kZ,fCorrectionRPhi->GetZaxis()->GetXmin()+0.001);
151   kZZ=TMath::Min(kZZ,fCorrectionRPhi->GetZaxis()->GetXmax()-0.001);
152
153   if (fCorrectionRPhi) {  
154     //    dlRPhi= -fCorrectionRPhi->Interpolate(sector,rr,kZZ);
155     dlRPhi= -fCorrectionRPhi->GetBinContent(fCorrectionRPhi->FindBin(sector,rr,kZZ));
156   }
157   if (fCorrectionR)    {
158     //    dlR= -fCorrectionR->Interpolate(sector,rr,kZZ);
159     dlR= -fCorrectionR->GetBinContent(fCorrectionR->FindBin(sector,rr,kZZ));
160   }
161   if (fCorrectionZ)    {
162     //    dlZ= -fCorrectionZ->Interpolate(sector,rr,kZZ);
163     dlZ= -fCorrectionZ->GetBinContent(fCorrectionZ->FindBin(sector,rr,kZZ));
164   }
165   Double_t dr    = fC0*dlR  + fC1*dlRPhi;
166   Double_t drphi = -fC1*dlR + fC0*dlRPhi;
167    // Calculate distorted position
168   if ( r > 0.0 ) {
169     r   =  r   + dr;
170     phi =  phi + drphi/r;
171   }
172   // Calculate correction in cartesian coordinates
173   dx[0] = r * TMath::Cos(phi) - x[0];
174   dx[1] = r * TMath::Sin(phi) - x[1];
175   dx[2] = dlZ; 
176
177 }
178
179 void AliTPCExBEffectiveSector::Print(const Option_t* option) const {
180   //
181   // Print function to check the settings (e.g. the twist in the X direction)
182   // option=="a" prints the C0 and C1 coefficents for calibration purposes
183   //
184
185   TString opt = option; opt.ToLower();
186   printf("%s\t%s\n",GetName(),GetTitle());  
187   if (opt.Contains("a")) { // Print all details
188     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
189     printf(" - C0: %1.4f, C1: %1.4f \n",fC0,fC1);
190   }    
191 }
192
193 void  AliTPCExBEffectiveSector::MakeResidualMap(THnSparse * hisInput, const char *sname, Int_t ptype, Int_t dtype){
194   //
195   // Make cluster residual map from the n-dimensional histogram
196   // hisInput supposed to have given format:
197   //          - 4 Dim  - delta, sector, localX, kZ
198   // Vertex position assumed to be at (0,0,0)          
199   TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
200   //
201   Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
202   Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
203   Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
204   TF1 *fgaus=0;
205   TH3F * hisResMap3D = 
206     new TH3F("his3D","his3D",
207              nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
208              nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
209              nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
210   hisResMap3D->GetXaxis()->SetTitle("sector");
211   hisResMap3D->GetYaxis()->SetTitle("localX");
212   hisResMap3D->GetZaxis()->SetTitle("kZ");
213
214   TH2F * hisResMap2D[4] ={0,0,0,0};
215   for (Int_t i=0; i<4; i++){
216     hisResMap2D[i]=
217       new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
218                nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
219                nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
220     hisResMap2D[i]->GetXaxis()->SetTitle("sector");
221     hisResMap2D[i]->GetYaxis()->SetTitle("localX");
222   }
223   //
224   //
225   //
226   TF1 * f1= 0;
227   Int_t axis0[4]={0,1,2,3};
228   Int_t axis1[4]={0,1,2,3};
229   for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
230     // phi- sector  range
231     hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1);
232     THnSparse *his1=hisInput->Projection(4,axis0); 
233     Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1);
234     //
235     for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
236       // local x range
237       // kz fits
238       his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1);
239       THnSparse *his2=his1->Projection(4,axis1); 
240       Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2);
241       //      
242       //A side
243       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
244       TH1 * hisA = his2->Projection(0);
245       Double_t meanA= hisA->GetMean();
246       Double_t rmsA= hisA->GetRMS();
247       Double_t entriesA= hisA->GetEntries();
248       delete hisA;
249       //C side
250       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
251       TH1 * hisC = his2->Projection(0);
252       Double_t meanC= hisC->GetMean();
253       Double_t rmsC= hisC->GetRMS();
254       Double_t entriesC= hisC->GetEntries();
255       delete hisC;
256       his2->GetAxis(3)->SetRangeUser(-1.2,1.2);      
257       TH2 * hisAC       = his2->Projection(0,3);
258       TProfile *profAC  = hisAC->ProfileX(); 
259       delete hisAC;
260       profAC->Fit("pol1","QNR","QNR",0.05,1);
261       if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
262       Double_t offsetA=f1->GetParameter(0);
263       Double_t slopeA=f1->GetParameter(1);
264       Double_t offsetAE=f1->GetParError(0);
265       Double_t slopeAE=f1->GetParError(1); 
266       Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters();
267       profAC->Fit("pol1","QNR","QNR",-1.1,-0.1);
268       f1=(TF1*)gROOT->FindObject("pol1");
269       Double_t offsetC=f1->GetParameter(0);
270       Double_t slopeC=f1->GetParameter(1); 
271       Double_t offsetCE=f1->GetParError(0);
272       Double_t slopeCE=f1->GetParError(1); 
273       Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters();
274       printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C);
275
276       (*pcstream)<<"deltaFit"<<
277         "sector="<<sector<<
278         "localX="<<localX<<
279         "meanA="<<meanA<<
280         "rmsA="<<rmsA<<
281         "entriesA="<<entriesA<<
282         "meanC="<<meanC<<
283         "rmsC="<<rmsC<<
284         "entriesC="<<entriesC<<
285         "offsetA="<<offsetA<<
286         "slopeA="<<slopeA<<
287         "offsetAE="<<offsetAE<<
288         "slopeAE="<<slopeAE<<
289         "chi2A="<<chi2A<<
290         "offsetC="<<offsetC<<
291         "slopeC="<<slopeC<<
292         "offsetCE="<<offsetCE<<
293         "slopeCE="<<slopeCE<<
294         "chi2C="<<chi2C<<
295         "\n";
296       //
297       hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA);
298       hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA);
299       hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC);
300       hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC);
301       
302       for (Int_t ibin3=1; ibin3<nbins3; ibin3++){
303         Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3);
304         if (TMath::Abs(kZ)<0.05) continue;  // crossing 
305         his2->GetAxis(3)->SetRange(ibin3,ibin3);
306         if (TMath::Abs(kZ)>0.15){
307           his2->GetAxis(3)->SetRange(ibin3,ibin3);
308         }
309         TH1 * his = his2->Projection(0);
310         Double_t mean= his->GetMean();
311         Double_t rms= his->GetRMS();
312         Double_t entries= his->GetEntries();
313         //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms);
314         hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
315         Double_t phi=TMath::Pi()*sector/9;
316         if (phi>TMath::Pi()) phi+=TMath::Pi();
317         Double_t meanG=0;
318         Double_t rmsG=0;
319         if (entries>50){
320           if (!fgaus) {     
321             his->Fit("gaus","Q","goff");
322             fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
323           }
324           if (fgaus) {
325             his->Fit(fgaus,"Q","goff");
326             meanG=fgaus->GetParameter(1);
327             rmsG=fgaus->GetParameter(2);
328           }
329         }
330         Double_t dsec=sector-Int_t(sector)-0.5;
331         Double_t snp=dsec*TMath::Pi()/9.;
332         (*pcstream)<<"delta"<<
333           "ptype="<<ptype<<
334           "dtype="<<dtype<<
335           "sector="<<sector<<
336           "dsec="<<dsec<<
337           "snp="<<snp<<
338           "phi="<<phi<<
339           "localX="<<localX<<
340           "kZ="<<kZ<<
341           "theta="<<kZ<<
342           "mean="<<mean<<
343           "rms="<<rms<<
344           "meanG="<<meanG<<
345           "rmsG="<<rmsG<<
346           "entries="<<entries<<
347           "meanA="<<meanA<<
348           "rmsA="<<rmsA<<
349           "entriesA="<<entriesA<<
350           "meanC="<<meanC<<
351           "rmsC="<<rmsC<<
352           "entriesC="<<entriesC<<
353           "offsetA="<<offsetA<<
354           "slopeA="<<slopeA<<
355           "chi2A="<<chi2A<<
356           "offsetC="<<offsetC<<
357           "slopeC="<<slopeC<<
358           "chi2C="<<chi2C<<
359           "\n";
360         delete his;
361       }
362       delete his2;
363     }
364     delete his1;
365   }
366   hisResMap3D->Write();
367   hisResMap2D[0]->Write();
368   hisResMap2D[1]->Write();
369   hisResMap2D[2]->Write();
370   hisResMap2D[3]->Write();
371   delete pcstream;
372 }