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 : Comments to be written here:
18 :
19 : 1. What do we calibrate.
20 :
21 : Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
22 :
23 : AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
24 :
25 : */
26 :
27 : #include "Riostream.h"
28 : #include "TDatabasePDG.h"
29 : #include "TGraphErrors.h"
30 : #include "TH1F.h"
31 : #include "THnSparse.h"
32 : #include "TList.h"
33 : #include "TMath.h"
34 : #include "TTimeStamp.h"
35 : #include "TTree.h"
36 : #include "TVectorD.h"
37 : //#include "TChain.h"
38 : //#include "TFile.h"
39 :
40 : #include "AliDCSSensor.h"
41 : #include "AliDCSSensorArray.h"
42 : #include "AliESDEvent.h"
43 : #include "AliESDInputHandler.h"
44 : #include "AliESDVertex.h"
45 : #include "AliESDfriend.h"
46 : #include "AliLog.h"
47 : #include "AliRelAlignerKalman.h"
48 : #include "AliTPCCalROC.h"
49 : #include "AliTPCParam.h"
50 : #include "AliTPCreco.h"
51 : #include "AliTPCTracklet.h"
52 : #include "AliTPCcalibDB.h"
53 : #include "AliTPCcalibLaser.h"
54 : #include "AliTPCcalibTime.h"
55 : #include "AliTPCclusterMI.h"
56 : #include "AliTPCseed.h"
57 : #include "AliTrackPointArray.h"
58 : #include "AliTracker.h"
59 : #include "AliKFVertex.h"
60 : #include <AliLog.h>
61 :
62 6 : ClassImp(AliTPCcalibTime)
63 :
64 : Double_t AliTPCcalibTime::fgResHistoMergeCut = 20000000.;
65 :
66 : AliTPCcalibTime::AliTPCcalibTime()
67 0 : :AliTPCcalibBase(),
68 0 : fMemoryMode(1), // 0 -do not fill THnSparse with residuals 1- fill only important QA THn 2 - Fill all THnsparse for calibration
69 0 : fLaser(0), // pointer to laser calibration
70 0 : fDz(0), // current delta z
71 0 : fCutMaxD(3), // maximal distance in rfi ditection
72 0 : fCutMaxDz(25), // maximal distance in rfi ditection
73 0 : fCutTheta(0.03), // maximal distan theta
74 0 : fCutMinDir(-0.99), // direction vector products
75 0 : fCutTracks(2500),
76 0 : fArrayLaserA(0), //laser fit parameters C
77 0 : fArrayLaserC(0), //laser fit parameters A
78 0 : fArrayDz(0), //NEW! Tmap of V drifts for different triggers
79 0 : fAlignITSTPC(0), //alignemnt array ITS TPC match
80 0 : fAlignTRDTPC(0), //alignemnt array TRD TPC match
81 0 : fAlignTOFTPC(0), //alignemnt array TOF TPC match
82 0 : fTimeKalmanBin(60*15), //time bin width for kalman - 15 minutes default
83 0 : fTimeBins(0),
84 0 : fTimeStart(0),
85 0 : fTimeEnd(0),
86 0 : fPtBins(0),
87 0 : fPtStart(0),
88 0 : fPtEnd(0),
89 0 : fVdriftBins(0),
90 0 : fVdriftStart(0),
91 0 : fVdriftEnd(0),
92 0 : fRunBins(0),
93 0 : fRunStart(0),
94 0 : fRunEnd(0)
95 0 : {
96 : //
97 : // default constructor
98 : //
99 0 : AliDebug(5,"Default Constructor");
100 0 : for (Int_t i=0;i<3;i++) {
101 0 : fHistVdriftLaserA[i]=0;
102 0 : fHistVdriftLaserC[i]=0;
103 : }
104 0 : for (Int_t i=0;i<10;i++) {
105 0 : fCosmiMatchingHisto[i]=0;
106 : }
107 : //
108 0 : for (Int_t i=0;i<5;i++) {
109 0 : fResHistoTPCCE[i]=0;
110 0 : fResHistoTPCITS[i]=0;
111 0 : fResHistoTPCTRD[i]=0;
112 0 : fResHistoTPCTOF[i]=0;
113 0 : fResHistoTPCvertex[i]=0;
114 0 : fTPCVertex[i]=0;
115 : }
116 0 : for (Int_t i=0;i<12;i++) {
117 0 : fTPCVertex[i]=0;
118 : }
119 0 : for (Int_t i=0;i<5;i++) {
120 0 : fTPCVertexCorrelation[i]=0;
121 : }
122 : static Int_t counter=0;
123 : if (1) {
124 0 : TTimeStamp s;
125 0 : Int_t time=s;
126 0 : AliDebug(5,Form("Counter Constructor\t%d\t%d",counter,time));
127 0 : counter++;
128 0 : }
129 :
130 0 : }
131 :
132 : AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeVdrift, Int_t memoryMode)
133 0 : :AliTPCcalibBase(),
134 0 : fMemoryMode(memoryMode), // 0 -do not fill THnSparse with residuals 1- fill only important QA THn 2 - Fill all THnsparse for calibration
135 0 : fLaser(0), // pointer to laser calibration
136 0 : fDz(0), // current delta z
137 0 : fCutMaxD(5*0.5356), // maximal distance in rfi ditection
138 0 : fCutMaxDz(40), // maximal distance in rfi ditection
139 0 : fCutTheta(5*0.004644),// maximal distan theta
140 0 : fCutMinDir(-0.99), // direction vector products
141 0 : fCutTracks(2500),
142 0 : fArrayLaserA(new TObjArray(1000)), //laser fit parameters C
143 0 : fArrayLaserC(new TObjArray(1000)), //laser fit parameters A
144 0 : fArrayDz(0), //Tmap of V drifts for different triggers
145 0 : fAlignITSTPC(0), //alignemnt array ITS TPC match
146 0 : fAlignTRDTPC(0), //alignemnt array TRD TPC match
147 0 : fAlignTOFTPC(0), //alignemnt array TOF TPC match
148 0 : fTimeKalmanBin(60*15), //time bin width for kalman - 15 minutes default
149 0 : fTimeBins(0),
150 0 : fTimeStart(0),
151 0 : fTimeEnd(0),
152 0 : fPtBins(0),
153 0 : fPtStart(0),
154 0 : fPtEnd(0),
155 0 : fVdriftBins(0),
156 0 : fVdriftStart(0),
157 0 : fVdriftEnd(0),
158 0 : fRunBins(0),
159 0 : fRunStart(0),
160 0 : fRunEnd(0)
161 0 : {
162 : //
163 : // Non deafaul constructor - to be used in the Calibration setups
164 : //
165 :
166 0 : SetName(name);
167 0 : SetTitle(title);
168 0 : for (Int_t i=0;i<3;i++) {
169 0 : fHistVdriftLaserA[i]=0;
170 0 : fHistVdriftLaserC[i]=0;
171 : }
172 :
173 0 : for (Int_t i=0;i<5;i++) {
174 0 : fResHistoTPCCE[i]=0;
175 0 : fResHistoTPCITS[i]=0;
176 0 : fResHistoTPCTRD[i]=0;
177 0 : fResHistoTPCTOF[i]=0;
178 0 : fResHistoTPCvertex[i]=0;
179 : }
180 :
181 :
182 0 : AliDebug(5,"Non Default Constructor");
183 0 : fTimeBins =(EndTime-StartTime)/deltaIntegrationTimeVdrift;
184 0 : fTimeStart =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi();
185 0 : fTimeEnd =EndTime; //(((TObjString*)(mapGRP->GetValue("fAliceStopTime")))->GetString()).Atoi();
186 0 : fPtBins = 400;
187 0 : fPtStart = -0.04;
188 0 : fPtEnd = 0.04;
189 0 : fVdriftBins = 500;
190 0 : fVdriftStart= -0.1;
191 0 : fVdriftEnd = 0.1;
192 0 : fRunBins = 1000001;
193 0 : fRunStart = -1.5;
194 0 : fRunEnd = 999999.5;
195 :
196 0 : Int_t binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins };
197 0 : Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart , fRunStart};
198 0 : Double_t xmaxVdriftLaser[4] = {fTimeEnd , fPtEnd , fVdriftEnd , fRunEnd };
199 0 : TString axisTitle[4]={
200 0 : "T",
201 0 : "#delta_{P/T}",
202 0 : "value",
203 0 : "run"
204 : };
205 0 : TString histoName[3]={
206 0 : "Loffset",
207 0 : "Lcorr",
208 0 : "Lgy"
209 : };
210 :
211 :
212 0 : for (Int_t i=0;i<3;i++) {
213 0 : fHistVdriftLaserA[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
214 0 : fHistVdriftLaserC[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
215 0 : fHistVdriftLaserA[i]->SetName(histoName[i]);
216 0 : fHistVdriftLaserC[i]->SetName(histoName[i]);
217 0 : for (Int_t iaxis=0; iaxis<4;iaxis++){
218 0 : fHistVdriftLaserA[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
219 0 : fHistVdriftLaserC[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
220 : }
221 : }
222 0 : fBinsVdrift[0] = fTimeBins;
223 0 : fBinsVdrift[1] = fPtBins;
224 0 : fBinsVdrift[2] = fVdriftBins;
225 0 : fBinsVdrift[3] = fRunBins;
226 0 : fXminVdrift[0] = fTimeStart;
227 0 : fXminVdrift[1] = fPtStart;
228 0 : fXminVdrift[2] = fVdriftStart;
229 0 : fXminVdrift[3] = fRunStart;
230 0 : fXmaxVdrift[0] = fTimeEnd;
231 0 : fXmaxVdrift[1] = fPtEnd;
232 0 : fXmaxVdrift[2] = fVdriftEnd;
233 0 : fXmaxVdrift[3] = fRunEnd;
234 :
235 0 : fArrayDz=new TObjArray();
236 0 : fAlignITSTPC = new TObjArray; //alignemnt array ITS TPC match
237 0 : fAlignTRDTPC = new TObjArray; //alignemnt array ITS TPC match
238 0 : fAlignTOFTPC = new TObjArray; //alignemnt array ITS TPC match
239 0 : fAlignITSTPC->SetOwner(kTRUE);
240 0 : fAlignTRDTPC->SetOwner(kTRUE);
241 0 : fAlignTOFTPC->SetOwner(kTRUE);
242 :
243 :
244 0 : fCosmiMatchingHisto[0]=new TH1F("Cosmics matching p0a","p0-all" ,100,-10*0.5356 ,10*0.5356 );
245 0 : fCosmiMatchingHisto[1]=new TH1F("Cosmics matching p1a","p1-all" ,100,-10*4.541 ,10*4.541 );
246 0 : fCosmiMatchingHisto[2]=new TH1F("Cosmics matching p2a","p2-all" ,100,-10*0.01134 ,10*0.01134 );
247 0 : fCosmiMatchingHisto[3]=new TH1F("Cosmics matching p3a","p3-all" ,100,-10*0.004644,10*0.004644);
248 0 : fCosmiMatchingHisto[4]=new TH1F("Cosmics matching p4a","p4-all" ,100,-10*0.03773 ,10*0.03773 );
249 0 : fCosmiMatchingHisto[5]=new TH1F("Cosmics matching p0p","p0-isPair",100,-10*0.5356 ,10*0.5356 );
250 0 : fCosmiMatchingHisto[6]=new TH1F("Cosmics matching p1p","p1-isPair",100,-10*4.541 ,10*4.541 );
251 0 : fCosmiMatchingHisto[7]=new TH1F("Cosmics matching p2p","p2-isPair",100,-10*0.01134 ,10*0.01134 );
252 0 : fCosmiMatchingHisto[8]=new TH1F("Cosmics matching p3p","p3-isPair",100,-10*0.004644,10*0.004644);
253 0 : fCosmiMatchingHisto[9]=new TH1F("Cosmics matching p4p","p4-isPair",100,-10*0.03773 ,10*0.03773 );
254 0 : for (Int_t i=0;i<12;i++) {
255 0 : fTPCVertex[i]=0;
256 : }
257 0 : for (Int_t i=0;i<5;i++) {
258 0 : fTPCVertexCorrelation[i]=0;
259 : }
260 0 : BookDistortionMaps();
261 :
262 0 : }
263 :
264 0 : AliTPCcalibTime::~AliTPCcalibTime(){
265 : //
266 : // Virtual Destructor
267 : //
268 : static Int_t counter=0;
269 : if (1) {
270 0 : TTimeStamp s;
271 0 : Int_t time=s;
272 0 : AliDebug(5,Form("Counter Destructor\t%s\t%d\t%d",GetName(),counter,time));
273 0 : counter++;
274 0 : }
275 0 : for(Int_t i=0;i<3;i++){
276 0 : if(fHistVdriftLaserA[i]){
277 0 : delete fHistVdriftLaserA[i];
278 0 : fHistVdriftLaserA[i]=NULL;
279 0 : }
280 0 : if(fHistVdriftLaserC[i]){
281 0 : delete fHistVdriftLaserC[i];
282 0 : fHistVdriftLaserC[i]=NULL;
283 0 : }
284 : }
285 0 : if(fArrayDz){
286 0 : fArrayDz->SetOwner();
287 0 : fArrayDz->Delete();
288 0 : delete fArrayDz;
289 0 : fArrayDz=NULL;
290 0 : }
291 0 : for(Int_t i=0;i<10;i++){
292 0 : if(fCosmiMatchingHisto[i]){
293 0 : delete fCosmiMatchingHisto[i];
294 0 : fCosmiMatchingHisto[i]=NULL;
295 0 : }
296 : }
297 :
298 0 : for (Int_t i=0;i<5;i++) {
299 0 : delete fResHistoTPCCE[i];
300 0 : delete fResHistoTPCITS[i];
301 0 : delete fResHistoTPCTRD[i];
302 0 : delete fResHistoTPCTOF[i];
303 0 : delete fResHistoTPCvertex[i];
304 0 : fResHistoTPCCE[i]=0;
305 0 : fResHistoTPCITS[i]=0;
306 0 : fResHistoTPCTRD[i]=0;
307 0 : fResHistoTPCTOF[i]=0;
308 0 : fResHistoTPCvertex[i]=0;
309 : }
310 :
311 0 : for (Int_t i=0;i<12;i++) if (fTPCVertex[i]) delete fTPCVertex[i];
312 0 : for (Int_t i=0;i<5;i++) if (fTPCVertexCorrelation[i]) delete fTPCVertexCorrelation[i];
313 :
314 0 : if (fAlignITSTPC){
315 0 : fAlignITSTPC->SetOwner(kTRUE);
316 0 : fAlignTRDTPC->SetOwner(kTRUE);
317 0 : fAlignTOFTPC->SetOwner(kTRUE);
318 :
319 0 : fAlignITSTPC->Delete();
320 0 : fAlignTRDTPC->Delete();
321 0 : fAlignTOFTPC->Delete();
322 0 : delete fAlignITSTPC;
323 0 : delete fAlignTRDTPC;
324 0 : delete fAlignTOFTPC;
325 : }
326 :
327 0 : if (fArrayLaserA) {
328 0 : fArrayLaserA->SetOwner();
329 0 : fArrayLaserA->Delete();
330 0 : delete fArrayLaserA;
331 : }
332 :
333 0 : if (fArrayLaserA) {
334 0 : fArrayLaserC->SetOwner();
335 0 : fArrayLaserC->Delete();
336 0 : delete fArrayLaserC;
337 : }
338 :
339 0 : }
340 :
341 : // Bool_t AliTPCcalibTime::IsLaser(const AliESDEvent *const /*event*/) const{
342 : // //
343 : // // Indicator is laser event not yet implemented - to be done using trigger info or event specie
344 : // //
345 : // return kTRUE; //More accurate creteria to be added
346 : // }
347 : // Bool_t AliTPCcalibTime::IsCosmics(const AliESDEvent *const /*event*/){
348 : // //
349 : // // Indicator is cosmic event not yet implemented - to be done using trigger info or event specie
350 : // //
351 :
352 : // return kTRUE; //More accurate creteria to be added
353 : // }
354 : // Bool_t AliTPCcalibTime::IsBeam(const AliESDEvent *const /*event*/) const{
355 : // //
356 : // // Indicator is physic event not yet implemented - to be done using trigger info or event specie
357 : // //
358 :
359 : // return kTRUE; //More accurate creteria to be added
360 : // }
361 : void AliTPCcalibTime::ResetCurrent(){
362 : //
363 : //ResetCurrent
364 : //
365 0 : fDz=0; //Reset current dz
366 0 : }
367 :
368 :
369 :
370 : void AliTPCcalibTime::Process(AliESDEvent *event){
371 : //
372 : // main function to make calibration
373 : //
374 0 : if(!event) return;
375 0 : if (event->GetNumberOfTracks()<2) return;
376 0 : AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
377 0 : if (!ESDfriend) {
378 0 : return;
379 : }
380 0 : if (ESDfriend->TestSkipBit()) return;
381 :
382 0 : ResetCurrent();
383 : //if(IsLaser (event))
384 0 : ProcessLaser (event);
385 : //if(IsCosmics(event))
386 0 : ProcessCosmic(event);
387 : //if(IsBeam (event))
388 0 : ProcessBeam (event);
389 0 : }
390 :
391 : void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){
392 : //
393 : // Fit drift velocity using laser
394 : //
395 : // 0. cuts
396 : const Int_t kMinTracks = 40; // minimal number of laser tracks
397 : const Int_t kMinTracksSide = 20; // minimal number of tracks per side
398 : const Float_t kMaxDeltaZ = 30.; // maximal trigger delay
399 : const Float_t kMaxDeltaV = 0.05; // maximal deltaV
400 : const Float_t kMaxRMS = 0.1; // maximal RMS of tracks
401 : //
402 : /*
403 : TCut cutRMS("sqrt(laserA.fElements[4])<0.1&&sqrt(laserC.fElements[4])<0.1");
404 : TCut cutZ("abs(laserA.fElements[0]-laserC.fElements[0])<3");
405 : TCut cutV("abs(laserA.fElements[1]-laserC.fElements[1])<0.01");
406 : TCut cutY("abs(laserA.fElements[2]-laserC.fElements[2])<2");
407 : TCut cutAll = cutRMS+cutZ+cutV+cutY;
408 : */
409 0 : if (event->GetNumberOfTracks()<kMinTracks) return;
410 : //
411 0 : if(!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
412 0 : fLaser->Process(event);
413 0 : if (fLaser->GetNtracks()<kMinTracks) return; // small amount of tracks cut
414 0 : if (fLaser->fFitAside->GetNrows()==0 && fLaser->fFitCside->GetNrows()==0) return; // no fit neither a or C side
415 : //
416 : // debug streamer - activate stream level
417 : // Use it for tuning of the cuts
418 : //
419 : // cuts to be applied
420 : //
421 0 : Int_t isReject[2]={0,0};
422 : //
423 : // not enough tracks
424 0 : if (TMath::Abs((*fLaser->fFitAside)[3]) < kMinTracksSide) isReject[0]|=1;
425 0 : if (TMath::Abs((*fLaser->fFitCside)[3]) < kMinTracksSide) isReject[1]|=1;
426 : // unreasonable z offset
427 0 : if (TMath::Abs((*fLaser->fFitAside)[0])>kMaxDeltaZ) isReject[0]|=2;
428 0 : if (TMath::Abs((*fLaser->fFitCside)[0])>kMaxDeltaZ) isReject[1]|=2;
429 : // unreasonable drift velocity
430 0 : if (TMath::Abs((*fLaser->fFitAside)[1]-1)>kMaxDeltaV) isReject[0]|=4;
431 0 : if (TMath::Abs((*fLaser->fFitCside)[1]-1)>kMaxDeltaV) isReject[1]|=4;
432 : // big chi2
433 0 : if (TMath::Sqrt(TMath::Abs((*fLaser->fFitAside)[4]))>kMaxRMS ) isReject[0]|=8;
434 0 : if (TMath::Sqrt(TMath::Abs((*fLaser->fFitCside)[4]))>kMaxRMS ) isReject[1]|=8;
435 :
436 :
437 :
438 :
439 0 : if (fStreamLevel>0){
440 0 : printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
441 :
442 0 : TTreeSRedirector *cstream = GetDebugStreamer();
443 0 : if (cstream){
444 0 : TTimeStamp tstamp(fTime);
445 0 : (*cstream)<<"laserInfo"<<
446 0 : "run="<<fRun<< // run number
447 0 : "event="<<fEvent<< // event number
448 0 : "time="<<fTime<< // time stamp of event
449 0 : "trigger="<<fTrigger<< // trigger
450 0 : "mag="<<fMagF<< // magnetic field
451 : //laser
452 0 : "rejectA="<<isReject[0]<<
453 0 : "rejectC="<<isReject[1]<<
454 0 : "laserA.="<<fLaser->fFitAside<<
455 0 : "laserC.="<<fLaser->fFitCside<<
456 0 : "laserAC.="<<fLaser->fFitACside<<
457 0 : "trigger="<<event->GetFiredTriggerClasses()<<
458 : "\n";
459 0 : }
460 0 : }
461 : //
462 : // fill histos
463 : //
464 0 : TVectorD vdriftA(5), vdriftC(5),vdriftAC(6);
465 0 : vdriftA=*(fLaser->fFitAside);
466 0 : vdriftC=*(fLaser->fFitCside);
467 0 : vdriftAC=*(fLaser->fFitACside);
468 : Int_t npointsA=0, npointsC=0;
469 : Float_t chi2A=0, chi2C=0;
470 0 : npointsA= TMath::Nint(vdriftA[3]);
471 0 : chi2A= vdriftA[4];
472 0 : npointsC= TMath::Nint(vdriftC[3]);
473 0 : chi2C= vdriftC[4];
474 :
475 0 : if (npointsA>kMinTracksSide || npointsC>kMinTracksSide){
476 0 : TVectorD *fitA = new TVectorD(6);
477 0 : TVectorD *fitC = new TVectorD(6);
478 0 : for (Int_t ipar=0; ipar<5; ipar++){
479 0 : (*fitA)[ipar]=vdriftA[ipar];
480 0 : (*fitC)[ipar]=vdriftC[ipar];
481 : }
482 0 : (*fitA)[5]=fTime;
483 0 : (*fitC)[5]=fTime;
484 0 : fArrayLaserA->AddLast(fitA);
485 0 : fArrayLaserC->AddLast(fitC);
486 0 : }
487 : //
488 :
489 0 : TTimeStamp tstamp(fTime);
490 0 : Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
491 0 : Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
492 : Double_t driftA=0, driftC=0;
493 0 : if (vdriftA[1]>1.-kMaxDeltaV) driftA = 1./vdriftA[1]-1.;
494 0 : if (vdriftC[1]>1.-kMaxDeltaV) driftC = 1./vdriftC[1]-1.;
495 : //
496 0 : Double_t vecDriftLaserA[4]={static_cast<Double_t>(fTime),(ptrelative0+ptrelative1)/2.0,driftA,static_cast<Double_t>(event->GetRunNumber())};
497 0 : Double_t vecDriftLaserC[4]={static_cast<Double_t>(fTime),(ptrelative0+ptrelative1)/2.0,driftC,static_cast<Double_t>(event->GetRunNumber())};
498 : // Double_t vecDrift[4] ={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitACside))[1])-1,event->GetRunNumber()};
499 :
500 0 : for (Int_t icalib=0;icalib<3;icalib++){
501 0 : if (icalib==0){ //z0 shift
502 0 : vecDriftLaserA[2]=vdriftA[0]/250.;
503 0 : vecDriftLaserC[2]=vdriftC[0]/250.;
504 0 : }
505 0 : if (icalib==1){ //vdrel shift
506 0 : vecDriftLaserA[2]=driftA;
507 0 : vecDriftLaserC[2]=driftC;
508 0 : }
509 0 : if (icalib==2){ //gy shift - full gy - full drift
510 0 : vecDriftLaserA[2]=vdriftA[2]/250.;
511 0 : vecDriftLaserC[2]=vdriftC[2]/250.;
512 0 : }
513 : //if (isReject[0]==0) fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
514 : //if (isReject[1]==0) fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
515 0 : fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
516 0 : fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
517 : }
518 0 : }
519 :
520 : void AliTPCcalibTime::ProcessCosmic(const AliESDEvent *const event){
521 : //
522 : // process Cosmic event - track matching A side C side
523 : //
524 0 : if (!event) {
525 0 : Printf("ERROR: ESD not available");
526 0 : return;
527 : }
528 0 : if (event->GetTimeStamp() == 0 ) {
529 0 : Printf("no time stamp!");
530 0 : return;
531 : }
532 :
533 : //fd
534 : // Find cosmic pairs
535 : //
536 : // Track0 is choosen in upper TPC part
537 : // Track1 is choosen in lower TPC part
538 : //
539 : const Int_t kMinClustersCross =30;
540 : const Int_t kMinClusters =80;
541 0 : Int_t ntracks=event->GetNumberOfTracks();
542 0 : if (ntracks==0) return;
543 0 : if (ntracks > fCutTracks) return;
544 :
545 0 : if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
546 0 : AliESDfriend *esdFriend=(AliESDfriend*)(((AliESDEvent*)event)->FindListObject("AliESDfriend"));
547 :
548 0 : TObjArray tpcSeeds(ntracks);
549 0 : Double_t vtxx[3]={0,0,0};
550 0 : Double_t svtxx[3]={0.000001,0.000001,100.};
551 0 : AliESDVertex vtx(vtxx,svtxx);
552 : //
553 : // track loop
554 : //
555 0 : TArrayI clusterSideA(ntracks);
556 0 : TArrayI clusterSideC(ntracks);
557 0 : for (Int_t i=0;i<ntracks;++i) {
558 0 : clusterSideA[i]=0;
559 0 : clusterSideC[i]=0;
560 0 : AliESDtrack *track = event->GetTrack(i);
561 :
562 0 : const AliExternalTrackParam * trackIn = track->GetInnerParam();
563 0 : const AliExternalTrackParam * trackOut = track->GetOuterParam();
564 0 : if (!trackIn) continue;
565 0 : if (!trackOut) continue;
566 :
567 0 : AliESDfriendTrack *friendTrack = (AliESDfriendTrack*) track->GetFriendTrack();
568 0 : if (!friendTrack) continue;
569 0 : if (friendTrack) ProcessSame(track,friendTrack,event);
570 0 : if (friendTrack) ProcessAlignITS(track,friendTrack,event,esdFriend);
571 0 : if (friendTrack) ProcessAlignTRD(track,friendTrack);
572 0 : if (friendTrack) ProcessAlignTOF(track,friendTrack);
573 : TObject *calibObject;
574 : AliTPCseed *seed = 0;
575 0 : for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
576 0 : if (seed) {
577 0 : tpcSeeds.AddAt(seed,i);
578 : Int_t nA=0, nC=0;
579 0 : for (Int_t irow=kMaxRow;irow--;) {
580 0 : AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
581 0 : if (!cl) continue;
582 0 : if ((cl->GetDetector()%36)<18) nA++;
583 0 : if ((cl->GetDetector()%36)>=18) nC++;
584 0 : }
585 0 : clusterSideA[i]=nA;
586 0 : clusterSideC[i]=nC;
587 0 : }
588 0 : }
589 0 : if (ntracks<2) return;
590 : //
591 : // Find pairs
592 : //
593 :
594 0 : for (Int_t i=0;i<ntracks;++i) {
595 0 : AliESDtrack *track0 = event->GetTrack(i);
596 : // track0 - choosen upper part
597 0 : if (!track0) continue;
598 0 : if (!track0->GetOuterParam()) continue;
599 0 : if (track0->GetOuterParam()->GetAlpha()<0) continue;
600 0 : Double_t d1[3];
601 0 : track0->GetDirection(d1);
602 0 : for (Int_t j=0;j<ntracks;++j) {
603 0 : if (i==j) continue;
604 0 : AliESDtrack *track1 = event->GetTrack(j);
605 : //track 1 lower part
606 0 : if (!track1) continue;
607 0 : if (!track1->GetOuterParam()) continue;
608 0 : if (track0->GetTPCNcls()+ track1->GetTPCNcls()< kMinClusters) continue;
609 0 : Int_t nAC = TMath::Max( TMath::Min(clusterSideA[i], clusterSideC[j]),
610 0 : TMath::Min(clusterSideC[i], clusterSideA[j]));
611 0 : if (nAC<kMinClustersCross) continue;
612 0 : Int_t nA0=clusterSideA[i];
613 0 : Int_t nC0=clusterSideC[i];
614 0 : Int_t nA1=clusterSideA[j];
615 0 : Int_t nC1=clusterSideC[j];
616 : // if (track1->GetOuterParam()->GetAlpha()>0) continue;
617 : //
618 0 : Double_t d2[3];
619 0 : track1->GetDirection(d2);
620 :
621 0 : AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
622 0 : AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
623 0 : if (! seed0) continue;
624 0 : if (! seed1) continue;
625 0 : Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
626 0 : Float_t dist0 = track0->GetLinearD(0,0);
627 0 : Float_t dist1 = track1->GetLinearD(0,0);
628 : //
629 : // conservative cuts - convergence to be guarantied
630 : // applying before track propagation
631 0 : if (TMath::Abs(TMath::Abs(dist0)-TMath::Abs(dist1))>fCutMaxD) continue; // distance to the 0,0
632 0 : if (TMath::Abs(dir)<TMath::Abs(fCutMinDir)) continue; // direction vector product
633 0 : Float_t bz = AliTracker::GetBz();
634 0 : Float_t dvertex0[2]; //distance to 0,0
635 0 : Float_t dvertex1[2]; //distance to 0,0
636 0 : track0->GetDZ(0,0,0,bz,dvertex0);
637 0 : track1->GetDZ(0,0,0,bz,dvertex1);
638 0 : if (TMath::Abs(dvertex0[1])>250) continue;
639 0 : if (TMath::Abs(dvertex1[1])>250) continue;
640 : //
641 : //
642 : //
643 0 : Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
644 0 : AliExternalTrackParam param0(*track0);
645 0 : AliExternalTrackParam param1(*track1);
646 : //
647 : // Propagate using Magnetic field and correct fo material budget
648 : //
649 0 : AliTracker::PropagateTrackTo(¶m0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE);
650 0 : AliTracker::PropagateTrackTo(¶m1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE);
651 : //
652 : // Propagate rest to the 0,0 DCA - z should be ignored
653 : //
654 : //Bool_t b0 = ;
655 0 : param0.PropagateToDCA(&vtx,bz,1000);
656 : //Bool_t b1 =
657 0 : param1.PropagateToDCA(&vtx,bz,1000);
658 0 : param0.GetDZ(0,0,0,bz,dvertex0);
659 0 : param1.GetDZ(0,0,0,bz,dvertex1);
660 0 : Double_t xyz0[3];
661 0 : Double_t xyz1[3];
662 0 : param0.GetXYZ(xyz0);
663 0 : param1.GetXYZ(xyz1);
664 0 : Bool_t isPair = IsPair(¶m0,¶m1);
665 0 : Bool_t isCross = IsCross(track0, track1);
666 0 : Bool_t isSame = IsSame(track0, track1);
667 :
668 0 : THnSparse* hist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
669 0 : TString shortName=hist->ClassName();
670 0 : shortName+="_MEAN_VDRIFT_COSMICS_";
671 0 : delete hist;
672 : hist=NULL;
673 :
674 0 : if((isSame) || (isCross && isPair)){
675 0 : if (track0->GetTPCNcls()+ track1->GetTPCNcls()> 80) {
676 0 : fDz = param0.GetZ() - param1.GetZ();
677 0 : Double_t sign=(nA0>nA1)? 1:-1;
678 0 : fDz*=sign;
679 0 : TTimeStamp tstamp(fTime);
680 0 : Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
681 0 : Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
682 0 : Double_t vecDrift[4]={static_cast<Double_t>(fTime),(ptrelative0+ptrelative1)/2.0,fDz/500.0,static_cast<Double_t>(event->GetRunNumber())};
683 : THnSparse* curHist=NULL;
684 0 : TString name="";
685 :
686 0 : name=shortName;
687 0 : name+=event->GetFiredTriggerClasses();
688 0 : name.ToUpper();
689 0 : curHist=(THnSparseF*)fArrayDz->FindObject(name);
690 0 : if(!curHist){
691 0 : curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
692 0 : fArrayDz->AddLast(curHist);
693 : }
694 : // curHist=(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
695 : // if(!curHist){
696 : // curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
697 : // fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
698 : // }
699 0 : curHist->Fill(vecDrift);
700 :
701 0 : name=shortName;
702 0 : name+="ALL";
703 0 : name.ToUpper();
704 0 : curHist=(THnSparseF*)fArrayDz->FindObject(name);
705 0 : if(!curHist){
706 0 : curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
707 0 : fArrayDz->AddLast(curHist);
708 : }
709 : // curHist=(THnSparseF*)(fMapDz->GetValue("all"));
710 : // if(!curHist){
711 : // curHist=new THnSparseF("all","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
712 : // fMapDz->Add(new TObjString("all"),curHist);
713 : // }
714 0 : curHist->Fill(vecDrift);
715 0 : }
716 : }
717 0 : TTreeSRedirector *cstream = GetDebugStreamer();
718 0 : if (fStreamLevel>0){
719 0 : if (cstream){
720 0 : (*cstream)<<"trackInfo"<<
721 0 : "tr0.="<<track0<<
722 0 : "tr1.="<<track1<<
723 0 : "p0.="<<¶m0<<
724 0 : "p1.="<<¶m1<<
725 0 : "nAC="<<nAC<<
726 0 : "nA0="<<nA0<<
727 0 : "nA1="<<nA1<<
728 0 : "nC0="<<nC0<<
729 0 : "nC1="<<nC1<<
730 0 : "isPair="<<isPair<<
731 0 : "isCross="<<isCross<<
732 0 : "isSame="<<isSame<<
733 0 : "fDz="<<fDz<<
734 0 : "fRun="<<fRun<<
735 0 : "fTime="<<fTime<<
736 : "\n";
737 : }
738 : }
739 0 : } // end 2nd order loop
740 0 : } // end 1st order loop
741 :
742 0 : if (fStreamLevel>0){
743 0 : TTreeSRedirector *cstream = GetDebugStreamer();
744 0 : if (cstream){
745 0 : (*cstream)<<"timeInfo"<<
746 0 : "run="<<fRun<< // run number
747 0 : "event="<<fEvent<< // event number
748 0 : "time="<<fTime<< // time stamp of event
749 0 : "trigger="<<fTrigger<< // trigger
750 0 : "mag="<<fMagF<< // magnetic field
751 : // Environment values
752 : //
753 : // accumulated values
754 : //
755 0 : "fDz="<<fDz<< //! current delta z
756 0 : "trigger="<<event->GetFiredTriggerClasses()<<
757 : "\n";
758 0 : }
759 0 : }
760 0 : if (GetDebugLevel()>20) printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
761 0 : }
762 :
763 : void AliTPCcalibTime::ProcessBeam(const AliESDEvent *const event){
764 : //
765 : // Process beam data - calculates vartex
766 : // from A side and C side
767 : // Histogram the differences
768 : //
769 : const Int_t kMinClusters =80;
770 : const Int_t kMinTracks =2; // minimal number of tracks to define the vertex
771 : const Int_t kMinTracksVertex=30; // minimal number of tracks to define the cumulative vertex
772 : const Double_t kMaxTgl =1.2; // maximal Tgl (z angle)
773 : const Double_t kMinPt =0.2; // minimal pt
774 : const Double_t kMaxD0 =5.; // cut on distance to the primary vertex first guess
775 : const Double_t kMaxZ0 =20;
776 : const Double_t kMaxD =2.5; // cut on distance to the primary vertex
777 : const Double_t kMaxZ =4; // maximal z distance between tracks form the same side
778 : const Double_t kMaxChi2 =15; // maximal chi2 of the TPCvertex
779 : const Double_t kCumulCovarXY=0.003; //increase the error of cumul vertex 30 microns profile
780 : const Double_t kCumulCovarZ=250.; //increase the error of cumul vertex
781 : const Double_t kMaxDvertex = 1.0; // cut to accept the vertex;
782 : //
783 0 : Int_t flags=0;
784 : const Int_t kBuffSize=100;
785 : static Double_t deltaZ[kBuffSize]={0};
786 : static Int_t counterZ=0;
787 0 : static AliKFVertex cumulVertexA, cumulVertexC, cumulVertexAC; // cumulative vertex
788 0 : AliKFVertex vertexA, vertexC;
789 :
790 0 : Float_t dca0[2]={0,0};
791 0 : Double_t dcaVertex[2]={0,0};
792 0 : Int_t ntracks=event->GetNumberOfTracks();
793 0 : if (ntracks==0) return;
794 0 : if (ntracks > fCutTracks) return;
795 : //
796 0 : AliESDfriend *esdFriend=(AliESDfriend*)(((AliESDEvent*)event)->FindListObject("AliESDfriend"));
797 : //
798 : // Divide tracks to A and C side tracks - using the cluster indexes
799 0 : TObjArray tracksA(ntracks);
800 0 : TObjArray tracksC(ntracks);
801 : //
802 0 : AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
803 0 : AliESDVertex *vertex = (AliESDVertex *)event->GetPrimaryVertex();
804 0 : AliESDVertex *vertexTracks = (AliESDVertex *)event->GetPrimaryVertexTracks();
805 0 : Double_t vertexZA[10000], vertexZC[10000];
806 : //
807 0 : Int_t ntracksA= 0;
808 0 : Int_t ntracksC= 0;
809 : //
810 0 : for (Int_t itrack=0;itrack<ntracks;itrack++) {
811 0 : AliESDtrack *track = event->GetTrack(itrack);
812 0 : AliESDfriendTrack *friendTrack = (AliESDfriendTrack*)track->GetFriendTrack();
813 0 : if (!friendTrack) continue;
814 0 : if (TMath::Abs(track->GetTgl())>kMaxTgl) continue;
815 0 : if (TMath::Abs(track->Pt())<kMinPt) continue;
816 0 : const AliExternalTrackParam * trackIn = track->GetInnerParam();
817 : TObject *calibObject=0;
818 : AliTPCseed *seed = 0;
819 : Int_t nA=0, nC=0;
820 0 : for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
821 0 : if (seed) {
822 0 : for (Int_t irow=kMaxRow;irow--;) {
823 0 : AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
824 0 : if (!cl) continue;
825 0 : if ((cl->GetDetector()%36)<18) nA++;
826 0 : if ((cl->GetDetector()%36)>=18) nC++;
827 0 : }
828 0 : if ((nA>kMinClusters || nC>kMinClusters) && (nA*nC==0) ){
829 0 : track->GetImpactParameters(dca0[0],dca0[1]);
830 0 : if (TMath::Abs(dca0[0])>kMaxD0) continue;
831 0 : if (TMath::Abs(dca0[1])>kMaxZ0) continue;
832 0 : AliExternalTrackParam pTPCvertex(*trackIn);
833 0 : if (!AliTracker::PropagateTrackToBxByBz(&pTPCvertex,4.+4.*TMath::Abs(dca0[0]),0.1,2,kTRUE)) continue;
834 0 : pTPCvertex.PropagateToDCA(vertex,AliTracker::GetBz(), kMaxD, dcaVertex,0);
835 0 : if (TMath::Abs(dcaVertex[0])>kMaxD) continue;
836 0 : if (nA>kMinClusters &&nC==0) { tracksA.AddLast(pTPCvertex.Clone()); vertexZA[ntracksA++] = pTPCvertex.GetZ();}
837 0 : if (nC>kMinClusters &&nA==0) {tracksC.AddLast(pTPCvertex.Clone()); vertexZC[ntracksC++] = pTPCvertex.GetZ();}
838 0 : }
839 : }
840 0 : }
841 0 : Double_t medianZA=TMath::Median(ntracksA, vertexZA); // tracks median
842 0 : Double_t medianZC=TMath::Median(ntracksC, vertexZC); // tracks median
843 : //
844 0 : ntracksA= tracksA.GetEntriesFast();
845 0 : ntracksC= tracksC.GetEntriesFast();
846 0 : if (ntracksA>kMinTracks && ntracksC>kMinTracks){
847 0 : deltaZ[counterZ%kBuffSize]=medianZA-medianZC;
848 0 : counterZ+=1;
849 0 : Double_t medianDelta=(counterZ>=kBuffSize)? TMath::Median(kBuffSize, deltaZ): TMath::Median(counterZ, deltaZ);
850 0 : if (TMath::Abs(medianDelta-(medianZA-medianZC))>kMaxZ) flags+=16;
851 : // increse the error of cumulative vertex at the beginning of event
852 0 : cumulVertexA.Covariance(0,0)+=kCumulCovarXY*kCumulCovarXY;
853 0 : cumulVertexA.Covariance(1,1)+=kCumulCovarXY*kCumulCovarXY;
854 0 : cumulVertexA.Covariance(2,2)+=kCumulCovarZ*kCumulCovarZ;
855 0 : cumulVertexC.Covariance(0,0)+=kCumulCovarXY*kCumulCovarXY;
856 0 : cumulVertexC.Covariance(1,1)+=kCumulCovarXY*kCumulCovarXY;
857 0 : cumulVertexC.Covariance(2,2)+=kCumulCovarZ*kCumulCovarZ;
858 0 : cumulVertexAC.Covariance(0,0)+=kCumulCovarXY*kCumulCovarXY;
859 0 : cumulVertexAC.Covariance(1,1)+=kCumulCovarXY*kCumulCovarXY;
860 0 : cumulVertexAC.Covariance(2,2)+=kCumulCovarZ*kCumulCovarZ;
861 : //
862 0 : for (Int_t iA=0; iA<ntracksA; iA++){
863 0 : if (flags!=0) continue;
864 0 : AliExternalTrackParam *aliTrack = (AliExternalTrackParam *)tracksA.At(iA);
865 0 : if (TMath::Abs(aliTrack->GetZ()-medianZA)>kMaxZ) continue;
866 0 : AliKFParticle part(*aliTrack,211);
867 0 : vertexA+=part;
868 0 : }
869 0 : for (Int_t iC=0; iC<ntracksC; iC++){
870 0 : if (flags!=0) continue;
871 0 : AliExternalTrackParam *aliTrack = (AliExternalTrackParam *)tracksC.At(iC);
872 0 : if (TMath::Abs(aliTrack->GetZ()-medianZC)>kMaxZ) continue;
873 0 : AliKFParticle part(*aliTrack,211);
874 0 : vertexC+=part;
875 0 : }
876 : //
877 0 : if (vertexA.GetNDF()<kMinTracks) flags+=32;
878 0 : if (vertexC.GetNDF()<kMinTracks) flags+=32;
879 0 : if (TMath::Abs(vertexA.Z()-medianZA)>kMaxZ) flags+=1; //apply cuts
880 0 : if (TMath::Abs(vertexC.Z()-medianZC)>kMaxZ) flags+=2;
881 0 : if (TMath::Abs(vertexA.GetChi2()/vertexA.GetNDF()+vertexC.GetChi2()/vertexC.GetNDF())> kMaxChi2) flags+=4;
882 : //
883 0 : if (flags==0){
884 0 : for (Int_t iA=0; iA<ntracksA; iA++){
885 0 : if (flags!=0) continue;
886 0 : AliExternalTrackParam *aliTrack = (AliExternalTrackParam *)tracksA.At(iA);
887 0 : if (TMath::Abs(aliTrack->GetZ()-medianZA)>kMaxZ) continue;
888 0 : AliKFParticle part(*aliTrack,211);
889 0 : cumulVertexA+=part;
890 0 : cumulVertexAC+=part;
891 0 : }
892 0 : for (Int_t iC=0; iC<ntracksC; iC++){
893 0 : if (flags!=0) continue;
894 0 : AliExternalTrackParam *aliTrack = (AliExternalTrackParam *)tracksC.At(iC);
895 0 : if (TMath::Abs(aliTrack->GetZ()-medianZC)>kMaxZ) continue;
896 0 : AliKFParticle part(*aliTrack,211);
897 0 : cumulVertexC+=part;
898 0 : cumulVertexAC+=part;
899 0 : }
900 : //
901 0 : if (TMath::Abs(cumulVertexA.X()-vertexA.X())>kMaxDvertex) flags+=64;
902 0 : if (TMath::Abs(cumulVertexA.Y()-vertexA.Y())>kMaxDvertex) flags+=64;
903 0 : if (TMath::Abs(cumulVertexA.Z()-vertexA.Z())>kMaxDvertex) flags+=64;
904 : //
905 0 : if (TMath::Abs(cumulVertexC.X()-vertexC.X())>kMaxDvertex) flags+=64;
906 0 : if (TMath::Abs(cumulVertexC.Y()-vertexC.Y())>kMaxDvertex) flags+=64;
907 0 : if (TMath::Abs(cumulVertexC.Z()-vertexC.Z())>kMaxDvertex) flags+=64;
908 :
909 :
910 0 : if ( flags==0 && cumulVertexC.GetNDF()>kMinTracksVertex&&cumulVertexA.GetNDF()>kMinTracksVertex){
911 0 : Double_t cont[2]={0,static_cast<Double_t>(fTime)};
912 : //
913 0 : cont[0]= cumulVertexA.X();
914 0 : fTPCVertex[0]->Fill(cont);
915 0 : cont[0]= cumulVertexC.X();
916 0 : fTPCVertex[1]->Fill(cont);
917 0 : cont[0]= 0.5*(cumulVertexA.X()-cumulVertexC.X());
918 0 : fTPCVertex[2]->Fill(cont);
919 0 : cont[0]= 0.5*(cumulVertexA.X()+cumulVertexC.X())-vertexSPD->GetX();
920 0 : fTPCVertex[3]->Fill(cont);
921 : //
922 0 : cont[0]= cumulVertexA.Y();
923 0 : fTPCVertex[4]->Fill(cont);
924 0 : cont[0]= cumulVertexC.Y();
925 0 : fTPCVertex[5]->Fill(cont);
926 0 : cont[0]= 0.5*(cumulVertexA.Y()-cumulVertexC.Y());
927 0 : fTPCVertex[6]->Fill(cont);
928 0 : cont[0]= 0.5*(cumulVertexA.Y()+cumulVertexC.Y())-vertexSPD->GetY();
929 0 : fTPCVertex[7]->Fill(cont);
930 : //
931 : //
932 0 : cont[0]= 0.5*(cumulVertexA.Z()+cumulVertexC.Z());
933 0 : fTPCVertex[8]->Fill(cont);
934 0 : cont[0]= 0.5*(cumulVertexA.Z()-cumulVertexC.Z());
935 0 : fTPCVertex[9]->Fill(cont);
936 0 : cont[0]= 0.5*(cumulVertexA.Z()-cumulVertexC.Z());
937 0 : fTPCVertex[10]->Fill(cont);
938 0 : cont[0]= 0.5*(cumulVertexA.Z()+cumulVertexC.Z())-vertexSPD->GetZ();
939 0 : fTPCVertex[11]->Fill(cont);
940 : //
941 0 : Double_t correl[2]={0,0};
942 : //
943 0 : correl[0]=cumulVertexC.Z();
944 0 : correl[1]=cumulVertexA.Z();
945 0 : fTPCVertexCorrelation[0]->Fill(correl); // fill A side :TPC
946 0 : correl[0]=cumulVertexA.Z();
947 0 : correl[1]=cumulVertexC.Z();
948 0 : fTPCVertexCorrelation[1]->Fill(correl); // fill C side :TPC
949 : //
950 0 : correl[0]=vertexSPD->GetZ();
951 0 : correl[1]=cumulVertexA.Z()-correl[0];
952 0 : fTPCVertexCorrelation[2]->Fill(correl); // fill A side :ITS
953 0 : correl[1]=cumulVertexC.Z()-correl[0];
954 0 : fTPCVertexCorrelation[3]->Fill(correl); // fill C side :ITS
955 0 : correl[1]=0.5*(cumulVertexA.Z()+cumulVertexC.Z())-correl[0];
956 0 : fTPCVertexCorrelation[4]->Fill(correl); // fill C side :ITS
957 0 : }
958 : }
959 0 : TTreeSRedirector *cstream = GetDebugStreamer();
960 0 : if (cstream){
961 : /*
962 : TCut cutChi2= "sqrt(vA.fChi2/vA.fNDF+vC.fChi2/vC.fNDF)<10"; // chi2 Cut e.g 10
963 : TCut cutXY= "sqrt((vA.fP[0]-vC.fP[0])^2+(vA.fP[0]-vC.fP[1])^2)<5"; // vertex Cut
964 : TCut cutZ= "abs(vA.fP[2]-mZA)<3&&abs(vC.fP[2]-mZC)<5"; // vertex Cut
965 : tree->Draw("sqrt(vA.fChi2/vA.fNDF)","sqrt(vA.fChi2/vA.fNDF)<100","")
966 :
967 : */
968 : //vertexA.Print();
969 : //vertexC.Print();
970 0 : (*cstream)<<"vertexTPC"<<
971 0 : "flags="<<flags<< // rejection flags
972 0 : "vSPD.="<<vertexSPD<< // SPD vertex
973 0 : "vT.="<<vertexTracks<< // track vertex
974 0 : "v.="<<vertex<< // esd vertex
975 0 : "mZA="<<medianZA<< // median Z position at vertex A side
976 0 : "mZC="<<medianZC<< // median Z position at vertex C side
977 0 : "mDelta="<<medianDelta<< // median delta A side -C side
978 0 : "counter="<<counterZ<< // counter Z
979 : //
980 0 : "vA.="<<&vertexA<< // vertex A side
981 0 : "vC.="<<&vertexC<< // vertex C side
982 0 : "cvA.="<<&cumulVertexA<< // cumulative vertex A side
983 0 : "cvC.="<<&cumulVertexC<< // cumulative vertex C side
984 0 : "cvAC.="<<&cumulVertexAC<< // cumulative vertex A+C side
985 0 : "nA="<<ntracksA<< // contributors
986 0 : "nC="<<ntracksC<< // contributors
987 : "\n";
988 : }
989 0 : }
990 0 : tracksA.Delete();
991 0 : tracksC.Delete();
992 0 : }
993 :
994 : void AliTPCcalibTime::Analyze(){
995 : //
996 : // Special macro to analyze result of calibration and extract calibration entries
997 : // Not yet ported to the Analyze function yet
998 : //
999 0 : }
1000 :
1001 : THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name) const
1002 : {
1003 : //
1004 : // Get histogram for given trigger mask
1005 : //
1006 0 : TIterator* iterator = fArrayDz->MakeIterator();
1007 0 : iterator->Reset();
1008 0 : TString newName=name;
1009 0 : newName.ToUpper();
1010 0 : THnSparse* newHist=new THnSparseF(newName,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
1011 : THnSparse* addHist=NULL;
1012 0 : while((addHist=(THnSparseF*)iterator->Next())){
1013 : // if(!addHist) continue;
1014 0 : TString histName=addHist->GetName();
1015 0 : if(!histName.Contains(newName)) continue;
1016 0 : addHist->Print();
1017 0 : newHist->Add(addHist);
1018 0 : }
1019 : return newHist;
1020 0 : }
1021 :
1022 : TObjArray* AliTPCcalibTime::GetHistoDrift() const
1023 : {
1024 : //
1025 : // return array of histograms
1026 : //
1027 0 : return fArrayDz;
1028 : }
1029 :
1030 : TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){
1031 : //
1032 : // Make a drift velocity (delta Z) graph
1033 : //
1034 0 : THnSparse* histoDrift=GetHistoDrift(name);
1035 : TGraphErrors* graphDrift=NULL;
1036 0 : if(histoDrift){
1037 0 : graphDrift=FitSlices(histoDrift,2,0,400,100,0.05,0.95, kTRUE);
1038 0 : TString end=histoDrift->GetName();
1039 0 : Int_t pos=end.Index("_");
1040 0 : end=end(pos,end.Capacity()-pos);
1041 0 : TString graphName=graphDrift->ClassName();
1042 0 : graphName+=end;
1043 0 : graphName.ToUpper();
1044 0 : graphDrift->SetName(graphName);
1045 0 : }
1046 0 : return graphDrift;
1047 0 : }
1048 :
1049 : TObjArray* AliTPCcalibTime::GetGraphDrift(){
1050 : //
1051 : // make a array of drift graphs
1052 : //
1053 0 : TObjArray* arrayGraphDrift=new TObjArray();
1054 0 : TIterator* iterator=fArrayDz->MakeIterator();
1055 0 : iterator->Reset();
1056 : THnSparse* addHist=NULL;
1057 0 : while((addHist=(THnSparseF*)iterator->Next())) arrayGraphDrift->AddLast(GetGraphDrift(addHist->GetName()));
1058 0 : return arrayGraphDrift;
1059 0 : }
1060 :
1061 : AliSplineFit* AliTPCcalibTime::GetFitDrift(const char* name){
1062 : //
1063 : // Make a fit AliSplinefit of drift velocity
1064 : //
1065 0 : TGraph* graphDrift=GetGraphDrift(name);
1066 : AliSplineFit* fitDrift=NULL;
1067 0 : if(graphDrift && graphDrift->GetN()){
1068 0 : fitDrift=new AliSplineFit();
1069 0 : fitDrift->SetGraph(graphDrift);
1070 0 : fitDrift->SetMinPoints(graphDrift->GetN()+1);
1071 0 : fitDrift->InitKnots(graphDrift,2,0,0.001);
1072 0 : fitDrift->SplineFit(0);
1073 0 : TString end=graphDrift->GetName();
1074 0 : Int_t pos=end.Index("_");
1075 0 : end=end(pos,end.Capacity()-pos);
1076 0 : TString fitName=fitDrift->ClassName();
1077 0 : fitName+=end;
1078 0 : fitName.ToUpper();
1079 : //fitDrift->SetName(fitName);
1080 0 : delete graphDrift;
1081 : graphDrift=NULL;
1082 0 : }
1083 0 : return fitDrift;
1084 0 : }
1085 :
1086 :
1087 : Long64_t AliTPCcalibTime::Merge(TCollection *const li) {
1088 : //
1089 : // Object specific merging procedure
1090 : //
1091 0 : TIterator* iter = li->MakeIterator();
1092 : AliTPCcalibTime* cal = 0;
1093 : //
1094 0 : while ((cal = (AliTPCcalibTime*)iter->Next())) {
1095 0 : if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
1096 0 : Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1097 0 : return -1;
1098 : }
1099 0 : for (Int_t imeas=0; imeas<3; imeas++){
1100 0 : if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
1101 0 : fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
1102 0 : fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
1103 0 : }
1104 : }
1105 : //
1106 0 : if (fTPCVertexCorrelation[0] && cal->fTPCVertexCorrelation[0]){
1107 0 : for (Int_t imeas=0; imeas<5; imeas++){
1108 0 : if (fTPCVertexCorrelation[imeas] && cal->fTPCVertexCorrelation[imeas]) fTPCVertexCorrelation[imeas]->Add(cal->fTPCVertexCorrelation[imeas]);
1109 : }
1110 0 : }
1111 :
1112 0 : if (fTPCVertex[0] && cal->fTPCVertex[0])
1113 0 : for (Int_t imeas=0; imeas<12; imeas++){
1114 0 : if (fTPCVertex[imeas] && cal->fTPCVertex[imeas]) fTPCVertex[imeas]->Add(cal->fTPCVertex[imeas]);
1115 0 : }
1116 :
1117 0 : if (fMemoryMode>0) for (Int_t imeas=0; imeas<5; imeas++){
1118 0 : if (fMemoryMode>1){
1119 0 : if ( GetResHistoTPCCE(imeas) && cal->GetResHistoTPCCE(imeas)){
1120 0 : if ((cal->GetResHistoTPCCE(imeas)->GetEntries()+GetResHistoTPCCE(imeas)->GetEntries()) < fgResHistoMergeCut)
1121 0 : fResHistoTPCCE[imeas]->Add(cal->fResHistoTPCCE[imeas]);
1122 : }else{
1123 0 : fResHistoTPCCE[imeas]=(THnSparse*)cal->fResHistoTPCCE[imeas]->Clone();
1124 : }
1125 : }
1126 : //
1127 0 : if ((fMemoryMode>0) &&cal->GetResHistoTPCITS(imeas) && cal->GetResHistoTPCITS(imeas)){
1128 0 : if (fMemoryMode>1 || (imeas%2)==1)
1129 : {
1130 0 : if (fResHistoTPCITS[imeas]->GetEntries()+(cal->fResHistoTPCITS[imeas])->GetEntries() < fgResHistoMergeCut)
1131 0 : fResHistoTPCITS[imeas]->Add(cal->fResHistoTPCITS[imeas]);
1132 : else
1133 0 : AliInfo(Form("fResHistoTPCITS[%i] full (has %.0f merged tracks, trying to add %.0f, max allowed: %.0f)",imeas, fResHistoTPCITS[imeas]->GetEntries(), cal->fResHistoTPCITS[imeas]->GetEntries(), fgResHistoMergeCut));
1134 : }
1135 0 : if (fMemoryMode>1)
1136 : {
1137 0 : if (fResHistoTPCvertex[imeas]->GetEntries()+(cal->fResHistoTPCvertex[imeas])->GetEntries() < fgResHistoMergeCut)
1138 0 : fResHistoTPCvertex[imeas]->Add(cal->fResHistoTPCvertex[imeas]);
1139 : }
1140 : }
1141 : //
1142 0 : if ((fMemoryMode>1) && cal->fResHistoTPCTRD[imeas]){
1143 0 : if (fResHistoTPCTRD[imeas] && cal->fResHistoTPCTRD[imeas])
1144 : {
1145 0 : if (fResHistoTPCTRD[imeas]->GetEntries()+(cal->fResHistoTPCTRD[imeas])->GetEntries() < fgResHistoMergeCut)
1146 0 : fResHistoTPCTRD[imeas]->Add(cal->fResHistoTPCTRD[imeas]);
1147 : }
1148 : else
1149 0 : fResHistoTPCTRD[imeas]=(THnSparse*)cal->fResHistoTPCTRD[imeas]->Clone();
1150 : }
1151 : //
1152 0 : if ((fMemoryMode>1) && cal->fResHistoTPCTOF[imeas]){
1153 0 : if (fResHistoTPCTOF[imeas] && cal->fResHistoTPCTOF[imeas])
1154 : {
1155 0 : if (fResHistoTPCTOF[imeas]->GetEntries()+(cal->fResHistoTPCTOF[imeas])->GetEntries() < fgResHistoMergeCut)
1156 0 : fResHistoTPCTOF[imeas]->Add(cal->fResHistoTPCTOF[imeas]);
1157 : }
1158 : else
1159 0 : fResHistoTPCTOF[imeas]=(THnSparse*)cal->fResHistoTPCTOF[imeas]->Clone();
1160 : }
1161 : //
1162 0 : if (cal->fArrayLaserA){
1163 0 : fArrayLaserA->Expand(fArrayLaserA->GetEntriesFast()+cal->fArrayLaserA->GetEntriesFast());
1164 0 : fArrayLaserC->Expand(fArrayLaserC->GetEntriesFast()+cal->fArrayLaserC->GetEntriesFast());
1165 0 : for (Int_t ical=0; ical<cal->fArrayLaserA->GetEntriesFast(); ical++){
1166 0 : if (cal->fArrayLaserA->UncheckedAt(ical)) fArrayLaserA->AddLast(cal->fArrayLaserA->UncheckedAt(ical)->Clone());
1167 0 : if (cal->fArrayLaserC->UncheckedAt(ical)) fArrayLaserC->AddLast(cal->fArrayLaserC->UncheckedAt(ical)->Clone());
1168 : }
1169 0 : }
1170 :
1171 0 : }
1172 : // TObjArray* addArray=cal->GetHistoDrift();
1173 : // if(!addArray) return 0;
1174 : // TIterator* iterator = addArray->MakeIterator();
1175 : // iterator->Reset();
1176 : // THnSparse* addHist=NULL;
1177 : // if ((fMemoryMode>1)) while((addHist=(THnSparseF*)iterator->Next())){
1178 : // // if(!addHist) continue;
1179 : // addHist->Print();
1180 : // THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName());
1181 : // if(!localHist){
1182 : // localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
1183 : // fArrayDz->AddLast(localHist);
1184 : // }
1185 : // localHist->Add(addHist);
1186 : // }
1187 : // delete iterator;
1188 0 : for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
1189 : //
1190 : // Merge alignment
1191 : //
1192 : const Int_t kMinUpdates=10;
1193 : const Float_t kMaxOut=0.1;
1194 0 : for (Int_t itype=0; itype<3; itype++){
1195 : //
1196 : //
1197 : TObjArray *arr0= 0;
1198 : TObjArray *arr1= 0;
1199 0 : if (itype==0) {arr0=fAlignITSTPC; arr1=cal->fAlignITSTPC;}
1200 0 : if (itype==1) {arr0=fAlignTRDTPC; arr1=cal->fAlignTRDTPC;}
1201 0 : if (itype==2) {arr0=fAlignTOFTPC; arr1=cal->fAlignTOFTPC;}
1202 0 : if (!arr1) continue;
1203 0 : if (!arr0) arr0=new TObjArray(arr1->GetEntriesFast());
1204 0 : if (arr1->GetEntriesFast()>arr0->GetEntriesFast()){
1205 0 : arr0->Expand(arr1->GetEntriesFast());
1206 0 : }
1207 0 : for (Int_t i=0;i<arr1->GetEntriesFast(); i++){
1208 0 : AliRelAlignerKalman *kalman1 = (AliRelAlignerKalman *)arr1->UncheckedAt(i);
1209 0 : AliRelAlignerKalman *kalman0 = (AliRelAlignerKalman *)arr0->UncheckedAt(i);
1210 0 : if (!kalman1) continue;
1211 0 : if (kalman1->GetNUpdates()<kMinUpdates) continue;
1212 0 : if (kalman1->GetNOutliers()>(kalman1->GetNUpdates()*kMaxOut)) continue;
1213 0 : if (!kalman0) {arr0->AddAt(new AliRelAlignerKalman(*kalman1),i); continue;}
1214 0 : kalman0->SetRejectOutliers(kFALSE);
1215 0 : kalman0->Merge(kalman1);
1216 0 : }
1217 0 : }
1218 :
1219 : }
1220 0 : delete iter;
1221 0 : return 0;
1222 0 : }
1223 :
1224 : Bool_t AliTPCcalibTime::IsPair(const AliExternalTrackParam *tr0, const AliExternalTrackParam *tr1){
1225 : /*
1226 : // 0. Same direction - OPOSITE - cutDir +cutT
1227 : TCut cutDir("cutDir","dir<-0.99")
1228 : // 1.
1229 : TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
1230 : //
1231 : // 2. The same rphi
1232 : TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
1233 : //
1234 : TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
1235 : // 1/Pt diff cut
1236 : */
1237 0 : const Double_t *p0 = tr0->GetParameter();
1238 0 : const Double_t *p1 = tr1->GetParameter();
1239 0 : fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
1240 0 : fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
1241 0 : fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
1242 0 : fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
1243 0 : fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
1244 :
1245 0 : if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
1246 0 : if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
1247 0 : if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
1248 0 : Double_t d0[3], d1[3];
1249 0 : tr0->GetDirection(d0);
1250 0 : tr1->GetDirection(d1);
1251 0 : if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
1252 :
1253 0 : fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
1254 0 : fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
1255 0 : fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
1256 0 : fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
1257 0 : fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
1258 :
1259 0 : return kTRUE;
1260 0 : }
1261 : Bool_t AliTPCcalibTime::IsCross(const AliESDtrack *const tr0, const AliESDtrack *const tr1){
1262 : //
1263 : // check if the cosmic pair of tracks crossed A/C side
1264 : //
1265 0 : Bool_t result= tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0;
1266 0 : if (result==kFALSE) return result;
1267 : result=kTRUE;
1268 0 : return result;
1269 0 : }
1270 :
1271 : Bool_t AliTPCcalibTime::IsSame(const AliESDtrack *const tr0, const AliESDtrack *const tr1){
1272 : //
1273 : // track crossing the CE
1274 : // 0. minimal number of clusters
1275 : // 1. Same sector +-1
1276 : // 2. Inner and outer track param on opposite side
1277 : // 3. Outer and inner track parameter close each to other
1278 : // 3.
1279 : Bool_t result=kTRUE;
1280 : //
1281 : // inner and outer on opposite sides in z
1282 : //
1283 : const Int_t knclCut0 = 30;
1284 : const Double_t kalphaCut = 0.4;
1285 : //
1286 : // 0. minimal number of clusters
1287 : //
1288 0 : if (tr0->GetTPCNcls()<knclCut0) return kFALSE;
1289 0 : if (tr1->GetTPCNcls()<knclCut0) return kFALSE;
1290 : //
1291 : // 1. alpha cut - sector+-1
1292 : //
1293 0 : if (TMath::Abs(tr0->GetOuterParam()->GetAlpha()-tr1->GetOuterParam()->GetAlpha())>kalphaCut) return kFALSE;
1294 : //
1295 : // 2. Z crossing
1296 : //
1297 0 : if (tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0) result&=kFALSE;
1298 0 : if (tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0) result&=kFALSE;
1299 0 : if (result==kFALSE){
1300 0 : return result;
1301 : }
1302 : //
1303 : //
1304 0 : const Double_t *p0I = tr0->GetInnerParam()->GetParameter();
1305 0 : const Double_t *p1I = tr1->GetInnerParam()->GetParameter();
1306 0 : const Double_t *p0O = tr0->GetOuterParam()->GetParameter();
1307 0 : const Double_t *p1O = tr1->GetOuterParam()->GetParameter();
1308 : //
1309 0 : if (TMath::Abs(p0I[0]-p1I[0])>fCutMaxD) result&=kFALSE;
1310 0 : if (TMath::Abs(p0I[1]-p1I[1])>fCutMaxDz) result&=kFALSE;
1311 0 : if (TMath::Abs(p0I[2]-p1I[2])>fCutTheta) result&=kFALSE;
1312 0 : if (TMath::Abs(p0I[3]-p1I[3])>fCutTheta) result&=kFALSE;
1313 0 : if (TMath::Abs(p0O[0]-p1O[0])>fCutMaxD) result&=kFALSE;
1314 0 : if (TMath::Abs(p0O[1]-p1O[1])>fCutMaxDz) result&=kFALSE;
1315 0 : if (TMath::Abs(p0O[2]-p1O[2])>fCutTheta) result&=kFALSE;
1316 0 : if (TMath::Abs(p0O[3]-p1O[3])>fCutTheta) result&=kFALSE;
1317 0 : if (result==kTRUE){
1318 : result=kTRUE; // just to put break point here
1319 0 : }
1320 0 : return result;
1321 0 : }
1322 :
1323 :
1324 : void AliTPCcalibTime::ProcessSame(const AliESDtrack *const track, AliESDfriendTrack *const friendTrack, const AliESDEvent *const event){
1325 : //
1326 : // Process TPC tracks crossing CE
1327 : //
1328 : // 0. Select only track crossing the CE
1329 : // 1. Cut on the track length
1330 : // 2. Refit the the track on A and C side separatelly
1331 : // 3. Fill time histograms
1332 : const Int_t kMinNcl=100;
1333 : const Int_t kMinNclS=25; // minimul number of clusters on the sides
1334 0 : const Double_t pimass=TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1335 : const Double_t kMaxDy=1; // maximal distance in y
1336 : const Double_t kMaxDsnp=0.05; // maximal distance in snp
1337 : const Double_t kMaxDtheta=0.05; // maximal distance in theta
1338 :
1339 0 : if (!friendTrack->GetTPCOut()) return;
1340 : //
1341 : // 0. Select only track crossing the CE
1342 : //
1343 0 : if (track->GetInnerParam()->GetZ()*friendTrack->GetTPCOut()->GetZ()>0) return;
1344 : //
1345 : // 1. cut on track length
1346 : //
1347 0 : if (track->GetTPCNcls()<kMinNcl) return;
1348 : //
1349 : // 2. Refit track sepparatel on A and C side
1350 : //
1351 : TObject *calibObject;
1352 : AliTPCseed *seed = 0;
1353 0 : for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
1354 0 : if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
1355 : }
1356 0 : if (!seed) return;
1357 : //
1358 0 : AliExternalTrackParam trackIn(*track->GetInnerParam());
1359 0 : AliExternalTrackParam trackOut(*track->GetOuterParam());
1360 0 : Double_t cov[3]={0.01,0.,0.01}; //use the same errors
1361 0 : Double_t xyz[3]={0,0.,0.0};
1362 : Double_t bz =0;
1363 0 : Int_t nclIn=0,nclOut=0;
1364 0 : trackIn.ResetCovariance(1000.);
1365 0 : trackOut.ResetCovariance(1000.);
1366 : //
1367 : //2.a Refit inner
1368 : //
1369 0 : Int_t sideIn=0;
1370 0 : for (Int_t irow=0;irow<kMaxRow;irow++) {
1371 0 : AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
1372 0 : if (!cl) continue;
1373 0 : if (cl->GetX()<80) continue;
1374 0 : if (sideIn==0){
1375 0 : if (cl->GetDetector()%36<18) sideIn=1;
1376 0 : if (cl->GetDetector()%36>=18) sideIn=-1;
1377 : }
1378 0 : if (sideIn== -1 && (cl->GetDetector()%36)<18) break;
1379 0 : if (sideIn== 1 &&(cl->GetDetector()%36)>=18) break;
1380 0 : Int_t sector = cl->GetDetector();
1381 0 : Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
1382 0 : if (TMath::Abs(dalpha)>0.01){
1383 0 : if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
1384 : }
1385 0 : Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
1386 0 : trackIn.GetXYZ(xyz);
1387 0 : bz = AliTracker::GetBz(xyz);
1388 0 : AliTracker::PropagateTrackToBxByBz(&trackIn,r[0],pimass,1.,kFALSE);
1389 0 : if (!trackIn.PropagateTo(r[0],bz)) break;
1390 0 : nclIn++;
1391 0 : trackIn.Update(&r[1],cov);
1392 0 : }
1393 : //
1394 : //2.b Refit outer
1395 : //
1396 0 : Int_t sideOut=0;
1397 0 : for (Int_t irow=kMaxRow;irow--;) {
1398 0 : AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
1399 0 : if (!cl) continue;
1400 0 : if (cl->GetX()<80) continue;
1401 0 : if (sideOut==0){
1402 0 : if (cl->GetDetector()%36<18) sideOut=1;
1403 0 : if (cl->GetDetector()%36>=18) sideOut=-1;
1404 0 : if (sideIn==sideOut) break;
1405 : }
1406 0 : if (sideOut== -1 && (cl->GetDetector()%36)<18) break;
1407 0 : if (sideOut== 1 &&(cl->GetDetector()%36)>=18) break;
1408 : //
1409 0 : Int_t sector = cl->GetDetector();
1410 0 : Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
1411 0 : if (TMath::Abs(dalpha)>0.01){
1412 0 : if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
1413 : }
1414 0 : Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
1415 0 : trackOut.GetXYZ(xyz);
1416 0 : bz = AliTracker::GetBz(xyz);
1417 0 : AliTracker::PropagateTrackToBxByBz(&trackOut,r[0],pimass,1.,kFALSE);
1418 0 : if (!trackOut.PropagateTo(r[0],bz)) break;
1419 0 : nclOut++;
1420 0 : trackOut.Update(&r[1],cov);
1421 0 : }
1422 0 : trackOut.Rotate(trackIn.GetAlpha());
1423 0 : Double_t meanX = (trackIn.GetX()+trackOut.GetX())*0.5;
1424 0 : trackIn.PropagateTo(meanX,bz);
1425 0 : trackOut.PropagateTo(meanX,bz);
1426 0 : if (TMath::Abs(trackIn.GetY()-trackOut.GetY())>kMaxDy) return;
1427 0 : if (TMath::Abs(trackIn.GetSnp()-trackOut.GetSnp())>kMaxDsnp) return;
1428 0 : if (TMath::Abs(trackIn.GetTgl()-trackOut.GetTgl())>kMaxDtheta) return;
1429 0 : if (TMath::Min(nclIn,nclOut)>kMinNclS){
1430 0 : FillResHistoTPCCE(&trackIn,&trackOut);
1431 : }
1432 0 : TTreeSRedirector *cstream = GetDebugStreamer();
1433 0 : if (cstream){
1434 0 : TVectorD gxyz(3);
1435 0 : trackIn.GetXYZ(gxyz.GetMatrixArray());
1436 0 : TTimeStamp tstamp(fTime);
1437 0 : (*cstream)<<"tpctpc"<<
1438 0 : "run="<<fRun<< // run number
1439 0 : "event="<<fEvent<< // event number
1440 0 : "time="<<fTime<< // time stamp of event
1441 0 : "trigger="<<fTrigger<< // trigger
1442 0 : "mag="<<fMagF<< // magnetic field
1443 : //
1444 0 : "sideIn="<<sideIn<< // side at inner part
1445 0 : "sideOut="<<sideOut<< // side at puter part
1446 0 : "xyz.="<<&gxyz<< // global position
1447 0 : "tIn.="<<&trackIn<< // refitterd track in
1448 0 : "tOut.="<<&trackOut<< // refitter track out
1449 0 : "nclIn="<<nclIn<< //
1450 0 : "nclOut="<<nclOut<< //
1451 : "\n";
1452 0 : }
1453 : //
1454 : // 3. Fill time histograms
1455 : // Debug stremaer expression
1456 : // chainTPCTPC->Draw("(tIn.fP[1]-tOut.fP[1])*sign(-tIn.fP[3]):tIn.fP[3]","min(nclIn,nclOut)>30","")
1457 0 : if (TMath::Min(nclIn,nclOut)>kMinNclS){
1458 0 : fDz = trackOut.GetZ()-trackIn.GetZ();
1459 0 : if (trackOut.GetTgl()<0) fDz*=-1.;
1460 0 : TTimeStamp tstamp(fTime);
1461 0 : Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1462 0 : Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1463 0 : Double_t vecDrift[4]={static_cast<Double_t>(fTime),(ptrelative0+ptrelative1)/2.0,fDz/500.0,static_cast<Double_t>(event->GetRunNumber())};
1464 : //
1465 : // fill histograms per trigger class and itegrated
1466 : //
1467 : THnSparse* curHist=NULL;
1468 0 : for (Int_t itype=0; itype<2; itype++){
1469 0 : TString name="MEAN_VDRIFT_CROSS_";
1470 0 : if (itype==0){
1471 0 : name+=event->GetFiredTriggerClasses();
1472 0 : name.ToUpper();
1473 : }else{
1474 0 : name+="ALL";
1475 : }
1476 0 : curHist=(THnSparseF*)fArrayDz->FindObject(name);
1477 0 : if(!curHist){
1478 0 : curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
1479 0 : fArrayDz->AddLast(curHist);
1480 : }
1481 0 : curHist->Fill(vecDrift);
1482 0 : }
1483 0 : }
1484 :
1485 0 : }
1486 :
1487 : void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, const AliESDfriendTrack *const friendTrack, const AliESDEvent *const event, AliESDfriend *const esdFriend){
1488 : //
1489 : // Process track - Update TPC-ITS alignment
1490 : // Updates:
1491 : // 0. Apply standartd cuts
1492 : // 1. Recalucluate the current statistic median/RMS
1493 : // 2. Apply median+-rms cut
1494 : // 3. Update kalman filter
1495 : //
1496 : const Int_t kMinTPC = 80; // minimal number of TPC cluster
1497 : const Int_t kMinITS = 3; // minimal number of ITS cluster
1498 : const Double_t kMinZ = 10; // maximal dz distance
1499 : const Double_t kMaxDy = 2.; // maximal dy distance
1500 : const Double_t kMaxAngle= 0.07; // maximal angular distance
1501 : const Double_t kSigmaCut= 5; // maximal sigma distance to median
1502 : const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1503 : const Double_t kT0Err = 3.; // initial uncertainty of the T0 time
1504 : const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1505 : const Double_t kOutCut = 3.0; // outlyer cut in AliRelAlgnmentKalman
1506 : const Double_t kMinPt = 0.3; // minimal pt
1507 : const Double_t kMax1Pt=0.5; //maximal 1/pt distance
1508 : const Int_t kN=50; // deepnes of history
1509 : static Int_t kglast=0;
1510 0 : static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
1511 : //
1512 : // 0. Apply standard cuts
1513 : //
1514 0 : Int_t dummycl[1000];
1515 0 : if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
1516 0 : if (!track->IsOn(AliESDtrack::kTPCrefit)) return;
1517 0 : if (!track->GetInnerParam()) return;
1518 0 : if (!track->GetOuterParam()) return;
1519 0 : if (track->GetInnerParam()->Pt()<kMinPt) return;
1520 : // exclude crossing track
1521 0 : if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1522 0 : if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ/3.) return;
1523 0 : if (track->GetInnerParam()->GetX()>90) return;
1524 : //
1525 0 : AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetInnerParam()));
1526 : //
1527 0 : AliExternalTrackParam pITS; // ITS standalone if possible
1528 0 : AliExternalTrackParam pITS2; //TPC-ITS track
1529 0 : if (friendTrack->GetITSOut()){
1530 0 : pITS2=(*(friendTrack->GetITSOut())); //TPC-ITS track - snapshot ITS out
1531 0 : pITS2.Rotate(pTPC.GetAlpha());
1532 0 : AliTracker::PropagateTrackToBxByBz(&pITS2,pTPC.GetX(),0.1,0.1,kFALSE);
1533 : }
1534 :
1535 : AliESDfriendTrack *itsfriendTrack=0;
1536 : //
1537 : // try to find standalone ITS track corresponing to the TPC if possible
1538 : //
1539 0 : Bool_t hasAlone=kFALSE;
1540 0 : Int_t ntracks=event->GetNumberOfTracks();
1541 0 : for (Int_t i=0; i<ntracks; i++){
1542 0 : AliESDtrack * trackITS = event->GetTrack(i);
1543 0 : if (!trackITS) continue;
1544 0 : if (trackITS->GetITSclusters(dummycl)<kMinITS) continue; // minimal amount of clusters
1545 0 : itsfriendTrack = (AliESDfriendTrack*)trackITS->GetFriendTrack();
1546 0 : if (!itsfriendTrack) continue;
1547 0 : if (!itsfriendTrack->GetITSOut()) continue;
1548 :
1549 0 : if (TMath::Abs(pTPC.GetTgl()-itsfriendTrack->GetITSOut()->GetTgl())> kMaxAngle) continue;
1550 0 : if (TMath::Abs(pTPC.GetSigned1Pt()-itsfriendTrack->GetITSOut()->GetSigned1Pt())> kMax1Pt) continue;
1551 0 : pITS=(*(itsfriendTrack->GetITSOut()));
1552 : //
1553 0 : pITS.Rotate(pTPC.GetAlpha());
1554 0 : AliTracker::PropagateTrackToBxByBz(&pITS,pTPC.GetX(),0.1,0.1,kFALSE);
1555 0 : if (TMath::Abs(pTPC.GetY()-pITS.GetY())> kMaxDy) continue;
1556 0 : if (TMath::Abs(pTPC.GetSnp()-pITS.GetSnp())> kMaxAngle) continue;
1557 0 : hasAlone=kTRUE;
1558 0 : }
1559 0 : if (!hasAlone) {
1560 0 : if (track->GetITSclusters(dummycl)<kMinITS) return;
1561 0 : pITS=pITS2; // use combined track if it has ITS
1562 : }
1563 : //
1564 0 : if (TMath::Abs(pITS.GetY()-pTPC.GetY()) >kMaxDy) return;
1565 0 : if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1566 0 : if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1567 : //
1568 : // 1. Update median and RMS info
1569 : //
1570 0 : TVectorD vecDelta(5),vecMedian(5), vecRMS(5);
1571 0 : TVectorD vecDeltaN(5);
1572 0 : Double_t sign=(pITS.GetParameter()[1]>0)? 1.:-1.;
1573 0 : vecDelta[4]=0;
1574 0 : for (Int_t i=0;i<4;i++){
1575 0 : vecDelta[i]=(pITS.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1576 0 : kgdP[i][kglast%kN]=vecDelta[i];
1577 : }
1578 0 : kglast=(kglast+1);
1579 0 : Int_t entries=(kglast<kN)?kglast:kN;
1580 0 : for (Int_t i=0;i<4;i++){
1581 0 : vecMedian[i] = TMath::Median(entries,kgdP[i]);
1582 0 : vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1583 0 : vecDeltaN[i] = 0;
1584 0 : if (vecRMS[i]>0.){
1585 0 : vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1586 0 : vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1587 0 : }
1588 : }
1589 : //
1590 : // 2. Apply median+-rms cut
1591 : //
1592 0 : if (kglast<3) return; //median and RMS to be defined
1593 0 : if ( vecDeltaN[4]/4.>kSigmaCut) return;
1594 : //
1595 : // 3. Update alignment
1596 : //
1597 0 : Int_t htime = (fTime-fTimeKalmanBin/2)/fTimeKalmanBin; //time bins number
1598 0 : if (fAlignITSTPC->GetEntriesFast()<htime){
1599 0 : fAlignITSTPC->Expand(htime*2+20);
1600 : }
1601 0 : AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignITSTPC->At(htime);
1602 0 : if (!align){
1603 : // make Alignment object if doesn't exist
1604 0 : align=new AliRelAlignerKalman();
1605 0 : align->SetRunNumber(fRun);
1606 0 : (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1607 0 : (*align->GetStateCov())(7,7)=kT0Err*kT0Err;
1608 0 : (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1609 0 : align->SetOutRejSigma(kOutCut+kOutCut*kN);
1610 : // align->SetRejectOutliers(kFALSE);
1611 0 : align->SetRejectOutliers(kTRUE);
1612 0 : align->SetRejectOutliersSigma2Median(kTRUE);
1613 :
1614 0 : align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
1615 0 : align->SetMagField(fMagF);
1616 0 : fAlignITSTPC->AddAt(align,htime);
1617 : }
1618 0 : align->AddTrackParams(&pITS,&pTPC);
1619 0 : Double_t averageTime = fTime;
1620 0 : if (align->GetTimeStamp()>0&&align->GetNUpdates()>0){
1621 0 : averageTime=((Double_t(align->GetTimeStamp())*Double_t(align->GetNUpdates())+Double_t(fTime)))/(Double_t(align->GetNUpdates())+1.);
1622 0 : }
1623 0 : align->SetTimeStamp(Int_t(averageTime));
1624 :
1625 0 : align->SetRunNumber(fRun );
1626 0 : Float_t dca[2],cov[3];
1627 0 : track->GetImpactParameters(dca,cov);
1628 0 : if (TMath::Abs(dca[0])<kMaxDy){
1629 0 : FillResHistoTPCITS(&pTPC,&pITS);
1630 0 : FillResHistoTPC(track);
1631 : }
1632 : //
1633 0 : Int_t nupdates=align->GetNUpdates();
1634 0 : align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates+1));
1635 : // align->SetRejectOutliers(kFALSE);
1636 0 : align->SetRejectOutliers(kTRUE);
1637 0 : align->SetRejectOutliersSigma2Median(kTRUE);
1638 :
1639 0 : TTreeSRedirector *cstream = GetDebugStreamer();
1640 0 : if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1641 0 : TVectorD gpTPC(3), gdTPC(3);
1642 0 : TVectorD gpITS(3), gdITS(3);
1643 0 : pTPC.GetXYZ(gpTPC.GetMatrixArray());
1644 0 : pTPC.GetDirection(gdTPC.GetMatrixArray());
1645 0 : pITS.GetXYZ(gpITS.GetMatrixArray());
1646 0 : pITS.GetDirection(gdITS.GetMatrixArray());
1647 0 : (*cstream)<<"itstpc"<<
1648 0 : "run="<<fRun<< // run number
1649 0 : "event="<<fEvent<< // event number
1650 0 : "time="<<fTime<< // time stamp of event
1651 0 : "trigger="<<fTrigger<< // trigger
1652 0 : "mag="<<fMagF<< // magnetic field
1653 : //
1654 0 : "hasAlone="<<hasAlone<< // has ITS standalone ?
1655 0 : "track.="<<track<< // track info
1656 0 : "nmed="<<kglast<< // number of entries to define median and RMS
1657 0 : "vMed.="<<&vecMedian<< // median of deltas
1658 0 : "vRMS.="<<&vecRMS<< // rms of deltas
1659 0 : "vDelta.="<<&vecDelta<< // delta in respect to median
1660 0 : "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1661 0 : "a.="<<align<< // current alignment
1662 0 : "pITS.="<<&pITS<< // track param ITS
1663 0 : "pITS2.="<<&pITS2<< // track param ITS+TPC
1664 0 : "pTPC.="<<&pTPC<< // track param TPC
1665 0 : "gpTPC.="<<&gpTPC<< // global position TPC
1666 0 : "gdTPC.="<<&gdTPC<< // global direction TPC
1667 0 : "gpITS.="<<&gpITS<< // global position ITS
1668 0 : "gdITS.="<<&gdITS<< // global position ITS
1669 : "\n";
1670 0 : }
1671 0 : }
1672 :
1673 :
1674 :
1675 :
1676 : void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack *const track, const AliESDfriendTrack *const friendTrack){
1677 : //
1678 : // Process track - Update TPC-TRD alignment
1679 : // Updates:
1680 : // 0. Apply standartd cuts
1681 : // 1. Recalucluate the current statistic median/RMS
1682 : // 2. Apply median+-rms cut
1683 : // 3. Update kalman filter
1684 : //
1685 : const Int_t kMinTPC = 80; // minimal number of TPC cluster
1686 : const Int_t kMinTRD = 60; // minimal number of TRD cluster
1687 : // const Double_t kMinZ = 20; // maximal dz distance
1688 : const Double_t kMaxDy = 5.; // maximal dy distance
1689 : const Double_t kMaxAngle= 0.1; // maximal angular distance
1690 : const Double_t kSigmaCut= 10; // maximal sigma distance to median
1691 : const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1692 : const Double_t kT0Err = 3.; // initial uncertainty of the T0 time
1693 : const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1694 : const Double_t kOutCut = 3.0; // outlyer cut in AliRelAlgnmentKalman
1695 : const Double_t kRefX = 330; // reference X
1696 : const Int_t kN=50; // deepnes of history
1697 : static Int_t kglast=0;
1698 0 : static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
1699 : //
1700 : // 0. Apply standard cuts
1701 : //
1702 0 : Int_t dummycl[1000];
1703 0 : if (track->GetTRDclusters(dummycl)<kMinTRD) return; // minimal amount of clusters
1704 0 : if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
1705 : // if (!friendTrack->GetTRDIn()) return;
1706 : // if (!track->IsOn(AliESDtrack::kTRDrefit)) return;
1707 0 : if (!track->IsOn(AliESDtrack::kTRDout)) return;
1708 0 : if (!track->GetInnerParam()) return;
1709 0 : if (!friendTrack->GetTPCOut()) return;
1710 : // exclude crossing track
1711 0 : if (friendTrack->GetTPCOut()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1712 : //
1713 0 : AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(friendTrack->GetTPCOut()));
1714 0 : AliTracker::PropagateTrackToBxByBz(&pTPC,kRefX,0.1,0.1,kFALSE);
1715 : AliExternalTrackParam *pTRDtrack = 0;
1716 : TObject *calibObject=0;
1717 0 : for (Int_t l=0;(calibObject=((AliESDfriendTrack*)friendTrack)->GetCalibObject(l));++l) {
1718 0 : if ((dynamic_cast< AliTPCseed*>(calibObject))) continue;
1719 0 : if ((pTRDtrack=dynamic_cast< AliExternalTrackParam*>(calibObject))) break;
1720 : }
1721 0 : if (!pTRDtrack) return;
1722 : // AliExternalTrackParam pTRD(*(friendTrack->GetTRDIn()));
1723 0 : AliExternalTrackParam pTRD(*(pTRDtrack));
1724 0 : pTRD.Rotate(pTPC.GetAlpha());
1725 : // pTRD.PropagateTo(pTPC.GetX(),fMagF);
1726 0 : AliTracker::PropagateTrackToBxByBz(&pTRD,pTPC.GetX(),0.1,0.1,kFALSE);
1727 :
1728 0 : ((Double_t*)pTRD.GetCovariance())[2]+=3.*3.; // increas sys errors
1729 0 : ((Double_t*)pTRD.GetCovariance())[9]+=0.1*0.1; // increse sys errors
1730 :
1731 0 : if (TMath::Abs(pTRD.GetY()-pTPC.GetY()) >kMaxDy) return;
1732 0 : if (TMath::Abs(pTRD.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1733 : // if (TMath::Abs(pTRD.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1734 : //
1735 : // 1. Update median and RMS info
1736 : //
1737 0 : TVectorD vecDelta(5),vecMedian(5), vecRMS(5);
1738 0 : TVectorD vecDeltaN(5);
1739 0 : Double_t sign=(pTRD.GetParameter()[1]>0)? 1.:-1.;
1740 0 : vecDelta[4]=0;
1741 0 : for (Int_t i=0;i<4;i++){
1742 0 : vecDelta[i]=(pTRD.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1743 0 : kgdP[i][kglast%kN]=vecDelta[i];
1744 : }
1745 0 : kglast=(kglast+1);
1746 0 : Int_t entries=(kglast<kN)?kglast:kN;
1747 0 : for (Int_t i=0;i<4;i++){
1748 0 : vecMedian[i] = TMath::Median(entries,kgdP[i]);
1749 :
1750 0 : vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1751 0 : vecDeltaN[i] = 0;
1752 0 : if (vecRMS[i]>0.){
1753 0 : vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1754 0 : vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1755 0 : }
1756 : }
1757 : //
1758 : // 2. Apply median+-rms cut
1759 : //
1760 0 : if (kglast<3) return; //median and RMS to be defined
1761 0 : if ( vecDeltaN[4]/4.>kSigmaCut) return;
1762 : //
1763 : // 3. Update alignment
1764 : //
1765 : //Int_t htime = fTime/3600; //time in hours
1766 0 : Int_t htime = (Int_t)(fTime-fTimeKalmanBin/2)/fTimeKalmanBin; //time in half hour
1767 0 : if (fAlignTRDTPC->GetEntriesFast()<htime){
1768 0 : fAlignTRDTPC->Expand(htime*2+20);
1769 : }
1770 0 : AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTRDTPC->At(htime);
1771 0 : if (!align){
1772 : // make Alignment object if doesn't exist
1773 0 : align=new AliRelAlignerKalman();
1774 0 : align->SetRunNumber(fRun);
1775 0 : (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1776 0 : (*align->GetStateCov())(7,7)=kT0Err*kT0Err;
1777 0 : (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1778 0 : align->SetOutRejSigma(kOutCut+kOutCut*kN);
1779 : // align->SetRejectOutliers(kFALSE);
1780 0 : align->SetRejectOutliers(kTRUE);
1781 0 : align->SetRejectOutliersSigma2Median(kTRUE);
1782 :
1783 0 : align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
1784 0 : align->SetMagField(fMagF);
1785 0 : fAlignTRDTPC->AddAt(align,htime);
1786 : }
1787 0 : align->AddTrackParams(&pTRD,&pTPC);
1788 : //align->SetTimeStamp(fTime);
1789 0 : Double_t averageTime = fTime;
1790 0 : if (align->GetTimeStamp()>0 && align->GetNUpdates()>0) {
1791 0 : averageTime = (((Double_t)fTime) + ((Double_t)align->GetTimeStamp())*align->GetNUpdates()) / (align->GetNUpdates() + 1.);
1792 : //printf("align->GetTimeStamp() %d, align->GetNUpdates() %d \n", align->GetTimeStamp(), align->GetNUpdates());
1793 0 : }
1794 0 : align->SetTimeStamp((Int_t)averageTime);
1795 :
1796 : //printf("fTime %d, averageTime %d \n", fTime, (Int_t)averageTime);
1797 :
1798 0 : align->SetRunNumber(fRun );
1799 0 : Float_t dca[2],cov[3];
1800 0 : track->GetImpactParameters(dca,cov);
1801 0 : if (TMath::Abs(dca[0])<kMaxDy){
1802 0 : FillResHistoTPCTRD(&pTPC,&pTRD); //only primaries
1803 : }
1804 : //
1805 0 : Int_t nupdates=align->GetNUpdates();
1806 0 : align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates+1));
1807 : // align->SetRejectOutliers(kFALSE);
1808 0 : align->SetRejectOutliers(kTRUE);
1809 0 : align->SetRejectOutliersSigma2Median(kTRUE);
1810 :
1811 0 : TTreeSRedirector *cstream = GetDebugStreamer();
1812 0 : if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1813 0 : TVectorD gpTPC(3), gdTPC(3);
1814 0 : TVectorD gpTRD(3), gdTRD(3);
1815 0 : pTPC.GetXYZ(gpTPC.GetMatrixArray());
1816 0 : pTPC.GetDirection(gdTPC.GetMatrixArray());
1817 0 : pTRD.GetXYZ(gpTRD.GetMatrixArray());
1818 0 : pTRD.GetDirection(gdTRD.GetMatrixArray());
1819 0 : (*cstream)<<"trdtpc"<<
1820 0 : "run="<<fRun<< // run number
1821 0 : "event="<<fEvent<< // event number
1822 0 : "time="<<fTime<< // time stamp of event
1823 0 : "trigger="<<fTrigger<< // trigger
1824 0 : "mag="<<fMagF<< // magnetic field
1825 : //
1826 0 : "nmed="<<kglast<< // number of entries to define median and RMS
1827 0 : "vMed.="<<&vecMedian<< // median of deltas
1828 0 : "vRMS.="<<&vecRMS<< // rms of deltas
1829 0 : "vDelta.="<<&vecDelta<< // delta in respect to median
1830 0 : "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1831 0 : "t.="<<track<< // ful track - find proper cuts
1832 0 : "a.="<<align<< // current alignment
1833 0 : "pTRD.="<<&pTRD<< // track param TRD
1834 0 : "pTPC.="<<&pTPC<< // track param TPC
1835 0 : "gpTPC.="<<&gpTPC<< // global position TPC
1836 0 : "gdTPC.="<<&gdTPC<< // global direction TPC
1837 0 : "gpTRD.="<<&gpTRD<< // global position TRD
1838 0 : "gdTRD.="<<&gdTRD<< // global position TRD
1839 : "\n";
1840 0 : }
1841 0 : }
1842 :
1843 :
1844 : void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack *const track, const AliESDfriendTrack *const friendTrack){
1845 : //
1846 : //
1847 : // Process track - Update TPC-TOF alignment
1848 : // Updates:
1849 : // -1. Make a TOF "track"
1850 : // 0. Apply standartd cuts
1851 : // 1. Recalucluate the current statistic median/RMS
1852 : // 2. Apply median+-rms cut
1853 : // 3. Update kalman filter
1854 : //
1855 : const Int_t kMinTPC = 80; // minimal number of TPC cluster
1856 : // const Double_t kMinZ = 10; // maximal dz distance
1857 : const Double_t kMaxDy = 5.; // maximal dy distance
1858 : const Double_t kMaxAngle= 0.05; // maximal angular distance
1859 : const Double_t kSigmaCut= 5; // maximal sigma distance to median
1860 : const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1861 : const Double_t kT0Err = 3.; // initial uncertainty of the T0 time
1862 : const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1863 :
1864 : const Double_t kOutCut = 3.0; // outlyer cut in AliRelAlgnmentKalman
1865 : const Int_t kN=50; // deepnes of history
1866 : static Int_t kglast=0;
1867 0 : static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
1868 : //
1869 : // -1. Make a TOF track-
1870 : // Clusters are not in friends - use alingment points
1871 : //
1872 0 : if (track->GetTOFsignal()<=0) return;
1873 0 : if (!friendTrack->GetTPCOut()) return;
1874 0 : if (!track->GetInnerParam()) return;
1875 0 : if (!friendTrack->GetTPCOut()) return;
1876 0 : const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
1877 0 : if (!points) return;
1878 0 : AliExternalTrackParam pTPC(*(friendTrack->GetTPCOut()));
1879 0 : AliExternalTrackParam pTOF(pTPC);
1880 0 : Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
1881 0 : Int_t npoints = points->GetNPoints();
1882 0 : AliTrackPoint point;
1883 : Int_t naccept=0;
1884 : //
1885 0 : for (Int_t ipoint=0;ipoint<npoints;ipoint++){
1886 0 : points->GetPoint(point,ipoint);
1887 0 : Float_t xyz[3];
1888 0 : point.GetXYZ(xyz);
1889 0 : Double_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1890 0 : if (r<370) continue; // I will be happy to use flag in case it will be implemented in AliTrackPoint - should be requested
1891 0 : if (r>400) continue;
1892 0 : AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,2.,kTRUE);
1893 0 : AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,0.1,kTRUE);
1894 0 : AliTrackPoint lpoint = point.Rotate(pTPC.GetAlpha());
1895 0 : pTPC.PropagateTo(lpoint.GetX(),fMagF);
1896 0 : pTOF=pTPC;
1897 0 : ((Double_t*)pTOF.GetParameter())[0] =lpoint.GetY();
1898 0 : ((Double_t*)pTOF.GetParameter())[1] =lpoint.GetZ();
1899 0 : ((Double_t*)pTOF.GetCovariance())[0]+=3.*3./12.;
1900 0 : ((Double_t*)pTOF.GetCovariance())[2]+=3.*3./12.;
1901 0 : ((Double_t*)pTOF.GetCovariance())[5]+=0.1*0.1;
1902 0 : ((Double_t*)pTOF.GetCovariance())[9]+=0.1*0.1;
1903 0 : naccept++;
1904 0 : }
1905 0 : if (naccept==0) return; // no tof match clusters
1906 : //
1907 : // 0. Apply standard cuts
1908 : //
1909 0 : if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
1910 : // exclude crossing track
1911 0 : if (friendTrack->GetTPCOut()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1912 : //
1913 0 : if (TMath::Abs(pTOF.GetY()-pTPC.GetY()) >kMaxDy) return;
1914 0 : if (TMath::Abs(pTOF.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1915 0 : if (TMath::Abs(pTOF.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1916 : //
1917 : // 1. Update median and RMS info
1918 : //
1919 0 : TVectorD vecDelta(5),vecMedian(5), vecRMS(5);
1920 0 : TVectorD vecDeltaN(5);
1921 0 : Double_t sign=(pTOF.GetParameter()[1]>0)? 1.:-1.;
1922 0 : vecDelta[4]=0;
1923 0 : for (Int_t i=0;i<4;i++){
1924 0 : vecDelta[i]=(pTOF.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1925 0 : kgdP[i][kglast%kN]=vecDelta[i];
1926 : }
1927 0 : kglast=(kglast+1);
1928 0 : Int_t entries=(kglast<kN)?kglast:kN;
1929 : Bool_t isOK=kTRUE;
1930 0 : for (Int_t i=0;i<4;i++){
1931 0 : vecMedian[i] = TMath::Median(entries,kgdP[i]);
1932 0 : vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1933 0 : vecDeltaN[i] = 0;
1934 0 : if (vecRMS[i]>0.){
1935 0 : vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/(vecRMS[i]+1.);
1936 0 : vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1937 0 : if (TMath::Abs(vecDeltaN[i])>kSigmaCut) isOK=kFALSE;
1938 : }
1939 : }
1940 : //
1941 : // 2. Apply median+-rms cut
1942 : //
1943 0 : if (kglast<10) return; //median and RMS to be defined
1944 0 : if (!isOK) return;
1945 : //
1946 : // 3. Update alignment
1947 : //
1948 : //Int_t htime = fTime/3600; //time in hours
1949 0 : Int_t htime = (Int_t)(fTime-fTimeKalmanBin)/fTimeKalmanBin; //time bin
1950 0 : if (fAlignTOFTPC->GetEntriesFast()<htime){
1951 0 : fAlignTOFTPC->Expand(htime*2+20);
1952 : }
1953 0 : AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTOFTPC->At(htime);
1954 0 : if (!align){
1955 : // make Alignment object if doesn't exist
1956 0 : align=new AliRelAlignerKalman();
1957 0 : align->SetRunNumber(fRun);
1958 0 : (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1959 0 : (*align->GetStateCov())(7,7)=kT0Err*kT0Err;
1960 0 : (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1961 0 : align->SetOutRejSigma(kOutCut+kOutCut*kN);
1962 : // align->SetRejectOutliers(kFALSE);
1963 0 : align->SetRejectOutliers(kTRUE);
1964 0 : align->SetRejectOutliersSigma2Median(kTRUE);
1965 :
1966 0 : align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
1967 0 : align->SetMagField(fMagF);
1968 0 : fAlignTOFTPC->AddAt(align,htime);
1969 : }
1970 0 : align->AddTrackParams(&pTOF,&pTPC);
1971 0 : Float_t dca[2],cov[3];
1972 0 : track->GetImpactParameters(dca,cov);
1973 0 : if (TMath::Abs(dca[0])<kMaxDy){
1974 0 : FillResHistoTPCTOF(&pTPC,&pTOF);
1975 : }
1976 : //align->SetTimeStamp(fTime);
1977 0 : Double_t averageTime = fTime;
1978 0 : if (align->GetTimeStamp()>0 && align->GetNUpdates()>0) {
1979 0 : averageTime = (((Double_t)fTime) + ((Double_t)align->GetTimeStamp())*align->GetNUpdates()) / (align->GetNUpdates() + 1.);
1980 : //printf("align->GetTimeStamp() %d, align->GetNUpdates() %d \n", align->GetTimeStamp(), align->GetNUpdates());
1981 0 : }
1982 0 : align->SetTimeStamp((Int_t)averageTime);
1983 :
1984 : //printf("fTime %d, averageTime %d \n", fTime, (Int_t)averageTime);
1985 :
1986 0 : align->SetRunNumber(fRun );
1987 : //
1988 0 : Int_t nupdates=align->GetNUpdates();
1989 0 : align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates+1));
1990 : // align->SetRejectOutliers(kFALSE);
1991 0 : align->SetRejectOutliers(kTRUE);
1992 0 : align->SetRejectOutliersSigma2Median(kTRUE);
1993 :
1994 0 : TTreeSRedirector *cstream = GetDebugStreamer();
1995 0 : if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1996 0 : TVectorD gpTPC(3), gdTPC(3);
1997 0 : TVectorD gpTOF(3), gdTOF(3);
1998 0 : pTPC.GetXYZ(gpTPC.GetMatrixArray());
1999 0 : pTPC.GetDirection(gdTPC.GetMatrixArray());
2000 0 : pTOF.GetXYZ(gpTOF.GetMatrixArray());
2001 0 : pTOF.GetDirection(gdTOF.GetMatrixArray());
2002 0 : (*cstream)<<"toftpc"<<
2003 0 : "run="<<fRun<< // run number
2004 0 : "event="<<fEvent<< // event number
2005 0 : "time="<<fTime<< // time stamp of event
2006 0 : "trigger="<<fTrigger<< // trigger
2007 0 : "mag="<<fMagF<< // magnetic field
2008 : //
2009 0 : "nmed="<<kglast<< // number of entries to define median and RMS
2010 0 : "vMed.="<<&vecMedian<< // median of deltas
2011 0 : "vRMS.="<<&vecRMS<< // rms of deltas
2012 0 : "vDelta.="<<&vecDelta<< // delta in respect to median
2013 0 : "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
2014 0 : "t.="<<track<< // ful track - find proper cuts
2015 0 : "a.="<<align<< // current alignment
2016 0 : "pTOF.="<<&pTOF<< // track param TOF
2017 0 : "pTPC.="<<&pTPC<< // track param TPC
2018 0 : "gpTPC.="<<&gpTPC<< // global position TPC
2019 0 : "gdTPC.="<<&gdTPC<< // global direction TPC
2020 0 : "gpTOF.="<<&gpTOF<< // global position TOF
2021 0 : "gdTOF.="<<&gdTOF<< // global position TOF
2022 : "\n";
2023 0 : }
2024 0 : }
2025 :
2026 :
2027 : void AliTPCcalibTime::BookDistortionMaps(){
2028 : //
2029 : // Book ndimensional histograms of distortions/residuals
2030 : // Only primary tracks are selected for analysis
2031 : //
2032 :
2033 0 : Double_t xminTrack[5], xmaxTrack[5];
2034 0 : Int_t binsTrack[5];
2035 0 : TString axisName[5];
2036 0 : TString axisTitle[5];
2037 : //
2038 0 : binsTrack[0] =50;
2039 0 : axisName[0] ="#Delta";
2040 0 : axisTitle[0] ="#Delta";
2041 : //
2042 0 : binsTrack[1] =44;
2043 0 : xminTrack[1] =-1.1; xmaxTrack[1]=1.1;
2044 0 : axisName[1] ="tanTheta";
2045 0 : axisTitle[1] ="tan(#Theta)";
2046 : //
2047 0 : binsTrack[2] =180;
2048 0 : xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi();
2049 0 : axisName[2] ="phi";
2050 0 : axisTitle[2] ="#phi";
2051 : //
2052 0 : binsTrack[3] =20;
2053 0 : xminTrack[3] =-1.; xmaxTrack[3]=1.; // 0.33 GeV cut
2054 0 : axisName[3] ="snp";
2055 0 : axisTitle[3] ="snp";
2056 : //
2057 0 : binsTrack[4] =10;
2058 0 : xminTrack[4] =120.; xmaxTrack[4]=215.; // crossing radius for CE only
2059 0 : axisName[4] ="r";
2060 0 : axisTitle[4] ="r(cm)";
2061 : //
2062 : // delta y
2063 0 : xminTrack[0] =-1.5; xmaxTrack[0]=1.5; //
2064 0 : fResHistoTPCCE[0] = new THnSparseS("TPCCE#Delta_{Y} (cm)","#Delta_{Y} (cm)", 5, binsTrack,xminTrack, xmaxTrack);
2065 0 : fResHistoTPCITS[0] = new THnSparseS("TPCITS#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
2066 0 : fResHistoTPCvertex[0] = new THnSparseS("TPCVertex#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
2067 0 : xminTrack[0] =-1.5; xmaxTrack[0]=1.5; //
2068 0 : fResHistoTPCTRD[0] = new THnSparseS("TPCTRD#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
2069 0 : xminTrack[0] =-5; xmaxTrack[0]=5; //
2070 0 : fResHistoTPCTOF[0] = new THnSparseS("TPCTOF#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
2071 : //
2072 : // delta z
2073 0 : xminTrack[0] =-6.; xmaxTrack[0]=6.; //
2074 0 : fResHistoTPCCE[1] = new THnSparseS("TPCCE#Delta_{Z} (cm)","#Delta_{Z} (cm)", 5, binsTrack,xminTrack, xmaxTrack);
2075 0 : fResHistoTPCITS[1] = new THnSparseS("TPCITS#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
2076 0 : fResHistoTPCvertex[1] = new THnSparseS("TPCVertex#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
2077 0 : fResHistoTPCTRD[1] = new THnSparseS("TPCTRD#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
2078 0 : xminTrack[0] =-5.; xmaxTrack[0]=5.; //
2079 0 : fResHistoTPCTOF[1] = new THnSparseS("TPCTOF#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
2080 : //
2081 : // delta snp-P2
2082 0 : xminTrack[0] =-0.015; xmaxTrack[0]=0.015; //
2083 0 : fResHistoTPCCE[2] = new THnSparseS("TPCCE#Delta_{#phi}","#Delta_{#phi}", 5, binsTrack,xminTrack, xmaxTrack);
2084 0 : fResHistoTPCITS[2] = new THnSparseS("TPCITS#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
2085 0 : fResHistoTPCvertex[2] = new THnSparseS("TPCITSv#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
2086 0 : fResHistoTPCTRD[2] = new THnSparseS("TPCTRD#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
2087 0 : fResHistoTPCTOF[2] = new THnSparseS("TPCTOF#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
2088 : //
2089 : // delta theta-P3
2090 0 : xminTrack[0] =-0.05; xmaxTrack[0]=0.05; //
2091 0 : fResHistoTPCCE[3] = new THnSparseS("TPCCE#Delta_{#theta}","#Delta_{#theta}", 5, binsTrack,xminTrack, xmaxTrack);
2092 0 : fResHistoTPCITS[3] = new THnSparseS("TPCITS#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
2093 0 : fResHistoTPCvertex[3] = new THnSparseS("TPCITSv#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
2094 0 : fResHistoTPCTRD[3] = new THnSparseS("TPCTRD#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
2095 0 : fResHistoTPCTOF[3] = new THnSparseS("TPCTOF#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
2096 : //
2097 : // delta theta-P4
2098 0 : xminTrack[0] =-0.2; xmaxTrack[0]=0.2; //
2099 0 : fResHistoTPCCE[4] = new THnSparseS("TPCCE#Delta_{1/pt}","#Delta_{1/pt}", 5, binsTrack,xminTrack, xmaxTrack);
2100 0 : fResHistoTPCITS[4] = new THnSparseS("TPCITS#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
2101 0 : fResHistoTPCvertex[4] = new THnSparseS("TPCITSv#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
2102 0 : fResHistoTPCTRD[4] = new THnSparseS("TPCTRD#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
2103 0 : fResHistoTPCTOF[4] = new THnSparseS("TPCTOF#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
2104 : //
2105 0 : for (Int_t ivar=0;ivar<4;ivar++){
2106 0 : for (Int_t ivar2=0;ivar2<5;ivar2++){
2107 0 : fResHistoTPCCE[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
2108 0 : fResHistoTPCCE[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
2109 0 : if (ivar2<4){
2110 0 : fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
2111 0 : fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
2112 0 : fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
2113 0 : fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
2114 0 : fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
2115 0 : fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
2116 : }
2117 : }
2118 : }
2119 : //
2120 : // Book vertex: time histograms
2121 : //
2122 0 : Int_t binsVertex[2]={500, fTimeBins};
2123 0 : Double_t aminVertex[2]={-5,fTimeStart};
2124 0 : Double_t amaxVertex[2]={5, fTimeEnd};
2125 0 : const char* hnames[12]={"TPCXAside", "TPCXCside","TPCXACdiff","TPCXAPCdiff",
2126 : "TPCYAside", "TPCYCside","TPCYACdiff","TPCYAPCdiff",
2127 : "TPCZAPCside", "TPCZAMCside","TPCZACdiff","TPCZAPCdiff"};
2128 0 : const char* anames[12]={"x (cm) - A side ", "x (cm) - C side","#Delta_{x} (cm) - TPC-A-C","#Delta_{x} (cm) - TPC-Common",
2129 : "y (cm) - A side ", "y (cm) - C side","#Delta_{x} (cm) - TPC-A-C","#Delta_{y} (cm) - TPC-Common",
2130 : "z (cm)", "#Delta_{Z} (cm) A-C side","#Delta_{x} (cm) - TPC-A-C","#Delta_{Z} (cm) TPC-common"};
2131 0 : for (Int_t ihis=0; ihis<12; ihis++) {
2132 0 : if (ihis>=8) aminVertex[0]=-20.;
2133 0 : if (ihis>=8) amaxVertex[0]=20.;
2134 0 : fTPCVertex[ihis]=new THnSparseF(hnames[ihis],hnames[ihis],2,binsVertex,aminVertex,amaxVertex);
2135 0 : fTPCVertex[ihis]->GetAxis(1)->SetTitle("Time");
2136 0 : fTPCVertex[ihis]->GetAxis(0)->SetTitle(anames[ihis]);
2137 : }
2138 :
2139 0 : Int_t binsVertexC[2]={40, 300};
2140 0 : Double_t aminVertexC[2]={-20,-30};
2141 0 : Double_t amaxVertexC[2]={20,30};
2142 0 : const char* hnamesC[5]={"TPCA_TPC","TPCC_TPC","TPCA_ITS","TPCC_ITS","TPC_ITS"};
2143 0 : for (Int_t ihis=0; ihis<5; ihis++) {
2144 0 : fTPCVertexCorrelation[ihis]=new THnSparseF(hnamesC[ihis],hnamesC[ihis],2,binsVertexC,aminVertexC,amaxVertexC);
2145 0 : fTPCVertexCorrelation[ihis]->GetAxis(1)->SetTitle("z (cm)");
2146 0 : fTPCVertexCorrelation[ihis]->GetAxis(0)->SetTitle("z (cm)");
2147 : }
2148 0 : }
2149 :
2150 :
2151 : void AliTPCcalibTime::FillResHistoTPCCE(const AliExternalTrackParam * pTPCIn, const AliExternalTrackParam * pTPCOut ){
2152 : //
2153 : // fill residual histograms pTPCOut-pTPCin - trac crossing CE
2154 : // Histogram
2155 : //
2156 0 : if (fMemoryMode<2) return;
2157 0 : Double_t histoX[5];
2158 0 : Double_t xyz[3];
2159 0 : pTPCIn->GetXYZ(xyz);
2160 0 : Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
2161 0 : histoX[1]= pTPCIn->GetTgl();
2162 0 : histoX[2]= phi;
2163 0 : histoX[3]= pTPCIn->GetSnp();
2164 0 : histoX[4]= pTPCIn->GetX();
2165 0 : AliExternalTrackParam lout(*pTPCOut);
2166 0 : lout.Rotate(pTPCIn->GetAlpha());
2167 0 : lout.PropagateTo(pTPCIn->GetX(),fMagF);
2168 : //
2169 0 : for (Int_t ihisto=0; ihisto<5; ihisto++){
2170 0 : histoX[0]=lout.GetParameter()[ihisto]-pTPCIn->GetParameter()[ihisto];
2171 0 : fResHistoTPCCE[ihisto]->Fill(histoX);
2172 : }
2173 0 : }
2174 : void AliTPCcalibTime::FillResHistoTPCITS(const AliExternalTrackParam * pTPCIn, const AliExternalTrackParam * pITSOut ){
2175 : //
2176 : // fill residual histograms pTPCIn-pITSOut
2177 : // Histogram is filled only for primary tracks
2178 : //
2179 0 : Double_t histoX[4];
2180 0 : Double_t xyz[3];
2181 0 : pTPCIn->GetXYZ(xyz);
2182 0 : Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
2183 0 : histoX[1]= pTPCIn->GetTgl();
2184 0 : histoX[2]= phi;
2185 0 : histoX[3]= pTPCIn->GetSnp();
2186 0 : AliExternalTrackParam lits(*pITSOut);
2187 0 : lits.Rotate(pTPCIn->GetAlpha());
2188 0 : lits.PropagateTo(pTPCIn->GetX(),fMagF);
2189 : //
2190 0 : for (Int_t ihisto=0; ihisto<5; ihisto++){
2191 0 : histoX[0]=pTPCIn->GetParameter()[ihisto]-lits.GetParameter()[ihisto];
2192 0 : fResHistoTPCITS[ihisto]->Fill(histoX);
2193 : }
2194 0 : }
2195 :
2196 :
2197 : void AliTPCcalibTime::FillResHistoTPC(const AliESDtrack * pTrack){
2198 : //
2199 : // fill residual histograms pTPC - vertex
2200 : // Histogram is filled only for primary tracks
2201 : //
2202 0 : if (fMemoryMode<2) return;
2203 0 : Double_t histoX[4];
2204 0 : const AliExternalTrackParam * pTPCIn = pTrack->GetInnerParam();
2205 0 : AliExternalTrackParam pTPCvertex(*(pTrack->GetInnerParam()));
2206 : //
2207 0 : if (!(pTrack->GetConstrainedParam())) return;
2208 0 : AliExternalTrackParam lits(*(pTrack->GetConstrainedParam()));
2209 0 : if (TMath::Abs(pTrack->GetY())>3) return; // beam pipe
2210 0 : pTPCvertex.Rotate(lits.GetAlpha());
2211 : //pTPCvertex.PropagateTo(pTPCvertex->GetX(),fMagF);
2212 0 : AliTracker::PropagateTrackToBxByBz(&pTPCvertex,lits.GetX(),0.1,2,kFALSE);
2213 0 : AliTracker::PropagateTrackToBxByBz(&pTPCvertex,lits.GetX(),0.1,0.1,kFALSE);
2214 0 : Double_t xyz[3];
2215 0 : pTPCIn->GetXYZ(xyz);
2216 0 : Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
2217 0 : histoX[1]= pTPCIn->GetTgl();
2218 0 : histoX[2]= phi;
2219 0 : histoX[3]= pTPCIn->GetSnp();
2220 : //
2221 0 : Float_t dca[2], cov[3];
2222 0 : pTrack->GetImpactParametersTPC(dca,cov);
2223 0 : for (Int_t ihisto=0; ihisto<5; ihisto++){
2224 0 : histoX[0]=pTPCvertex.GetParameter()[ihisto]-lits.GetParameter()[ihisto];
2225 : // if (ihisto<2) histoX[0]=dca[ihisto];
2226 0 : fResHistoTPCvertex[ihisto]->Fill(histoX);
2227 : }
2228 0 : }
2229 :
2230 :
2231 : void AliTPCcalibTime::FillResHistoTPCTRD(const AliExternalTrackParam * pTPCOut, const AliExternalTrackParam * pTRDIn ){
2232 : //
2233 : // fill resuidual histogram TPCout-TRDin
2234 : //
2235 0 : if (fMemoryMode<2) return;
2236 0 : Double_t histoX[4];
2237 0 : Double_t xyz[3];
2238 0 : pTPCOut->GetXYZ(xyz);
2239 0 : Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
2240 0 : histoX[1]= pTPCOut->GetTgl();
2241 0 : histoX[2]= phi;
2242 0 : histoX[3]= pTPCOut->GetSnp();
2243 : //
2244 0 : AliExternalTrackParam ltrd(*pTRDIn);
2245 0 : ltrd.Rotate(pTPCOut->GetAlpha());
2246 : // ltrd.PropagateTo(pTPCOut->GetX(),fMagF);
2247 0 : AliTracker::PropagateTrackToBxByBz(<rd,pTPCOut->GetX(),0.1,0.1,kFALSE);
2248 :
2249 0 : for (Int_t ihisto=0; ihisto<5; ihisto++){
2250 0 : histoX[0]=pTPCOut->GetParameter()[ihisto]-ltrd.GetParameter()[ihisto];
2251 0 : fResHistoTPCTRD[ihisto]->Fill(histoX);
2252 : }
2253 :
2254 0 : }
2255 :
2256 : void AliTPCcalibTime::FillResHistoTPCTOF(const AliExternalTrackParam * pTPCOut, const AliExternalTrackParam * pTOFIn ){
2257 : //
2258 : // fill resuidual histogram TPCout-TOFin
2259 : // track propagated to the TOF position
2260 0 : if (fMemoryMode<2) return;
2261 0 : Double_t histoX[4];
2262 0 : Double_t xyz[3];
2263 :
2264 0 : AliExternalTrackParam ltpc(*pTPCOut);
2265 0 : ltpc.Rotate(pTOFIn->GetAlpha());
2266 0 : AliTracker::PropagateTrackToBxByBz(<pc,pTOFIn->GetX(),0.1,0.1,kFALSE);
2267 : //
2268 0 : ltpc.GetXYZ(xyz);
2269 0 : Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
2270 0 : histoX[1]= ltpc.GetTgl();
2271 0 : histoX[2]= phi;
2272 0 : histoX[3]= ltpc.GetSnp();
2273 : //
2274 0 : for (Int_t ihisto=0; ihisto<2; ihisto++){
2275 0 : histoX[0]=ltpc.GetParameter()[ihisto]-pTOFIn->GetParameter()[ihisto];
2276 0 : fResHistoTPCTOF[ihisto]->Fill(histoX);
2277 : }
2278 :
2279 0 : }
|
