Line data Source code
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 : /// \class AliTPCcalibAlignInterpolation
18 : /// Class to produce TPC time dependent space point distortion maps using the ITS, TRD and TOF
19 : /// as a reference detector
20 : ///
21 : /// Related to task https://alice.its.cern.ch/jira/browse/ATO-108
22 : /// - code created addopting compiled macro for open gating grid analysis form TPC git:
23 : /// $NOTES/SpaceChargeDistortion/code/spaceChargeDistortions.C
24 : ///
25 : /// \author Marian Ivanov, marian.ivanov@cern.ch
26 :
27 : #include "TROOT.h"
28 : #include "TMath.h"
29 : #include "TFile.h"
30 : #include "TTree.h"
31 : #include "AliESDEvent.h"
32 : #include "AliESDfriend.h"
33 : #include "TTreeStream.h"
34 : #include "AliESDfriendTrack.h"
35 : #include "AliExternalTrackParam.h"
36 : #include "AliTrackPointArray.h"
37 : #include "TChain.h"
38 : #include "AliXRDPROOFtoolkit.h"
39 : #include "AliTrackerBase.h"
40 : #include "AliGeomManager.h"
41 : #include "TVectorF.h"
42 : #include "TVectorD.h"
43 : #include "TStopwatch.h"
44 : #include "TProfile.h"
45 : #include "TGraphErrors.h"
46 : //#include "THnBase.h"
47 : #include "THn.h"
48 : #include "AliSysInfo.h"
49 : #include "TMatrixD.h"
50 : #include "TF1.h"
51 : #include "TDatabasePDG.h"
52 : #include "TTreeStream.h"
53 : #include "TStatToolkit.h"
54 : #include "AliTPCclusterMI.h"
55 : #include "AliTPCseed.h"
56 : #include "AliTPCcalibDB.h"
57 : #include "AliTPCTransform.h"
58 : #include "AliTPCRecoParam.h"
59 : #include "AliTPCreco.h"
60 : #include "AliTPCcalibAlignInterpolation.h"
61 : #include "AliPID.h"
62 : #include "AliCDBManager.h"
63 : #include "AliMagF.h"
64 : #include "AliGRPManager.h"
65 : #include <TGeoGlobalMagField.h>
66 : #include "TSystem.h"
67 : #include "TGrid.h"
68 : #include "TCut.h"
69 : #include "AliNDLocalRegression.h"
70 : #include "AliMathBase.h"
71 : #include "TStyle.h"
72 : #include "TCanvas.h"
73 : #include "AliCDBStorage.h"
74 : #include "AliCDBEntry.h"
75 : #include "AliCDBId.h"
76 : #include "AliTPCParam.h"
77 : #include "AliTPCClusterParam.h"
78 : #include "AliDAQ.h"
79 : #include "AliGRPObject.h"
80 :
81 : const Int_t AliTPCcalibAlignInterpolation_kMaxPoints=500;
82 :
83 6 : ClassImp(AliTPCcalibAlignInterpolation)
84 :
85 :
86 : AliTPCcalibAlignInterpolation::AliTPCcalibAlignInterpolation() :
87 0 : AliTPCcalibBase(),
88 0 : fOnTheFlyFill(0), // flag - on the fly filling of histograms
89 0 : fHisITSDRPhi(0), // TPC-ITS residual histograms
90 0 : fHisITSTRDDRPhi(0), // TPC-ITS+TRD residual histograms
91 0 : fHisITSTOFDRPhi(0), // TPC-ITS_TOF residual histograms
92 0 : fHisITSDZ(0), // TPC-ITS residual histograms
93 0 : fHisITSTRDDZ(0), // TPC-ITS+TRD residual histograms
94 0 : fHisITSTOFDZ(0), // TPC-ITS_TOF residual histograms
95 0 : fRhoTPC(0.9e-3),
96 0 : fX0TPC(28.94),
97 0 : fStreamer(0), // calibration streamer
98 0 : fStreamLevelTrack(0), // stream level - In mode 0 only basic information needed for calibration stored (see EStream
99 0 : fSyswatchStep(100), // dump system resource information after fSyswatchStep tracks
100 0 : fTrackCounter(0) // processed track counter
101 0 : {
102 :
103 0 : }
104 : AliTPCcalibAlignInterpolation::AliTPCcalibAlignInterpolation(const Text_t *name, const Text_t *title, Bool_t onTheFlyFill):
105 0 : AliTPCcalibBase(),
106 0 : fOnTheFlyFill(onTheFlyFill), // flag - on the fly filling of histograms
107 0 : fHisITSDRPhi(0), // TPC-ITS residual histograms
108 0 : fHisITSTRDDRPhi(0), // TPC-ITS+TRD residual histograms
109 0 : fHisITSTOFDRPhi(0), // TPC-ITS_TOF residual histograms
110 0 : fHisITSDZ(0), // TPC-ITS residual histograms
111 0 : fHisITSTRDDZ(0), // TPC-ITS+TRD residual histograms
112 0 : fHisITSTOFDZ(0), // TPC-ITS_TOF residual histograms
113 0 : fRhoTPC(0.9e-3),
114 0 : fX0TPC(28.94),
115 0 : fStreamer(0), // calibration streamer
116 0 : fStreamLevelTrack(0), // stream level - In mode 0 only basic information needed for calibration stored (see EStream
117 0 : fSyswatchStep(100), // dump system resource information after fSyswatchStep tracks
118 0 : fTrackCounter(0) // processed track counter
119 0 : {
120 : // create output histograms
121 0 : SetName(name);
122 0 : SetTitle(title);
123 0 : if (onTheFlyFill) CreateResidualHistosInterpolation();
124 0 : }
125 :
126 0 : AliTPCcalibAlignInterpolation::~AliTPCcalibAlignInterpolation(){
127 : //
128 : //
129 : //
130 0 : if (fStreamer){
131 : // fStreamer->GetFile()->Close();
132 0 : fStreamer->GetFile()->cd();
133 0 : if (fHisITSDRPhi) fHisITSDRPhi->Write();
134 0 : if (fHisITSTRDDRPhi) fHisITSTRDDRPhi->Write();
135 0 : if (fHisITSTOFDRPhi) fHisITSTOFDRPhi->Write();
136 : }
137 0 : delete fStreamer;
138 0 : fStreamer=0;
139 0 : delete fHisITSDRPhi;
140 0 : delete fHisITSTRDDRPhi;
141 0 : delete fHisITSTOFDRPhi;
142 0 : }
143 :
144 :
145 : void AliTPCcalibAlignInterpolation::Terminate(){
146 : //
147 : // Terminate function
148 : // call base terminate + Eval of fitters
149 : //
150 0 : Info("AliTPCcalibAlignInterpolation","Terminate");
151 0 : if (fStreamer){
152 : // fStreamer->GetFile()->Close();
153 0 : fStreamer->GetFile()->cd();
154 0 : if (fHisITSDRPhi) fHisITSDRPhi->Write();
155 0 : if (fHisITSTRDDRPhi) fHisITSTRDDRPhi->Write();
156 0 : if (fHisITSTOFDRPhi) fHisITSTOFDRPhi->Write();
157 : }
158 0 : delete fStreamer;
159 0 : fStreamer=0;
160 :
161 0 : AliTPCcalibBase::Terminate();
162 0 : }
163 :
164 :
165 : Bool_t AliTPCcalibAlignInterpolation::RefitITStrack(AliESDfriendTrack *friendTrack, Double_t mass, AliExternalTrackParam &trackITS,
166 : Double_t &chi2, Double_t &npoints, int* sortInd){
167 : //
168 : // Makes a refit of the ITS track
169 : // Input: AliESDfriendTrack, particle mass, outer ITS TrackParam
170 : // Output: AliExternalTrackParam of the ITS track refit at the last layer of ITS
171 : //
172 : const Double_t sigma2[6][3] = {
173 : {0.002*0.002, 0, 0.013*0.013},
174 : {0.002*0.002, 0, 0.013*0.013},
175 : {0.050*0.050, 0, 0.050*0.050},
176 : {0.050*0.050, 0, 0.050*0.050},
177 : {0.003*0.003, 0, 0.100*0.100},
178 : {0.003*0.003, 0, 0.100*0.100},
179 : }; // ITS intrincsic resolution in (y,z) - error from the points can be used SD layer 2-3 sighnificantly bigger error
180 : // !!!! We should set ITS error parameterization form outside !!!!
181 : const Double_t kMaxRadius=50;
182 : static Int_t sortedIndexLoc[AliTPCcalibAlignInterpolation_kMaxPoints]={0};
183 0 : chi2=0;
184 0 : npoints=0;
185 : //
186 0 : if (friendTrack->GetITSOut()==NULL) return kFALSE;
187 : //
188 0 : trackITS = *((AliExternalTrackParam*)friendTrack->GetITSOut());
189 : // Reset track to the vertex
190 : //if (!AliTrackerBase::PropagateTrackToBxByBz(&trackITS,0,mass,1,kFALSE)) return kFALSE;
191 0 : if (!AliTrackerBase::PropagateTrackParamOnlyTo(&trackITS,0.,5,kFALSE)) return kFALSE;
192 0 : trackITS.ResetCovariance(1000.);
193 :
194 : // Get space points
195 0 : AliTrackPointArray *pointarray = (AliTrackPointArray*)friendTrack->GetTrackPointArray();
196 0 : if (!pointarray){
197 0 : printf("Space points are not stored in the friendTrack!\n");
198 0 : return kFALSE;
199 : };
200 0 : Int_t nPoints = pointarray->GetNPoints(); // # space points of all available detectors
201 : // Sort space points first
202 : int *sortedIndex = sortInd;
203 0 : if (!sortedIndex) {
204 0 : SortPointArray(pointarray, sortedIndexLoc); // space point indices sorted by radius in increasing order
205 : sortedIndex = sortedIndexLoc;
206 0 : }
207 : //
208 : // Propagate track through ITS space points
209 0 : AliTrackPoint spacepoint;
210 : Int_t volId=0,modId=0,layerId=0;
211 :
212 0 : for (Int_t iPoint=0;iPoint<nPoints;iPoint++){
213 0 : pointarray->GetPoint(spacepoint,sortedIndex[iPoint]);
214 0 : int lr = AliGeomManager::VolUIDToLayer(spacepoint.GetVolumeID())-1;
215 0 : if (lr<0||lr>=6) continue;
216 0 : Double_t xyz[3] = {(Double_t)spacepoint.GetX(),(Double_t)spacepoint.GetY(),(Double_t)spacepoint.GetZ()};
217 0 : Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
218 0 : trackITS.Global2LocalPosition(xyz,alpha);
219 0 : if (xyz[0]>kMaxRadius) break; // use only ITS points - maybe we should indexes of elements
220 0 : if (!trackITS.Rotate(alpha)) return kFALSE;
221 0 : if (!AliTrackerBase::PropagateTrackToBxByBz(&trackITS,xyz[0],mass,1,kFALSE)) return kFALSE;
222 0 : Double_t pos[2] = {xyz[1], xyz[2]};
223 0 : const Double_t* cov = sigma2[lr];
224 0 : volId = spacepoint.GetVolumeID();
225 : // layerId = AliGeomManager::VolUIDToLayer(volId,modId);
226 : // if (layerId ==AliGeomManager::kSDD1 || layerId ==AliGeomManager::kSDD2) cov[0]*=16.; cov[2]*=16.;}
227 0 : double chi2cl = trackITS.GetPredictedChi2(pos,cov);
228 0 : chi2 += chi2cl;
229 0 : npoints++;
230 0 : if (!trackITS.Update(pos,cov)) return kFALSE;
231 0 : }
232 0 : return npoints>0;
233 0 : }
234 :
235 :
236 : Bool_t AliTPCcalibAlignInterpolation::RefitTRDtrack(AliESDfriendTrack *friendTrack, Double_t mass, AliExternalTrackParam &trackTRD,
237 : Double_t &chi2, Double_t &npoints, Int_t* sortInd){
238 : //
239 : // Makes a refit of the TRD track using TOF and TRD points
240 : // Input: AliESDfriendTrack, particle mass, inner TRD TrackParam
241 : // Output: AliExternalTrackParam of the TRD track refit - in the first layer of TRD
242 : // Here we forgot about the tiliting pads of TRD - we assume delta Z and delta y are uncorelated
243 : // given approximation is in average tru - in case avearaging of significantly bigger than pad length
244 : //
245 : const Double_t sigmaTRD2[2] = {0.04*0.04, 5*5};
246 : const Double_t sigmaTOF2[2] = {1, 1};
247 : static Int_t sortedIndexLoc[AliTPCcalibAlignInterpolation_kMaxPoints]={0};
248 : const Double_t kMaxRadius=390;
249 : const Double_t kMinRadius=280;
250 : //
251 0 : chi2=0;
252 0 : npoints=0;
253 : //
254 0 : if (friendTrack->GetTRDIn() == NULL) return kFALSE;
255 0 : trackTRD = *((AliExternalTrackParam*)friendTrack->GetTRDIn());
256 :
257 :
258 : // Reset track outside TRD
259 0 : if (!AliTrackerBase::PropagateTrackParamOnlyTo(&trackTRD,kMaxRadius,5,kFALSE)) return kFALSE;
260 : //if (!AliTrackerBase::PropagateTrackToBxByBz(&trackTRD,kMaxRadius,mass,1,kFALSE)) return kFALSE;
261 0 : trackTRD.ResetCovariance(1000.);
262 :
263 : // Get space points
264 0 : AliTrackPointArray *pointarray = (AliTrackPointArray*)friendTrack->GetTrackPointArray();
265 0 : if (!pointarray){
266 0 : printf("Space points are not stored in the friendTrack!\n");
267 0 : return kFALSE;
268 : };
269 0 : Int_t nPoints = pointarray->GetNPoints(); // # space points of all available detectors
270 : // Sort space points first
271 : int *sortedIndex = sortInd;
272 0 : if (!sortedIndex) {
273 0 : SortPointArray(pointarray, sortedIndexLoc); // space point indices sorted by radius in increasing order
274 : sortedIndex = sortedIndexLoc;
275 0 : }
276 :
277 : // Propagate track through TRD space points
278 0 : AliTrackPoint spacepoint;
279 0 : Int_t volId,modId,layerId, npfit=0;
280 0 : for (Int_t iPoint=nPoints;iPoint--;) {
281 0 : pointarray->GetPoint(spacepoint,sortedIndex[iPoint]);
282 0 : volId = spacepoint.GetVolumeID();
283 0 : layerId = AliGeomManager::VolUIDToLayer(volId,modId);
284 0 : if (layerId <AliGeomManager::kTRD1) break;
285 0 : if (layerId>AliGeomManager::kTOF) continue;
286 0 : Double_t xyz[3] = {(Double_t)spacepoint.GetX(),(Double_t)spacepoint.GetY(),(Double_t)spacepoint.GetZ()};
287 0 : Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
288 0 : trackTRD.Global2LocalPosition(xyz,alpha);
289 0 : if (xyz[0]<kMinRadius) break; // use only TRD points
290 0 : if (!trackTRD.Rotate(alpha)) return kFALSE;
291 0 : if (!AliTrackerBase::PropagateTrackToBxByBz(&trackTRD,xyz[0],mass,1,kFALSE)) return kFALSE;
292 0 : Double_t pos[2] = {xyz[1], xyz[2]};
293 0 : Double_t cov[3] = {sigmaTRD2[0],0,sigmaTRD2[1]};
294 0 : if (layerId==AliGeomManager::kTOF) {cov[0]=sigmaTOF2[0]; cov[2]=sigmaTOF2[1];};
295 0 : double chi2cl = trackTRD.GetPredictedChi2(pos,cov);
296 0 : chi2 += chi2cl;
297 0 : if (!trackTRD.Update(pos,cov)) return kFALSE;
298 0 : npfit++;
299 0 : }
300 0 : npoints = npfit;
301 0 : return npoints>0;
302 0 : }
303 :
304 :
305 : Bool_t AliTPCcalibAlignInterpolation::RefitTOFtrack(AliESDfriendTrack *friendTrack, Double_t mass, AliExternalTrackParam &trackTOF,
306 : Double_t &chi2, Double_t &npoints, Int_t* sortInd){
307 : //
308 : // Makes a refit of the TRD track
309 : // Input: AliESDfriendTrack, particle mass, OUTER ITS track - propagated to the TOF point and updated by TOF point
310 : // Output: AliExternalTrackParam of the TOF track refit - at the TOF point
311 : Double_t sigma2[2] = {1., 1.}; // should be parameterized
312 : const Double_t kTOFRadius = 370;
313 : //
314 0 : chi2=0;
315 0 : npoints=0;
316 : //
317 : static Int_t sortedIndexLoc[AliTPCcalibAlignInterpolation_kMaxPoints]={0};
318 : // if (!AliTrackerBase::PropagateTrackParamOnlyTo(&trackTOF,kTOFRadius,15,kTRUE)) return kFALSE;
319 0 : if (!AliTrackerBase::PropagateTrackToBxByBz(&trackTOF,kTOFRadius,mass,10,kTRUE)) return kFALSE;
320 : // RS why don't we reset the cov. matrix here?
321 0 : Int_t volId,modId,layerId;
322 :
323 : // Get space points
324 0 : AliTrackPointArray *pointarray = (AliTrackPointArray*)friendTrack->GetTrackPointArray();
325 0 : if (!pointarray){
326 0 : printf("Space points are not stored in the friendTrack!\n");
327 0 : return kFALSE;
328 : };
329 0 : Int_t nPoints = pointarray->GetNPoints(); // # space points of all available detectors
330 : // Sort space points first
331 : int *sortedIndex = sortInd;
332 0 : if (!sortedIndex) {
333 0 : SortPointArray(pointarray, sortedIndexLoc); // space point indices sorted by radius in increasing order
334 : sortedIndex = sortedIndexLoc;
335 0 : }
336 :
337 : // Propagate track through TRD space points
338 0 : AliTrackPoint spacepoint;
339 : int npfit = 0;
340 0 : for (Int_t iPoint=nPoints;iPoint--;){
341 0 : pointarray->GetPoint(spacepoint,sortedIndex[iPoint]);
342 0 : volId = spacepoint.GetVolumeID();
343 0 : layerId = AliGeomManager::VolUIDToLayer(volId,modId);
344 0 : if (layerId !=AliGeomManager::kTOF) continue;
345 :
346 0 : Double_t xyz[3] = {(Double_t)spacepoint.GetX(),(Double_t)spacepoint.GetY(),(Double_t)spacepoint.GetZ()};
347 0 : Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
348 0 : trackTOF.Global2LocalPosition(xyz,alpha);
349 0 : if (!trackTOF.Rotate(alpha)) return kFALSE;
350 0 : if (!AliTrackerBase::PropagateTrackToBxByBz(&trackTOF,xyz[0],mass,1,kFALSE)) return kFALSE;
351 0 : Double_t pos[2] = {xyz[1], xyz[2]};
352 0 : Double_t cov[3] = {sigma2[0],0,sigma2[1]};
353 0 : double chi2cl = trackTOF.GetPredictedChi2(pos,cov);
354 0 : chi2 += chi2cl;
355 0 : if (!trackTOF.Update(pos,cov)) return kFALSE;
356 : npfit++;
357 0 : break; // there is just 1 TOF poitn
358 0 : }
359 0 : npoints = npfit;
360 0 : return npoints>0;
361 0 : }
362 :
363 : Bool_t AliTPCcalibAlignInterpolation::SortPointArray(AliTrackPointArray *pointarray, Int_t * sortedIndex){
364 : //
365 : // Fill array of indexes to the pointArray (array sorted in increasing order)
366 : //
367 0 : if (sortedIndex==NULL) return kFALSE;
368 0 : Int_t nPoints = pointarray->GetNPoints();
369 0 : if (!nPoints) return kFALSE;
370 0 : Double_t rp[nPoints];
371 0 : const float* x = pointarray->GetX();
372 0 : const float* y = pointarray->GetY();
373 0 : for (Int_t iPoint=nPoints;iPoint--;) rp[iPoint] = x[iPoint]*x[iPoint]+y[iPoint]*y[iPoint];
374 0 : TMath::Sort(nPoints,rp,sortedIndex,kFALSE);
375 : return kTRUE;
376 0 : }
377 :
378 :
379 :
380 : void AliTPCcalibAlignInterpolation::ProcessStandalone(const char * inputList){
381 : //
382 : // Process ESD information standalone without full colibration train
383 : // Input:
384 : // inputList - list of the input ESD files
385 : //
386 : // code from test macro to be set here
387 :
388 0 : }
389 :
390 :
391 :
392 : void AliTPCcalibAlignInterpolation::Process(AliESDEvent *esdEvent){
393 : //
394 : // Create distortion maps out of residual histograms of ITS-TRD interpolation and TPC space points
395 : // JIRA ticket: https://alice.its.cern.ch/jira/browse/ATO-108
396 : //
397 : const Double_t kMaxChi2=10; // max track/track chi2
398 : const Double_t kMaxAlignTolerance=0.1; // max track/track chi2
399 : const Double_t kMaxSNP = 0.8; // max snp tolerated
400 : //
401 : static Bool_t firstCall = kTRUE;
402 0 : if (firstCall) {
403 0 : firstCall = kFALSE;
404 0 : ExtractTPCGasData();
405 0 : }
406 : //
407 0 : AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esdEvent->FindListObject("AliESDfriend"));
408 0 : if (!esdFriend) return;
409 0 : if (esdFriend->TestSkipBit()) return;
410 0 : Int_t nPrimTracks= (esdEvent->GetPrimaryVertex()!=NULL)? esdEvent->GetPrimaryVertex()->GetNContributors():0;
411 0 : Int_t nPrimTracksSPD= (esdEvent->GetPrimaryVertexSPD()!=NULL)? esdEvent->GetPrimaryVertexSPD()->GetNContributors():0;
412 0 : Int_t nTracks = esdEvent->GetNumberOfTracks(); // Get number of tracks in ESD
413 0 : Int_t nSPD= esdEvent->GetMultiplicity()->GetNumberOfITSClusters(0,1);
414 0 : Int_t nSDD= esdEvent->GetMultiplicity()->GetNumberOfITSClusters(2,3);
415 0 : Int_t nSSD= esdEvent->GetMultiplicity()->GetNumberOfITSClusters(4,5);
416 0 : if (nTracks==0) return;
417 0 : if (!fStreamer) fStreamer = new TTreeSRedirector("ResidualHistos.root","recreate");
418 0 : AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
419 0 : TVectorF vecNClTPC(72);
420 0 : TVectorF vecNClTPCused(72);
421 0 : for (Int_t isec=0; isec<72;isec++){
422 0 : vecNClTPC[isec]=esdFriend->GetNclustersTPC(isec);
423 0 : vecNClTPCused[isec]=esdFriend->GetNclustersTPCused(isec);
424 : }
425 0 : ULong64_t gid = esdEvent->GetHeader()->GetEventIdAsLong();
426 0 : Int_t timeStamp= esdEvent->GetTimeStamp();
427 0 : (*fStreamer)<<"eventInfo"<< // store event info - used to calculate per sector currents
428 0 : "gid="<<gid<<
429 0 : "timeStamp="<<timeStamp<<
430 0 : "nSPD="<<nSPD<<
431 0 : "nSDD="<<nSDD<<
432 0 : "nSSD="<<nSSD<<
433 0 : "nPrimTrack="<<nPrimTracks<<
434 0 : "nPrimTrackSPD="<<nPrimTracksSPD<<
435 0 : "nTracks="<<nTracks<<
436 0 : "vecNClTPC.="<<&vecNClTPC<<
437 0 : "vecNClTPCused.="<<&vecNClTPCused<<
438 : "\n";
439 :
440 :
441 : //
442 : const Int_t nPointsAlloc=AliTPCcalibAlignInterpolation_kMaxPoints;
443 : const Int_t kMaxLayer=kMaxRow;
444 0 : AliExternalTrackParam trackArrayITS[kMaxLayer];
445 0 : AliExternalTrackParam trackArrayTRD[kMaxLayer];
446 0 : AliExternalTrackParam trackArrayTOF[kMaxLayer];
447 0 : AliExternalTrackParam trackArrayITSTRD[kMaxLayer];
448 0 : AliExternalTrackParam trackArrayITSTOF[kMaxLayer];
449 0 : AliTPCclusterMI clusterArray[kMaxLayer];
450 : //
451 : //MakeResidualHistosInterpolation();
452 : //
453 0 : Int_t sortedIndex[AliTPCcalibAlignInterpolation_kMaxPoints];
454 0 : TVectorF deltaITS0(kMaxLayer), deltaITS1(kMaxLayer);
455 0 : TVectorF deltaTRD0(kMaxLayer), deltaTRD1(kMaxLayer);
456 0 : TVectorF deltaTOF0(kMaxLayer), deltaTOF1(kMaxLayer);
457 0 : TVectorF vecR(kMaxLayer), vecPhi(kMaxLayer), vecZ(kMaxLayer), vecSec(kMaxLayer);
458 : static int evCnt=0;
459 0 : Bool_t backupUseComposedCorrection = transform->GetCurrentRecoParamNonConst()->GetUseComposedCorrection();
460 0 : transform->GetCurrentRecoParamNonConst()->SetUseComposedCorrection(kFALSE);
461 0 : Bool_t backupUseCorrectionMaps = transform->GetCurrentRecoParamNonConst()->GetUseCorrectionMap();
462 0 : transform->GetCurrentRecoParamNonConst()->SetUseCorrectionMap(kFALSE);
463 0 : Bool_t backupAccountDistortions = transform->GetCurrentRecoParamNonConst()->GetAccountDistortions();
464 0 : transform->GetCurrentRecoParamNonConst()->SetAccountDistortions(kFALSE);
465 :
466 0 : for (Int_t iTrack=0;iTrack<nTracks;iTrack++){ // Track loop
467 : // 0.) For each track in each event, get the AliESDfriendTrack
468 0 : AliESDtrack *esdTrack = esdEvent->GetTrack(iTrack);
469 0 : AliESDfriendTrack *friendTrack = (AliESDfriendTrack*)esdTrack->GetFriendTrack();
470 0 : if (!friendTrack) continue;
471 0 : if (esdTrack->GetITSNcls()<4 || esdTrack->GetTPCNcls()<15) continue;
472 0 : Double_t mass = esdTrack->GetMass(); // particle mass
473 0 : Double_t tofDiff=esdTrack->GetTOFExpTDiffSpec(AliPID::kPion);
474 : // Get TPC seed
475 : TObject *calibObject=0;
476 0 : AliTPCseed *seed = (AliTPCseed*)friendTrack->GetTPCseed();
477 0 : if (!seed) continue;
478 : //
479 : // 1.) Start with AliExternalTrackParam *ITSOut and *TRDIn
480 : //
481 0 : AliExternalTrackParam paramITS;
482 0 : Double_t itsChi2=0, itsNCl=0;
483 0 : AliExternalTrackParam paramTRD;
484 0 : Double_t trdChi2=0, trdNCl=0;
485 0 : AliExternalTrackParam paramTOF;
486 0 : Double_t tofChi2=0, tofNCl=0;
487 : //
488 : // prepare sorted points
489 0 : AliTrackPointArray *pointarray = (AliTrackPointArray*)friendTrack->GetTrackPointArray();
490 0 : if (!pointarray) continue;
491 0 : Int_t nPointsAll = pointarray->GetNPoints(); // # space points of all available detectors
492 0 : SortPointArray(pointarray, sortedIndex); // space point indices sorted by radius in increasing order
493 :
494 : // 2.) ITS, TRD and ITS-TRD refits
495 : //
496 0 : if (!RefitITStrack(friendTrack,mass,paramITS, itsChi2, itsNCl, sortedIndex)) continue;
497 0 : if (itsNCl<4) continue;
498 : //
499 0 : RefitTRDtrack(friendTrack,mass,paramTRD, trdChi2, trdNCl, sortedIndex);
500 0 : paramTOF=paramITS;
501 0 : RefitTOFtrack(friendTrack,mass,paramTOF, tofChi2, tofNCl, sortedIndex);
502 0 : if (fTrackCounter%fSyswatchStep==0) AliSysInfo::AddStamp("Refitting",fTrackCounter,1,0,0);
503 0 : if ((trdNCl+tofNCl+itsNCl)<5) continue; // - use ITS only tracks also -should it be option?
504 : //
505 : // 3.) Propagate to TPC volume, histogram residuals to TPC clusters and dump all information to TTree
506 : //
507 0 : AliTrackPoint spacepoint;
508 : Int_t volId,modId,layerId;
509 0 : fTrackCounter++; // increase total track number
510 : //
511 : // 3.a) Make a local copy of clusters and apply transformation
512 : //
513 : //
514 0 : for (Int_t iPoint=kMaxLayer;iPoint--;){
515 : //
516 : // reset track interpolation statuses
517 0 : trackArrayITS[iPoint].SetUniqueID(0);
518 0 : trackArrayTRD[iPoint].SetUniqueID(0);
519 0 : trackArrayITSTRD[iPoint].SetUniqueID(0);
520 0 : trackArrayTOF[iPoint].SetUniqueID(0);
521 0 : trackArrayITSTOF[iPoint].SetUniqueID(0);
522 : //
523 0 : const AliTPCclusterMI *cluster=seed->GetClusterPointer(iPoint);
524 0 : if (!cluster){
525 0 : clusterArray[iPoint].SetVolumeId(0);
526 0 : continue;
527 : }
528 0 : clusterArray[iPoint]=*cluster;
529 0 : Int_t i[1]={cluster->GetDetector()};
530 0 : Double_t x[3]={static_cast<Double_t>(cluster->GetRow()),cluster->GetPad(),cluster->GetTimeBin()};
531 0 : transform->Transform(x,i,0,1);
532 0 : clusterArray[iPoint].SetX(x[0]);
533 0 : clusterArray[iPoint].SetY(x[1]);
534 0 : clusterArray[iPoint].SetZ(x[2]);
535 0 : }
536 : //
537 : // 4.) Propagate ITS tracks outward
538 : //
539 0 : Bool_t itsOK=kTRUE;
540 : int npUpdITS = 0;
541 0 : for (Int_t iPoint=0;iPoint<kMaxLayer;iPoint++) {
542 : //trackArrayITS[iPoint].SetUniqueID(0);
543 0 : AliTPCclusterMI &cluster=clusterArray[iPoint];
544 0 : if (cluster.GetVolumeId()==0) continue;
545 0 : double alpSect = ((cluster.GetDetector()%18)+0.5)*20*TMath::DegToRad();
546 0 : double xyz[3] = {cluster.GetX(),cluster.GetY(),cluster.GetZ()}; // sector tracking frame
547 0 : paramITS.Local2GlobalPosition(xyz,alpSect);
548 0 : Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
549 0 : paramITS.Global2LocalPosition(xyz,alpha); // cluster frame
550 0 : if (!(itsOK=paramITS.Rotate(alpha))) break;
551 : // full material correction makes sense only when crossing the boundary of the TPC
552 0 : itsOK = (++npUpdITS)==1 ?
553 0 : AliTrackerBase::PropagateTrackToBxByBz(¶mITS,xyz[0],mass,1,kFALSE) :
554 0 : PropagateInTPCTo(¶mITS,xyz[0],fRhoTPC,fX0TPC,mass) &&
555 0 : TMath::Abs(paramITS.GetSnp())<kMaxSNP;
556 0 : if (itsOK){
557 0 : trackArrayITS[iPoint]=paramITS;
558 0 : trackArrayITS[iPoint].SetUniqueID(1);
559 : }
560 0 : else break; // no sense to propagate farther
561 0 : }
562 0 : if (!itsOK) continue; // no point in continuing if ITS failed
563 0 : if (fTrackCounter%fSyswatchStep==0) AliSysInfo::AddStamp("ExtrapolateITS",fTrackCounter,2,0,0);
564 : //
565 : // 5.) Propagate TRD/TOF tracks inwards
566 : //
567 0 : Bool_t trdOK=(trdNCl>0);
568 0 : Bool_t tofOK=(tofNCl>0);
569 : int npUpdTRD = 0, npUpdTOF = 0;
570 : //
571 0 : for (Int_t iPoint=kMaxLayer;iPoint--;){
572 0 : if (!trdOK && !tofOK) break; // no sense to continue;
573 0 : AliTPCclusterMI &cluster=clusterArray[iPoint];
574 : // if (cluster==NULL) continue;
575 0 : if (cluster.GetVolumeId()==0) continue;
576 0 : double alpSect = ((cluster.GetDetector()%18)+0.5)*20*TMath::DegToRad();
577 0 : double xyz[3] = {cluster.GetX(),cluster.GetY(),cluster.GetZ()}; // sector tracking frame
578 0 : paramTRD.Local2GlobalPosition(xyz,alpSect);
579 0 : Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
580 0 : paramTRD.Global2LocalPosition(xyz,alpha); // cluster frame
581 :
582 0 : if (trdOK){
583 : // material correction makes sense only when crossing the boundary of the TPC
584 0 : trdOK = paramTRD.Rotate(alpha) && ((++npUpdTRD)==1 ?
585 0 : AliTrackerBase::PropagateTrackToBxByBz(¶mTRD,xyz[0],mass,1,kFALSE) :
586 0 : PropagateInTPCTo(¶mTRD,xyz[0],fRhoTPC,fX0TPC,mass))
587 0 : && TMath::Abs(paramTRD.GetSnp())<kMaxSNP;
588 0 : if (trdOK){
589 0 : trackArrayTRD[iPoint]=paramTRD;
590 0 : trackArrayTRD[iPoint].SetUniqueID(1);
591 : //
592 0 : trackArrayITSTRD[iPoint]=paramTRD;
593 0 : if (trackArrayITS[iPoint].GetUniqueID()) { // update by ITS only if the latter is OK
594 0 : AliTrackerBase::UpdateTrack(trackArrayITSTRD[iPoint], trackArrayITS[iPoint]);
595 0 : Double_t chi2=trackArrayITSTRD[iPoint].GetY()-trackArrayITS[iPoint].GetY();
596 0 : chi2*=chi2;
597 0 : chi2/=trackArrayITSTRD[iPoint].GetSigmaY2()+trackArrayITS[iPoint].GetSigmaY2()+kMaxAlignTolerance;
598 0 : if (chi2<kMaxChi2) trackArrayITSTRD[iPoint].SetUniqueID(1);
599 0 : }
600 : }
601 : }
602 0 : if (tofOK){
603 : // material correction makes sense only when crossing the boundary of the TPC
604 0 : tofOK = paramTOF.Rotate(alpha) && ((++npUpdTOF)==1 ?
605 0 : AliTrackerBase::PropagateTrackToBxByBz(¶mTOF,xyz[0],mass,1,kFALSE) :
606 0 : PropagateInTPCTo(¶mTOF,xyz[0],fRhoTPC,fX0TPC,mass))
607 0 : && TMath::Abs(paramTOF.GetSnp())<kMaxSNP;
608 0 : if (tofOK){
609 0 : trackArrayTOF[iPoint]=paramTOF;
610 0 : trackArrayTOF[iPoint].SetUniqueID(1);
611 :
612 0 : trackArrayITSTOF[iPoint]=paramTOF;
613 0 : if (trackArrayITS[iPoint].GetUniqueID()) { // update by ITS only if the latter is OK
614 0 : AliTrackerBase::UpdateTrack(trackArrayITSTOF[iPoint], trackArrayITS[iPoint]);
615 0 : Double_t chi2=trackArrayITSTOF[iPoint].GetY()-trackArrayITS[iPoint].GetY();
616 0 : chi2*=chi2;
617 0 : chi2/=trackArrayITSTOF[iPoint].GetSigmaY2()+trackArrayITS[iPoint].GetSigmaY2()+kMaxAlignTolerance;
618 0 : if (chi2<kMaxChi2) trackArrayITSTOF[iPoint].SetUniqueID(1);
619 0 : }
620 : }
621 : }
622 : //
623 0 : }
624 0 : if (fTrackCounter%fSyswatchStep==0) AliSysInfo::AddStamp("InterpolateTRD",fTrackCounter,3,0,0);
625 :
626 0 : if ( ((fStreamLevelTrack&kStremInterpolation)>0)&&(fTrackCounter%fSyswatchStep==0)){
627 : //if ((fTrackCounter%fSyswatchStep==0)){
628 0 : for (Int_t iPoint=0;iPoint<kMaxLayer;iPoint++){
629 0 : AliTPCclusterMI &cluster=clusterArray[iPoint];
630 : //if (cluster==NULL) continue;
631 0 : if (cluster.GetVolumeId()==0) continue;
632 :
633 0 : (*fStreamer)<<"interpolation"<<
634 0 : "itrack="<<fTrackCounter<< // total track #
635 0 : "cluster.="<<&cluster<< // space points //
636 0 : "itsNCl="<<itsNCl<<
637 0 : "trdNCl="<<trdNCl<<
638 0 : "tofNCl="<<tofNCl<<
639 0 : "itsOK="<<itsOK<<
640 0 : "trdOK="<<trdOK<<
641 0 : "tofOK="<<tofOK<<
642 : //
643 0 : "itsChi2="<<itsChi2<<
644 0 : "trdChi2="<<trdChi2<<
645 0 : "tofChi2="<<tofChi2<<
646 0 : "tofBC="<<tofDiff<<
647 : //
648 0 : "trackITS.="<<&trackArrayITS[iPoint]<< // ITS fit
649 0 : "trackTRD.="<<&trackArrayTRD[iPoint]<< // TRD fit
650 0 : "trackTOF.="<<&trackArrayTOF[iPoint]<< // TOF fit
651 0 : "trackITSTRD.="<<&trackArrayITSTRD[iPoint]<< // ITS-TRD fit
652 0 : "trackITSTOF.="<<&trackArrayITSTOF[iPoint]<< // ITS-TOF fit
653 : "\n";
654 0 : }
655 0 : }
656 0 : UShort_t counter=0;
657 : Double_t rounding=200;
658 : //
659 0 : memset( deltaITS0.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
660 0 : memset( deltaITS1.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
661 0 : memset( deltaTRD0.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
662 0 : memset( deltaTRD1.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
663 0 : memset( deltaTOF0.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
664 0 : memset( deltaTOF1.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
665 : //
666 0 : memset( vecR.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
667 0 : memset( vecPhi.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
668 0 : memset( vecZ.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
669 0 : memset( vecSec.GetMatrixArray(), 0,kMaxLayer*sizeof(Float_t));
670 : //
671 0 : for (Int_t iPoint=0;iPoint<kMaxLayer;iPoint++){
672 : //
673 0 : deltaITS0[counter]=deltaITS1[counter]=deltaTRD0[counter]=deltaTRD1[counter]=deltaTOF0[counter]=deltaTOF1[counter]=-999;
674 0 : vecR[counter] = -999;
675 : //
676 0 : AliTPCclusterMI &cluster=clusterArray[iPoint];
677 0 : if (cluster.GetVolumeId()==0) continue;
678 0 : Double_t zsignSector=((cluster.GetDetector()%36)<18) ? 1.:-1.;
679 : //if (zsignSector*cluster.GetZ()<0.) continue;
680 : //
681 0 : if (trackArrayITS[iPoint].GetUniqueID()>0) { // deltas make sense only if ITS was ok
682 0 : deltaITS0[counter]= TMath::Nint(trackArrayITS[iPoint].GetY()*rounding)/rounding;
683 0 : deltaITS1[counter]= TMath::Nint((trackArrayITS[iPoint].GetZ()-cluster.GetZ())*rounding)/rounding;
684 : //
685 0 : if (trackArrayITSTRD[iPoint].GetUniqueID()>0){
686 0 : deltaTRD0[counter]= TMath::Nint(trackArrayITSTRD[iPoint].GetY()*rounding)/rounding;
687 0 : deltaTRD1[counter]= TMath::Nint((trackArrayITSTRD[iPoint].GetZ()-cluster.GetZ())*rounding)/rounding;
688 0 : }
689 0 : if (trackArrayITSTOF[iPoint].GetUniqueID()>0){
690 0 : deltaTOF0[counter]= TMath::Nint(trackArrayITSTOF[iPoint].GetY()*rounding)/rounding;
691 0 : deltaTOF1[counter]= TMath::Nint((trackArrayITSTOF[iPoint].GetZ()-cluster.GetZ())*rounding)/rounding;
692 0 : }
693 : // vecR(kMaxLayer), vecPhi(kMaxLayer), vecZ(kMaxLayer);
694 0 : vecR[counter]=trackArrayITS[iPoint].GetX();
695 0 : vecPhi[counter]=trackArrayITS[iPoint].GetAlpha();
696 0 : vecZ[counter]=trackArrayITS[iPoint].GetZ();
697 0 : vecSec[counter]=cluster.GetDetector();
698 0 : counter++;
699 0 : }
700 0 : }
701 0 : AliExternalTrackParam * ip = (AliExternalTrackParam *)esdTrack->GetInnerParam();
702 0 : Bool_t saveBit = ip->TestBit(kAlignmentBugFixedBit);
703 0 : ip->SetBit(kAlignmentBugFixedBit,kTRUE); // Flag that the alignment bug is fixed
704 0 : Int_t timeStamp= esdEvent->GetTimeStamp();
705 0 : (*fStreamer)<<"delta"<<
706 0 : "nTracks="<<nTracks<< // number of tracks in event (pileup indicator)
707 0 : "nPrimTracks="<<nPrimTracks<< // number of tracks pointed to primary vertes of selected event
708 0 : "timeStamp="<<timeStamp<< // time stamp
709 0 : "itrack="<<fTrackCounter<< // total track #
710 0 : "gid="<<gid<< // global ID of the event
711 0 : "itsNCl="<<itsNCl<<
712 0 : "trdNCl="<<trdNCl<<
713 0 : "tofNCl="<<tofNCl<<
714 0 : "itsOK="<<itsOK<<
715 0 : "trdOK="<<trdOK<<
716 0 : "tofOK="<<tofOK<<
717 0 : "itsChi2="<<itsChi2<<
718 0 : "trdChi2="<<trdChi2<<
719 0 : "tofChi2="<<tofChi2<<
720 0 : "tofBC="<<tofDiff<<
721 : //
722 0 : "track.="<<ip<< // track parameters at inner wal of TPC
723 0 : "npValid="<<counter<<
724 0 : "vecR.="<<&vecR<<
725 0 : "vecPhi.="<<&vecPhi<<
726 0 : "vecSec.="<<&vecSec<< // sector number
727 0 : "vecZ.="<<&vecZ<<
728 0 : "its0.="<<&deltaITS0<<
729 0 : "its1.="<<&deltaITS1<<
730 0 : "trd0.="<<&deltaTRD0<<
731 0 : "trd1.="<<&deltaTRD1<<
732 0 : "tof0.="<<&deltaTOF0<<
733 0 : "tof1.="<<&deltaTOF1<<
734 : "\n";
735 0 : if (fTrackCounter%fSyswatchStep==0) AliSysInfo::AddStamp("FittTree",fTrackCounter,4,0,0);
736 0 : if (fTrackCounter%fSyswatchStep==0) AliSysInfo::AddStamp("FillHistos",fTrackCounter,5,0,0);
737 : //
738 0 : ip->SetBit(kAlignmentBugFixedBit,saveBit);
739 0 : }
740 0 : transform->GetCurrentRecoParamNonConst()->SetUseComposedCorrection( backupUseComposedCorrection);
741 0 : transform->GetCurrentRecoParamNonConst()->SetUseCorrectionMap(backupUseCorrectionMaps);
742 0 : transform->GetCurrentRecoParamNonConst()->SetAccountDistortions(backupAccountDistortions);
743 :
744 : //
745 : // end of track loop
746 0 : }
747 :
748 : void AliTPCcalibAlignInterpolation::CreateResidualHistosInterpolation(Double_t dy, Double_t dz, Int_t selHis){
749 : //
750 : // Make cluster residual histograms
751 : //
752 0 : Double_t xminTrack[9], xmaxTrack[9];
753 0 : Double_t xminTrackITS[9], xmaxTrackITS[9];
754 0 : Int_t binsTrack[9], binsTrackITS[9];
755 0 : TString axisName[9],axisTitle[9];
756 : //
757 : // 0 - local q/pt
758 : // 1 - global phi in sector number as float
759 : // 2 - local r
760 : // 3 - local kz
761 : // 4 - delta of interest
762 :
763 : //
764 : // gx,gy,gz - will be taken from the TPC
765 : //
766 : //
767 0 : axisName[kQ2PT]="qpt"; axisTitle[kQ2PT]="q/pt (c/GeV)"; // to fill : track.GetSigned1Pt()
768 0 : binsTrack[kQ2PT]=7; xminTrack[kQ2PT]=-3; xmaxTrack[kQ2PT]=3;
769 0 : binsTrackITS[kQ2PT]=7; xminTrackITS[kQ2PT]=-3; xmaxTrackITS[kQ2PT]=3;
770 : //
771 0 : axisName[kSect]="sector"; axisTitle[kSect]="Sector Number"; // to fill: 9*atan2(gy,gx)/pi+ if (sector>0) sector+18
772 0 : binsTrack[kSect]=181; xminTrack[kSect]=-0.05; xmaxTrack[kSect]=18.05;
773 0 : binsTrackITS[kSect]=181; xminTrackITS[kSect]=-0.05; xmaxTrackITS[kSect]=18.05;
774 : //
775 0 : axisName[kLocX]="R"; axisTitle[kLocX]="r (cm)"; // to fill: gr=sqrt(gy**2+gx**2)
776 0 : binsTrack[kLocX]=53; xminTrack[kLocX]=85.; xmaxTrack[kLocX]=245.;
777 0 : binsTrackITS[kLocX]=53; xminTrackITS[kLocX]=85.; xmaxTrackITS[kLocX]=245.;
778 : //
779 0 : axisName[kZ2X]="kZ"; axisTitle[kZ2X]="z/r"; // to fill : gz/gr
780 0 : binsTrack[kZ2X]=20; xminTrack[kZ2X]=-1.0; xmaxTrack[kZ2X]=1.0; // +-1 for ITS+TRD and ITS+TOF
781 0 : binsTrackITS[kZ2X]=20; xminTrackITS[kZ2X]=-1.8; xmaxTrackITS[kZ2X]=1.8; // +-1.8 for the ITS
782 : //
783 0 : axisName[kDelt]="delta"; axisTitle[kDelt]="#Delta (cm)"; // to fill local(clusterY-track.y)
784 0 : binsTrack[kDelt]=100; xminTrack[kDelt]=-dy; xmaxTrack[kDelt]=dy;
785 0 : binsTrackITS[kDelt]=100; xminTrackITS[kDelt]=-dy; xmaxTrackITS[kDelt]=dy;
786 : //
787 0 : binsTrack[kDelt]=TMath::Min(Int_t(20.+2.*dy/0.05),100); // buffer should be smaller than 1 GBy
788 0 : if (selHis==0 ||selHis<0) fHisITSDRPhi = new THnF("deltaRPhiTPCITS","#Delta_{Y} (cm)", kNDim, binsTrackITS,xminTrackITS, xmaxTrackITS);
789 0 : if (selHis==1 ||selHis<0) fHisITSTRDDRPhi = new THnF("deltaRPhiTPCITSTRD","#Delta_{Y} (cm) TPC-(ITS+TRD)", kNDim, binsTrack,xminTrack, xmaxTrack);
790 0 : if (selHis==2 ||selHis<0) fHisITSTOFDRPhi = new THnF("deltaRPhiTPCITSTOF","#Delta_{Y} (cm) TPC-(ITS+TOF)", kNDim, binsTrack,xminTrack, xmaxTrack);
791 : //
792 0 : binsTrack[kDelt]=TMath::Min(Int_t(20.+2.*dz/0.05),100); // buffer should be smaller than 1 GBy
793 0 : xminTrack[kDelt]=-dz; xmaxTrack[kDelt]=dz;
794 0 : xminTrackITS[kDelt]=-dz; xmaxTrackITS[kDelt]=dz;
795 0 : if (selHis==3 ||selHis<0) fHisITSDZ = new THnF("deltaZTPCITS","#Delta_{Z} (cm)", kNDim, binsTrackITS,xminTrackITS, xmaxTrackITS);
796 0 : if (selHis==4 ||selHis<0) fHisITSTRDDZ = new THnF("deltaZTPCITSTRD","#Delta_{Z} (cm) TPC-(ITS+TRD)", kNDim, binsTrack,xminTrack, xmaxTrack);
797 0 : if (selHis==5 ||selHis<0) fHisITSTOFDZ = new THnF("deltaZTPCITSTOF","#Delta_{Z} (cm) TPC-(ITS+TOF)", kNDim, binsTrack,xminTrack, xmaxTrack);
798 : //
799 : //
800 : //
801 0 : THn *hisToFill[6]={GetHisITSDRPhi(), GetHisITSTRDDRPhi(), GetHisITSTOFDRPhi(), GetHisITSDZ(), GetHisITSTRDDZ(), GetHisITSTOFDZ()};
802 0 : for (Int_t ihis=0; ihis<6; ihis++){
803 0 : if (hisToFill[ihis]) for (Int_t ivar2=0;ivar2<kNDim;ivar2++){
804 0 : hisToFill[ihis]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
805 0 : hisToFill[ihis]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
806 0 : }
807 : }
808 :
809 0 : }
810 :
811 :
812 :
813 : void AliTPCcalibAlignInterpolation::CreateDistortionMapsFromFile(const char * inputFile, const char *outputFile){
814 : //
815 : // Create distortion maps from residual histograms
816 : // TPC cluster to ITS, ITS-TRD and ITS-TOF track fits
817 : //
818 0 : TFile *fHistos = TFile::Open(inputFile);
819 :
820 0 : THnF *histoITS = (THnF*) fHistos->Get("deltaRPhiTPCITS");
821 0 : THnF *histoITSTRD= (THnF*) fHistos->Get("deltaRPhiTPCITSTRD");
822 0 : THnF *histoITSTOF = (THnF*) fHistos->Get("deltaRPhiTPCITSTOF");
823 0 : THnF *histoITSZ = (THnF*) fHistos->Get("deltaZTPCITS");
824 0 : THnF *histoITSTRDZ= (THnF*) fHistos->Get("deltaZTPCITSTRD");
825 0 : THnF *histoITSTOFZ = (THnF*) fHistos->Get("deltaZTPCITSTOF");
826 :
827 0 : TTreeSRedirector * pcstream = new TTreeSRedirector(outputFile,"recreate");
828 :
829 0 : TMatrixD projectionInfo(5,5);
830 0 : projectionInfo(0,0)=0; projectionInfo(0,1)=0; projectionInfo(0,2)=0;
831 0 : projectionInfo(1,0)=4; projectionInfo(1,1)=0; projectionInfo(1,2)=1;
832 0 : projectionInfo(2,0)=3; projectionInfo(2,1)=0; projectionInfo(2,2)=1;
833 0 : projectionInfo(3,0)=2; projectionInfo(3,1)=0; projectionInfo(3,2)=1;
834 0 : projectionInfo(4,0)=1; projectionInfo(4,1)=0; projectionInfo(4,2)=1;
835 :
836 0 : TStatToolkit::MakeDistortionMap(4, histoITS, pcstream, projectionInfo);
837 0 : TStatToolkit::MakeDistortionMap(4, histoITSTRD, pcstream, projectionInfo);
838 0 : TStatToolkit::MakeDistortionMap(4, histoITSTOF, pcstream, projectionInfo);
839 0 : TStatToolkit::MakeDistortionMap(4, histoITSZ, pcstream, projectionInfo);
840 0 : TStatToolkit::MakeDistortionMap(4, histoITSTRDZ, pcstream, projectionInfo);
841 0 : TStatToolkit::MakeDistortionMap(4, histoITSTOFZ, pcstream, projectionInfo);
842 0 : delete pcstream;
843 : //
844 0 : }
845 :
846 : void AliTPCcalibAlignInterpolation::FillHistogramsFromChain(const char * residualList, Double_t dy, Double_t dz,
847 : Int_t startTime, Int_t stopTime, Int_t maxStat,
848 : Int_t selHis,const char * residualInfoFile,
849 : Bool_t fixAlignmentBug){
850 : /**
851 : * Trees with point-track residuals to residual histogram
852 : * @param residualList text file with tree list
853 : * Output - ResidualHistograms.root file with hitogram within AliTPCcalibAlignInterpolation object
854 : */
855 : //
856 : //
857 : //
858 0 : ::Info(" AliTPCcalibAlignInterpolation::FillHistogramsFromChain","Start %s\n", residualList);
859 : Int_t cacheSize= 200000000;
860 0 : if (gSystem->Getenv("treeCacheSize")) cacheSize=TString(gSystem->Getenv("treeCacheSize")).Atoi();
861 : Bool_t autoCache = kFALSE;
862 0 : if (gSystem->Getenv("autoCacheSize")) {
863 0 : autoCache = Bool_t(TString(gSystem->Getenv("autoCacheSize")).Atoi());
864 0 : }
865 : Int_t cacheLearnEntries = 1;
866 0 : if (gSystem->Getenv("cacheLearnEntriesProjection")) {
867 0 : cacheLearnEntries = TString(gSystem->Getenv("cacheLearnEntriesProjection")).Atoi();
868 0 : }
869 0 : Printf("************* cacheSize = %d, autoCache = %d, cacheLearnEntries = %d", cacheSize, (Int_t)autoCache, cacheLearnEntries);
870 :
871 : const Int_t kNDim1 = kNDim-1;
872 : const Double_t kernelSigma2I[4]={1./0.25,1./0.25,1./0.25,1./0.25}; // inverse kernel sigma in bin width units
873 : const Double_t kFillGap=0.02 ; // weight for the "non primary distortion info" -
874 : // used to fill the gap without measurement (PHOS hole)
875 : const Double_t kFillGapITS=0.01;
876 : // track and cluster quality cuts - see also AliTPCcalibAlignInterpolation::CalculateDistance
877 : const Int_t kMaxSkippedCluster=10; // 10 cluster
878 : const Float_t kMaxRMSTrackCut=2.0; // maximal RMS (cm) between the tracks
879 : const Float_t kMaxRMSClusterCut=0.3; // maximal RMS (cm) between the cluster and local mean
880 : const Float_t kMaxDeltaClusterCut=0.5; // maximal delta(cm) between the cluster and local mean
881 : //
882 : // gap weight is kFillGap + Exp(-dist): don't calculate exponent if dist is >kMaxExpArg
883 0 : const Double_t kMaxExpArg = -TMath::Log(TMath::Max(kFillGap*0.1, 1e-3));
884 : //
885 0 : Int_t runNumber=TString(gSystem->Getenv("runNumber")).Atoi();
886 0 : float bz = fixAlignmentBug ? InitForAlignmentBugFix(runNumber) : 0;
887 :
888 : // 0.) Load current information file and bookd variables
889 : //
890 : const Int_t nSec=81; // 72 sector +5 sumarry info+ 4 medians +
891 : const Double_t kMaxZSect[2]={2.49725e+02,2.49698e+02};
892 0 : TVectorF meanNcl(nSec); // mean current estator ncl per sector
893 0 : TVectorF meanNclUsed(nSec); // mean current estator ncl per sector
894 0 : Double_t meanTime=0, maxTime=startTime, minTime=stopTime;
895 0 : Int_t currentTrack=0;
896 :
897 0 : TFile *finfo = TFile::Open(residualInfoFile);
898 : TTree *treeInfo=0;
899 0 : if (finfo) {
900 0 : treeInfo=(TTree*)finfo->Get("summaryTime");
901 0 : }else{
902 0 : ::Fatal("AliTPCcalibAlignInterpolation::FillHistogramsFromChain","residualInfoFile %s does not exist",residualInfoFile);
903 : }
904 0 : TGraphErrors * nclArray[nSec]={0};
905 0 : TGraphErrors * nclArrayUsed[nSec]={0};
906 :
907 0 : if (treeInfo) {
908 0 : for (Int_t iSec=0; iSec<nSec; iSec++){
909 0 : nclArray[iSec]=0;
910 0 : nclArrayUsed[iSec]=0;
911 0 : treeInfo->SetBranchAddress(TString::Format("grNcl%d.",iSec).Data(),&nclArray[iSec]);
912 0 : treeInfo->SetBranchAddress(TString::Format("grNclUsed%d.",iSec).Data(),&nclArrayUsed[iSec]);
913 : }
914 0 : treeInfo->GetEntry(0);
915 : }else{
916 0 : ::Fatal("AliTPCcalibAlignInterpolation::FillHistogramsFromChain","residualInfoFile %s does not contain tree summaryTime",residualInfoFile);
917 : }
918 : //
919 : // 0.a) Load drift velocity calibration in case availbel
920 : //
921 0 : TVectorD *vdriftParam=0;
922 0 : TGraphErrors *vdriftGraph=0;
923 0 : TFile *fdrift = TFile::Open("fitDrift.root");
924 0 : AliSysInfo::AddStamp("FillHistogramsFromChain.LoadDriftBegin",1,0);
925 0 : if (fdrift){
926 0 : TTree * tree = (TTree*)fdrift->Get("fitTimeStat");
927 0 : if (tree==NULL){
928 0 : ::Error("LoadDriftCalibration FAILED", "tree fitTimeStat not avaliable in file fitDrift.root");
929 : }else{
930 0 : tree->SetBranchAddress("grTRDReg.",&vdriftGraph);
931 0 : tree->SetBranchAddress("paramRobust.",&vdriftParam);
932 0 : tree->GetEntry(0);
933 0 : if (vdriftGraph==NULL || vdriftGraph->GetN()<=0){
934 0 : ::Info("LoadDriftCalibration FAILED", "ITS/TRD drift calibration not availalble. Trying ITS/TOF");
935 0 : tree->SetBranchAddress("grTOFReg.",&vdriftGraph);
936 0 : tree->GetEntry(0);
937 : }else{
938 0 : ::Info("LoadDriftCalibration", "tree fitTimeStat loaded from the tree");
939 : }
940 : }
941 :
942 0 : }else{
943 0 : ::Error("LoadDriftCalibration FAILED", "fitDrift.root not present");
944 : }
945 0 : AliSysInfo::AddStamp("FillHistogramsFromChain.LoadDriftEND",1,1);
946 : //
947 : // 1.) Fill histograms and mean informations
948 : //
949 : const Int_t knPoints=kMaxRow;
950 0 : AliTPCcalibAlignInterpolation * calibInterpolation = new AliTPCcalibAlignInterpolation("calibInterpolation","calibInterpolation",kFALSE);
951 0 : calibInterpolation->CreateResidualHistosInterpolation(dy,dz,selHis);
952 0 : TString branches[6]={"its0.","trd0.","tof0.", "its1.","trd1.","tof1."};
953 : //
954 0 : TVectorF *vecDeltaOther= 0;
955 0 : TVectorF *vecDeltaOtherITS= 0;
956 0 : TVectorF *vecDelta= 0;
957 0 : TVectorF *vecDeltaITS= 0;
958 0 : TVectorF *vecR=0;
959 0 : TVectorF *vecSec=0;
960 0 : TVectorF *vecPhi=0;
961 0 : TVectorF *vecZ=0;
962 0 : TVectorF *vecLocalDelta = new TVectorF(knPoints);
963 0 : Int_t timeStamp=0;
964 0 : AliExternalTrackParam *param = 0;
965 : //
966 0 : TString esdList0 = gSystem->GetFromPipe(TString::Format("cat %s",residualList).Data());
967 0 : TObjArray *esdArray= esdList0.Tokenize("\n");
968 0 : Int_t nesd = esdArray->GetEntriesFast();
969 : //
970 0 : THn *hisToFill[6]={calibInterpolation->GetHisITSDRPhi(), calibInterpolation->GetHisITSTRDDRPhi(), calibInterpolation->GetHisITSTOFDRPhi(), calibInterpolation->GetHisITSDZ(), calibInterpolation->GetHisITSTRDDZ(), calibInterpolation->GetHisITSTOFDZ()};
971 : TTreeSRedirector * fout = 0;
972 0 : if (selHis<0) {
973 0 : if (startTime<=0) fout=new TTreeSRedirector("ResidualHistograms.root","recreate");
974 0 : if (startTime>0) fout=new TTreeSRedirector(TString::Format("ResidualHistograms_Time%d.root",startTime).Data(),"recreate");
975 : }
976 0 : if (selHis>=0) {
977 0 : if (startTime<=0) fout=new TTreeSRedirector(TString::Format("ResidualHistograms_His%d.root",selHis).Data(),"recreate");
978 0 : if (startTime>0) fout=new TTreeSRedirector(TString::Format("ResidualHistograms_His%d_Time%d.root",selHis,startTime).Data(),"recreate");
979 : }
980 : TH1 * hisTime=0;
981 0 : if (startTime>0) hisTime=new TH1F("hisTrackTime","hisTrackTime",(stopTime-startTime)/20,startTime,stopTime);
982 0 : TStopwatch timerAll;
983 0 : UShort_t npValid=knPoints;
984 0 : Long64_t fillCounter=0;
985 0 : Long64_t clusterCounter=0;
986 0 : AliSysInfo::AddStamp("FillHistogramsFromChain.BEGIN",2,0);
987 0 : for (Int_t ihis=0; ihis<6; ihis++){
988 0 : if (selHis>=0 && ihis!=selHis) continue;
989 0 : Double_t binWidth[4]={0};
990 0 : for (Int_t idim=0; idim<4; idim++) binWidth[idim]=hisToFill[ihis]->GetAxis(idim)->GetBinWidth(1);
991 0 : for (Int_t iesd=0; iesd<nesd; iesd++){
992 0 : TStopwatch timerFile;
993 0 : TString fileNameString(esdArray->At(iesd)->GetName());
994 0 : if (fileNameString.Contains("alien://") && (!gGrid || (gGrid && !gGrid->IsConnected()))) TGrid::Connect("alien://");
995 0 : TFile *esdFile = TFile::Open(fileNameString.Data(),"read");
996 0 : if (!esdFile) continue;
997 0 : TTree *tree = (TTree*)esdFile->Get("delta");
998 0 : if (!autoCache) {
999 0 : tree->SetCacheSize(cacheSize);
1000 0 : tree->SetCacheLearnEntries(cacheLearnEntries);
1001 : }
1002 0 : tree->SetBranchStatus("*",kFALSE);
1003 0 : if (!tree) continue;
1004 0 : ::Info(" AliTPCcalibAlignInterpolation::FillHistogramsFromChain", "Processing file \t %s\n",esdArray->At(iesd)->GetName());
1005 0 : AliSysInfo::AddStamp(esdArray->At(iesd)->GetName(),ihis,iesd,currentTrack);
1006 0 : tree->SetBranchStatus("timeStamp",kTRUE);
1007 0 : TBranch *br = tree->GetBranch("timeStamp");
1008 0 : tree->SetBranchStatus("vecR.",kTRUE);
1009 0 : tree->SetBranchStatus("vecSec.",kTRUE);
1010 0 : tree->SetBranchStatus("vecPhi.",kTRUE);
1011 0 : tree->SetBranchStatus("vecZ.",kTRUE);
1012 0 : tree->SetBranchStatus("track.*",kTRUE);
1013 0 : tree->SetBranchAddress("vecR.",&vecR);
1014 0 : tree->SetBranchAddress("vecSec.",&vecSec);
1015 0 : tree->SetBranchAddress("vecPhi.",&vecPhi);
1016 0 : tree->SetBranchAddress("vecZ.",&vecZ);
1017 0 : tree->SetBranchAddress("track.",¶m);
1018 : //
1019 : // before doing anything else, check if alignment bug indeed must be fixed
1020 0 : Bool_t fixAlignmentBugForFile = fixAlignmentBug;
1021 0 : if (fixAlignmentBugForFile) {
1022 0 : TBranch *brParam = tree->GetBranch("track.");
1023 0 : brParam->GetEntry(0);
1024 0 : if (param->TestBit(kAlignmentBugFixedBit)) {
1025 0 : ::Info(" AliTPCcalibAlignInterpolation::FillHistogramsFromChain",
1026 : "AlignmentBugFix is requested but is not needed for this chunk\n");
1027 : fixAlignmentBugForFile = kFALSE;
1028 0 : }
1029 0 : }
1030 : //
1031 0 : br->SetAddress(&timeStamp);
1032 0 : if (tree->GetBranch("npValid")!=NULL) {
1033 0 : tree->SetBranchStatus("npValid",kTRUE);
1034 0 : tree->SetBranchAddress("npValid",&npValid);
1035 : }
1036 0 : tree->SetBranchStatus(branches[ihis],kTRUE);
1037 0 : tree->SetBranchAddress(branches[ihis],&vecDelta);
1038 : //
1039 : // if aligment bug fix is needed, we need also other delta
1040 0 : if (fixAlignmentBugForFile) {
1041 0 : int ihisOther = ihis<=2 ? ihis+3 : ihis-3;
1042 0 : tree->SetBranchStatus(branches[ihisOther],kTRUE);
1043 0 : tree->SetBranchAddress(branches[ihisOther],&vecDeltaOther);
1044 0 : }
1045 :
1046 0 : if (ihis<=2 &&ihis!=0){
1047 0 : tree->SetBranchStatus(branches[0],kTRUE);
1048 0 : tree->SetBranchAddress(branches[0],&vecDeltaITS);
1049 0 : if (fixAlignmentBugForFile) {
1050 0 : tree->SetBranchStatus(branches[3],kTRUE);
1051 0 : tree->SetBranchAddress(branches[3],&vecDeltaOtherITS);
1052 : }
1053 : }
1054 0 : else if (ihis>2 && ihis!=3){
1055 0 : tree->SetBranchStatus(branches[3],kTRUE);
1056 0 : tree->SetBranchAddress(branches[3],&vecDeltaITS);
1057 0 : if (fixAlignmentBugForFile) {
1058 0 : tree->SetBranchStatus(branches[0],kTRUE);
1059 0 : tree->SetBranchAddress(branches[0],&vecDeltaOtherITS);
1060 : }
1061 : }
1062 :
1063 : // prepare aux info for histo bin calculation
1064 0 : Long64_t nBProd[kNDim] = {0};
1065 0 : int nBinDim[kNDim];
1066 0 : double bsize[kNDim],bsizeI[kNDim],limMin[kNDim],limMax[kNDim];
1067 0 : nBProd[kNDim1] = 1;
1068 0 : THn* curHis = hisToFill[ihis];
1069 0 : TNDArrayT<float>& arrND = (TNDArrayT<float>&)curHis->GetArray(); // for direct access
1070 0 : for (int i=kNDim;i--;) {
1071 0 : TAxis* ax = curHis->GetAxis(i);
1072 0 : limMin[i] = ax->GetXmin();
1073 0 : limMax[i] = ax->GetXmax();
1074 0 : nBinDim[i] = ax->GetNbins();
1075 0 : bsize[i] = (limMax[i]-limMin[i])/nBinDim[i];
1076 0 : bsizeI[i] = 1./bsize[i];
1077 0 : if (i<kNDim1) nBProd[i] = nBProd[i+1]*(nBinDim[i+1]+2); // +2 to account for under/over-flows
1078 : }
1079 : //
1080 0 : Int_t ntracks=tree->GetEntries();
1081 0 : if (!autoCache) {
1082 0 : tree->SetCacheSize(0);
1083 0 : tree->SetCacheSize(cacheSize);
1084 0 : tree->SetCacheLearnEntries(cacheLearnEntries);
1085 0 : tree->AddBranchToCache("*", kTRUE);
1086 : }
1087 : //
1088 0 : for (Int_t itrack=0; itrack<ntracks; itrack++){
1089 0 : if (startTime>0){
1090 0 : br->GetEntry(itrack);
1091 0 : if (timeStamp<startTime || timeStamp>stopTime) continue;
1092 0 : hisTime->Fill(timeStamp);
1093 : }
1094 0 : tree->GetEntry(itrack);
1095 0 : Double_t corrTime = (vdriftGraph!=NULL) ? vdriftGraph->Eval(timeStamp):0;
1096 0 : const Float_t *vSec= vecSec->GetMatrixArray();
1097 0 : const Float_t *vPhi= vecPhi->GetMatrixArray();
1098 0 : const Float_t *vR = vecR->GetMatrixArray();
1099 0 : const Float_t *vZ = vecZ->GetMatrixArray();
1100 0 : const Float_t *vDelta = vecDelta->GetMatrixArray();
1101 0 : const Float_t *vDeltaITS = (vecDeltaITS!=NULL) ? vecDeltaITS->GetMatrixArray():0;
1102 0 : const Float_t *vDeltaOther = vecDeltaOther ? vecDeltaOther->GetMatrixArray():0;
1103 0 : const Float_t *vDeltaOtherITS = vecDeltaOtherITS ? vecDeltaOtherITS->GetMatrixArray():0;
1104 0 : const Float_t *vLocalDelta=vecLocalDelta->GetMatrixArray();
1105 : //
1106 : // calculate distance and aplly track and cluster quality cut
1107 0 : Float_t rmsTrack=3, rmsCluster=1;
1108 0 : Int_t nSkippedCluster=AliTPCcalibAlignInterpolation::CalculateDistance(*vecDelta,*vecDeltaITS, *vecSec, *vecLocalDelta, npValid, rmsTrack, rmsCluster,1.5);
1109 0 : if (nSkippedCluster>kMaxSkippedCluster) continue;
1110 0 : if (rmsTrack>kMaxRMSTrackCut) continue;
1111 0 : if (rmsCluster>kMaxRMSClusterCut) continue;
1112 : //
1113 0 : currentTrack++;
1114 :
1115 0 : if (timeStamp<minTime) minTime=timeStamp;
1116 0 : if (timeStamp>maxTime) maxTime=timeStamp;
1117 0 : meanTime+=timeStamp;
1118 0 : if (treeInfo) for (Int_t iSec=0; iSec<nSec; iSec++){
1119 0 : meanNcl[iSec]+=nclArray[iSec]->Eval(timeStamp);
1120 0 : meanNclUsed[iSec]+=nclArrayUsed[iSec]->Eval(timeStamp);
1121 0 : }
1122 :
1123 0 : if (maxStat>0 &¤tTrack>maxStat) break;
1124 0 : double xxx[kNDim] = {0.};
1125 : //for (Int_t ipoint=0; ipoint<knPoints; ipoint++){
1126 0 : for (Int_t ipoint=0; ipoint<npValid; ipoint++){
1127 0 : if (vR[ipoint]<=0 || vDelta[ipoint]<-990.) continue;
1128 0 : if (TMath::Abs(vDelta[ipoint])<0.000001) continue; // RS Do we need this?
1129 0 : if (TMath::Abs(vLocalDelta[ipoint])> kMaxDeltaClusterCut) continue;
1130 0 : float phiUse = vPhi[ipoint], rUse = vR[ipoint], zUse = vZ[ipoint];
1131 0 : float deltaITSUse = (vDeltaITS) ? vDeltaITS[ipoint]:0;
1132 0 : float deltaRefUse = vDelta[ipoint];
1133 0 : int rocID = TMath::Nint(vSec[ipoint]);
1134 :
1135 0 : if (fixAlignmentBugForFile) { // was it produced w/o bug fix
1136 0 : float dAux = vDeltaOther[ipoint];
1137 0 : if (ihis<3) FixAlignmentBug(rocID, param->GetParameter()[4], bz, phiUse, rUse, zUse, deltaRefUse, dAux);
1138 0 : else FixAlignmentBug(rocID, param->GetParameter()[4], bz, phiUse, rUse, zUse, dAux, deltaRefUse);
1139 : //
1140 0 : if (vDeltaITS) {
1141 0 : dAux = vDeltaOtherITS[ipoint];
1142 0 : float phiAux = vPhi[ipoint], rAux = vR[ipoint], zAux = vZ[ipoint];
1143 0 : if (ihis<3) FixAlignmentBug(rocID, param->GetParameter()[4], bz, phiAux, rAux, zAux, deltaITSUse, dAux);
1144 0 : else FixAlignmentBug(rocID, param->GetParameter()[4], bz, phiAux, rAux, zAux, dAux, deltaITSUse);
1145 0 : }
1146 : //
1147 0 : }
1148 : //
1149 0 : Double_t sector=9.*phiUse/TMath::Pi();
1150 0 : if (sector<0) sector+=18;
1151 0 : Double_t deltaPhi=phiUse-TMath::Pi()*(Int_t(sector)+0.5)/9.;
1152 0 : Double_t localX = TMath::Cos(deltaPhi)*rUse;
1153 :
1154 0 : xxx[kQ2PT] = param->GetParameter()[4];
1155 0 : xxx[kSect] = sector;
1156 0 : xxx[kLocX] = localX;
1157 0 : Double_t side=-1.+2.*((rocID%36)<18);
1158 0 : double maxZ = kMaxZSect[side<0];
1159 0 : xxx[kZ2X] = (zUse*side)<-1 ? side*0.001 : zUse/localX; // do not mix z on A side and C side ?? RS
1160 : // apply drift velocity calibration if available
1161 :
1162 0 : if (ihis>2){ // if z residuals and vdrift calibration existing
1163 0 : Double_t drift = (side>0) ? maxZ-zUse : zUse+maxZ;
1164 0 : Double_t gy = TMath::Sin(phiUse)*localX;
1165 : Double_t pvecFit[3];
1166 : pvecFit[0]= side; // z shift (cm)
1167 0 : pvecFit[1]= drift*gy/maxZ; // global y gradient
1168 : pvecFit[2]= drift; // drift length
1169 0 : Double_t expected = (vdriftParam!=NULL) ? (*vdriftParam)[0]+(*vdriftParam)[1]*pvecFit[0]+(*vdriftParam)[2]*pvecFit[1]+(*vdriftParam)[3]*pvecFit[2]:0;
1170 0 : deltaRefUse= side*(deltaRefUse*side-(expected+corrTime*drift));
1171 0 : deltaITSUse= side*(deltaITSUse*side-(expected+corrTime*drift));
1172 0 : }
1173 0 : xxx[kDelt] = deltaRefUse;
1174 0 : clusterCounter++;
1175 : //
1176 : // calculate axis bins and global bin
1177 0 : Int_t binIndex[kNDim]={0};
1178 : //
1179 0 : if (xxx[kNDim1]<limMin[kNDim1]) binIndex[kNDim1] = 0; // underflow
1180 0 : else if (xxx[kNDim1]<limMax[kNDim1]) binIndex[kNDim1] = 1+int((xxx[kNDim1] - limMin[kNDim1])*bsizeI[kNDim1]); // range
1181 0 : else binIndex[kNDim1] = nBinDim[kNDim1]+1; // oveflow
1182 : //
1183 0 : ULong64_t binToFill = binIndex[kNDim1]; // global bin
1184 0 : for (int id=kNDim1;id--;) {
1185 0 : if (xxx[id]<limMin[id]) binIndex[id] = 0; // underflow
1186 0 : else if (xxx[id]<limMax[id]) binIndex[id] = 1+int((xxx[id] - limMin[id])*bsizeI[id]); // range
1187 0 : else binIndex[id] = nBinDim[id]+1; // oveflow
1188 0 : binToFill += binIndex[id]*nBProd[id];
1189 : }
1190 : //
1191 0 : if (vDeltaITS){
1192 0 : xxx[kDelt] = deltaITSUse;
1193 : int binDeltITS;
1194 0 : if (deltaITSUse<limMin[kDelt]) binDeltITS = 0; // underflow
1195 0 : else if (deltaITSUse<limMax[kDelt]) binDeltITS = 1+int((deltaITSUse - limMin[kDelt])*bsizeI[kDelt]); // range
1196 0 : else binDeltITS = nBinDim[kDelt]+1; // oveflow
1197 0 : Long64_t binToFillITS = binToFill + (binDeltITS-binIndex[kDelt])*nBProd[kDelt]; // global bin for ITS
1198 0 : arrND.At(binToFillITS) += kFillGapITS; // curHis->Fill(xxx,kFillGapITS);
1199 0 : xxx[kDelt] = deltaRefUse;
1200 0 : }
1201 0 : arrND.At(binToFill) += 1.; // curHis->Fill(xxx,1.);
1202 :
1203 0 : Double_t dbinCenter[kNDim];
1204 0 : for (Int_t idim=kNDim;idim--;) {
1205 : // fractional distance to the center of the bin
1206 : // double bincenter = limMin[idim]+bsize[idim]*(binIndex[idim]-0.5); // binindex=0 is underflow!
1207 : // dbinCenter[idim] = (xxx[idim]-bincenter)*bsizeI[idim];
1208 0 : dbinCenter[idim] = (xxx[idim]-limMin[idim])*bsize[idim] - (binIndex[idim]-0.5); // binindex=0 is underflow!
1209 : }
1210 : // dbinCenter[kDelt] = 0;
1211 : //
1212 0 : for (Int_t iqpt=-1; iqpt<=1; iqpt++){ //qpt
1213 0 : int qptBin = binIndex[kQ2PT]+iqpt;
1214 0 : if (qptBin<0||qptBin>nBinDim[kQ2PT]) continue;
1215 0 : double dqpt = dbinCenter[kQ2PT] + iqpt;
1216 0 : dqpt *= dqpt*kernelSigma2I[kQ2PT];
1217 0 : binToFill += iqpt*nBProd[kQ2PT];
1218 : //for (Int_t isec=0; isec<=0; isec++){ //sector - (Not defined yet if we should make bin respone functio and unfold later)
1219 : // int secBin = binIndex[kSect]+isec;
1220 : // if (secBin<0||secBin>nBinDim[kSect]) continue;
1221 : // double dsec = dbinCenter[kSect] + isec;
1222 : // dsec *= dsec*kernelSigma2I[kSect];
1223 : // binToFill += isec*nBProd[kSect];
1224 : // for (Int_t ilocx=0; ilocx<=0; ilocx++){ //local X
1225 : // int locxBin = binIndex[kLocX]+ilocx;
1226 : // if (locxBin<0||locxBin>nBinDim[kLocX]) continue;
1227 : // double dlocx = dbinCenter[kLocX] + ilocx;
1228 : // dlocx *= dlocx*kernelSigma2I[kLocX];
1229 : // binToFill += ilocx*nBProd[kLocX];
1230 0 : for (Int_t iz2x=-2; iz2x<=2; iz2x++){ // Z/x
1231 0 : if ( xxx[kZ2X]*(xxx[kZ2X]+bsize[kZ2X]*iz2x) < 0 ) continue; // do not mix a and C side
1232 0 : int z2xBin = binIndex[kZ2X]+iz2x;
1233 0 : if (z2xBin<0||z2xBin>nBinDim[kZ2X]) continue;
1234 0 : double dz2x = dbinCenter[kZ2X] + iz2x;
1235 0 : dz2x *= dz2x*kernelSigma2I[kZ2X];
1236 0 : binToFill += iz2x*nBProd[kZ2X];
1237 : //
1238 : {
1239 : // Looping is over, fill histo
1240 0 : Double_t weightAll= 2.*(dz2x+dqpt); // 2.*(dz2x+dqpt+dsec+dlocx);
1241 0 : weightAll = weightAll>kMaxExpArg ? kFillGap : kFillGap+TMath::Exp(-weightAll);
1242 0 : arrND.At(binToFill) += weightAll; // curHis->Fill(...)
1243 0 : if (fillCounter==0) printf("Start to Fill\n");
1244 0 : fillCounter++;
1245 : }
1246 : //
1247 0 : binToFill -= iz2x*nBProd[kZ2X];
1248 0 : } // z2x fill loop
1249 : // binToFill -= ilocx*nBProd[kLocX];
1250 : // } // xloc fill loop
1251 : // binToFill -= isec*nBProd[kSect];
1252 : // } // sector fill loop
1253 0 : binToFill -= iqpt*nBProd[kQ2PT]; // restore
1254 0 : } // qpt fill loop
1255 0 : }
1256 0 : }
1257 0 : tree->PrintCacheStats();
1258 :
1259 0 : timerFile.Print();
1260 0 : delete tree;
1261 0 : delete esdFile;
1262 :
1263 0 : }
1264 0 : fout->GetFile()->cd();
1265 0 : hisToFill[ihis]->Write();
1266 0 : }
1267 0 : AliSysInfo::AddStamp("FillHistogramsFromChain.END",2,1);
1268 0 : if (hisTime) hisTime->Write();
1269 0 : ::Info(" AliTPCcalibAlignInterpolation::FillHistogramsFromChain","End of processing\n");
1270 0 : timerAll.Print();
1271 0 : printf("StatInfo.fillCounter:\t%lld\n",fillCounter);
1272 0 : printf("StatInfo.clusterCounter:\t%lld\n",clusterCounter);
1273 0 : printf("StatInfo.trackCounter:\t%d\n",currentTrack);
1274 : //
1275 : // 2.) Fill metadata information
1276 : //
1277 0 : if (currentTrack>0){
1278 0 : meanTime/=currentTrack;
1279 0 : if (treeInfo) for (Int_t iSec=0; iSec<nSec; iSec++){
1280 0 : meanNcl[iSec]/=currentTrack;
1281 0 : meanNclUsed[iSec]/=currentTrack;
1282 0 : }
1283 : }
1284 0 : (*fout)<<"metaData"<<
1285 0 : "runNumber="<<runNumber<< // runNumber
1286 0 : "selHis="<<selHis<< // selected histogram type
1287 0 : "fillCounter="<<fillCounter<< // number of histogram fills
1288 0 : "clusterCounter="<<clusterCounter<< // number of clusters used for fill
1289 0 : "startTime="<<startTime<< // start time as requested
1290 0 : "stopTime="<<stopTime<< // stop time as requested
1291 0 : "meanTime="<<meanTime<< // mean time
1292 0 : "minTime="<<minTime<< // minimal time stamp in data sample
1293 0 : "maxTime="<<maxTime<< // maximal time stamp in data sample
1294 0 : "ntracksUsed="<<currentTrack<< // number of tracks acumulated in time interval
1295 0 : "meanNcl.="<<&meanNcl<< // current estimator - mean number of clusters
1296 0 : "meanNclUsed.="<<&meanNclUsed; // current estimator - mean number of clusters
1297 :
1298 0 : for (Int_t iSec=0; iSec<nSec; iSec++){
1299 0 : (*fout)<<"metaData"<<
1300 0 : TString::Format("grNcl%d.=",iSec).Data()<< nclArray[iSec]<<
1301 0 : TString::Format("grNclUsed%d.=",iSec).Data()<< nclArrayUsed[iSec];
1302 : }
1303 0 : (*fout)<<"metaData"<<"\n";
1304 :
1305 0 : delete fout;
1306 0 : }
1307 :
1308 :
1309 : void AliTPCcalibAlignInterpolation::FillHistogramsFromStreamers(const char * residualList, Double_t dy, Double_t dz, Int_t downscale){
1310 : /**
1311 : * Input list of ErrParam trees as defined in the AliTPCtracker in debug mode
1312 : * @param residualList text file with tree list
1313 : * Output - ResidualHistograms.root file with hitogram within AliTPCcalibAlignInterpolation object
1314 : residualList="residual.list"
1315 : dy=1; dz=1
1316 : */
1317 : //
1318 : //
1319 : //
1320 0 : AliTPCcalibAlignInterpolation * calibInterpolation = new AliTPCcalibAlignInterpolation("calibInterpolation","calibInterpolation",kFALSE);
1321 0 : calibInterpolation->CreateResidualHistosInterpolation(dy,dz);
1322 0 : TString esdList0 = gSystem->GetFromPipe(TString::Format("cat %s",residualList).Data());
1323 0 : TObjArray *esdArray= esdList0.Tokenize("\n");
1324 0 : Int_t nesd = esdArray->GetEntriesFast();
1325 : //
1326 0 : THn *hisToFill[6]={calibInterpolation->GetHisITSDRPhi(), calibInterpolation->GetHisITSTRDDRPhi(), calibInterpolation->GetHisITSTOFDRPhi(), calibInterpolation->GetHisITSDZ(), calibInterpolation->GetHisITSTRDDZ(), calibInterpolation->GetHisITSTOFDZ()};
1327 : //
1328 : //
1329 0 : AliExternalTrackParam * param=0;
1330 0 : AliTPCclusterMI * cl=0;
1331 0 : Int_t iter=0;
1332 : Int_t currentCl=0;
1333 0 : for (Int_t iesd=0; iesd<nesd; iesd++){
1334 0 : TString fileNameString(esdArray->At(iesd)->GetName());
1335 0 : if (fileNameString.Contains("alien://") && (!gGrid || (gGrid && !gGrid->IsConnected()))) TGrid::Connect("alien://");
1336 0 : TFile *esdFile = TFile::Open(fileNameString.Data(),"read");
1337 0 : if (!esdFile) continue;
1338 0 : TTree *tree = (TTree*)esdFile->Get("ErrParam");
1339 0 : if (!tree) continue;
1340 0 : tree->SetBranchAddress("Cl.",&cl);
1341 0 : tree->SetBranchAddress("T.",¶m);
1342 0 : tree->SetBranchAddress("iter",&iter);
1343 0 : Int_t nCl=tree->GetEntries();
1344 0 : for (Int_t iCl=0; iCl<nCl; iCl+=downscale){
1345 0 : tree->GetEntry(iCl);
1346 0 : if (iCl%100000==0) printf("%d\n",iCl);
1347 0 : currentCl++;
1348 0 : double alpSect = ((cl->GetDetector()%18)+0.5)*20*TMath::DegToRad();
1349 0 : double xyz[3] = {cl->GetX(),cl->GetY(),cl->GetZ()}; // sector tracking frame
1350 0 : param->Local2GlobalPosition(xyz,alpSect);
1351 0 : Double_t phi = TMath::ATan2(xyz[1],xyz[0]);
1352 0 : Double_t radius=TMath::Sqrt(xyz[1]*xyz[1]+xyz[0]*xyz[0]);
1353 0 : param->Rotate(phi);
1354 0 : param->PropagateTo(radius,0.); // for big distortion we should query field, for small deltas we are using straight approximtion
1355 0 : Double_t sector=9*phi/TMath::Pi();
1356 0 : if (sector<0) sector+=18;
1357 0 : Double_t deltaY=param->GetY();
1358 0 : Double_t deltaZ=param->GetZ()-cl->GetZ();
1359 0 : Double_t localX = cl->GetX();
1360 0 : Double_t zsignSector=((cl->GetDetector()%36)<18) ? 1.:-1.;
1361 0 : if (zsignSector*cl->GetZ()<0.) continue;
1362 0 : Double_t xxx[5]={ param->GetParameter()[4], sector, localX, cl->GetZ()/cl->GetX(), deltaY};
1363 0 : hisToFill[iter]->Fill(xxx);
1364 0 : xxx[4]=deltaZ;
1365 0 : hisToFill[3+iter]->Fill(xxx);
1366 0 : }
1367 0 : }
1368 0 : TFile * fout = TFile::Open("ResidualHistograms.root","recreate");
1369 0 : calibInterpolation->GetHisITSDRPhi()->Write("deltaYIter0");
1370 0 : calibInterpolation->GetHisITSTRDDRPhi()->Write("deltaYIter1");
1371 0 : calibInterpolation->GetHisITSTOFDRPhi()->Write("deltaYIter2");
1372 0 : calibInterpolation->GetHisITSDZ()->Write("deltaZIter0");
1373 0 : calibInterpolation->GetHisITSTRDDZ()->Write("deltaZIter1");
1374 0 : calibInterpolation->GetHisITSTOFDZ()->Write("deltaZIter2");
1375 0 : delete fout;
1376 0 : }
1377 :
1378 :
1379 :
1380 :
1381 : TTree* AliTPCcalibAlignInterpolation::AddFriendDistortionTree(TTree * tree, const char * fname, const char *treeName, const char *friendAlias){
1382 : //
1383 : //
1384 : //
1385 0 : TFile * fin = TFile::Open(fname);
1386 0 : if (fin==NULL) {
1387 0 : ::Error("AliTPCcalibAlignInterpolation::AddFriendDistortionTree", "file %s not readable", fname);
1388 0 : return 0;
1389 : }
1390 0 : TTree * treeFriend = (TTree*) fin->Get(treeName);
1391 :
1392 0 : if (treeFriend==NULL){
1393 0 : ::Error("AliTPCcalibAlignInterpolation::AddFriendDistortionTree", "file %s not readable", fname);
1394 0 : return 0;
1395 : }
1396 0 : if (tree==NULL) {
1397 : tree = treeFriend;
1398 0 : tree->SetName("Default");
1399 0 : tree->SetTitle("Default");
1400 0 : treeFriend = (TTree*) fin->Get(treeName);
1401 0 : }
1402 : {
1403 0 : tree->AddFriend(treeFriend,TString::Format("%s",friendAlias).Data());
1404 0 : tree->SetAlias(TString::Format("%sOK",friendAlias).Data(),TString::Format("%s.rms>0&&abs(%s.mean-%s.meanG)<2&&%s.chi2G>0&&%s.rmsG<2&&%s.rmsG/%s.rms<2",friendAlias,friendAlias,friendAlias,friendAlias,friendAlias,friendAlias,friendAlias).Data());
1405 0 : tree->SetAlias(TString::Format("%sDrawOK",friendAlias).Data(),TString::Format("%s.rms>0&&abs(%s.mean-%s.meanG)<4&&%s.chi2G>0",friendAlias,friendAlias,friendAlias,friendAlias).Data());
1406 : }
1407 0 : return tree;
1408 0 : }
1409 :
1410 : //_____________________________________________________________________________
1411 : Bool_t AliTPCcalibAlignInterpolation::PropagateInTPCTo(AliExternalTrackParam* t, Double_t xk, Double_t rho,Double_t x0, Double_t mass)
1412 : {
1413 : //-----------------------------------------------------------------
1414 : // This function propagates a track to a reference plane x=xk.
1415 : // rho - density of the crossed matrial (g/cm^3)
1416 : // x0 - radiation length of the crossed material (g/cm^2)
1417 : //-----------------------------------------------------------------
1418 : //
1419 0 : Double_t old[3]={t->GetX(),t->GetY(),t->GetZ()};
1420 0 : Double_t b[3]; AliTrackerBase::GetBxByBz(old,b);
1421 0 : if (!t->PropagateToBxByBz(xk,b)) return kFALSE;
1422 :
1423 0 : Double_t d = TMath::Sqrt((t->GetX()-old[0])*(t->GetX()-old[0]) +
1424 0 : (t->GetY()-old[1])*(t->GetY()-old[1]) +
1425 0 : (t->GetZ()-old[2])*(t->GetZ()-old[2]));
1426 0 : if (old[0] < xk) d = -d;
1427 0 : if (!t->CorrectForMeanMaterial(d*rho/x0,d*rho,mass,
1428 0 : kFALSE,AliExternalTrackParam::BetheBlochGas)) return kFALSE;
1429 :
1430 0 : return kTRUE;
1431 0 : }
1432 :
1433 : //_____________________________________________________________________________
1434 : void AliTPCcalibAlignInterpolation::ExtractTPCGasData()
1435 : {
1436 : // get TPC gas rho and X0
1437 0 : double p0[3] = {90,1,45};
1438 0 : double p1[3] = {240,1,120};
1439 0 : double par[10];
1440 0 : AliTrackerBase::MeanMaterialBudget(p0,p1,par);
1441 0 : fRhoTPC = par[0]>0 ? par[0] : 0.9e-3;
1442 0 : double l = par[4];
1443 0 : fX0TPC = par[1]>0 ? par[4]/par[1] : 28.94;
1444 : //
1445 0 : AliInfoF("Propagation in TPC will use rho=%.2f X0=%.2f",fRhoTPC,fX0TPC);
1446 0 : }
1447 :
1448 :
1449 : void AliTPCcalibAlignInterpolation::MakeEventStatInfo(const char * inputList, Int_t timeInterval, Int_t id, Int_t skip){
1450 : //
1451 : /// Code to query statistical event information from the ResidualTrees.root file
1452 : /// output written to file residualInfo.root
1453 : /// \param const char * inputList - ascii file with input list
1454 : /// \param Int_t timeInterval - length of time interval (beginning of time intervals rounded)
1455 : /// \param id - additional ID added to the tree
1456 : /// \param skip - parameter skip file
1457 : /// Algorithm:
1458 : /// 1.) Cache information per files - beginTime and endTime for file
1459 : /// 2.) Cache information per time interval
1460 :
1461 : /*
1462 : run=240204;
1463 : GetResidualStatInfo("cat residual.list",300,run,1);
1464 : */
1465 0 : TObjArray *array = TString(gSystem->GetFromPipe(TString::Format("%s",inputList).Data())).Tokenize("\n");
1466 0 : Int_t nFiles=array->GetEntries();
1467 0 : if (nFiles<=0) {
1468 0 : ::Error("GetResidualStatInfo", "Wrong input list %s", inputList);
1469 0 : return;
1470 : }
1471 0 : TStopwatch timer;
1472 : //
1473 : // 1.) Cache information per files - beginTime and endTime for file
1474 : //
1475 0 : TStopwatch timer1;
1476 0 : TTreeSRedirector * pcstream = new TTreeSRedirector("residualInfo.root", "recreate");
1477 : Int_t cacheSize=100000000; // 100 MBy cache
1478 0 : if (gSystem->Getenv("treeCacheSize")) cacheSize=TString(gSystem->Getenv("treeCacheSize")).Atoi();
1479 : Bool_t autoCache = kFALSE;
1480 0 : if (gSystem->Getenv("autoCacheSize")) autoCache = Bool_t(TString(gSystem->Getenv("autoCacheSize")).Atoi());
1481 : Int_t cacheLearnEntries = 1;
1482 0 : if (gSystem->Getenv("cacheLearnEntriesStatInfo")) {
1483 0 : cacheLearnEntries = TString(gSystem->Getenv("cacheLearnEntriesStatInfo")).Atoi();
1484 0 : }
1485 0 : Printf("************* cacheSize = %d, autoCache = %d, cacheLearnEntries = %d", cacheSize, (Int_t)autoCache, cacheLearnEntries);
1486 :
1487 0 : TChain * chainInfo = AliXRDPROOFtoolkit::MakeChain("residual.list","eventInfo",0,-1);
1488 0 : if (!autoCache) {
1489 0 : chainInfo->SetCacheSize(cacheSize);
1490 0 : chainInfo->SetCacheLearnEntries(cacheLearnEntries);
1491 : }
1492 0 : TChain * chainTracks = AliXRDPROOFtoolkit::MakeChain("residual.list","delta",0,-1);
1493 0 : if (!autoCache) {
1494 0 : chainTracks->SetCacheSize(cacheSize);
1495 0 : chainTracks->SetCacheLearnEntries(cacheLearnEntries);
1496 : }
1497 : //
1498 : Int_t gidRounding=128; // git has to be rounded
1499 0 : Int_t neventsAll=chainInfo->GetEntries(); // total amount of events
1500 0 : Int_t ntracksAll=chainTracks->GetEntries(); // total amount of tracks
1501 :
1502 : //chainInfo->SetEstimate(-1); // using -1 make a problem - too much memory allocate with autimatic switch . crash in case of 4MBy limits in the next Draw
1503 0 : chainInfo->SetEstimate(neventsAll);
1504 0 : chainInfo->Draw("timeStamp:gid/128","timeStamp>0","goff");
1505 : //
1506 : //
1507 0 : Long64_t minTime=0,maxTime=0;
1508 0 : double minGID=0,maxGID=0,meanGID=0,meanTime=0;
1509 0 : if (neventsAll) {
1510 0 : double minTimeD=0,maxTimeD=0;
1511 0 : TStatToolkit::GetMinMaxMean(chainInfo->GetV1(),neventsAll,minTimeD,maxTimeD, meanTime);
1512 0 : minTime = minTimeD;
1513 0 : maxTime = maxTimeD;
1514 0 : TStatToolkit::GetMinMaxMean(chainInfo->GetV2(),neventsAll,minGID,maxGID, meanGID);
1515 0 : minGID*=128;
1516 0 : maxGID*=128;
1517 0 : meanGID*=128;
1518 0 : }
1519 0 : (*pcstream)<<"summary1"<<
1520 0 : "id="<<id<< // chain id - usually should be run number
1521 0 : "nevents="<<neventsAll<< // total number of events
1522 0 : "ntracks="<<ntracksAll<< // total number of tracks
1523 0 : "minTime="<<minTime<< // minimal time stamp in sample
1524 0 : "maxTime="<<maxTime<< // max time stamp in sample
1525 0 : "meanTime="<<meanTime<< // mean time
1526 0 : "minGID="<<minGID<< // minimal event gid in sample (rounded to 128)
1527 0 : "maxGID="<<maxGID<< // max event gid in sample (rounded to 128)
1528 0 : "meanGID="<<meanGID<< // mean event gid in sample (rounded to 128)
1529 : "\n";
1530 0 : delete pcstream;
1531 0 : ::Info("GetResidualStatInfo","Total time");
1532 0 : timer1.Print();
1533 : //
1534 : // 2.) Cache information per time interval
1535 : //
1536 0 : TStopwatch timer2;
1537 0 : pcstream = new TTreeSRedirector("residualInfo.root", "update");
1538 : // Int_t entries = neventsAll;
1539 :
1540 0 : Long64_t minTimeQA = timeInterval*(minTime/timeInterval);
1541 0 : Long64_t maxTimeQA = timeInterval*(1+(maxTime/timeInterval));
1542 0 : Int_t nIntervals=(maxTimeQA-minTimeQA)/timeInterval;
1543 0 : Int_t nIntervalsQA=(maxTimeQA-minTimeQA)/15;
1544 : //
1545 0 : TH1F * hisEvent= new TH1F("hisEvent","hisEvent",nIntervalsQA,minTimeQA,maxTimeQA);
1546 : const Int_t nSec=81; // 72 sector +5 sumarry info+ 4 medians
1547 0 : TProfile * profArrayNcl[nSec]={0};
1548 0 : TProfile * profArrayNclUsed[nSec]={0};
1549 0 : TGraphErrors * grArrayNcl[nSec]={0};
1550 0 : TGraphErrors * grArrayNclUsed[nSec]={0};
1551 0 : TProfile * profArrayITSNcl[3]={0};
1552 0 : TGraphErrors * grArrayITSNcl[3]={0};
1553 :
1554 0 : for (Int_t isec=0; isec<nSec; isec++){
1555 0 : profArrayNcl[isec]=new TProfile(TString::Format("TPCnclSec%d",isec).Data(), TString::Format("TPCnclSec%d",isec).Data(), nIntervalsQA,minTimeQA,maxTimeQA);
1556 0 : profArrayNclUsed[isec]=new TProfile(TString::Format("TPCnclUsedSec%d",isec).Data(), TString::Format("TPCnclUsedSec%d",isec).Data(), nIntervalsQA,minTimeQA,maxTimeQA);
1557 : }
1558 0 : for (Int_t iits=0; iits<3; iits++){
1559 0 : profArrayITSNcl[iits]=new TProfile(TString::Format("ITSnclSec%d",iits).Data(), TString::Format("ITSnclSec%d",iits).Data(), nIntervalsQA,minTimeQA,maxTimeQA);
1560 : }
1561 :
1562 0 : TVectorF *vecNClTPC=0;
1563 0 : TVectorF *vecNClTPCused=0;
1564 0 : Int_t nITS[3]={0};
1565 0 : Int_t timeStamp=0;
1566 0 : for (Int_t iFile=0; iFile<nFiles; iFile+=skip){
1567 0 : timer.Start();
1568 0 : TString fileName = array->At(iFile)->GetName();
1569 0 : if (fileName.Contains("alien://") && (!gGrid || (gGrid && !gGrid->IsConnected()))) TGrid::Connect("alien://");
1570 0 : printf("%d\t%s\n",iFile,fileName.Data());
1571 0 : AliSysInfo::AddStamp(fileName.Data(),1,iFile);
1572 0 : TFile * f = TFile::Open(fileName.Data());
1573 0 : if (f==NULL) continue;
1574 0 : TTree * treeInfo = (TTree*)f->Get("eventInfo");
1575 0 : if (treeInfo==NULL) continue;
1576 0 : if (!autoCache) {
1577 0 : treeInfo->SetCacheSize(cacheSize);
1578 0 : treeInfo->SetCacheLearnEntries(cacheLearnEntries);
1579 : }
1580 0 : treeInfo->SetBranchAddress("vecNClTPC.",&vecNClTPC);
1581 0 : treeInfo->SetBranchAddress("vecNClTPCused.",&vecNClTPCused);
1582 0 : treeInfo->SetBranchAddress("nSPD",&nITS[0]);
1583 0 : treeInfo->SetBranchAddress("nSDD",&nITS[1]);
1584 0 : treeInfo->SetBranchAddress("nSSD",&nITS[2]);
1585 0 : treeInfo->AddBranchToCache(treeInfo->GetBranch("vecNClTPC."), kTRUE);
1586 0 : treeInfo->AddBranchToCache(treeInfo->GetBranch("vecNClTPCused."), kTRUE);
1587 0 : treeInfo->AddBranchToCache(treeInfo->GetBranch("nSPD"), kTRUE);
1588 0 : treeInfo->AddBranchToCache(treeInfo->GetBranch("nSDD"), kTRUE);
1589 0 : treeInfo->AddBranchToCache(treeInfo->GetBranch("nSSD"), kTRUE);
1590 0 : Bool_t hasTimeStamp=(treeInfo->GetBranch("timeStamp")!=NULL);
1591 0 : if (hasTimeStamp) treeInfo->SetBranchAddress("timeStamp",&timeStamp);
1592 0 : if (!hasTimeStamp) ((TBranch*)(treeInfo->GetListOfBranches()->At(1)))->SetAddress(&timeStamp);
1593 0 : treeInfo->AddBranchToCache(treeInfo->GetBranch("timeStamp"), kTRUE);
1594 0 : Int_t treeEntries=treeInfo->GetEntries();
1595 0 : for (Int_t iEntry=0; iEntry<treeEntries; iEntry++){
1596 0 : treeInfo->GetEntry(iEntry);
1597 0 : hisEvent->Fill(timeStamp);
1598 0 : for (Int_t isec=0; isec<72; isec++){
1599 0 : profArrayNcl[isec]->Fill(timeStamp, (*vecNClTPC)[isec]);
1600 0 : profArrayNclUsed[isec]->Fill(timeStamp, (*vecNClTPC)[isec]);
1601 0 : if (isec<36){
1602 0 : if (isec<18) profArrayNcl[72]->Fill(timeStamp, (*vecNClTPC)[isec]);
1603 0 : if (isec>=18) profArrayNcl[73]->Fill(timeStamp, (*vecNClTPC)[isec]);
1604 0 : if (isec<18) profArrayNclUsed[72]->Fill(timeStamp, (*vecNClTPCused)[isec]);
1605 0 : if (isec>=18) profArrayNclUsed[73]->Fill(timeStamp, (*vecNClTPCused)[isec]);
1606 : }else{
1607 0 : if ((isec%36)<18) profArrayNcl[74]->Fill(timeStamp, (*vecNClTPC)[isec]);
1608 0 : if ((isec%36)>=18) profArrayNcl[75]->Fill(timeStamp, (*vecNClTPC)[isec]);
1609 0 : if ((isec%36)<18) profArrayNclUsed[74]->Fill(timeStamp, (*vecNClTPCused)[isec]);
1610 0 : if ((isec%36)>=18) profArrayNclUsed[75]->Fill(timeStamp, (*vecNClTPCused)[isec]);
1611 : }
1612 0 : profArrayNcl[76]->Fill(timeStamp, (*vecNClTPC)[isec]);
1613 0 : profArrayNclUsed[76]->Fill(timeStamp, (*vecNClTPCused)[isec]);
1614 : }
1615 0 : profArrayNcl[77]->Fill(timeStamp, TMath::Median(18, &((vecNClTPC->GetMatrixArray())[0])));
1616 0 : profArrayNcl[78]->Fill(timeStamp, TMath::Median(18, &((vecNClTPC->GetMatrixArray())[18])));
1617 0 : profArrayNcl[79]->Fill(timeStamp, TMath::Median(18, &((vecNClTPC->GetMatrixArray())[36])));
1618 0 : profArrayNcl[80]->Fill(timeStamp, TMath::Median(18, &((vecNClTPC->GetMatrixArray())[54])));
1619 : //
1620 0 : profArrayNclUsed[77]->Fill(timeStamp, TMath::Median(18, &((vecNClTPCused->GetMatrixArray())[0])));
1621 0 : profArrayNclUsed[78]->Fill(timeStamp, TMath::Median(18, &((vecNClTPCused->GetMatrixArray())[18])));
1622 0 : profArrayNclUsed[79]->Fill(timeStamp, TMath::Median(18, &((vecNClTPCused->GetMatrixArray())[36])));
1623 0 : profArrayNclUsed[80]->Fill(timeStamp, TMath::Median(18, &((vecNClTPCused->GetMatrixArray())[54])));
1624 0 : for (Int_t iits=0; iits<3; iits++){
1625 0 : profArrayITSNcl[iits]->Fill(timeStamp,nITS[iits]);
1626 : }
1627 : }
1628 0 : timer.Print();
1629 0 : treeInfo->PrintCacheStats();
1630 0 : AliSysInfo::AddStamp(fileName.Data(),2,iFile);
1631 0 : delete treeInfo;
1632 0 : delete f;
1633 0 : AliSysInfo::AddStamp(fileName.Data(),3,iFile);
1634 :
1635 0 : }
1636 0 : timer2.Print();
1637 0 : TGraphErrors grEvent(hisEvent);
1638 0 : (*pcstream)<<"summaryTime"<<
1639 0 : "id="<<id<<
1640 0 : "grEvent.="<<&grEvent;
1641 0 : for (Int_t isec=0; isec<nSec; isec++){
1642 0 : grArrayNcl[isec] = new TGraphErrors((profArrayNcl[isec]));
1643 0 : grArrayNclUsed[isec] = new TGraphErrors((profArrayNclUsed[isec]));
1644 0 : (*pcstream)<<"summaryTime"<<
1645 0 : TString::Format("grNcl%d.=",isec).Data()<<grArrayNcl[isec]<<
1646 0 : TString::Format("grNclUsed%d.=",isec).Data()<<grArrayNclUsed[isec];
1647 : }
1648 0 : for (Int_t iits=0; iits<3; iits++){
1649 0 : grArrayITSNcl[iits] = new TGraphErrors((profArrayITSNcl[iits]));
1650 0 : (*pcstream)<<"summaryTime"<<
1651 0 : TString::Format("grITSNcl%d.=",iits).Data()<<grArrayITSNcl[iits];
1652 : }
1653 :
1654 :
1655 0 : (*pcstream)<<"summaryTime"<<"\n";
1656 0 : for (Int_t isec=0; isec<nSec; isec++){
1657 0 : delete profArrayNcl[isec];
1658 0 : delete profArrayNclUsed[isec];
1659 0 : delete grArrayNcl[isec];
1660 0 : delete grArrayNclUsed[isec];
1661 : }
1662 0 : delete hisEvent;
1663 0 : delete pcstream;
1664 :
1665 0 : printf("StatInfo.minTime\t%lld\n",minTime); //this formatting does not work on my Debian why it was changed ?
1666 0 : printf("StatInfo.maxTime\t%lld\n",maxTime);
1667 : //printf("StatInfo.minTime\t%f\n",Double_t(minTime)); //this formatting does not work on my Debian why it was changed ?
1668 : //printf("StatInfo.maxTime\t%f\n",Double_t(maxTime));
1669 :
1670 :
1671 0 : delete array;
1672 0 : }
1673 :
1674 :
1675 :
1676 : Bool_t AliTPCcalibAlignInterpolation::FitDrift(double deltaT, double sigmaT, double time0, double_t time1,Bool_t fixAlignmentBug, Bool_t tofBCValidation){
1677 : //
1678 : // Fit time dependence of the drift velocity for ITS-TRD and ITS-TOF scenario
1679 : /* Intput:
1680 : "residual.list" - ascii file with files assumed to be in working directory
1681 : double deltaT - time binning for drift velocity
1682 : double sigmaT - kernel width for time smoothing
1683 : double time0 - starting time for drift velocity calculation
1684 : double time1 - stop time for drift velocty calculation
1685 : * in case time0 and time1 not specified - full time range in the selected sample used (time0=minTime, time1=maxTime)
1686 : Output:
1687 : "fitDrift.root" - small file with the drift velocity calibration created
1688 : robustFit tree - drift vlocity, time0, z shift and gy gradiend calibration using TLinearFitter::EvalRobust estimator
1689 : - 20 statistically independent values to QA procedure
1690 : - resulting values choosen using robust median estimator
1691 : fitTime - set of graphs drift velocity as function of time
1692 : - n different graph values
1693 : - resulting time dependent graph used as an logal regression with kernel sigmaT per interval
1694 : fitTimeStat - comparing median and local regression statistic
1695 : */
1696 : /*
1697 : To do:
1698 : - 1.) make version:
1699 : - 1.a) with- without bunch 0 crossing
1700 : - 1.b) disentagle fit for the TRD and for the TOF
1701 : - 2.) use time bin sigma estimate for the outlier rejection in the AliNDLocalFit fit
1702 :
1703 : */
1704 :
1705 : /*
1706 : fileList="residual.list"
1707 : time0=0; time1=0;
1708 : deltaT=60; sigmaT=600
1709 : fitDrift(60,00,0);
1710 : */
1711 : const Double_t kMinEntries=1000;
1712 : const Double_t kInvalidR = 70;
1713 : const Double_t kInvalidRes = -900;
1714 : const Double_t kMaxDist0=20;
1715 : const Double_t kMaxDist1=5;
1716 : const Double_t kDumpSample=0.01;
1717 : const Double_t kBCcutMin=-5;
1718 : const Double_t kBCcutMax=20;
1719 : const Double_t robFraction=0.95;
1720 : //const Double_t regRobust=0.0000000001;
1721 :
1722 : Int_t maxEntries=1000000;
1723 : Int_t maxPointsRobust=4000000;
1724 : //
1725 : const Double_t kMaxZSect[2]={2.49725e+02,2.49698e+02};
1726 0 : TCut selection="";
1727 : Int_t entriesAll=0;
1728 0 : Int_t runNumber=TString(gSystem->Getenv("runNumber")).Atoi();
1729 :
1730 : Int_t cacheSize=100000000; // 100 MBy cache
1731 0 : if (gSystem->Getenv("treeCacheSize")) cacheSize=TString(gSystem->Getenv("treeCacheSize")).Atoi();
1732 : Bool_t autoCache = kFALSE;
1733 0 : if (gSystem->Getenv("autoCacheSize")) autoCache = Bool_t(TString(gSystem->Getenv("autoCacheSize")).Atoi());
1734 : Int_t cacheLearnEntries = 1;
1735 0 : if (gSystem->Getenv("cacheLearnEntriesVDrift")) {
1736 0 : cacheLearnEntries = TString(gSystem->Getenv("cacheLearnEntriesVDrift")).Atoi();
1737 0 : }
1738 0 : Printf("************* cacheSize = %d, autoCache = %d, cacheLearnEntries = %d", cacheSize, (Int_t)autoCache, cacheLearnEntries);
1739 :
1740 : //
1741 : //
1742 0 : if (deltaT<=0 || sigmaT<=0){
1743 0 : ::Error("AliTPCcalibAlignInterpolation::FitDrift FAILED ","Invalid parameter value for the deltaT %.1f and sigmaT %.1f", deltaT, sigmaT);
1744 0 : return kFALSE;
1745 : }
1746 0 : if (TString(gSystem->GetFromPipe("cat residual.list | grep -c alien://")).Atoi()>0) TGrid::Connect("alien");
1747 :
1748 0 : TChain * chainDelta = AliXRDPROOFtoolkit::MakeChain("residual.list","delta",0,-1);
1749 : //
1750 : // before doing anything else, check if alignment bug indeed must be fixed
1751 0 : Bool_t fixAlignmentBugForFile = fixAlignmentBug;
1752 0 : AliExternalTrackParam* param = 0;
1753 0 : if (fixAlignmentBugForFile) {
1754 0 : TBranch *brParam = chainDelta->GetBranch("track.");
1755 0 : brParam->SetAddress(¶m);
1756 0 : brParam->GetEntry(0);
1757 0 : if (param->TestBit(kAlignmentBugFixedBit)) {
1758 0 : ::Info(" AliTPCcalibAlignInterpolation::FitDrift",
1759 : "AlignmentBugFix is requested but is not needed for this chunk\n");
1760 : fixAlignmentBugForFile = kFALSE;
1761 0 : }
1762 0 : }
1763 : //
1764 0 : if (!autoCache) {
1765 0 : chainDelta->SetCacheSize(0);
1766 0 : chainDelta->SetCacheSize(cacheSize);
1767 0 : chainDelta->SetCacheLearnEntries(cacheLearnEntries);
1768 : }
1769 0 : AliSysInfo::AddStamp("FitDrift.chainDeltaGetEntriesBegin",1,0,0);
1770 0 : entriesAll = chainDelta->GetEntries();
1771 0 : AliSysInfo::AddStamp("FitDrift.chainDeltaGetEntriesEnd",1,1,0);
1772 :
1773 0 : if (entriesAll<kMinEntries) {
1774 0 : ::Error("fitDrift FAILED","Not enough tracks in the chain. Ntracks=%d",entriesAll);
1775 0 : return kFALSE;
1776 : }
1777 0 : maxEntries=TMath::Min(maxEntries, entriesAll);
1778 0 : TTreeSRedirector *pcstream = new TTreeSRedirector("fitDrift.root","recreate");
1779 0 : if (time0==time1){
1780 0 : AliSysInfo::AddStamp("FitDrift.chainInfoGetTimeBegin",2,0,0);
1781 0 : TChain * chainInfo= AliXRDPROOFtoolkit::MakeChain("residual.list","eventInfo",0,-1);
1782 0 : if (!autoCache) {
1783 0 : chainInfo->SetCacheSize(cacheSize);
1784 0 : chainInfo->SetCacheLearnEntries(cacheLearnEntries);
1785 : }
1786 : // chainInfo->SetEstimate(-1); // i cashed here -1 does not work
1787 0 : chainInfo->SetEstimate(maxEntries); // i cashed here -1 does not work
1788 0 : Int_t entries = chainInfo->Draw("timeStamp","","goff",maxEntries);
1789 0 : chainInfo->PrintCacheStats();
1790 0 : AliSysInfo::AddStamp("FitDrift.chainInfoGetTimeEND",2,1,0);
1791 0 : if (entries) TStatToolkit::GetMinMax(chainInfo->GetV1(),entries,time0,time1);
1792 0 : }
1793 : // 0.) Cache variables: to be done using loop
1794 : // Variables to cache:
1795 : // tof1.fElements
1796 : // trd1.fElements
1797 : // vecZ.fElements
1798 : // vecR.fElements
1799 : // vecSec.fElements
1800 : // vecPhi.fElements
1801 : // timeStamp
1802 : // npValid
1803 : //
1804 0 : AliSysInfo::AddStamp("FitDrift.StartCache",3,0,0);
1805 :
1806 0 : float bz = fixAlignmentBugForFile ? InitForAlignmentBugFix(runNumber) : 0;
1807 :
1808 : // check if TOF was present
1809 : const float tofBCMin = -25.f;
1810 : const float tofBCMax = 50.f;
1811 : //
1812 0 : if (tofBCValidation) {
1813 0 : AliCDBManager* man = AliCDBManager::Instance();
1814 0 : if (!man->IsDefaultStorageSet()) man->SetDefaultStorage("raw://");
1815 0 : if (man->GetRun()<0) man->SetRun(runNumber);
1816 0 : AliGRPObject* grp = (AliGRPObject*)man->Get(AliCDBPath("GRP/GRP/Data"))->GetObject();
1817 0 : Int_t activeDetectors = grp->GetDetectorMask();
1818 0 : if (!(activeDetectors&AliDAQ::DetectorPattern("TOF"))) {
1819 0 : AliWarningClass("Disabling TOF BC validation since TOF is not in the GRP");
1820 : tofBCValidation = kFALSE;
1821 0 : }
1822 : else {
1823 0 : AliInfoClassF("TOF BC validation is enabled with window %.2f : %.2f ns",tofBCMin,tofBCMax);
1824 : }
1825 0 : }
1826 : //
1827 0 : chainDelta->SetEstimate(maxEntries*160/5.);
1828 0 : TString toRead = "tof1.fElements:trd1.fElements:vecZ.fElements:vecR.fElements:vecSec.fElements:vecPhi.fElements:timeStamp:tofBC:its1.fElements";
1829 0 : TString cutEv = "Entry$%5==0 && vecR.fElements>10";
1830 0 : if (tofBCValidation) cutEv += Form(" && tofOK && !(tofBC<%.3f || tofBC>%.3f)",tofBCMin,tofBCMax);
1831 : //
1832 0 : if (fixAlignmentBugForFile) toRead += ":tof0.fElements:trd0.fElements:track.fP[4]"; // needed only in case we need to fix alignment bug
1833 : //
1834 0 : Int_t entriesFit0 = chainDelta->Draw(toRead.Data(),cutEv.Data(),"goffpara",maxEntries);
1835 0 : chainDelta->PrintCacheStats();
1836 0 : AliInfoClassF("Selected %d points from first %d entries",entriesFit0,maxEntries);
1837 :
1838 : float lossFactorOK = 0.5f; // allowed loss factor
1839 0 : float lossFactor = entriesFit0*5./160./maxEntries;
1840 : int prescaling = 10;
1841 0 : if (lossFactor < lossFactorOK) {
1842 0 : prescaling = TMath::Max(lossFactor/lossFactorOK,2.0f);
1843 0 : }
1844 0 : Int_t entriesFit=entriesFit0/prescaling;
1845 0 : AliInfoClassF("Selection loss factor: %f -> random prescaling: %d -> entriesFit: %d",lossFactor,prescaling,entriesFit);
1846 :
1847 0 : AliSysInfo::AddStamp("FitDrift.EndCache",3,1,0);
1848 :
1849 0 : AliSysInfo::AddStamp("FitDrift.BeginFill",4,1,0);
1850 :
1851 0 : TVectorD * deltaZTOF = new TVectorD(entriesFit); //
1852 0 : TVectorD * deltaZTRD = new TVectorD(entriesFit); //
1853 0 : TVectorD * deltaYTOF = new TVectorD(entriesFit); //
1854 0 : TVectorD * deltaYTRD = new TVectorD(entriesFit); //
1855 0 : TVectorD * vecZ = new TVectorD(entriesFit); //
1856 0 : TVectorD * vecR = new TVectorD(entriesFit); //
1857 0 : TVectorD * vecSec = new TVectorD(entriesFit); //
1858 0 : TVectorD * vecPhi = new TVectorD(entriesFit); //
1859 0 : TVectorD * vecTime = new TVectorD(entriesFit); //
1860 0 : TVectorD * vecTOFBC = new TVectorD(entriesFit); //
1861 0 : TVectorD * deltaZITS = new TVectorD(entriesFit); //
1862 :
1863 0 : for (Int_t i=0; i<entriesFit; i++){
1864 0 : Int_t index=gRandom->Rndm()*entriesFit0;
1865 0 : (*deltaZTOF)[i]= chainDelta->GetVal(0)[index];
1866 0 : (*deltaZTRD)[i]= chainDelta->GetVal(1)[index];
1867 0 : (*vecZ)[i]= chainDelta->GetVal(2)[index];
1868 0 : (*vecR)[i]= chainDelta->GetVal(3)[index];
1869 0 : (*vecSec)[i]= chainDelta->GetVal(4)[index];
1870 0 : (*vecPhi)[i]= chainDelta->GetVal(5)[index];
1871 0 : (*vecTime)[i]= chainDelta->GetVal(6)[index];
1872 0 : (*vecTOFBC)[i]= chainDelta->GetVal(7)[index];
1873 0 : (*deltaZITS)[i]= chainDelta->GetVal(8)[index];
1874 : //
1875 0 : if ( (*vecR)[i] < kInvalidR ) continue; // there was no cluster at this point
1876 : // if needed, fix alignment bug
1877 0 : if (fixAlignmentBugForFile) {
1878 0 : float dyTOF = chainDelta->GetVal(9)[index];
1879 0 : float dyTRD = chainDelta->GetVal(10)[index];
1880 0 : float q2pt = chainDelta->GetVal(11)[index];
1881 : //
1882 0 : float dzITS = (*deltaZITS)[i];
1883 0 : float dzTOF = (*deltaZTOF)[i];
1884 0 : float dzTRD = (*deltaZTRD)[i];
1885 0 : float rTOF = (*vecR)[i], rTRD = rTOF;
1886 0 : float zTOF = (*vecZ)[i] + dzTOF-dzITS , zTRD = (*vecZ)[i] + dzTRD-dzITS;
1887 0 : float phiTOF= (*vecPhi)[i], phiTRD = phiTOF;
1888 :
1889 0 : int rocID = TMath::Nint((*vecSec)[i]);
1890 :
1891 0 : if (dyTOF>kInvalidRes) {
1892 0 : FixAlignmentBug(rocID, q2pt, bz, phiTOF, rTOF, zTOF, dyTOF, dzTOF);
1893 0 : (*deltaZTOF)[i] = dzTOF;
1894 0 : (*vecR)[i] = rTOF;
1895 0 : (*vecZ)[i] = zTOF;
1896 0 : (*vecPhi)[i] = phiTOF;
1897 0 : }
1898 0 : if (dyTRD>kInvalidRes) {
1899 0 : FixAlignmentBug(rocID, q2pt, bz, phiTRD, rTRD, zTRD, dyTRD, dzTRD);
1900 : //
1901 0 : (*deltaZTRD)[i] = dzTRD;
1902 0 : (*vecR)[i] = rTRD;
1903 0 : (*vecZ)[i] = zTRD;
1904 0 : (*vecPhi)[i] = phiTRD;
1905 0 : }
1906 : //
1907 0 : }
1908 0 : }
1909 0 : AliSysInfo::AddStamp("FitDrift.EndFill",4,1,0);
1910 :
1911 :
1912 : //
1913 : // 1.) Make robust first estimate
1914 : //
1915 0 : TVectorD paramRobust(5);
1916 0 : TVectorD paramTRD(5);
1917 0 : TVectorD paramTOF(5);
1918 0 : TVectorD paramRobustBC(5);
1919 0 : TVectorD paramTRDBC(5);
1920 0 : TVectorD paramTOFBC(5);
1921 : //
1922 0 : AliSysInfo::AddStamp("FitDrift.StartRobust",2,0,0);
1923 0 : TF1 * fpol1 = new TF1("f1","[0]+[1]*x",0,250);
1924 : //
1925 0 : if (pcstream->GetFile()->Get("robustFit")==0){
1926 0 : for (Int_t iter=0; iter<20; iter++){
1927 0 : AliSysInfo::AddStamp("FitDrift.RobustIter",3,iter,0);
1928 0 : TLinearFitter *fitterRobust= new TLinearFitter(4,TString::Format("hyp%d",3).Data());
1929 0 : TLinearFitter *fitterTRD= new TLinearFitter(4,TString::Format("hyp%d",3).Data());
1930 0 : TLinearFitter *fitterTOF= new TLinearFitter(4,TString::Format("hyp%d",3).Data());
1931 0 : TLinearFitter *fitterRobustBC= new TLinearFitter(4,TString::Format("hyp%d",3).Data());
1932 0 : TLinearFitter *fitterTRDBC= new TLinearFitter(4,TString::Format("hyp%d",3).Data());
1933 0 : TLinearFitter *fitterTOFBC= new TLinearFitter(4,TString::Format("hyp%d",3).Data());
1934 0 : TH2F * hisDeltaZ[12];
1935 0 : for (Int_t ihis=0; ihis<12; ihis++) {
1936 0 : hisDeltaZ[ihis] = new TH2F("hisZ","hisZ",21,40,250,100,-5,5);
1937 : }
1938 :
1939 0 : Int_t maxPoints= TMath::Min(maxPointsRobust,entriesFit);
1940 0 : for (Int_t ipoint=0; ipoint<maxPoints; ipoint++){
1941 0 : Double_t radius= (*vecR)[ipoint];
1942 :
1943 0 : if (radius<kInvalidR) continue; // bad point
1944 :
1945 0 : Double_t pvecFit[10]={0};
1946 0 : if (gRandom->Rndm()>0.1) continue;
1947 0 : Int_t sector = TMath::Nint((*vecSec)[ipoint]);
1948 0 : Double_t side = -1.+((sector%36)<18)*2.;
1949 0 : Double_t z= (*vecZ)[ipoint];
1950 0 : double maxZ = kMaxZSect[side<0];
1951 0 : Double_t drift = (side>0) ? maxZ-(*vecZ)[ipoint] : (*vecZ)[ipoint]+maxZ;
1952 0 : if (drift>maxZ) continue;
1953 0 : Double_t phi = (*vecPhi)[ipoint];
1954 0 : Double_t gy = TMath::Sin(phi)*radius;
1955 0 : Float_t tofBC=(*vecTOFBC)[ipoint];
1956 0 : pvecFit[0]= side; // z shift (cm)
1957 0 : pvecFit[1]= drift*gy/maxZ; // global y gradient
1958 0 : pvecFit[2]= drift; // drift length
1959 0 : Double_t dZTOF=(*deltaZTOF)[ipoint];
1960 0 : Double_t dZTRD=(*deltaZTRD)[ipoint];
1961 0 : Int_t hisIndex=(((1-side)/2)*6);
1962 :
1963 0 : if (iter==0){
1964 0 : if ( dZTOF>kInvalidRes &&TMath::Abs(dZTOF)<kMaxDist0) fitterRobust->AddPoint(pvecFit,dZTOF*side,1);
1965 0 : if ( dZTRD>kInvalidRes &&TMath::Abs(dZTRD)<kMaxDist0) fitterRobust->AddPoint(pvecFit,dZTRD*side,1);
1966 : }else{
1967 0 : Double_t expected = paramRobust[0]+paramRobust[1]*pvecFit[0]+paramRobust[2]*pvecFit[1]+paramRobust[3]*pvecFit[2];
1968 0 : if ( dZTRD>kInvalidRes &&TMath::Abs(dZTRD*side-expected)<kMaxDist1) {
1969 0 : fitterRobust->AddPoint(pvecFit,dZTRD*side,1);
1970 0 : fitterTRD->AddPoint(pvecFit,dZTRD*side,1);
1971 0 : hisDeltaZ[hisIndex+0]->Fill(drift, dZTRD*side-expected);
1972 0 : hisDeltaZ[hisIndex+1]->Fill(drift, dZTRD*side-expected);
1973 0 : if (tofBC>kBCcutMin&&tofBC<kBCcutMax){
1974 0 : fitterRobustBC->AddPoint(pvecFit,dZTRD*side,1);
1975 0 : fitterTRDBC->AddPoint(pvecFit,dZTRD*side,1);
1976 0 : hisDeltaZ[hisIndex+3]->Fill(drift, dZTRD*side-expected);
1977 0 : hisDeltaZ[hisIndex+4]->Fill(drift, dZTRD*side-expected);
1978 : }
1979 : }
1980 0 : if ( dZTOF>kInvalidRes &&TMath::Abs(dZTOF*side-expected)<kMaxDist1) {
1981 0 : hisDeltaZ[hisIndex+0]->Fill(drift, dZTOF*side-expected);
1982 0 : hisDeltaZ[hisIndex+2]->Fill(drift, dZTOF*side-expected);
1983 0 : fitterRobust->AddPoint(pvecFit,dZTOF*side,1);
1984 0 : fitterTOF->AddPoint(pvecFit,dZTOF*side,1);
1985 0 : if (tofBC>kBCcutMin&&tofBC<kBCcutMax){
1986 0 : fitterRobustBC->AddPoint(pvecFit,dZTOF*side,1);
1987 0 : fitterTOFBC->AddPoint(pvecFit,dZTOF*side,1);
1988 0 : hisDeltaZ[hisIndex+3]->Fill(drift, dZTOF*side-expected);
1989 0 : hisDeltaZ[hisIndex+5]->Fill(drift, dZTOF*side-expected);
1990 : }
1991 : }
1992 :
1993 0 : if (gRandom->Rndm()<kDumpSample){
1994 0 : Double_t cTime = (*vecTime)[ipoint];
1995 0 : (*pcstream)<<"dumpSample"<<
1996 0 : "iter="<<iter<<
1997 0 : "run="<<runNumber<<
1998 0 : "ctime="<<cTime<<
1999 0 : "sector="<<sector<<
2000 0 : "side="<<side<<
2001 0 : "drift="<<drift<<
2002 0 : "radius="<<radius<<
2003 0 : "phi="<<phi<<
2004 0 : "tofBC="<<tofBC<<
2005 0 : "gy="<<gy<<
2006 0 : "z="<<z<<
2007 0 : "expected="<<expected<<
2008 0 : "paramRobust.="<<¶mRobust<< // drift fit using all tracks
2009 0 : "dZTOF="<<dZTOF<<
2010 0 : "dZTRD="<<dZTRD<<
2011 : "\n";
2012 0 : }
2013 0 : }
2014 0 : }
2015 0 : printf("iter=%d\n",iter);
2016 0 : if (fitterRobust->GetNpoints()<kMinEntries*0.1){
2017 0 : ::Error("fitDrift FAILED","Not enough points in the chain. Iter=%d\tN=%d" , iter, fitterRobust->GetNpoints());
2018 0 : delete pcstream;
2019 0 : return kFALSE;
2020 : }
2021 : //fitterRobust->EvalRobust(robFraction, regRobust); // ROOT has to be patched to use regularization
2022 0 : fitterRobust->EvalRobust(robFraction);
2023 : //fitterRobust->Eval(); // EvalRobust sometimes failed - looks like related to the random selection of subset of data - can be only one side
2024 0 : fitterRobust->GetParameters(paramRobust);
2025 0 : Double_t npoints= fitterRobust->GetNpoints();
2026 0 : Double_t chi2= TMath::Sqrt(fitterRobust->GetChisquare()/npoints);
2027 0 : paramRobust.Print();
2028 :
2029 0 : Int_t nTRD = fitterTRD->GetNpoints();
2030 0 : Int_t nTOF = fitterTOF->GetNpoints();
2031 0 : Int_t nRobustBC = fitterRobustBC->GetNpoints();
2032 0 : Int_t nTRDBC = fitterTRDBC->GetNpoints();
2033 0 : Int_t nTOFBC = fitterTOFBC->GetNpoints();
2034 :
2035 0 : if (nRobustBC>kMinEntries*0.1){ fitterRobustBC->EvalRobust(robFraction); fitterRobustBC->GetParameters(paramRobustBC);}
2036 0 : if (nTRD>kMinEntries*0.1){ fitterTRD->EvalRobust(robFraction); fitterTRD->GetParameters(paramTRD);}
2037 0 : if (nTOF>kMinEntries*0.1){ fitterTOF->EvalRobust(robFraction); fitterTOF->GetParameters(paramTOF);}
2038 0 : if (nTRDBC>kMinEntries*0.1){ fitterTRDBC->EvalRobust(robFraction); fitterTRDBC->GetParameters(paramTRDBC);}
2039 0 : if (nTOFBC>kMinEntries*0.1){ fitterTOFBC->EvalRobust(robFraction); fitterTOFBC->GetParameters(paramTOFBC); }
2040 :
2041 : // if (nRobustBC>kMinEntries*0.1){ fitterRobustBC->EvalRobust(robFraction, regRobust); fitterRobustBC->GetParameters(paramRobustBC);}
2042 : // if (nTRD>kMinEntries*0.1){ fitterTRD->EvalRobust(robFraction, regRobust); fitterTRD->GetParameters(paramTRD);}
2043 : // if (nTOF>kMinEntries*0.1){ fitterTOF->EvalRobust(robFraction, regRobust); fitterTOF->GetParameters(paramTOF);}
2044 : // if (nTRDBC>kMinEntries*0.1){ fitterTRDBC->EvalRobust(robFraction, regRobust); fitterTRDBC->GetParameters(paramTRDBC);}
2045 : // if (nTOFBC>kMinEntries*0.1){ fitterTOFBC->EvalRobust(robFraction, regRobust); fitterTOFBC->GetParameters(paramTOFBC); }
2046 : //
2047 0 : TGraphErrors * grDeltaZ[12] = {0};
2048 0 : TGraphErrors * grRMSZ[12] = {0};
2049 0 : TVectorD * fitDeltaZ[12] = {0};
2050 0 : TObjArray fitArray(3);
2051 0 : for (Int_t ihis=0; ihis<12; ihis++){
2052 0 : hisDeltaZ[ihis]->FitSlicesY(0,0,-1,0,"QNR",&fitArray);
2053 0 : grDeltaZ[ihis] = new TGraphErrors(((TH1D*)fitArray.At(1)));
2054 0 : grRMSZ[ihis] = new TGraphErrors(((TH1D*)fitArray.At(2)));
2055 0 : fitDeltaZ[ihis] = new TVectorD(2);
2056 0 : grDeltaZ[ihis]->Fit(fpol1,"q","q");
2057 0 : fpol1->GetParameters(fitDeltaZ[ihis]->GetMatrixArray());
2058 0 : fitArray.Delete();
2059 0 : delete hisDeltaZ[ihis];
2060 0 : hisDeltaZ[ihis] = 0;
2061 : }
2062 :
2063 0 : delete fitterRobust;
2064 0 : delete fitterTRD;
2065 0 : delete fitterTOF;
2066 0 : delete fitterRobustBC;
2067 0 : delete fitterTRDBC;
2068 0 : delete fitterTOFBC;
2069 :
2070 0 : (*pcstream)<<"robustFit"<<
2071 0 : "iter="<<iter<<
2072 0 : "time0="<<time0<<
2073 0 : "time1="<<time1<<
2074 0 : "run="<<runNumber<<
2075 0 : "npoints="<<npoints<<
2076 0 : "chi2="<<chi2<<
2077 0 : "nTRD="<<nTRD<< // number of TRD points
2078 0 : "nTOF="<<nTOF<< // number of TRD points
2079 0 : "nTRDBC="<<nTRDBC<< // number of TRD points BC
2080 0 : "nTOFBC="<<nTOFBC<< // number of TRD points BC
2081 : //
2082 0 : "paramRobust.="<<¶mRobust<< // drift fit using all tracks
2083 0 : "paramTRD.="<<¶mTRD<< // drift fit using TRD tracks
2084 0 : "paramTOF.="<<¶mTOF<< // drift fit using TOF tracks
2085 0 : "paramRobustBC.="<<¶mRobustBC<< // drift fit using all tracks with BC
2086 0 : "paramTRDBC.="<<¶mTRDBC<< // drift fit using TRD tracks with BC
2087 0 : "paramTOFBC.="<<¶mTOFBC; // drift fit using TOF tracks with BC
2088 :
2089 0 : for (Int_t ihis=0; ihis<12; ihis++){
2090 0 : (*pcstream)<<"robustFit"<<
2091 0 : TString::Format("grDeltaZ%d.=",ihis).Data()<<grDeltaZ[ihis]<< // residual histogram drift fit - mean
2092 0 : TString::Format("grRMSZ%d.=",ihis).Data()<<grDeltaZ[ihis]<< // residual histogram drift fit - rms
2093 0 : TString::Format("fitDeltaZ%d.=",ihis).Data()<<fitDeltaZ[ihis]; // residual histogram drift fit - linear fit
2094 : }
2095 0 : (*pcstream)<<"robustFit"<<
2096 : "\n";
2097 0 : for (Int_t ihis=0; ihis<12; ihis++){
2098 0 : delete grDeltaZ[ihis];
2099 0 : delete grRMSZ[ihis];
2100 0 : delete fitDeltaZ[ihis];
2101 : }
2102 0 : }
2103 : // delete pcstream;
2104 : //pcstream = new TTreeSRedirector("fitDrift.root","update");
2105 : }
2106 0 : delete fpol1;
2107 : //
2108 0 : TTree * treeRobust= (TTree*)(pcstream->GetFile()->Get("robustFit"));
2109 0 : Int_t entriesR= treeRobust->Draw("paramRobust.fElements[0]:paramRobust.fElements[1]:paramRobust.fElements[2]:paramRobust.fElements[3]","","goffPara");
2110 0 : for (Int_t ipar=0; ipar<4; ipar++) {paramRobust[ipar]=TMath::Median(entriesR-5, &(treeRobust->GetVal(ipar)[5]));}
2111 : //
2112 0 : treeRobust->Draw("paramRobustBC.fElements[0]:paramRobustBC.fElements[1]:paramRobustBC.fElements[2]:paramRobustBC.fElements[3]","","goffPara");
2113 0 : for (Int_t ipar=0; ipar<4; ipar++) {paramRobustBC[ipar]=TMath::Median(entriesR-5, &(treeRobust->GetVal(ipar)[5]));}
2114 0 : treeRobust->Draw("paramTRD.fElements[0]:paramTRD.fElements[1]:paramTRD.fElements[2]:paramTRD.fElements[3]","","goffPara");
2115 0 : for (Int_t ipar=0; ipar<4; ipar++) {paramTRD[ipar]=TMath::Median(entriesR-5, &(treeRobust->GetVal(ipar)[5]));}
2116 0 : treeRobust->Draw("paramTOF.fElements[0]:paramTOF.fElements[1]:paramTOF.fElements[2]:paramTOF.fElements[3]","","goffPara");
2117 0 : for (Int_t ipar=0; ipar<4; ipar++) {paramTOF[ipar]=TMath::Median(entriesR-5, &(treeRobust->GetVal(ipar)[5]));}
2118 0 : treeRobust->Draw("paramTRDBC.fElements[0]:paramTRDBC.fElements[1]:paramTRDBC.fElements[2]:paramTRDBC.fElements[3]","","goffPara");
2119 0 : for (Int_t ipar=0; ipar<4; ipar++) {paramTRDBC[ipar]=TMath::Median(entriesR-5, &(treeRobust->GetVal(ipar)[5]));}
2120 0 : treeRobust->Draw("paramTOFBC.fElements[0]:paramTOFBC.fElements[1]:paramTOFBC.fElements[2]:paramTOFBC.fElements[3]","","goffPara");
2121 0 : for (Int_t ipar=0; ipar<4; ipar++) {paramTOFBC[ipar]=TMath::Median(entriesR-5, &(treeRobust->GetVal(ipar)[5]));}
2122 0 : paramRobust.Print();
2123 : //
2124 : //
2125 : // 3.) Make drift fit per time interval
2126 : //
2127 0 : Int_t nTimeBins= Int_t((time1-time0)/deltaT)+1;
2128 0 : Int_t nParam = nTimeBins+1;
2129 0 : TVectorD vecFit(nParam);
2130 0 : Double_t *pvecFit=vecFit.GetMatrixArray();
2131 : //TLinearFitter *fitterTRD= new TLinearFitter(nParam,TString::Format("hyp%d",nParam+1).Data());
2132 : //TLinearFitter *fitterTOF= new TLinearFitter(nParam,TString::Format("hyp%d",nParam+1).Data());
2133 : //
2134 0 : TObjArray arrayFit(3);
2135 0 : TH2F hisTRD("hisTRD","hisTRD",nTimeBins,time0,time1,100,-0.02,0.02);
2136 0 : TH2F hisTOF("hisTOF","hisTOF",nTimeBins,time0,time1,100,-0.02,0.02);
2137 :
2138 0 : for (Int_t iter=0; iter<10; iter++){
2139 0 : hisTRD.Reset();
2140 0 : hisTOF.Reset();
2141 0 : for (Int_t ipoint=0; ipoint<entriesFit; ipoint++){
2142 0 : if (ipoint%10!=iter) continue; // points correlated - can be skipped
2143 0 : Double_t radius= (*vecR)[ipoint];
2144 :
2145 0 : if (radius<kInvalidR) continue; // no cluster
2146 :
2147 0 : Int_t sector = TMath::Nint((*vecSec)[ipoint]);
2148 0 : Double_t side = -1.+((sector%36)<18)*2.;
2149 0 : Double_t z= (*vecZ)[ipoint];
2150 0 : double maxZ = kMaxZSect[side<0];
2151 0 : Double_t drift = (side>0) ? maxZ-(*vecZ)[ipoint] : (*vecZ)[ipoint]+maxZ;
2152 0 : if (drift>maxZ) continue;
2153 0 : Double_t phi = (*vecPhi)[ipoint];
2154 0 : Double_t gy = TMath::Sin(phi)*radius;
2155 0 : pvecFit[0]= side; // z shift (cm)
2156 0 : pvecFit[1]= drift*gy/maxZ; // global y gradient
2157 0 : pvecFit[2]= drift; // drift length
2158 0 : Double_t expected = paramRobust[0]+paramRobust[1]*pvecFit[0]+paramRobust[2]*pvecFit[1]+paramRobust[3]*pvecFit[2];
2159 0 : Int_t time=(*vecTime)[ipoint];
2160 : //
2161 0 : Double_t dZTOF=(*deltaZTOF)[ipoint];
2162 0 : if (dZTOF>kInvalidRes) {
2163 0 : if (TMath::Abs(dZTOF*side-expected)<kMaxDist1) {
2164 : dZTOF *= side;
2165 : // fitterTOF->AddPoint(pvecFit,dZTOF,1);
2166 0 : hisTOF.Fill(time,(dZTOF-expected)/drift,drift/maxZ);
2167 : }
2168 : }
2169 0 : Double_t dZTRD=(*deltaZTRD)[ipoint];
2170 0 : if ( dZTRD>kInvalidRes) {
2171 0 : if (TMath::Abs(dZTRD*side-expected)<kMaxDist1) {
2172 : dZTRD*=side;
2173 : //fitterTOF->AddPoint(pvecFit,dZTRD,1);
2174 0 : hisTRD.Fill(time,(dZTRD-expected)/drift,drift/maxZ);
2175 : }
2176 : }
2177 :
2178 0 : }
2179 0 : printf("iter=%d\n",iter);
2180 : TGraphErrors *grTRD=NULL, *grTOF=NULL;
2181 : TGraphErrors *grTRD2=NULL, *grTOF2=NULL;
2182 0 : Int_t nclTRD=hisTRD.GetEntries();
2183 0 : Int_t nclTOF=hisTOF.GetEntries();
2184 0 : if (nclTRD>kMinEntries){
2185 0 : hisTRD.FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2186 0 : grTRD = new TGraphErrors((TH1D*)arrayFit.At(1));
2187 0 : grTRD2 = new TGraphErrors((TH1D*)arrayFit.At(2));
2188 0 : arrayFit.Delete();
2189 : }
2190 0 : if (nclTOF>kMinEntries){
2191 0 : hisTOF.FitSlicesY(0,0,-1,0,"QNR",&arrayFit);
2192 0 : grTOF = new TGraphErrors((TH1D*)arrayFit.At(1));
2193 0 : grTOF2 = new TGraphErrors((TH1D*)arrayFit.At(2));
2194 0 : arrayFit.Delete();
2195 : }
2196 0 : if (grTRD==NULL) {
2197 : grTRD=grTOF; // we should have at minimum one of the histograms not empty
2198 : grTRD2=grTOF2;
2199 0 : }
2200 0 : if (grTOF==NULL) {
2201 : grTOF=grTRD;
2202 : grTOF2=grTRD2;
2203 0 : }
2204 0 : (*pcstream)<<"fitTime"<<
2205 0 : "iter="<<iter<<
2206 0 : "nclTRD="<<nclTRD<< //
2207 0 : "nclTOF="<<nclTOF<< //
2208 0 : "grTRD.="<<grTRD<< // time dependent drift correction TRD - mean in time bin
2209 0 : "grTOF.="<<grTOF<< // time dependent drift correction TOF - mean in time bin
2210 0 : "grTRD2.="<<grTRD2<< // time dependent drift correction TRD - sigma in time bin
2211 0 : "grTOF2.="<<grTOF2<< // time dependent drift correction TOF - sigma in time bin
2212 0 : "time0="<<time0<<
2213 0 : "time1="<<time1<<
2214 0 : "run="<<runNumber<<
2215 0 : "paramRobust.="<<¶mRobust<< // time independent parameters
2216 0 : "paramTRD.="<<¶mTRD<< // drift fit using TRD tracks
2217 0 : "paramTOF.="<<¶mTOF<< // drift fit using TOF tracks
2218 0 : "paramRobustBC.="<<¶mRobustBC<< // drift fit using all tracks with BC
2219 0 : "paramTRDBC.="<<¶mTRDBC<< // drift fit using TRD tracks with BC
2220 0 : "paramTOFBC.="<<¶mTOFBC<< // drift fit using TOF tracks with BC
2221 : "\n";
2222 0 : }
2223 0 : delete pcstream;
2224 0 : pcstream = new TTreeSRedirector("fitDrift.root","update");
2225 : //
2226 : // 3.) Make local regression and median fit
2227 : //
2228 0 : TTree * treeFit= (TTree*)(pcstream->GetFile()->Get("fitTime"));
2229 0 : Int_t entriesGr = treeFit->Draw("grTRD.fY:grTOF.fY:grTRD.fX:Iteration$","1","goffpara");
2230 0 : Int_t nbins = TMath::MaxElement(entriesGr, treeFit->GetV4())+1;
2231 :
2232 0 : Double_t dtime0=0,dtime1=0;
2233 0 : if (entriesGr) TStatToolkit::GetMinMax(treeFit->GetV3(),entriesGr,dtime0,dtime1);
2234 0 : Int_t ngraphs =entriesGr/nbins;
2235 : //
2236 : // 3.a) local regression fit
2237 : //
2238 0 : treeFit->Draw("grTRD.fY-grTRD.fY[Iteration$+1]:grTRD.fX","1","goff");
2239 0 : Double_t deltaRMS,deltaMean;
2240 0 : AliMathBase::EvaluateUni(entriesGr,treeFit->GetV1(),deltaMean, deltaRMS,entriesGr*0.8);
2241 0 : TCut cutTRD=TString::Format("abs(grTRD.fY-grTOF.fY)<0.01&&abs(grTRD.fY-grTRD.fY[Iteration$+1]-%f)<%f",deltaMean, 3*deltaRMS).Data();
2242 0 : TCut cutTOF=TString::Format("abs(grTRD.fY-grTOF.fY)<0.01&&abs(grTOF.fY-grTOF.fY[Iteration$+1]-%f)<%f",deltaMean, 3*deltaRMS).Data();
2243 :
2244 0 : THnD *hN = new THnD("hN","hN", 1, &nbins, &dtime0, &dtime1);
2245 0 : AliNDLocalRegression * pfitTRD = new AliNDLocalRegression("pfitTRD","pfitTRD");
2246 0 : AliNDLocalRegression * pfitTOF = new AliNDLocalRegression("pfitTOF","pfitTOF");
2247 0 : pfitTRD->SetHistogram((THn*)(hN->Clone()));
2248 0 : pfitTOF->SetHistogram((THn*)(hN->Clone()));
2249 0 : pfitTRD->MakeFit(treeFit, "grTRD.fY:grTRD.fEY+0.01", "grTRD.fX",cutTRD, TString::Format("(grTRD.fX/grTRD.fX)+%f",sigmaT),"2:2",0.0001);
2250 0 : pfitTOF->MakeFit(treeFit, "grTOF.fY:grTOF.fEY+0.01", "grTOF.fX",cutTOF, TString::Format("(grTRD.fX/grTRD.fX)+%f",sigmaT),"2:2",0.0001);
2251 0 : AliNDLocalRegression::AddVisualCorrection(pfitTRD,104);
2252 0 : AliNDLocalRegression::AddVisualCorrection(pfitTOF,204);
2253 : //
2254 : //
2255 0 : TVectorD medianTRD(nbins);
2256 0 : TVectorD medianTOF(nbins);
2257 0 : TVectorD medianQA(nbins);
2258 0 : TVectorD regTRD(nbins);
2259 0 : TVectorD regTOF(nbins);
2260 0 : TVectorD regQA(nbins);
2261 : //
2262 0 : TVectorD rmsTRD(nbins);
2263 0 : TVectorD rmsTOF(nbins);
2264 0 : TVectorD rmsQA(nbins);
2265 0 : TVectorD vecTimeg(nbins);
2266 0 : TVectorD vecWorking(entriesGr);
2267 0 : treeFit->Draw("grTRD.fY:grTOF.fY:grTRD.fX:Iteration$","1","goffpara");
2268 0 : for (Int_t itype=0; itype<2; itype++){
2269 0 : for (Int_t itime=0; itime<nbins; itime++){
2270 0 : for (Int_t igr=0; igr<ngraphs; igr++){
2271 0 : Int_t index=igr*nbins+itime;
2272 0 : vecWorking[igr]=treeFit->GetVal(itype)[index];
2273 : }
2274 0 : Double_t grtime = treeFit->GetVal(2)[itime];
2275 0 : vecTimeg[itime]=treeFit->GetVal(2)[itime];
2276 0 : if (itype==0) {
2277 0 : medianTRD[itime]=TMath::Median(ngraphs,vecWorking.GetMatrixArray());
2278 0 : regTRD[itime]= pfitTRD->Eval(&grtime);
2279 0 : rmsTRD[itime]=TMath::RMS(ngraphs,vecWorking.GetMatrixArray())/TMath::Sqrt(ngraphs);
2280 0 : }
2281 0 : if (itype==1) {
2282 0 : medianTOF[itime]=TMath::Median(ngraphs,vecWorking.GetMatrixArray());
2283 0 : regTOF[itime]= pfitTOF->Eval(&grtime);
2284 0 : rmsTOF[itime]=TMath::RMS(ngraphs,vecWorking.GetMatrixArray())/TMath::Sqrt(ngraphs);
2285 0 : }
2286 0 : }
2287 : }
2288 0 : for (Int_t itime=0; itime<nbins; itime++){
2289 0 : medianQA[itime]= medianTRD[itime]-medianTOF[itime];
2290 0 : regQA[itime]= regTRD[itime]-regTOF[itime];
2291 0 : rmsQA[itime]=TMath::Sqrt(rmsTRD[itime]*rmsTRD[itime]+rmsTOF[itime]*rmsTOF[itime]);
2292 : }
2293 0 : TGraphErrors *grmedTRD = new TGraphErrors(nbins, vecTimeg.GetMatrixArray(), medianTRD.GetMatrixArray(),0, rmsTRD.GetMatrixArray());
2294 0 : TGraphErrors *grmedTOF = new TGraphErrors(nbins, vecTimeg.GetMatrixArray(), medianTOF.GetMatrixArray(),0, rmsTOF.GetMatrixArray());
2295 0 : TGraphErrors *grmedQA = new TGraphErrors(nbins, vecTimeg.GetMatrixArray(), medianQA.GetMatrixArray(),0, rmsQA.GetMatrixArray());
2296 0 : TGraphErrors *grregTRD = new TGraphErrors(nbins, vecTimeg.GetMatrixArray(), regTRD.GetMatrixArray(),0, rmsTRD.GetMatrixArray());
2297 0 : TGraphErrors *grregTOF = new TGraphErrors(nbins, vecTimeg.GetMatrixArray(), regTOF.GetMatrixArray(),0, rmsTOF.GetMatrixArray());
2298 0 : TGraphErrors *grregQA = new TGraphErrors(nbins, vecTimeg.GetMatrixArray(), regQA.GetMatrixArray(),0, rmsQA.GetMatrixArray());
2299 : //
2300 0 : (*pcstream)<<"fitTimeStat"<<
2301 0 : "grTRDReg.="<<grregTRD<< // time dependent drift correction TRD - regression estimator
2302 0 : "grTOFReg.="<<grregTOF<< // time dependent drift correction TOF - regression estimator
2303 0 : "grQAReg.="<<grregQA<< // time dependent drift correction TOF - regression estimator
2304 0 : "grTRDMed.="<<grmedTRD<< // time dependent drift correction TRD - median estimator
2305 0 : "grTOFMed.="<<grmedTOF<< // time dependent drift correction TOF - median estimator
2306 0 : "grQAMed.="<<grmedQA<< // time dependent drift correction TOF - median estimator
2307 0 : "time0="<<time0<<
2308 0 : "time1="<<time1<<
2309 0 : "run="<<runNumber<<
2310 0 : "paramRobust.="<<¶mRobust<< // time independent parameters
2311 : "\n";
2312 0 : delete grmedTRD;
2313 0 : delete grmedTOF;
2314 0 : delete grmedQA;
2315 0 : delete grregTRD;
2316 0 : delete grregTOF;
2317 0 : delete grregQA;
2318 :
2319 0 : delete deltaZTOF;
2320 0 : delete deltaZTRD;
2321 0 : delete vecZ;
2322 0 : delete vecR;
2323 0 : delete vecSec;
2324 0 : delete vecPhi;
2325 0 : delete vecTime;
2326 0 : delete vecTOFBC;
2327 :
2328 0 : delete pcstream;
2329 : return kTRUE;
2330 0 : }
2331 :
2332 :
2333 :
2334 : void AliTPCcalibAlignInterpolation::MakeNDFit(const char * inputFile, const char * inputTree, Float_t sector0, Float_t sector1, Float_t theta0, Float_t theta1){
2335 : //
2336 : ///
2337 : /// Make ND local regression, QA for later usage
2338 : /// Parameters:
2339 :
2340 : // Algorithm:
2341 : // 1.) Make NDLocal regression fits
2342 : // 2.) Make regularization (smoothing)
2343 : // 3.) Make QA trending variable
2344 : // 4.) Make QA plots
2345 : // 5.) Export QA trending variables into trending tree
2346 : //
2347 : /*
2348 : Example usage:
2349 : const char * inputFile="ResidualMapFull_1.root"
2350 : const char * inputTree="deltaRPhiTPCITSTRDDist"
2351 : Float_t sector0=3, sector1 =5;
2352 : Float_t theta0=0, theta1=1;
2353 : AliTPCcalibAlignInterpolation::MakeNDFit(inputFile,inputTree, sector0,sector1, theta0,theta1);
2354 : */
2355 :
2356 0 : TTreeSRedirector * pcstream = new TTreeSRedirector(TString::Format("%sFit_sec%d_%d_theta%d_%d.root",inputTree,Int_t(sector0),Int_t(sector1),Int_t(theta0),Int_t(theta1)).Data(),"recreate");
2357 0 : TTreeSRedirector * pcstreamFit = new TTreeSRedirector(TString::Format("fitTree_%sFit_sec%d_%d_theta%d_%d.root",inputTree,Int_t(sector0),Int_t(sector1),Int_t(theta0),Int_t(theta1)).Data(),"recreate");
2358 0 : Int_t runNumber=TString(gSystem->Getenv("runNumber")).Atoi();
2359 0 : TFile * fdist = TFile::Open(inputFile);
2360 0 : if (!fdist){
2361 0 : ::Error("MakeNDFit","Intput file %s not accessible\n",inputFile);
2362 0 : return;
2363 : }
2364 0 : TTree *treeDist = (TTree*)fdist->Get(inputTree);
2365 0 : if (!treeDist){
2366 0 : ::Error("MakeNDFit","Intput tree %s not accessible\n",inputTree);
2367 0 : return;
2368 : }
2369 0 : TTree *treeMeta = (TTree*)fdist->Get("metaData");
2370 0 : pcstream->GetFile()->cd();
2371 0 : TTree *treeMetaCopy = treeMeta->CopyTree("1");
2372 0 : treeMetaCopy->Write("metaData");
2373 0 : delete treeMetaCopy;
2374 : //
2375 : // 1.) Make NDLocal regression fits
2376 : //
2377 : const Double_t pxmin=8.48499984741210938e+01; //param.GetPadRowRadii(0,0)-param.GetPadPitchLength(0,0)/2
2378 : const Double_t pxmax=2.46600006103515625e+02; //2.46600006103515625e+02param.GetPadRowRadii(36,param.GetNRow(36)-1)+param.GetPadPitchLength(36,param.GetNRow(36)-1)/2.
2379 : Int_t ndim=4;
2380 0 : Int_t nbins[4]= {30, (Int_t)((sector1-sector0-0.1)*15), int(abs(theta1-theta0)*10), 3}; // {radius, phi bin, }
2381 0 : Double_t xmin[4] = {pxmin, sector0+0.05, theta0, -2.0};
2382 0 : Double_t xmax[4] = {pxmax, sector1-0.05, theta1, 2.0};
2383 : //
2384 :
2385 0 : THnF* hN= new THnF("exampleFit","exampleFit", ndim, nbins, xmin,xmax);
2386 0 : treeDist->SetAlias("isOK","rms>0&&vecLTM.fElements[1]*binMedian>0");
2387 0 : treeDist->SetAlias("delta","vecLTM.fElements[1]");
2388 0 : TCut cutFit="isOK";
2389 0 : TCut cutAcceptFit=TString::Format("sectorCenter>%f&§orCenter<%f&&kZCenter>%f&&kZCenter<%f", sector0-0.5,sector1+0.5,theta0,theta1).Data();
2390 0 : TCut cutAcceptDraw=TString::Format("sectorCenter>%f&§orCenter<%f&&kZCenter>%f&&kZCenter<%f", sector0,sector1,theta0,theta1).Data();
2391 :
2392 0 : AliNDLocalRegression *fitCorrs[6]={0};
2393 0 : for (Int_t icorr=0; icorr<1; icorr++){
2394 0 : fitCorrs[icorr]= new AliNDLocalRegression();
2395 0 : fitCorrs[icorr]->SetName(TString::Format("%sFit%d_sec%d_%d_theta%d_%d",inputTree,icorr, Int_t(sector0),Int_t(sector1),Int_t(theta0),Int_t(theta1)).Data());
2396 0 : Int_t hashIndex=fitCorrs[icorr]->GetVisualCorrectionIndex();
2397 0 : fitCorrs[icorr]->SetHistogram((THn*)(hN->Clone()));
2398 0 : TStopwatch timer;
2399 0 : fitCorrs[0]->SetStreamer(pcstream);
2400 0 : if (icorr==0) fitCorrs[icorr]->MakeFit(treeDist,"delta:(1+sqrt(abs(qptCenter)))/sqrt(entries)", "RCenter:sectorCenter:kZCenter:qptCenter",cutFit+cutAcceptFit,"3:0.05:0.1:18","2:2:2:2",0.0001);
2401 0 : timer.Print();
2402 0 : AliNDLocalRegression::AddVisualCorrection(fitCorrs[icorr]);
2403 0 : treeDist->SetAlias(TString::Format("delta_Fit%d",icorr).Data(),TString::Format("AliNDLocalRegression::GetCorrND(%d,RCenter,sectorCenter,kZCenter,qptCenter+0)",hashIndex).Data());
2404 0 : }
2405 : //
2406 : // 2.) Make regularization (smoothing)
2407 : //
2408 : Int_t nDims=4;
2409 0 : Int_t indexes[4]={0,1,2,3};
2410 0 : Double_t relWeight0[12]={1,4,16, 1,4,16, 1,4,16, 1,4,16};
2411 0 : Double_t relWeightC[12]={0.5,4,16, 0.5,4,16, 0.5,4,16, 0.5,4,16};
2412 0 : fitCorrs[1]=(AliNDLocalRegression *)fitCorrs[0]->Clone();
2413 0 : fitCorrs[1]->AddWeekConstrainsAtBoundaries(nDims, indexes,relWeight0, 0);
2414 0 : fitCorrs[2]=(AliNDLocalRegression *)fitCorrs[1]->Clone();
2415 0 : fitCorrs[2]->AddWeekConstrainsAtBoundaries(nDims, indexes,relWeight0, 0);
2416 0 : fitCorrs[3]=(AliNDLocalRegression *)fitCorrs[1]->Clone();
2417 0 : fitCorrs[4]=(AliNDLocalRegression *)fitCorrs[2]->Clone();
2418 0 : fitCorrs[3]->AddWeekConstrainsAtBoundaries(nDims, indexes,relWeightC, 0, kTRUE);
2419 0 : fitCorrs[4]->AddWeekConstrainsAtBoundaries(nDims, indexes,relWeightC, 0, kTRUE);
2420 : //
2421 0 : fitCorrs[1]->SetName(TString::Format("%s_Smooth1",fitCorrs[0]->GetName()).Data());
2422 0 : fitCorrs[2]->SetName(TString::Format("%s_Smooth2",fitCorrs[1]->GetName()).Data());
2423 0 : fitCorrs[3]->SetName(TString::Format("%s_SmoothConst1",fitCorrs[0]->GetName()).Data());
2424 0 : fitCorrs[4]->SetName(TString::Format("%s_SmoothConst2",fitCorrs[1]->GetName()).Data());
2425 0 : for (Int_t icorr=1; icorr<5; icorr++){
2426 0 : Int_t hashIndex=fitCorrs[icorr]->GetVisualCorrectionIndex();
2427 0 : AliNDLocalRegression::AddVisualCorrection(fitCorrs[icorr]);
2428 0 : treeDist->SetAlias(TString::Format("delta_Fit%d",icorr).Data(),TString::Format("AliNDLocalRegression::GetCorrND(%d,RCenter,sectorCenter,kZCenter,qptCenter+0)",hashIndex).Data());
2429 : }
2430 : //
2431 : // 3.) Make QA variables and store fits
2432 : //
2433 0 : gStyle->SetOptFit(1);
2434 0 : treeDist->Draw(">>drawlist3",cutFit+cutAcceptDraw,"entrylist");
2435 0 : TEntryList *elist = (TEntryList*)gDirectory->Get("drawlist3");
2436 0 : treeDist->SetEntryList(elist);
2437 :
2438 0 : Double_t quantiles[10]={0.001,0.01,0.05,0.1,0.2, 0.8,0.9,0.99,0.999};
2439 0 : Double_t deltas[10]={0};
2440 0 : const char * pchVar[9]={"delta-delta_Fit0",
2441 : "delta-delta_Fit1",
2442 : "delta-delta_Fit2",
2443 : "delta-delta_Fit3",
2444 : "delta-delta_Fit4",
2445 : "delta_Fit0-delta_Fit1",
2446 : "delta_Fit0-delta_Fit2",
2447 : "delta_Fit0-delta_Fit3",
2448 : "delta_Fit0-delta_Fit4"};
2449 0 : const char * pchTittle[9]={"map-fit_{0} (cm)",
2450 : "map-fit_{0Reg1} (cm)",
2451 : "map-fit_{0Reg2} (cm)",
2452 : "map-fit_{0Reg1Const} (cm)",
2453 : "map-fit_{0Reg2Const} (cm)",
2454 : "fit_{O}-fit_{0Reg1} (cm)",
2455 : "fit_{0}-fit_{0Reg2} (cm)",
2456 : "fit_{0}-fit_{0Reg1Const} (cm)",
2457 : "fit_{0}-fit_{0Reg2Const} (cm)"};
2458 0 : TGraph * grQuantiles[18]={0};
2459 0 : TVectorD vecRMS1(18);
2460 0 : for (Int_t idiff=0; idiff<18; idiff++){
2461 : Int_t chEntries=0;
2462 0 : if (idiff<9){
2463 0 : chEntries=treeDist->Draw(pchVar[idiff%9],"1","goff");
2464 0 : }
2465 0 : if (idiff>=9){
2466 0 : chEntries=treeDist->Draw(TString::Format("(%s)/(rms/sqrt(entries))",pchVar[idiff%9]).Data(),"1","goff");
2467 0 : }
2468 0 : for (Int_t iq=0; iq<10; iq++){
2469 0 : deltas[iq]=TMath::KOrdStat(chEntries,treeDist->GetV1(), Int_t(chEntries*quantiles[iq]));
2470 : }
2471 0 : grQuantiles[idiff]=new TGraph(10,quantiles,deltas);
2472 0 : grQuantiles[idiff]->SetTitle(pchTittle[idiff%9]);
2473 0 : Double_t mean, rms=0;
2474 0 : AliMathBase::EvaluateUni(chEntries,treeDist->GetV1(),mean, rms,0.80*chEntries);
2475 0 : vecRMS1[idiff]=rms;
2476 0 : }
2477 0 : vecRMS1.Print();
2478 : TH1* his=0;
2479 0 : TVectorD slopePols1(5);
2480 0 : for (Int_t i=0; i<5; i++){
2481 0 : treeDist->Draw(TString::Format("delta:delta_Fit%d>>hisQA2D%d",i,i).Data(),cutFit+cutAcceptDraw,"colzgoff");
2482 0 : his=treeDist->GetHistogram();
2483 0 : his->Fit("pol1");
2484 0 : slopePols1[i]= his->GetFunction("pol1")->GetParameter(1);
2485 : }
2486 :
2487 0 : TFile * fout = pcstream->GetFile();
2488 0 : pcstream->GetFile()->cd();
2489 0 : for (Int_t iter=0; iter<5; iter++){
2490 0 : if (TString(gSystem->Getenv("AliTPCcalibAlignInterpolation_MakeNDFit_KeepCovariance")).Atoi()==0){
2491 0 : if (iter!=0) fitCorrs[iter]->CleanCovariance();
2492 : }
2493 0 : fitCorrs[iter]->Write();
2494 0 : if (TString(gSystem->Getenv("AliTPCcalibAlignInterpolation_MakeNDFit_DumpFitTree")).Atoi()>0){
2495 0 : fitCorrs[iter]->DumpToTree(4, (*pcstreamFit)<<TString::Format("tree%s", fitCorrs[iter]->GetName()).Data());
2496 0 : }
2497 : }
2498 : //
2499 : // 4.) Make standard QA plot
2500 : //
2501 0 : gStyle->SetLabelSize(0.06,"XYZ");
2502 0 : gStyle->SetTitleSize(0.06,"XYZ");
2503 0 : TCanvas *canvasQA = new TCanvas("canvasQA","canvasQA",1200,1000);
2504 0 : canvasQA->Divide(1,3);
2505 : //
2506 : // 4.1) delta:fit correaltion
2507 0 : canvasQA->cd(1)->SetLogz();
2508 0 : treeDist->Draw("delta:delta_Fit0>>hisQA2D",cutFit+cutAcceptDraw,"colzgoff");
2509 0 : his=treeDist->GetHistogram();
2510 0 : his->GetXaxis()->SetTitle("#Delta_{fit} (cm)");
2511 0 : his->GetYaxis()->SetTitle("#Delta_{map} (cm)");
2512 0 : his->Fit("pol1");
2513 0 : Double_t slopePol1= his->GetFunction("pol1")->GetParameter(1);
2514 0 : his->Write("hisQA2D");
2515 0 : his->Draw("colz");
2516 : //
2517 : // 4.2) delta-fitReg0 and fitReg0-fitReg1
2518 0 : canvasQA->cd(2)->SetLogy();
2519 0 : treeDist->SetLineColor(1); treeDist->SetMarkerColor(1);
2520 0 : treeDist->Draw("delta-delta_Fit0>>hisQA1D(300)",cutFit+cutAcceptDraw,"goff");
2521 0 : his=treeDist->GetHistogram();
2522 0 : his->GetXaxis()->SetTitle("#Delta_{map}-#Delta_{fit} (cm)");
2523 0 : his->Fit("gaus");
2524 0 : Double_t mean=his->GetMean();
2525 0 : Double_t rms=his->GetRMS();
2526 0 : Double_t rmsG= his->GetFunction("gaus")->GetParameter(2);
2527 : //
2528 0 : treeDist->GetHistogram()->Write("hisQA1D");
2529 0 : treeDist->SetLineColor(2); treeDist->SetMarkerColor(2);
2530 0 : treeDist->Draw("delta_Fit0-delta_Fit1>>hisQA1DifFit(300)",cutFit+cutAcceptDraw,"same");
2531 0 : his=treeDist->GetHistogram();
2532 0 : Double_t meanDiffFit=his->GetMean();
2533 0 : Double_t rmsDiffFit=his->GetRMS();
2534 0 : treeDist->GetHistogram()->Write("hisQA1DifFit");
2535 : //
2536 : // 4.3) (delta-fitReg0) "pull"
2537 0 : treeDist->SetLineColor(1); treeDist->SetMarkerColor(1);
2538 0 : canvasQA->cd(3)->SetLogy();
2539 0 : treeDist->Draw("(delta_Fit0-delta)/(rms/sqrt(entries))>>hisQAPull",cutFit+cutAcceptDraw,"goff");
2540 0 : his=treeDist->GetHistogram();
2541 0 : his->Fit("gaus");
2542 0 : Double_t meanPull=his->GetMean();
2543 0 : Double_t rmsPull=his->GetRMS();
2544 0 : Double_t rmsPullG= his->GetFunction("gaus")->GetParameter(2);
2545 0 : treeDist->GetHistogram()->Write("hisQAPull");
2546 0 : treeDist->SetLineColor(2); treeDist->SetMarkerColor(2);
2547 0 : treeDist->Draw("(delta_Fit0-delta_Fit1)/(rms/sqrt(entries))>>hisQAPullDifFit",cutFit+cutAcceptDraw,"same");
2548 0 : his=treeDist->GetHistogram();
2549 0 : Double_t meanPullDiffFit=his->GetMean();
2550 0 : Double_t rmsPullDiffFit=his->GetRMS();
2551 0 : treeDist->GetHistogram()->Write("hisQAPullDiffFit");
2552 :
2553 0 : canvasQA->SaveAs((TString::Format("%sFit_sec%d_%d_theta%d_%dQA.png",inputTree,Int_t(sector0),Int_t(sector1),Int_t(theta0),Int_t(theta1)).Data()));
2554 :
2555 0 : TCanvas *canvasQAFit = new TCanvas("canvasQAFit","canvasQAFit",1200,1000);
2556 0 : canvasQAFit->SetRightMargin(0.01);
2557 0 : canvasQAFit->Divide(1,5,0,0);
2558 0 : treeDist->SetMarkerStyle(25);
2559 0 : treeDist->SetMarkerSize(0.5);
2560 0 : TCut defaultCut="abs(qptCenter)<0.1";
2561 : //
2562 : {
2563 0 : canvasQAFit->cd(1)->SetRightMargin(0.1);
2564 0 : treeDist->Draw("delta:sectorCenter:RCenter",defaultCut+"abs(abs(kZCenter)-0.1)<0.06","colz");
2565 0 : canvasQAFit->cd(2)->SetRightMargin(0.1);
2566 0 : treeDist->Draw("delta-delta_Fit0:sectorCenter:RCenter",defaultCut+"abs(abs(kZCenter)-0.1)<0.06","colz");
2567 0 : canvasQAFit->cd(3)->SetRightMargin(0.1);
2568 0 : treeDist->Draw("delta-delta_Fit1:sectorCenter:RCenter",defaultCut+"abs(abs(kZCenter)-0.1)<0.06","colz");
2569 0 : canvasQAFit->cd(4)->SetRightMargin(0.1);
2570 0 : treeDist->Draw("delta-delta_Fit2:sectorCenter:RCenter",defaultCut+"abs(abs(kZCenter)-0.1)<0.06","colz");
2571 0 : canvasQAFit->cd(5)->SetRightMargin(0.1);
2572 0 : treeDist->Draw("delta_Fit0-delta_Fit2:sectorCenter:RCenter",defaultCut+"abs(abs(kZCenter)-0.1)<0.06","colz");
2573 : }
2574 0 : canvasQAFit->SaveAs((TString::Format("%sFit_sec%d_%d_theta%d_%dQAFit.png",inputTree,Int_t(sector0),Int_t(sector1),Int_t(theta0),Int_t(theta1)).Data()));
2575 0 : TObjString input=inputTree;
2576 : //
2577 : // 5.) Export trending variables - used for validation of the fit
2578 : //
2579 0 : treeMeta->Draw("runNumber:selHis:fillCounter:clusterCounter:meanTime:ntracksUsed:startTime:stopTime","1","goffpara");
2580 0 : (*pcstream)<<"NDFitTrending"<< // cp of subset of info from meta data (rest accessible in the metadata tree also avaialble in file)
2581 0 : "run="<<treeMeta->GetVal(0)[0]<<
2582 0 : "selHis="<<treeMeta->GetVal(1)[0]<<
2583 0 : "fillCounter="<<treeMeta->GetVal(2)[0]<<
2584 0 : "clusterCounter="<<treeMeta->GetVal(3)[0]<<
2585 0 : "meanTime="<<treeMeta->GetVal(4)[0]<<
2586 0 : "ntracksUsed="<<treeMeta->GetVal(5)[0]<<
2587 0 : "startTime="<<treeMeta->GetVal(6)[0]<<
2588 0 : "stopTime="<<treeMeta->GetVal(7)[0];
2589 0 : Int_t entriesCl= treeMeta->Draw("grNcl72.fY:grNcl73.fY:grNcl74.fY:grNcl75.fY","(grNcl72.fX>startTime&&grNcl72.fX<stopTime&&grNcl72.fY!=0)","goffpara");
2590 0 : TVectorF vecNcl(8);
2591 0 : if (entriesCl>0) {
2592 0 : for (Int_t icl=0; icl<4; icl++){
2593 0 : vecNcl[icl]=TMath::Median(entriesCl, treeMeta->GetVal(icl));
2594 : }
2595 0 : entriesCl= treeMeta->Draw("grNclUsed72.fY:grNclUsed73.fY:grNclUsed74.fY:grNclUsed75.fY","(grNcl72.fX>startTime&&grNcl72.fX<stopTime&&grNcl72.fY!=0)","goffpara");
2596 0 : for (Int_t icl=0; icl<4; icl++){
2597 0 : vecNcl[icl+4]=TMath::Median(entriesCl, treeMeta->GetVal(icl));
2598 : }
2599 0 : }
2600 0 : (*pcstream)<<"NDFitTrending"<< // cp of subset of info from meta data (rest accessible in the metadata tree also avaialble in file)
2601 0 : "vecNclCounter.="<<&vecNcl;
2602 : //
2603 0 : (*pcstream)<<"NDFitTrending"<<
2604 0 : "input.="<<&input<< // name of the input file
2605 0 : "inputTree="<<inputTree<< // name of the input tree
2606 0 : "runNumber="<<runNumber<< // run number ID
2607 0 : "sec0="<<sector0<< // range: sector0
2608 0 : "sec1="<<sector1<< // range: sector1
2609 0 : "theta0="<<theta0<< // range: theta0
2610 0 : "theta1="<<theta1<< // range: theta1
2611 : // // QA plots statistical info
2612 0 : "slopePols1.="<<&slopePols1<< // data:fit - pol1 fit for differnt Regularization
2613 0 : "slopePol1="<<slopePol1<<
2614 : //
2615 0 : "mean="<<mean<< // mean <value-fitND0>
2616 0 : "rms="<<rms<< // rms (value-fitND0>
2617 0 : "rmsG="<<rmsG<< // gaus fit rms (value-fitND0>
2618 0 : "meanPull="<<meanPull<< // normalized mean <value-fitND0>
2619 0 : "rmsPull="<<rmsPull<< // normalized error<value-fitND0>
2620 0 : "rmsPullG="<<rmsPull<< // gaus fit normalized error<value-fitND0>
2621 0 : "meanDiffFit="<<meanDiffFit<< //
2622 0 : "rmsDiffFit="<<rmsDiffFit<<
2623 0 : "meanPullDiffFit="<<meanPullDiffFit<<
2624 0 : "rmsPullDiffFit="<<rmsPullDiffFit<<
2625 0 : "vecRMS1.="<<&vecRMS1; // rms of the diffs of different estimators (Kernles, Regularization)
2626 0 : for (Int_t iq=0; iq<18; iq++){
2627 0 : (*pcstream)<<"NDFitTrending"<<
2628 0 : TString::Format("grQuantiles%d.=",iq).Data()<<grQuantiles[iq];
2629 : }
2630 0 : (*pcstream)<<"NDFitTrending"<<"\n";
2631 0 : delete treeMeta;
2632 0 : delete pcstream;
2633 0 : delete pcstreamFit;
2634 0 : }
2635 :
2636 :
2637 : TTree* AliTPCcalibAlignInterpolation::LoadDistortionTrees(const char * maplist, Int_t cacheSize, Int_t markerStyle, Float_t markerSize){
2638 : //
2639 : // Load distortion trees specified in the maplist
2640 : // Loading distortion maps as a friend trees used for
2641 : // - correction for reference distortions (e.g map at low IR)
2642 : // - distortion maps correlation studies
2643 : // - distortion maps scaling fitting
2644 : //
2645 : // To obtain run number and TimeBin ID we assume naming convention as used in the calibration procedure.
2646 : // This naming convention is hardwired in the code.
2647 : //
2648 : TTree * treeReturn=0;
2649 : TTree * tree=0;
2650 0 : TObjArray* array = TString(gSystem->GetFromPipe(TString::Format("cat %s",maplist).Data())).Tokenize("\n");
2651 0 : TObjArray*arrayOK=new TObjArray(array->GetEntries());
2652 0 : for (Int_t i=0; i<array->GetEntries(); i++){
2653 0 : printf("%s\n",array->At(i)->GetName());
2654 0 : TString fname(array->At(i)->GetName());
2655 0 : Int_t index=fname.Index("/000");
2656 0 : TString runName(&(fname[index+1]),9);
2657 0 : index=fname.Index("/Time");
2658 0 : TString timeString(&(fname[index+9]),4);
2659 0 : if (TString(array->At(i)->GetName()).Contains("_0.root")){
2660 0 : tree = AliTPCcalibAlignInterpolation::AddFriendDistortionTree(treeReturn,array->At(i)->GetName(),"deltaRPhiTPCITSDist",runName+"_"+timeString+"ITSY");
2661 0 : }
2662 0 : if (TString(array->At(i)->GetName()).Contains("_1.root")){
2663 0 : tree = AliTPCcalibAlignInterpolation::AddFriendDistortionTree(treeReturn,array->At(i)->GetName(),"deltaRPhiTPCITSTRDDist",runName+"_"+timeString+"TRDY");
2664 0 : }
2665 0 : if (TString(array->At(i)->GetName()).Contains("_2.root")){
2666 0 : tree = AliTPCcalibAlignInterpolation::AddFriendDistortionTree(treeReturn,array->At(i)->GetName(),"deltaRPhiTPCITSTOFDist",runName+"_"+timeString+"TOFY");
2667 0 : }
2668 0 : if (TString(array->At(i)->GetName()).Contains("_3.root")){
2669 0 : tree = AliTPCcalibAlignInterpolation::AddFriendDistortionTree(treeReturn,array->At(i)->GetName(),"deltaZTPCITSDist",runName+"_"+timeString+"ITSZ");
2670 0 : }
2671 0 : if (TString(array->At(i)->GetName()).Contains("_4.root")){
2672 0 : tree = AliTPCcalibAlignInterpolation::AddFriendDistortionTree(treeReturn,array->At(i)->GetName(),"deltaZTPCITSTRDDist",runName+"_"+timeString+"TRDZ");
2673 0 : }
2674 0 : if (TString(array->At(i)->GetName()).Contains("_5.root")){
2675 0 : tree = AliTPCcalibAlignInterpolation::AddFriendDistortionTree(treeReturn,array->At(i)->GetName(),"deltaZTPCITSTOFDist",runName+"_"+timeString+"TOFZ");
2676 0 : }
2677 0 : if (tree) {
2678 0 : arrayOK->AddLast(array->At(i));
2679 : treeReturn=tree;
2680 0 : }
2681 0 : }
2682 0 : treeReturn->SetCacheSize(cacheSize);
2683 0 : treeReturn->SetMarkerStyle(markerStyle);
2684 0 : treeReturn->SetMarkerSize(markerSize);
2685 0 : return treeReturn;
2686 0 : }
2687 :
2688 :
2689 : Bool_t AliTPCcalibAlignInterpolation::LoadNDLocalFit(TTree * tree, const char *chTree){
2690 : ///
2691 : /// Load ND local fits. We assume data are organized in particular directory structure, in directories together with input maps
2692 : ///
2693 : ///
2694 : /// Input:
2695 : /// \param TTree * tree - input tree with distortion maps and "friend trees" per time bins
2696 : /// \param const char *chTree - name of the "distortion branch"
2697 :
2698 0 : if ( tree->GetListOfFriends()->FindObject(chTree)==NULL){
2699 0 : ::Error("AliTPCcalibAlignInterpolation::LoadNDLocal","Tree %s does not exist",chTree);
2700 0 : return kFALSE;
2701 : }
2702 0 : TString floc = tree->GetListOfFriends()->FindObject(chTree)->GetTitle();
2703 0 : TString fdir = gSystem->DirName(floc);
2704 : //
2705 0 : TObjArray *ndFileList = ( gSystem->GetFromPipe(TString::Format("ls %s/delta*root",fdir.Data()).Data())).Tokenize("\n");
2706 :
2707 0 : if (ndFileList->GetEntries()==0){
2708 0 : ::Error(" AliTPCcalibAlignInterpolation::LoadNDLocal","File with NDLocal %s",chTree);
2709 0 : return kFALSE;
2710 : }
2711 0 : for (Int_t ind=0; ind<ndFileList->GetEntries(); ind++){
2712 : //
2713 0 : TString fname=ndFileList->At(ind)->GetName();
2714 0 : TFile * f = TFile::Open(fname.Data());
2715 0 : TList * arrKey = f->GetListOfKeys();
2716 0 : for (Int_t ikey=0; ikey<arrKey->GetEntries(); ikey++){
2717 0 : TString keyName=arrKey->At(ikey)->GetName();
2718 0 : if (keyName.Contains("delta")==0) continue;
2719 0 : TObject * o = f->Get(keyName);
2720 0 : AliNDLocalRegression * reg = dynamic_cast<AliNDLocalRegression*>(o);
2721 0 : if (reg==NULL){
2722 0 : delete o;
2723 0 : continue;
2724 : }
2725 0 : TString aliasName = TString::Format("%s.%s",chTree,reg->GetName());
2726 0 : aliasName.ReplaceAll("-","Min");
2727 0 : ::Info("AliTPCcalibAlignInterpolation::LoadNDLocal","Loaded ND local regression %s/%s as alias %s", fname.Data(), reg->GetName(),aliasName.Data());
2728 0 : reg->SetName(aliasName);
2729 0 : Int_t hashIndex=reg->GetVisualCorrectionIndex();
2730 0 : reg->AddVisualCorrection(reg, hashIndex);
2731 0 : tree->SetAlias(aliasName, TString::Format("AliNDLocalRegression::GetCorrND(%d,RCenter,sectorCenter,kZCenter,qptCenter+0)",hashIndex).Data());
2732 0 : tree->SetAlias(aliasName+"L", TString::Format("AliNDLocalRegression::GetCorrND(%d,RCenter,sectorCenter,kZCenter,qptCenter+0)-(AliNDLocalRegression::GetCorrND(%d,RCenter,sectorCenter,kZCenter,qptCenter-2)+AliNDLocalRegression::GetCorrND(%d,RCenter,sectorCenter,kZCenter,qptCenter+2))*0.5",hashIndex,hashIndex,hashIndex).Data());
2733 0 : tree->SetAlias(aliasName+"M", TString::Format("(AliNDLocalRegression::GetCorrND(%d,RCenter,sectorCenter,kZCenter,qptCenter-2)+AliNDLocalRegression::GetCorrND(%d,RCenter,sectorCenter,kZCenter,qptCenter+2))*0.5",hashIndex, hashIndex).Data());
2734 0 : }
2735 0 : }
2736 0 : return kTRUE;
2737 :
2738 0 : }
2739 :
2740 :
2741 : void AliTPCcalibAlignInterpolation::DrawMapEstimatorComparison(TTree * tree, const char* chtree, Float_t radius, Float_t kZ, TCut & selection, const char *figType){
2742 : // Predefined plot:
2743 : // Draw distortion map comparison
2744 : // Compare median and LTM estimator of the mean value of distortion in the bin
2745 : //
2746 :
2747 : // Example usage:
2748 : /*
2749 : const char* chtree="000244918_his1TRDY"; // Low IR one polarity
2750 : const char* chtree="000246391_his1TRDY"; // Low IR another polarity
2751 : Float_t radius=100;
2752 : Float_t kZ=0.1;
2753 : figType="png";
2754 : AliTPCcalibAlignInterpolation::DrawMapEstimatorComparison(tree, chtree, radius,kZ,"abs(qptCenter)<0.1",figType);
2755 : */
2756 0 : if (!tree) {
2757 0 : ::Error("DrawEstimatorComparison","Tree not available");
2758 0 : return;
2759 : }
2760 0 : if (chtree && tree->GetListOfFriends()->FindObject(chtree)==NULL){
2761 0 : ::Error("DrawEstimatorComparison","Ttree %s not available",chtree);
2762 0 : return;
2763 : }
2764 0 : gStyle->SetLabelSize(0.07,"XYZ");
2765 0 : gStyle->SetTitleSize(0.06,"XYZ");
2766 0 : gStyle->SetTitleOffset(1.0,"X");
2767 0 : gStyle->SetTitleOffset(0.6,"Y");
2768 0 : gStyle->SetTitleOffset(0.4,"Z");
2769 0 : gStyle->SetOptTitle(1);
2770 :
2771 :
2772 0 : TCanvas * canvasC = new TCanvas("canvasC","canvasC",1400,1000);
2773 : TPad * pad=0;
2774 0 : TPad *pads[4]={0};
2775 0 : for (Int_t ipad=0; ipad<4; ipad++){
2776 0 : pads[ipad] = new TPad("pad1","This is pad1",0.00,ipad/4., 1,(ipad+1.)/4.);
2777 0 : pads[ipad]->SetBottomMargin(0);
2778 0 : pads[ipad]->SetTopMargin(0);
2779 : }
2780 0 : pads[0]->SetBottomMargin(0.15);
2781 0 : pads[3]->SetTopMargin(0.05);
2782 0 : for (Int_t ipad=0; ipad<4; ipad++){
2783 0 : pads[ipad]->Draw();
2784 0 : pads[ipad]->SetGrid(1,1);
2785 : }
2786 :
2787 0 : if (chtree){
2788 0 : pads[3]->cd();
2789 0 : tree->Draw(TString::Format("%s.binMedian:sectorCenter:RCenter",chtree).Data(),selection+TString::Format("abs(kZCenter+(%2.2f))<0.06&&abs(RCenter-%2.2f)<20&&%s.rms>0",kZ,radius,chtree).Data(),"colz");
2790 0 : pads[2]->cd();
2791 0 : tree->Draw(TString::Format("%s.vecLTM.fElements[1]:sectorCenter:RCenter",chtree).Data(),selection+TString::Format("abs(kZCenter+(%2.2f))<0.06&&abs(RCenter-%2.2f)<20&&%s.rms>0",kZ,radius,chtree).Data(),"colz");
2792 0 : pads[1]->cd();
2793 0 : tree->Draw(TString::Format("%s.meanG:sectorCenter:RCenter",chtree).Data(),selection+TString::Format("abs(kZCenter+(%2.2f))<0.06&&abs(RCenter-%2.2f)<20&&%s.rms>0&&abs(%s.meanG-%s.binMedian)<1",kZ,radius,chtree,chtree,chtree).Data(),"colz");
2794 0 : pads[0]->cd();
2795 0 : tree->Draw(TString::Format("%s.vecLTM.fElements[1]-%s.binMedian:sectorCenter:RCenter",chtree,chtree).Data(),selection+TString::Format("abs(kZCenter+(%2.2f))<0.06&&abs(RCenter-%2.2f)<20&&%s.rms>0",kZ,radius,chtree).Data(),"colz");
2796 0 : }else{
2797 : pads[3]->cd();
2798 0 : tree->Draw("binMedian:sectorCenter:RCenter",selection+TString::Format("abs(kZCenter+(%2.2f))<0.06&&abs(RCenter-%2.2f)<20&&rms>0",kZ,radius).Data(),"colz");
2799 0 : pads[2]->cd();
2800 0 : tree->Draw("vecLTM.fElements[1]:sectorCenter:RCenter",selection+TString::Format("abs(kZCenter+(%2.2f))<0.06&&abs(RCenter-%2.2f)<20&&rms>0",kZ,radius).Data(),"colz");
2801 0 : pads[1]->cd();
2802 0 : tree->Draw("meanG:sectorCenter:RCenter",selection+TString::Format("abs(kZCenter+(%2.2f))<0.06&&abs(RCenter-%2.2f)<20&&rms>0&&abs(meanG-binMedian)<1",kZ,radius).Data(),"colz");
2803 0 : pads[0]->cd();
2804 0 : tree->Draw("binMedian-vecLTM.fElements[1]:sectorCenter:RCenter",selection+TString::Format("abs(kZCenter+(%2.2f))<0.06&&abs(RCenter-%2.2f)<20&&rms>0",kZ,radius).Data(),"colz");
2805 : }
2806 :
2807 0 : if (figType) canvasC->SaveAs(TString::Format("distortionMap_%s_R%2.0f_kZ%2.2f.%s",chtree,radius,kZ,figType).Data());
2808 0 : }
2809 :
2810 :
2811 :
2812 : Bool_t AliTPCcalibAlignInterpolation::DrawScalingComparison(TTree * tree, const char* chRef, const char *chBin0, const char *chBin1, Float_t R0, Float_t R1, Float_t kZ, const char *figType){
2813 : ///
2814 : /// Make predefined plot:
2815 : /// Draw distortion maps comparison and save figure (figType) in current working directory.
2816 : /// Fig 1.): Map run/timeBin chBin0 corrected for reference map (chRef)offset
2817 : /// Fig 2.): Map run/timeBin chBin1 corrected for reference map (chRef)offset
2818 : /// Fig 3.): Map (chBin1-chRef)-scale*(chBin0-chRef)
2819 : ///
2820 : /// Input:
2821 : /// \param TTree * tree - input tree with distortion maps and "friend trees" per time bins
2822 : /// \param chRef - reference distortion map
2823 : /// \param chBin0 - run or time bin shown in firs row
2824 : /// \param chBin1 - run or time bin
2825 : /// \return kTRUE if comparison of distortion map possible and figure saved
2826 : /// TString::Format("distortionMapScaling_%s_%s_%s_R0%2.0f_R1%2.0f_kZ%2.2f.%s",chBin0,chBin1,chRef,R0,R1,kZ,figType).Data();
2827 : /// Example usage:
2828 : /// To be used for systematic studies of TPC distortion.
2829 : /// E.g:
2830 : /*
2831 : TTree * tree = AliTPCcalibAlignInterpolation::LoadDistortionTrees("map.list");
2832 : const char* chRef="000246391_his1TRDY"; // Low IR refernce run for given B-field polarity
2833 : const char* chBin0="000246980_his1TRDY"; // time bin at the beginnig of the fill
2834 : const char* chBin1="000246994_his1TRDY"; // time bin at the end of the fill
2835 : AliTPCcalibAlignInterpolation::DrawScalingComparison(tree, chRef, chBin0, chBin1, 85, 130, 0.1, "png");
2836 : */
2837 : //
2838 : // 0.) Check variables
2839 : //
2840 0 : TCut defaultCut="abs(qptCenter)<0.1";
2841 0 : if (tree==NULL) {
2842 0 : ::Error("AliTPCcalibAlignInterpolation::DrawScalingComparison","Tree not set");
2843 0 : return kFALSE;
2844 : }
2845 0 : Int_t tentries[3]={0};
2846 0 : const char *vars[3]={chRef, chBin0, chBin1};
2847 0 : for (Int_t i=0; i<3;i++){
2848 0 : tentries[i]=tree->Draw(TString::Format("%s.binMedian",vars[i]),defaultCut,"goff", 10000);
2849 0 : if (tentries[i]<=0){
2850 0 : ::Error("AliTPCcalibAlignInterpolation::DrawScalingComparison","Expression %s or tree %s not valid ", vars[i], tree->GetName());
2851 0 : return kFALSE;
2852 : }
2853 : }
2854 : //
2855 : // 1.) get scaling factor using A side scaling (OROC scaling on C side more problematic)
2856 : //
2857 0 : Int_t entries = tree->Draw(TString::Format("%s.binMedian-%s.binMedian:%s.binMedian-%s.binMedian", chBin0,chRef, chBin1,chRef).Data(),defaultCut+"kZCenter>0","goff");
2858 0 : TGraph * gr0 = new TGraph(entries,tree->GetV1(),tree->GetV2());
2859 0 : TGraph * gr1 = new TGraph(entries,tree->GetV2(),tree->GetV1());
2860 0 : gr0->Fit("pol1");
2861 0 : gr1->Fit("pol1");
2862 0 : Double_t slope=(gr0->GetFunction("pol1")->GetParameter(1)+ 1/gr1->GetFunction("pol1")->GetParameter(1))*0.5;
2863 0 : tree->SetAlias("norm0",TString::Format("%f*(%s.binMedian-%s.binMedian)",slope,chBin0,chRef).Data());
2864 : //
2865 : // 2. Make plots
2866 : //
2867 0 : gStyle->SetLabelSize(0.07,"XYZ");
2868 0 : gStyle->SetLabelSize(0.03,"Y");
2869 0 : gStyle->SetTitleSize(0.06,"XYZ");
2870 0 : gStyle->SetTitleSize(0.04,"Y");
2871 0 : gStyle->SetTitleOffset(1.0,"X");
2872 0 : gStyle->SetTitleOffset(0.5,"Y");
2873 0 : gStyle->SetTitleOffset(0.4,"Z");
2874 0 : gStyle->SetOptTitle(1);
2875 : //
2876 0 : TCanvas * canvasC = new TCanvas("canvasC","canvasC",1400,1000);
2877 : TPad * pad=0;
2878 0 : TPad *pad3 = new TPad("pad1","This is pad1",0.00,0.0, 1,0.33);
2879 0 : TPad *pad2 = new TPad("pad2","This is pad2",0.00,0.33, 1,0.66);
2880 0 : TPad *pad1 = new TPad("pad3","This is pad3",0.00,0.66, 1,1);
2881 0 : pad1->SetBottomMargin(0);
2882 0 : pad2->SetBottomMargin(0);
2883 0 : pad3->SetBottomMargin(0.15);
2884 0 : pad2->SetTopMargin(0);
2885 0 : pad3->SetTopMargin(0);
2886 0 : pad1->Draw();
2887 0 : pad2->Draw();
2888 0 : pad3->Draw();
2889 0 : pad1->SetGrid(1,1);
2890 0 : pad2->SetGrid(1,1);
2891 0 : pad3->SetGrid(1,1);
2892 0 : TCut cutAccept=defaultCut+TString::Format("abs(kZCenter+(%2.2f))<0.06&&RCenter>%2.0f&&RCenter<%2.0f&&%s.rms>0",kZ,R0,R1,chRef).Data();
2893 : Int_t isOK=0;
2894 : const Int_t kMinEntries=100;
2895 : {
2896 0 : pad1->cd();
2897 0 : isOK+=tree->Draw(TString::Format("(%s.binMedian-%s.binMedian):sectorCenter:RCenter", chBin0,chRef),cutAccept,"colz")>kMinEntries;
2898 0 : pad2->cd();
2899 0 : isOK+=tree->Draw(TString::Format("(%s.binMedian-%s.binMedian):sectorCenter:RCenter", chBin1,chRef),cutAccept,"colz")>kMinEntries;
2900 0 : pad3->cd();
2901 0 : isOK+=tree->Draw(TString::Format("(%s.binMedian-%s.binMedian)-norm0:sectorCenter:RCenter", chBin1,chRef),cutAccept,"colz")>kMinEntries;
2902 : }
2903 0 : if (isOK<3){
2904 0 : ::Error("AliTPCcalibAlignInterpolation::DrawScalingComparison","Not enough points in selected region R<%2.0f,%2.0f>", R0,R1);
2905 0 : return kFALSE;
2906 : }
2907 0 : if (figType) canvasC->SaveAs(TString::Format("distortionMapScaling_%s_%s_%s_R0%2.0f_R1%2.0f_kZ%2.2f.%s",chBin0,chBin1,chRef,R0,R1,kZ,figType).Data());
2908 0 : return kTRUE;
2909 0 : }
2910 :
2911 :
2912 :
2913 : //_______________________________________________
2914 : double AliTPCcalibAlignInterpolation::GetTgPhi(double x, double y2x, double q2p, double b)
2915 : {
2916 : // calculate tangent of primary track at any frame at given x,y
2917 0 : double y = y2x*x;
2918 0 : double c = q2p*b*(-0.299792458e-3);
2919 0 : if (TMath::Abs(c)<1e-9) return y2x;
2920 0 : double r2 = y*y+x*x;
2921 0 : double det = 4./r2 - c*c;
2922 : double snp;
2923 0 : if (det<0) {
2924 0 : snp = TMath::Sign(-0.8,c);
2925 : //printf("track of q2p=%f cannot reach x:%f y:%f, define snp as %f \n",q2p,x,y,snp);
2926 0 : }
2927 : else {
2928 0 : snp = 0.5*(y*TMath::Sqrt(det)-c*x); // snp at vertex
2929 0 : snp += x*c; // snp at x,y
2930 : }
2931 0 : return snp/TMath::Sqrt((1-snp)*(1+snp));
2932 0 : }
2933 :
2934 : //______________________________________________
2935 : void AliTPCcalibAlignInterpolation::FixAlignmentBug(int sect, float q2pt, float bz,
2936 : float& alp, float& x, float &z, float &deltaY, float &deltaZ)
2937 : {
2938 : // fix alignment bug: https://alice.its.cern.ch/jira/browse/ATO-339?focusedCommentId=170850&page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel#comment-170850
2939 : //
2940 : // NOTE: deltaZ in the buggy code is calculated as Ztrack_with_bug - Zcluster_w/o_bug
2941 : static TGeoHMatrix *mCache[72] = {0};
2942 0 : if (sect<0||sect>=72) {
2943 0 : AliErrorClassF("Invalid sector %d",sect);
2944 0 : return;
2945 : }
2946 0 : int lr = sect/36 ? (AliGeomManager::kTPC2) : (AliGeomManager::kTPC1);
2947 0 : TGeoHMatrix* mgt = mCache[sect];
2948 0 : if (!mgt) {
2949 0 : int volID = AliGeomManager::LayerToVolUIDSafe(lr,sect%36);
2950 0 : mgt = new TGeoHMatrix(*AliGeomManager::GetTracking2LocalMatrix(volID));
2951 0 : mgt->MultiplyLeft(AliGeomManager::GetMatrix(volID));
2952 0 : mCache[sect] = mgt;
2953 0 : printf("Caching matrix for sector %d\n",sect);
2954 0 : }
2955 0 : double alpSect = ((sect%18)+0.5)*20.*TMath::DegToRad();
2956 :
2957 : // cluster in its proper alpha frame with alignment bug
2958 0 : double xyzClUse[3] = {x,0,z}; // this is what we read from the residual tree
2959 0 : double xyzTrUse[3] = {x, deltaY, z}; // track in bad cluster frame
2960 : //
2961 : // recover cluster Z position by adding deltaZ
2962 0 : double zClSave = xyzClUse[2] -= deltaZ; // here the cluster is not affected by Z alignment component of the bug!
2963 0 : static AliExternalTrackParam trDummy;
2964 0 : trDummy.Local2GlobalPosition(xyzClUse,alp); // misaligned cluster in global frame
2965 0 : double xyz0[3]={xyzClUse[0],xyzClUse[1],xyzClUse[2]};
2966 0 : mgt->MasterToLocal(xyz0,xyzClUse);
2967 : // we got ideal cluster in the sector tracking frame, but now the Z is wrong, since it was not affected by the bug!!!
2968 : //
2969 0 : xyzClUse[2] = zClSave;
2970 :
2971 : // go to ideal cluster frame
2972 0 : trDummy.Local2GlobalPosition(xyzClUse,alpSect); // ideal global
2973 0 : double alpFix = TMath::ATan2(xyzClUse[1],xyzClUse[0]); // fixed cluster phi
2974 0 : trDummy.Global2LocalPosition(xyzClUse,alpFix); // fixed cluster in in its frame
2975 : //
2976 0 : trDummy.Local2GlobalPosition(xyzTrUse,alp); // track in global frame
2977 0 : trDummy.Global2LocalPosition(xyzTrUse,alpFix); // track in cluster frame
2978 0 : alp = alpFix;
2979 : //
2980 0 : double dx = xyzTrUse[0] - xyzClUse[0]; // x might not be the same after alignment fix
2981 : // deduce track slopes assuming it comes from the vertex
2982 0 : double tgphi = GetTgPhi(xyzClUse[0],xyzTrUse[1]/xyzClUse[0],q2pt,bz);
2983 0 : xyzTrUse[1] -= dx*tgphi;
2984 0 : xyzTrUse[2] -= dx*xyzClUse[2]/xyzClUse[0]; // z2x
2985 : //
2986 0 : x = xyzClUse[0];
2987 0 : z = xyzTrUse[2]; // we still use track Z as a reference ...
2988 0 : deltaY = xyzTrUse[1]-xyzClUse[1];
2989 0 : deltaZ = xyzTrUse[2]-xyzClUse[2];
2990 : //
2991 0 : }
2992 :
2993 :
2994 : void AliTPCcalibAlignInterpolation::MakeVDriftOCDB(const char *inputFile, Int_t run, TString targetOCDBstorage, const char * testDiffCDB/*=0*/){
2995 : ///
2996 : /// Make OCDB entry using information using the vdrift calibration obtained in the AliTPCcalibAlignInterpolation::FitDrift function (suing ResidualTrees.root)
2997 : ///
2998 : ///
2999 : /*
3000 : char * inputFile= "fitDrift.root"
3001 : char * testDiffCDB="/cvmfs/alice.cern.ch/calibration/data/2015/OCDB/TPC/Calib/TimeDrift/Run244918_244918_v3_s0.root";
3002 : Int_t run=244918;
3003 :
3004 : .x $ALICE_PHYSICS/PWGPP/CalibMacros/CPass0/ConfigCalibTrain.C(run,"local:///cvmfs/alice.cern.ch/calibration/data/2015/OCDB/")
3005 : AliTPCcalibAlignInterpolation::MakeVDriftOCDB( inputFile,run,"", testDiffCDB)
3006 : */
3007 :
3008 : /*
3009 : The same calibration was done previoulsy using the AliTPCcalibTime class
3010 : Many graphs in previously done RUN1 calibration used for monitoring purposes and do not have equivalent in the new calibration
3011 : using ResidualTrees
3012 : Old CPass0 calibration () - 114 graphs exported and monitored later using trending (O (100 kBy) per run)
3013 : * 108 graphs for Align+drift = 9 params x 3 detectors x 4 version )
3014 : ** vdrift calibration (3 param) and alingment graphs (6)
3015 : ** each pair of detector ITS->TPC, TRD->TPC, TOF->TPC
3016 : ** Kalman filter storted in 4 version
3017 : *** raw input per time bin, kalman forward propagation, Kalaman back propagation, smoothed Kalman
3018 : ** usage
3019 : *** ITS->TPC smoothed version in case exist
3020 : *** TOF->TPC smoothed in case ITS->TPC not availale
3021 :
3022 : * 6 graphs for the Laser CE backup graphs
3023 : * backup of ClusterParam and RecoParam as it was used in the process of calibration
3024 : ** calibration depends
3025 :
3026 : New calibration input data ( O(3 MBy) for 2 hours measurement):
3027 :
3028 : Export variables:
3029 : only vdrift calibration graphs (_DELTAZ,DRIFTVD, TO and VDGY) will be exported
3030 : smoothed version equivalent (grTRDReg) and raw calibration equivalent (grTRDmed)
3031 : ITS - mean (TRD+TOF)*0.5
3032 : TRD - TRD
3033 : TOF - TOF
3034 : */
3035 0 : TObjArray * driftArray = new TObjArray();
3036 :
3037 : //
3038 : // 0. Initialize OCDB if not done before
3039 : //
3040 : AliTPCParam *param =0;
3041 0 : if (AliCDBManager::Instance()->GetDefaultStorage()!=NULL){
3042 : //
3043 : //
3044 0 : AliTPCClusterParam *clParam = AliTPCcalibDB::Instance()->GetClusterParam();
3045 0 : param= AliTPCcalibDB::Instance()->GetParameters();
3046 0 : TObjArray *recoParams = new TObjArray(4) ;
3047 0 : for (Int_t i=0;i<4;i++) recoParams->AddAt(AliTPCcalibDB::Instance()->GetRecoParam(i),i);
3048 0 : driftArray->AddLast(clParam);
3049 0 : driftArray->AddLast(recoParams);
3050 0 : ::Info("AliTPCcalibAlignInterpolation::MakeVDriftOCDB","Residual calibration used. We have to trust - your OCDB is correct, as we can nt check.");
3051 0 : }else{
3052 0 : ::Fatal("AliTPCcalibAlignInterpolation::MakeVDriftOCDB","Im sorry. OCDB has to be intilaized before. We need to get Parameters objects, which were used in previous calibration. In ideal case the OCDB entry to be setup to the entries, indeed used to produce Residual trees ... But this we can not CHECK");
3053 : }
3054 :
3055 : //
3056 : // 1. Read calibratio data from the tree
3057 : //
3058 0 : TVectorD *vdriftParam=0;
3059 0 : TGraphErrors *grTRDReg=0;
3060 0 : TGraphErrors *grTOFReg=0;
3061 0 : TGraphErrors *grTRDMed=0;
3062 0 : TGraphErrors *grTOFMed=0;
3063 0 : TFile *fdrift = TFile::Open(inputFile);
3064 0 : if (fdrift){
3065 0 : TTree * tree = (TTree*)fdrift->Get("fitTimeStat");
3066 0 : if (tree==NULL){
3067 0 : ::Fatal("MakeVDriftOCDB.LoadDriftCalibration FAILED", "tree fitTimeStat not avaliable in file fitDrift.root");
3068 0 : }else{
3069 0 : if (tree->GetBranch("grTRDReg.")) tree->SetBranchAddress("grTRDReg.",&grTRDReg);
3070 0 : if (tree->GetBranch("grTOFReg.")) tree->SetBranchAddress("grTOFReg.",&grTOFReg);
3071 0 : if (tree->GetBranch("grTRDMed.")) tree->SetBranchAddress("grTRDMed.",&grTRDMed);
3072 0 : if (tree->GetBranch("grTOFMed.")) tree->SetBranchAddress("grTOFMed.",&grTOFMed);
3073 0 : if (tree->GetBranch("paramRobust.")) tree->SetBranchAddress("paramRobust.",&vdriftParam);
3074 0 : tree->GetEntry(0);
3075 : }
3076 0 : }else{
3077 0 : ::Fatal("MakeVDriftOCDB.LoadDriftCalibration FAILED", "Input file %s not accesible", inputFile);
3078 : }
3079 0 : if (grTRDReg==NULL && grTOFReg==NULL){
3080 0 : ::Fatal("MakeVDriftOCDB.LoadDriftCalibration FAILED", "drift calibration enot present");
3081 0 : }
3082 : //
3083 : // 2.) Transform v drift calibration object into the format used in previous calibration
3084 : //
3085 : TGraphErrors *graphDELTAZ=0;
3086 : TGraphErrors *graphT0=0;
3087 : TGraphErrors *graphVDGY=0;
3088 : TGraphErrors *graphDRIFTVD=0;
3089 : TGraphErrors *graph=0;
3090 :
3091 0 : const char *chDet[3]={"ITS","TRD","TOF"};
3092 0 : Double_t atime[2]={0,0};
3093 0 : Double_t deltaZ[2]={0,0};
3094 0 : Double_t t0[2]={0,0};
3095 0 : Double_t vdgy[2]={0,0};
3096 0 : graphDRIFTVD=(grTRDMed!=NULL) ? grTRDMed:grTOFMed;
3097 : //
3098 0 : atime[0]=graphDRIFTVD->GetX()[0];
3099 0 : atime[1]=graphDRIFTVD->GetX()[graphDRIFTVD->GetN()-1];
3100 0 : for (Int_t ipoint=0; ipoint<=1; ipoint++){
3101 0 : deltaZ[ipoint]=-(*vdriftParam)[1]; // unit OK
3102 0 : vdgy[ipoint]=-(*vdriftParam)[2]; // units OK
3103 0 : t0[ipoint]=-(*vdriftParam)[0]/(1+(*vdriftParam)[3]); // t0 to be normalized to the ms
3104 0 : t0[ipoint]/=(param->GetDriftV()/1000000.);
3105 : }
3106 :
3107 :
3108 0 : Double_t vdrifCorrRun=1./(1.+(*vdriftParam)[3]);
3109 :
3110 : //
3111 0 : for (Int_t idet=0; idet<3; idet++){
3112 0 : for (Int_t itype=0; itype<2; itype++){ //0. original version 1.)smoothed version
3113 0 : TString grPrefix="ALIGN_";
3114 0 : grPrefix+=chDet[idet];
3115 0 : if (itype==0) grPrefix+="_TPC_";
3116 0 : if (itype==1) grPrefix+="B_TPC_";
3117 : //
3118 0 : graphDELTAZ=new TGraphErrors(2, atime, deltaZ);
3119 0 : graphDELTAZ->SetName(grPrefix+"DELTAZ");
3120 0 : driftArray->AddLast(graphDELTAZ);
3121 0 : graphT0=new TGraphErrors(2, atime, t0);
3122 0 : graphT0->SetName(grPrefix+"T0");
3123 0 : driftArray->AddLast(graphT0);
3124 0 : graphVDGY=new TGraphErrors(2, atime, vdgy);
3125 0 : graphVDGY->SetName(grPrefix+"VDGY");
3126 0 : driftArray->AddLast(graphVDGY);
3127 : //
3128 : // drift velocity
3129 : //graph=(grTRDMed!=NULL) ? new TGraphErrors(*grTRDMed): new TGraphErrors(*grTOFMed); // normal constructor give us seg fault
3130 0 : graph=(TGraphErrors*)((grTRDMed!=NULL) ? grTRDMed->Clone(): grTOFMed->Clone()); // normal constructor give us 0
3131 0 : Int_t npoints = graph->GetN();
3132 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++){
3133 0 : graph->GetY()[ipoint]=-1;
3134 0 : if (idet==1 && grTRDMed) { // TRD
3135 0 : if (itype==0) graph->GetY()[ipoint]=(vdrifCorrRun*(1-grTRDMed->GetY()[ipoint]))-1.;
3136 0 : if (itype==1) graph->GetY()[ipoint]=(vdrifCorrRun*(1-grTRDReg->GetY()[ipoint]))-1.;
3137 : }
3138 0 : if (idet==2 && grTOFMed) { // TOF
3139 0 : if (itype==0) graph->GetY()[ipoint]=(vdrifCorrRun*(1-grTOFMed->GetY()[ipoint]))-1.;
3140 0 : if (itype==1) graph->GetY()[ipoint]=(vdrifCorrRun*(1-grTOFReg->GetY()[ipoint]))-1.;
3141 : }
3142 0 : if (idet==0 && (grTRDMed&&grTOFMed)) { // (TRD+TOF)*0.5
3143 0 : if (itype==0) graph->GetY()[ipoint]=(vdrifCorrRun*(1-(grTRDMed->GetY()[ipoint]+grTOFMed->GetY()[ipoint])*0.5))-1.;
3144 0 : if (itype==1) graph->GetY()[ipoint]=(vdrifCorrRun*(1-(grTRDReg->GetY()[ipoint]+grTOFReg->GetY()[ipoint])*0.5))-1.;
3145 : }
3146 : }
3147 0 : graph->SetName(grPrefix+"DRIFTVD");
3148 0 : driftArray->AddLast(graph);
3149 0 : }
3150 : }
3151 : //
3152 : // 3. Store selected graphs in OCDB
3153 : //
3154 : AliCDBStorage* targetStorage = 0x0;
3155 0 : if (targetOCDBstorage.Length()==0) {
3156 0 : targetOCDBstorage+="local://"+gSystem->GetFromPipe("pwd")+"/OCDB";
3157 0 : targetStorage = AliCDBManager::Instance()->GetStorage(targetOCDBstorage.Data());
3158 0 : }
3159 0 : else if (targetOCDBstorage.CompareTo("same",TString::kIgnoreCase) == 0 ){
3160 0 : targetStorage = AliCDBManager::Instance()->GetDefaultStorage();
3161 0 : }
3162 : else {
3163 0 : targetStorage = AliCDBManager::Instance()->GetStorage(targetOCDBstorage.Data());
3164 : }
3165 :
3166 0 : AliCDBMetaData* metaData = new AliCDBMetaData;
3167 0 : metaData->SetObjectClassName("TObjArray");
3168 0 : metaData->SetResponsible("Marian Ivanov");
3169 0 : metaData->SetBeamPeriod(1);
3170 0 : metaData->SetAliRootVersion(">v5-07-20"); //AliRoot version
3171 0 : metaData->SetComment("AliTPCcalibAlignInterpolation Calibration of the time dependence of the drift velocity using Residual trees");
3172 0 : AliCDBId id1("TPC/Calib/TimeDrift", run, run);
3173 :
3174 : //now the entry owns the metadata, but NOT the data
3175 0 : AliCDBEntry *driftCDBentry=new AliCDBEntry(driftArray,id1,metaData,kFALSE);
3176 0 : targetStorage->Put(driftCDBentry);
3177 : //
3178 : //
3179 : // 4. Make diff to reference CDB if sepecfied
3180 : //
3181 0 : if (testDiffCDB){
3182 0 : TFile *f = TFile::Open(testDiffCDB);
3183 0 : AliCDBEntry * entry=(AliCDBEntry*)f->Get("AliCDBEntry");
3184 0 : TObjArray * array = (TObjArray*)entry->GetObject();
3185 : //array->ls();
3186 0 : TGraphErrors * grRef= (TGraphErrors * )array->FindObject("ALIGN_TRDB_TPC_DRIFTVD");
3187 0 : TGraphErrors * gr= (TGraphErrors * )driftArray->FindObject("ALIGN_TRDB_TPC_DRIFTVD");
3188 0 : grRef->SetMarkerStyle(25);
3189 0 : gr->SetMarkerStyle(21);
3190 0 : Double_t minimum=TMath::Min(grRef->GetMinimum(), gr->GetMinimum());
3191 0 : Double_t maximum=TMath::Max(grRef->GetMaximum(), gr->GetMaximum());
3192 0 : grRef->SetMinimum(minimum);
3193 0 : grRef->SetMaximum(maximum);
3194 :
3195 0 : grRef->Draw("ap");
3196 0 : gr->Draw("p");
3197 0 : }
3198 :
3199 0 : delete driftCDBentry;
3200 :
3201 :
3202 0 : }
3203 :
3204 :
3205 : Float_t AliTPCcalibAlignInterpolation::CalculateDistance(const TVectorF &track0, const TVectorF &track1, const TVectorF &vecSec, TVectorF &vecDelta, Int_t npValid, Float_t &rmsTrack, Float_t &rmsCluster, Float_t lpNorm){
3206 : /// Calculate track and cluster RMS distance. RMS distance is than used in toder to
3207 : /// 1.) Calcute rms distance between the track0 and reference track1
3208 : /// 2.) Calculate rms distance between the cluster and local median cluster position in region +-2 padrows (wihing one sector)
3209 : /// and store delta in exported vecDelta array
3210 : /// \param[in] track0 vector: cluster-track residuals for track of interest
3211 : /// \param[in] track1 vector: cluster-track residuals for reference track
3212 : /// \param[in] vecSec vector: cluster sector index
3213 : /// \param[in] lpNorm the p value for the p-type norm (https://en.wikipedia.org/wiki/Norm_(mathematics))
3214 : /// \param[out] rmsTrack calculated 2 track distance (lpNorm)
3215 : /// \param[out] rmsCluster calculated cluster local mean distance (lpNorm)
3216 : ///
3217 : /// Return value - number of the clusters with too big RMS
3218 : ///
3219 : // In the following routine AliTPCcalibAlignInterpolation::FillHistogramsFromChain cut applied on quality of track and custer infromation
3220 : // Suggested cuts (based on the resolution parametrrization studies)
3221 : // - nOut <10 - less than 15 clusters rejected - kMaxSkippedCluster=10;
3222 : // - rmsTrack<2 cm kMaxRMSTrackCut=2.0;
3223 : // - rmsCluster<0.3 cm kMaxRMSClusterCut=0.3;
3224 : // - |vecDelta|<0.5 cm kMaxDeltaClusterCut=0.5;
3225 : //
3226 : // Parameters of algorithm for the moment set as a constant
3227 : const Int_t kMinFractionPoints=0.5;
3228 : const Float_t kMaxDist=20;
3229 : const Float_t kMaxDistTrack=5;
3230 0 : Float_t maxRMSTrack=rmsTrack;
3231 0 : Float_t maxRMSCluster=rmsCluster;
3232 :
3233 : //
3234 : Float_t distance2=0;
3235 : Int_t npoints=0;
3236 : // cache matrix pointers
3237 0 : const Float_t *delta0=track0.GetMatrixArray();
3238 0 : const Float_t *delta1=track1.GetMatrixArray();
3239 0 : const Float_t *sec=vecSec.GetMatrixArray();
3240 : //
3241 : // 1.) Calcute rms distance between the track0 and reference track1
3242 : //
3243 0 : for (Int_t irow=0; irow<npValid; irow++){
3244 0 : vecDelta[irow]=1000;
3245 0 : if (TMath::Abs(delta0[irow])==0) continue; //
3246 0 : if (TMath::Abs(delta1[irow])==0) continue; //
3247 0 : if (TMath::Abs(delta0[irow])>kMaxDist) continue; //
3248 0 : if (TMath::Abs(delta1[irow])>kMaxDist) continue; //
3249 0 : if (TMath::Abs(delta0[irow]-delta1[irow])>kMaxDistTrack) continue; //
3250 0 : npoints++;
3251 0 : distance2+=(delta0[irow]-delta1[irow])*(delta0[irow]-delta1[irow]);
3252 0 : }
3253 0 : if (npoints<=80) return -1;
3254 0 : if (npoints<npValid*kMinFractionPoints) return -1; // not enough points
3255 0 : rmsTrack = TMath::Sqrt(distance2/npoints);
3256 0 : if (rmsTrack>maxRMSTrack) {
3257 0 : return -1;
3258 : }
3259 : //
3260 : // 2.) Calculate rms distance between the cluster and local median cluster posotion in region +-2 padrows (wihing one sector)
3261 : // and store delta in exported vecDelta array
3262 : Int_t nOutCounter=0; // counter of clusters
3263 0 : rmsCluster=0;
3264 : Float_t rmsSum=0;
3265 0 : Float_t deltaArray[10]={0};
3266 0 : for (Int_t irow=0; irow<npValid; irow++){
3267 : Int_t idelta=1;
3268 0 : vecDelta[irow]=1000;
3269 0 : if (TMath::Abs(delta0[irow])>kMaxDist) continue;
3270 0 : deltaArray[0]=delta0[irow];
3271 : Int_t nforMedian=1;
3272 0 : Int_t sector36=Int_t(sec[irow])%36;
3273 0 : for (idelta=1; idelta<3; idelta++) {
3274 0 : if ((irow-idelta)<0) break;
3275 0 : if ((irow+idelta)>=npValid) break;
3276 0 : if (TMath::Abs(sec[irow-idelta])>72) break;
3277 0 : if (TMath::Abs(sec[irow+idelta])>72) break;
3278 0 : if (Int_t(sec[irow-idelta])%36!=sector36) break; // sometimes looks like sector not initialized
3279 0 : if (Int_t(sec[irow+idelta])%36!=sector36) break;
3280 0 : if (TMath::Abs(delta0[irow-idelta])<kMaxDist) deltaArray[nforMedian++]=delta0[irow-idelta];
3281 0 : if (TMath::Abs(delta0[irow+idelta])<kMaxDist) deltaArray[nforMedian++]=delta0[irow+idelta];
3282 : }
3283 0 : if (irow==0) vecDelta[irow]=delta0[0]-delta0[1];
3284 0 : if (irow==npValid-1) vecDelta[irow]=delta0[irow]-delta0[irow-1];
3285 0 : if (nforMedian>1){
3286 0 : vecDelta[irow]=delta0[irow]-TMath::Mean(nforMedian, deltaArray);
3287 0 : }
3288 0 : if (vecDelta[irow]>maxRMSCluster){
3289 0 : nOutCounter++;
3290 0 : }else{
3291 0 : rmsCluster+=TMath::Power(TMath::Abs(vecDelta[irow]),lpNorm);
3292 0 : rmsSum++;
3293 : }
3294 0 : }
3295 0 : if (rmsSum<=0) return 0;
3296 0 : rmsCluster=TMath::Power(rmsCluster/rmsSum,1./lpNorm);
3297 0 : return nOutCounter;
3298 0 : }
3299 :
3300 :
3301 : Float_t AliTPCcalibAlignInterpolation::InitForAlignmentBugFix(int run, const char* ocdb)
3302 : {
3303 0 : ::Info(" AliTPCcalibAlignInterpolation::InitForAlignmentBugFix","Alignment bug fix is requested\n");
3304 : //
3305 : // this requires the field and the geometry ...
3306 0 : if (run<1) ::Fatal("tstw","InitForBugFix: Run number is not provided");
3307 0 : Bool_t geomOK = AliGeomManager::GetGeometry() != 0;
3308 0 : AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
3309 0 : if (!geomOK || !fld) { // need to setup ocdb?
3310 0 : AliCDBManager* man = AliCDBManager::Instance();
3311 0 : if (!man->IsDefaultStorageSet()) man->SetDefaultStorage(ocdb);
3312 0 : if (man->GetRun()!=run) man->SetRun(run);
3313 0 : }
3314 0 : if (!geomOK) {
3315 0 : AliGeomManager::LoadGeometry();
3316 0 : AliGeomManager::ApplyAlignObjsFromCDB("TPC");
3317 0 : }
3318 0 : if (!fld) {
3319 0 : AliGRPManager grpMan;
3320 0 : grpMan.ReadGRPEntry();
3321 0 : grpMan.SetMagField();
3322 0 : fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
3323 0 : }
3324 0 : return fld->SolenoidField();
3325 0 : }
|
