Add computation of dE/dx (truncated mean) after inward fit for ITS standalone tracks...
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTime.cxx
CommitLineData
c74c5f6c 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/*
74235403 17Comments to be written here:
c74c5f6c 18
74235403 191. What do we calibrate.
20
21 Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
c74c5f6c 22
23 AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
24
74235403 252. How to interpret results
26
273. Simple example
c74c5f6c 28
74235403 29 a) determine the required time range:
c74c5f6c 30
74235403 31 AliXRDPROOFtoolkit tool;
32 TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
33 chain->Draw("GetTimeStamp()")
c74c5f6c 34
74235403 35 b) analyse calibration object on Proof in calibration train
c74c5f6c 36
74235403 37 AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift);
c74c5f6c 38
74235403 39 c) plot results
40 .x ~/NimStyle.C
2be25cc9 41 gSystem->Load("libANALYSIS");
42 gSystem->Load("libTPCcalib");
c74c5f6c 43
dde68d36 44 TFile f("CalibObjectsTrain1.root");
45 AliTPCcalibTime *calib = (AliTPCcalibTime *)f->Get("calibTime");
46 calib->GetHistoDrift("all")->Projection(2,0)->Draw()
47 calib->GetFitDrift("all")->Draw("lp")
da6c0bc9 48
74235403 494. Analysis using debug streamers.
da6c0bc9 50
74235403 51 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
52 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
53 AliXRDPROOFtoolkit tool;
3e55050f 54 TChain * chainTime = tool.MakeChainRandom("time.txt","trackInfo",0,10000);
817766d5 55
3e55050f 56 AliXRDPROOFtoolkit::FilterList("timetpctpc.txt","* tpctpc",1)
57 AliXRDPROOFtoolkit::FilterList("timetoftpc.txt","* toftpc",1)
58 AliXRDPROOFtoolkit::FilterList("timeitstpc.txt","* itstpc",1)
59 AliXRDPROOFtoolkit::FilterList("timelaser.txt","* laserInfo",1)
60 TChain * chainTPCTPC = tool.MakeChainRandom("timetpctpc.txt.Good","tpctpc",0,10000);
817766d5 61 TChain * chainTPCITS = tool.MakeChainRandom("timeitstpc.txt.Good","itstpc",0,10000);
3e55050f 62 TChain * chainTPCTOF = tool.MakeChainRandom("timetoftpc.txt.Good","toftpc",0,10000);
817766d5 63 TChain * chainLaser = tool.MakeChainRandom("timelaser.txt.Good","laserInfo",0,10000);
74235403 64 chainTime->Lookup();
dde68d36 65 chainLaser->Lookup();
c74c5f6c 66*/
67
c74c5f6c 68#include "Riostream.h"
2811495d 69#include "TDatabasePDG.h"
70#include "TGraphErrors.h"
c74c5f6c 71#include "TH1F.h"
c74c5f6c 72#include "THnSparse.h"
73#include "TList.h"
74#include "TMath.h"
2811495d 75#include "TTimeStamp.h"
76#include "TTree.h"
c74c5f6c 77#include "TVectorD.h"
2811495d 78//#include "TChain.h"
79//#include "TFile.h"
80
81#include "AliDCSSensor.h"
82#include "AliDCSSensorArray.h"
c74c5f6c 83#include "AliESDEvent.h"
c74c5f6c 84#include "AliESDInputHandler.h"
2811495d 85#include "AliESDVertex.h"
86#include "AliESDfriend.h"
c74c5f6c 87#include "AliLog.h"
817766d5 88#include "AliRelAlignerKalman.h"
2811495d 89#include "AliTPCCalROC.h"
90#include "AliTPCParam.h"
c74c5f6c 91#include "AliTPCTracklet.h"
da6c0bc9 92#include "AliTPCcalibDB.h"
93#include "AliTPCcalibLaser.h"
2811495d 94#include "AliTPCcalibTime.h"
95#include "AliTPCclusterMI.h"
96#include "AliTPCseed.h"
817766d5 97#include "AliTrackPointArray.h"
2811495d 98#include "AliTracker.h"
817766d5 99
c74c5f6c 100ClassImp(AliTPCcalibTime)
101
102
103AliTPCcalibTime::AliTPCcalibTime()
da6c0bc9 104 :AliTPCcalibBase(),
da6c0bc9 105 fLaser(0), // pointer to laser calibration
106 fDz(0), // current delta z
d3ce44cb 107 fCutMaxD(3), // maximal distance in rfi ditection
108 fCutMaxDz(25), // maximal distance in rfi ditection
c74c5f6c 109 fCutTheta(0.03), // maximal distan theta
2be25cc9 110 fCutMinDir(-0.99), // direction vector products
b842d904 111 fCutTracks(100),
dde68d36 112 fArrayDz(0), //NEW! Tmap of V drifts for different triggers
817766d5 113 fAlignITSTPC(0), //alignemnt array ITS TPC match
114 fAlignTRDTPC(0), //alignemnt array TRD TPC match
115 fAlignTOFTPC(0), //alignemnt array TOF TPC match
2be25cc9 116 fTimeBins(0),
117 fTimeStart(0),
118 fTimeEnd(0),
119 fPtBins(0),
120 fPtStart(0),
121 fPtEnd(0),
122 fVdriftBins(0),
123 fVdriftStart(0),
124 fVdriftEnd(0),
125 fRunBins(0),
126 fRunStart(0),
127 fRunEnd(0)
c74c5f6c 128{
6ba68367 129 //
130 // default constructor
131 //
c74c5f6c 132 AliInfo("Default Constructor");
2be25cc9 133 for (Int_t i=0;i<3;i++) {
134 fHistVdriftLaserA[i]=0;
135 fHistVdriftLaserC[i]=0;
136 }
d3ce44cb 137 for (Int_t i=0;i<10;i++) {
74235403 138 fCosmiMatchingHisto[i]=0;
139 }
391ffdb2 140 //
141 for (Int_t i=0;i<5;i++) {
142 fResHistoTPCITS[i]=0;
143 fResHistoTPCTRD[i]=0;
0b736a46 144 fResHistoTPCTOF[i]=0;
391ffdb2 145 fResHistoTPCvertex[i]=0;
146 }
147
c74c5f6c 148}
149
74235403 150AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeVdrift)
c74c5f6c 151 :AliTPCcalibBase(),
d3ce44cb 152 fLaser(0), // pointer to laser calibration
153 fDz(0), // current delta z
154 fCutMaxD(5*0.5356), // maximal distance in rfi ditection
dde68d36 155 fCutMaxDz(40), // maximal distance in rfi ditection
d3ce44cb 156 fCutTheta(5*0.004644),// maximal distan theta
157 fCutMinDir(-0.99), // direction vector products
b842d904 158 fCutTracks(100),
dde68d36 159 fArrayDz(0), //Tmap of V drifts for different triggers
817766d5 160 fAlignITSTPC(0), //alignemnt array ITS TPC match
161 fAlignTRDTPC(0), //alignemnt array TRD TPC match
162 fAlignTOFTPC(0), //alignemnt array TOF TPC match
2be25cc9 163 fTimeBins(0),
164 fTimeStart(0),
165 fTimeEnd(0),
166 fPtBins(0),
167 fPtStart(0),
168 fPtEnd(0),
169 fVdriftBins(0),
170 fVdriftStart(0),
171 fVdriftEnd(0),
172 fRunBins(0),
173 fRunStart(0),
174 fRunEnd(0)
c74c5f6c 175{
6ba68367 176 //
177 // Non deafaul constructor - to be used in the Calibration setups
178 //
179
c74c5f6c 180 SetName(name);
181 SetTitle(title);
2be25cc9 182 for (Int_t i=0;i<3;i++) {
183 fHistVdriftLaserA[i]=0;
184 fHistVdriftLaserC[i]=0;
185 }
c74c5f6c 186
391ffdb2 187 for (Int_t i=0;i<5;i++) {
188 fResHistoTPCITS[i]=0;
189 fResHistoTPCTRD[i]=0;
0b736a46 190 fResHistoTPCTOF[i]=0;
391ffdb2 191 fResHistoTPCvertex[i]=0;
192 }
193
194
c74c5f6c 195 AliInfo("Non Default Constructor");
74235403 196 fTimeBins =(EndTime-StartTime)/deltaIntegrationTimeVdrift;
2be25cc9 197 fTimeStart =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi();
198 fTimeEnd =EndTime; //(((TObjString*)(mapGRP->GetValue("fAliceStopTime")))->GetString()).Atoi();
dde68d36 199 fPtBins = 400;
2be25cc9 200 fPtStart = -0.04;
201 fPtEnd = 0.04;
dde68d36 202 fVdriftBins = 500;
203 fVdriftStart= -0.1;
204 fVdriftEnd = 0.1;
32100b1d 205 fRunBins = 1000001;
dde68d36 206 fRunStart = -1.5;
32100b1d 207 fRunEnd = 999999.5;
2be25cc9 208
209 Int_t binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins };
210 Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart , fRunStart};
211 Double_t xmaxVdriftLaser[4] = {fTimeEnd , fPtEnd , fVdriftEnd , fRunEnd };
dde68d36 212 TString axisTitle[4]={
213 "T",
214 "#delta_{P/T}",
215 "value",
216 "run"
217 };
218 TString histoName[3]={
219 "Loffset",
220 "Lcorr",
221 "Lgy"
222 };
2be25cc9 223
dde68d36 224
2be25cc9 225 for (Int_t i=0;i<3;i++) {
226 fHistVdriftLaserA[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
227 fHistVdriftLaserC[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
dde68d36 228 fHistVdriftLaserA[i]->SetName(histoName[i]);
229 fHistVdriftLaserC[i]->SetName(histoName[i]);
230 for (Int_t iaxis=0; iaxis<4;iaxis++){
231 fHistVdriftLaserA[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
232 fHistVdriftLaserC[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
233 }
2be25cc9 234 }
235 fBinsVdrift[0] = fTimeBins;
236 fBinsVdrift[1] = fPtBins;
237 fBinsVdrift[2] = fVdriftBins;
238 fBinsVdrift[3] = fRunBins;
239 fXminVdrift[0] = fTimeStart;
240 fXminVdrift[1] = fPtStart;
241 fXminVdrift[2] = fVdriftStart;
242 fXminVdrift[3] = fRunStart;
243 fXmaxVdrift[0] = fTimeEnd;
244 fXmaxVdrift[1] = fPtEnd;
245 fXmaxVdrift[2] = fVdriftEnd;
246 fXmaxVdrift[3] = fRunEnd;
247
dde68d36 248 fArrayDz=new TObjArray();
817766d5 249 fAlignITSTPC = new TObjArray; //alignemnt array ITS TPC match
250 fAlignTRDTPC = new TObjArray; //alignemnt array ITS TPC match
251 fAlignTOFTPC = new TObjArray; //alignemnt array ITS TPC match
3e55050f 252 fAlignITSTPC->SetOwner(kTRUE);
253 fAlignTRDTPC->SetOwner(kTRUE);
254 fAlignTOFTPC->SetOwner(kTRUE);
255
c74c5f6c 256
d3ce44cb 257 fCosmiMatchingHisto[0]=new TH1F("Cosmics matching","p0-all" ,100,-10*0.5356 ,10*0.5356 );
258 fCosmiMatchingHisto[1]=new TH1F("Cosmics matching","p1-all" ,100,-10*4.541 ,10*4.541 );
259 fCosmiMatchingHisto[2]=new TH1F("Cosmics matching","p2-all" ,100,-10*0.01134 ,10*0.01134 );
260 fCosmiMatchingHisto[3]=new TH1F("Cosmics matching","p3-all" ,100,-10*0.004644,10*0.004644);
261 fCosmiMatchingHisto[4]=new TH1F("Cosmics matching","p4-all" ,100,-10*0.03773 ,10*0.03773 );
262 fCosmiMatchingHisto[5]=new TH1F("Cosmics matching","p0-isPair",100,-10*0.5356 ,10*0.5356 );
263 fCosmiMatchingHisto[6]=new TH1F("Cosmics matching","p1-isPair",100,-10*4.541 ,10*4.541 );
264 fCosmiMatchingHisto[7]=new TH1F("Cosmics matching","p2-isPair",100,-10*0.01134 ,10*0.01134 );
265 fCosmiMatchingHisto[8]=new TH1F("Cosmics matching","p3-isPair",100,-10*0.004644,10*0.004644);
266 fCosmiMatchingHisto[9]=new TH1F("Cosmics matching","p4-isPair",100,-10*0.03773 ,10*0.03773 );
391ffdb2 267 BookDistortionMaps();
74235403 268}
c74c5f6c 269
270AliTPCcalibTime::~AliTPCcalibTime(){
271 //
6ba68367 272 // Virtual Destructor
c74c5f6c 273 //
d3ce44cb 274 for(Int_t i=0;i<3;i++){
275 if(fHistVdriftLaserA[i]){
276 delete fHistVdriftLaserA[i];
277 fHistVdriftLaserA[i]=NULL;
278 }
279 if(fHistVdriftLaserC[i]){
280 delete fHistVdriftLaserC[i];
281 fHistVdriftLaserC[i]=NULL;
282 }
283 }
dde68d36 284 if(fArrayDz){
285 fArrayDz->SetOwner();
286 fArrayDz->Delete();
287 delete fArrayDz;
288 fArrayDz=NULL;
2be25cc9 289 }
d3ce44cb 290 for(Int_t i=0;i<5;i++){
291 if(fCosmiMatchingHisto[i]){
292 delete fCosmiMatchingHisto[i];
293 fCosmiMatchingHisto[i]=NULL;
294 }
74235403 295 }
391ffdb2 296
297 for (Int_t i=0;i<5;i++) {
298 delete fResHistoTPCITS[i];
299 delete fResHistoTPCTRD[i];
0b736a46 300 delete fResHistoTPCTOF[i];
391ffdb2 301 delete fResHistoTPCvertex[i];
302 fResHistoTPCITS[i]=0;
303 fResHistoTPCTRD[i]=0;
0b736a46 304 fResHistoTPCTOF[i]=0;
391ffdb2 305 fResHistoTPCvertex[i]=0;
306 }
307
308
3e55050f 309 fAlignITSTPC->SetOwner(kTRUE);
310 fAlignTRDTPC->SetOwner(kTRUE);
311 fAlignTOFTPC->SetOwner(kTRUE);
312
817766d5 313 fAlignITSTPC->Delete();
314 fAlignTRDTPC->Delete();
315 fAlignTOFTPC->Delete();
3e55050f 316 delete fAlignITSTPC;
317 delete fAlignTRDTPC;
318 delete fAlignTOFTPC;
c74c5f6c 319}
2be25cc9 320
2811495d 321Bool_t AliTPCcalibTime::IsLaser(const AliESDEvent *const /*event*/){
6ba68367 322 //
323 // Indicator is laser event not yet implemented - to be done using trigger info or event specie
324 //
dde68d36 325 return kTRUE; //More accurate creteria to be added
da6c0bc9 326}
2811495d 327Bool_t AliTPCcalibTime::IsCosmics(const AliESDEvent *const /*event*/){
6ba68367 328 //
329 // Indicator is cosmic event not yet implemented - to be done using trigger info or event specie
330 //
331
dde68d36 332 return kTRUE; //More accurate creteria to be added
333}
2811495d 334Bool_t AliTPCcalibTime::IsBeam(const AliESDEvent *const /*event*/){
6ba68367 335 //
336 // Indicator is physic event not yet implemented - to be done using trigger info or event specie
337 //
338
dde68d36 339 return kTRUE; //More accurate creteria to be added
340}
341void AliTPCcalibTime::ResetCurrent(){
342 fDz=0; //Reset current dz
2be25cc9 343}
6ba68367 344
345
346
2be25cc9 347void AliTPCcalibTime::Process(AliESDEvent *event){
6ba68367 348 //
349 // main function to make calibration
350 //
2be25cc9 351 if(!event) return;
352 if (event->GetNumberOfTracks()<2) return;
353 ResetCurrent();
dde68d36 354 if(IsLaser (event)) ProcessLaser (event);
355 if(IsCosmics(event)) ProcessCosmic(event);
356 if(IsBeam (event)) ProcessBeam (event);
2be25cc9 357}
358
359void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){
c74c5f6c 360 //
2be25cc9 361 // Fit drift velocity using laser
362 //
363 // 0. cuts
dde68d36 364 const Int_t kMinTracks = 40; // minimal number of laser tracks
365 const Int_t kMinTracksSide = 20; // minimal number of tracks per side
366 const Float_t kMaxDeltaZ = 30.; // maximal trigger delay
367 const Float_t kMaxDeltaV = 0.05; // maximal deltaV
2be25cc9 368 const Float_t kMaxRMS = 0.1; // maximal RMS of tracks
dde68d36 369 //
2be25cc9 370 /*
371 TCut cutRMS("sqrt(laserA.fElements[4])<0.1&&sqrt(laserC.fElements[4])<0.1");
372 TCut cutZ("abs(laserA.fElements[0]-laserC.fElements[0])<3");
373 TCut cutV("abs(laserA.fElements[1]-laserC.fElements[1])<0.01");
374 TCut cutY("abs(laserA.fElements[2]-laserC.fElements[2])<2");
375 TCut cutAll = cutRMS+cutZ+cutV+cutY;
376 */
377 if (event->GetNumberOfTracks()<kMinTracks) return;
c74c5f6c 378 //
2be25cc9 379 if(!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
380 fLaser->Process(event);
381 if (fLaser->GetNtracks()<kMinTracks) return; // small amount of tracks cut
dde68d36 382 if (fLaser->fFitAside->GetNrows()==0 && fLaser->fFitCside->GetNrows()==0) return; // no fit neither a or C side
c74c5f6c 383 //
2be25cc9 384 // debug streamer - activate stream level
385 // Use it for tuning of the cuts
da6c0bc9 386 //
dde68d36 387 // cuts to be applied
388 //
389 Int_t isReject[2]={0,0};
390 //
391 // not enough tracks
392 if (TMath::Abs((*fLaser->fFitAside)[3]) < kMinTracksSide) isReject[0]|=1;
393 if (TMath::Abs((*fLaser->fFitCside)[3]) < kMinTracksSide) isReject[1]|=1;
394 // unreasonable z offset
395 if (TMath::Abs((*fLaser->fFitAside)[0])>kMaxDeltaZ) isReject[0]|=2;
396 if (TMath::Abs((*fLaser->fFitCside)[0])>kMaxDeltaZ) isReject[1]|=2;
397 // unreasonable drift velocity
398 if (TMath::Abs((*fLaser->fFitAside)[1]-1)>kMaxDeltaV) isReject[0]|=4;
399 if (TMath::Abs((*fLaser->fFitCside)[1]-1)>kMaxDeltaV) isReject[1]|=4;
400 // big chi2
401 if (TMath::Sqrt(TMath::Abs((*fLaser->fFitAside)[4]))>kMaxRMS ) isReject[0]|=8;
402 if (TMath::Sqrt(TMath::Abs((*fLaser->fFitCside)[4]))>kMaxRMS ) isReject[1]|=8;
403
404
405
406
2be25cc9 407 if (fStreamLevel>0){
408 printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
409
da6c0bc9 410 TTreeSRedirector *cstream = GetDebugStreamer();
411 if (cstream){
412 TTimeStamp tstamp(fTime);
2be25cc9 413 (*cstream)<<"laserInfo"<<
da6c0bc9 414 "run="<<fRun<< // run number
415 "event="<<fEvent<< // event number
416 "time="<<fTime<< // time stamp of event
417 "trigger="<<fTrigger<< // trigger
418 "mag="<<fMagF<< // magnetic field
da6c0bc9 419 //laser
dde68d36 420 "rejectA="<<isReject[0]<<
421 "rejectC="<<isReject[1]<<
2be25cc9 422 "laserA.="<<fLaser->fFitAside<<
423 "laserC.="<<fLaser->fFitCside<<
424 "laserAC.="<<fLaser->fFitACside<<
425 "trigger="<<event->GetFiredTriggerClasses()<<
da6c0bc9 426 "\n";
427 }
428 }
2be25cc9 429 //
2be25cc9 430 // fill histos
431 //
432 TVectorD vdriftA(5), vdriftC(5),vdriftAC(5);
433 vdriftA=*(fLaser->fFitAside);
434 vdriftC=*(fLaser->fFitCside);
435 vdriftAC=*(fLaser->fFitACside);
436 Int_t npointsA=0, npointsC=0;
437 Float_t chi2A=0, chi2C=0;
438 npointsA= TMath::Nint(vdriftA[3]);
439 chi2A= vdriftA[4];
440 npointsC= TMath::Nint(vdriftC[3]);
441 chi2C= vdriftC[4];
2be25cc9 442
2be25cc9 443 TTimeStamp tstamp(fTime);
444 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
445 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
dde68d36 446 Double_t driftA=0, driftC=0;
447 if (vdriftA[1]>1.-kMaxDeltaV) driftA = 1./vdriftA[1]-1.;
448 if (vdriftC[1]>1.-kMaxDeltaV) driftC = 1./vdriftC[1]-1.;
449 //
450 Double_t vecDriftLaserA[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftA,event->GetRunNumber()};
451 Double_t vecDriftLaserC[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftC,event->GetRunNumber()};
452 // Double_t vecDrift[4] ={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitACside))[1])-1,event->GetRunNumber()};
453
454 for (Int_t icalib=0;icalib<3;icalib++){
455 if (icalib==0){ //z0 shift
456 vecDriftLaserA[2]=vdriftA[0]/250.;
457 vecDriftLaserC[2]=vdriftC[0]/250.;
2be25cc9 458 }
dde68d36 459 if (icalib==1){ //vdrel shift
460 vecDriftLaserA[2]=driftA;
461 vecDriftLaserC[2]=driftC;
2be25cc9 462 }
dde68d36 463 if (icalib==2){ //gy shift - full gy - full drift
464 vecDriftLaserA[2]=vdriftA[2]/250.;
465 vecDriftLaserC[2]=vdriftC[2]/250.;
2be25cc9 466 }
391ffdb2 467 //if (isReject[0]==0) fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
468 //if (isReject[1]==0) fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
469 fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
470 fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
2be25cc9 471 }
dde68d36 472
473// THnSparse* curHist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
474// TString shortName=curHist->ClassName();
475// shortName+="_MEAN_DRIFT_LASER_";
476// delete curHist;
477// curHist=NULL;
478// TString name="";
479
480// name=shortName;
481// name+=event->GetFiredTriggerClasses();
482// name.ToUpper();
483// curHist=(THnSparseF*)fArrayDz->FindObject(name);
484// if(!curHist){
485// curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
486// fArrayDz->AddLast(curHist);
487// }
488// curHist->Fill(vecDrift);
489
490// name=shortName;
491// name+="ALL";
492// name.ToUpper();
493// curHist=(THnSparseF*)fArrayDz->FindObject(name);
494// if(!curHist){
495// curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
496// fArrayDz->AddLast(curHist);
497// }
498// curHist->Fill(vecDrift);
2be25cc9 499}
c74c5f6c 500
2811495d 501void AliTPCcalibTime::ProcessCosmic(const AliESDEvent *const event){
6ba68367 502 //
503 // process Cosmic event - track matching A side C side
504 //
c74c5f6c 505 if (!event) {
506 Printf("ERROR: ESD not available");
507 return;
508 }
509 if (event->GetTimeStamp() == 0 ) {
510 Printf("no time stamp!");
511 return;
512 }
74235403 513
2be25cc9 514 //fd
c74c5f6c 515 // Find cosmic pairs
516 //
517 // Track0 is choosen in upper TPC part
518 // Track1 is choosen in lower TPC part
519 //
cc65e4f5 520 const Int_t kMinClustersCross =30;
521 const Int_t kMinClusters =80;
c74c5f6c 522 Int_t ntracks=event->GetNumberOfTracks();
523 if (ntracks==0) return;
74235403 524 if (ntracks > fCutTracks) return;
525
a390f11f 526 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
3e0d9569 527 AliESDfriend *esdFriend=(AliESDfriend*)(((AliESDEvent*)event)->FindListObject("AliESDfriend"));
c74c5f6c 528
529 TObjArray tpcSeeds(ntracks);
530 Double_t vtxx[3]={0,0,0};
531 Double_t svtxx[3]={0.000001,0.000001,100.};
532 AliESDVertex vtx(vtxx,svtxx);
533 //
534 // track loop
535 //
cc65e4f5 536 TArrayI clusterSideA(ntracks);
537 TArrayI clusterSideC(ntracks);
c74c5f6c 538 for (Int_t i=0;i<ntracks;++i) {
cc65e4f5 539 clusterSideA[i]=0;
540 clusterSideC[i]=0;
74235403 541 AliESDtrack *track = event->GetTrack(i);
542
543 const AliExternalTrackParam * trackIn = track->GetInnerParam();
544 const AliExternalTrackParam * trackOut = track->GetOuterParam();
545 if (!trackIn) continue;
546 if (!trackOut) continue;
547
6ba68367 548 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
b1c27e4f 549 if (!friendTrack) continue;
3e55050f 550 if (friendTrack) ProcessSame(track,friendTrack,event);
6ba68367 551 if (friendTrack) ProcessAlignITS(track,friendTrack,event,esdFriend);
0dac7b3a 552 if (friendTrack) ProcessAlignTRD(track,friendTrack);
817766d5 553 if (friendTrack) ProcessAlignTOF(track,friendTrack);
74235403 554 TObject *calibObject;
555 AliTPCseed *seed = 0;
556 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
cc65e4f5 557 if (seed) {
558 tpcSeeds.AddAt(seed,i);
559 Int_t nA=0, nC=0;
560 for (Int_t irow=159;irow>0;irow--) {
561 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
562 if (!cl) continue;
563 if ((cl->GetDetector()%36)<18) nA++;
564 if ((cl->GetDetector()%36)>=18) nC++;
565 }
566 clusterSideA[i]=nA;
567 clusterSideC[i]=nC;
568 }
c74c5f6c 569 }
c74c5f6c 570 if (ntracks<2) return;
571 //
572 // Find pairs
573 //
d3ce44cb 574
c74c5f6c 575 for (Int_t i=0;i<ntracks;++i) {
74235403 576 AliESDtrack *track0 = event->GetTrack(i);
c74c5f6c 577 // track0 - choosen upper part
578 if (!track0) continue;
579 if (!track0->GetOuterParam()) continue;
580 if (track0->GetOuterParam()->GetAlpha()<0) continue;
581 Double_t d1[3];
582 track0->GetDirection(d1);
583 for (Int_t j=0;j<ntracks;++j) {
584 if (i==j) continue;
74235403 585 AliESDtrack *track1 = event->GetTrack(j);
586 //track 1 lower part
587 if (!track1) continue;
588 if (!track1->GetOuterParam()) continue;
cc65e4f5 589 if (track0->GetTPCNcls()+ track1->GetTPCNcls()< kMinClusters) continue;
590 Int_t nAC = TMath::Max( TMath::Min(clusterSideA[i], clusterSideC[j]),
591 TMath::Min(clusterSideC[i], clusterSideA[j]));
592 if (nAC<kMinClustersCross) continue;
593 Int_t nA0=clusterSideA[i];
594 Int_t nC0=clusterSideC[i];
595 Int_t nA1=clusterSideA[j];
596 Int_t nC1=clusterSideC[j];
3e55050f 597 // if (track1->GetOuterParam()->GetAlpha()>0) continue;
74235403 598 //
599 Double_t d2[3];
600 track1->GetDirection(d2);
601
602 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
603 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
604 if (! seed0) continue;
605 if (! seed1) continue;
606 Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
607 Float_t dist0 = track0->GetLinearD(0,0);
608 Float_t dist1 = track1->GetLinearD(0,0);
609 //
610 // conservative cuts - convergence to be guarantied
611 // applying before track propagation
3e55050f 612 if (TMath::Abs(TMath::Abs(dist0)-TMath::Abs(dist1))>fCutMaxD) continue; // distance to the 0,0
613 if (TMath::Abs(dir)<TMath::Abs(fCutMinDir)) continue; // direction vector product
74235403 614 Float_t bz = AliTracker::GetBz();
615 Float_t dvertex0[2]; //distance to 0,0
616 Float_t dvertex1[2]; //distance to 0,0
617 track0->GetDZ(0,0,0,bz,dvertex0);
618 track1->GetDZ(0,0,0,bz,dvertex1);
619 if (TMath::Abs(dvertex0[1])>250) continue;
620 if (TMath::Abs(dvertex1[1])>250) continue;
621 //
622 //
623 //
624 Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
625 AliExternalTrackParam param0(*track0);
626 AliExternalTrackParam param1(*track1);
627 //
628 // Propagate using Magnetic field and correct fo material budget
629 //
9963b5e2 630 AliTracker::PropagateTrackTo(&param0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE);
631 AliTracker::PropagateTrackTo(&param1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE);
74235403 632 //
633 // Propagate rest to the 0,0 DCA - z should be ignored
634 //
635 //Bool_t b0 = ;
636 param0.PropagateToDCA(&vtx,bz,1000);
637 //Bool_t b1 =
638 param1.PropagateToDCA(&vtx,bz,1000);
74235403 639 param0.GetDZ(0,0,0,bz,dvertex0);
640 param1.GetDZ(0,0,0,bz,dvertex1);
dde68d36 641 Double_t xyz0[3];
642 Double_t xyz1[3];
74235403 643 param0.GetXYZ(xyz0);
644 param1.GetXYZ(xyz1);
645 Bool_t isPair = IsPair(&param0,&param1);
d3ce44cb 646 Bool_t isCross = IsCross(track0, track1);
dde68d36 647 Bool_t isSame = IsSame(track0, track1);
d3ce44cb 648
dde68d36 649 THnSparse* hist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
650 TString shortName=hist->ClassName();
651 shortName+="_MEAN_VDRIFT_COSMICS_";
652 delete hist;
653 hist=NULL;
654
3e55050f 655 if((isSame) || (isCross && isPair)){
656 if (track0->GetTPCNcls()+ track1->GetTPCNcls()> 80) {
74235403 657 fDz = param0.GetZ() - param1.GetZ();
cc65e4f5 658 Double_t sign=(nA0>nA1)? 1:-1;
659 fDz*=sign;
74235403 660 TTimeStamp tstamp(fTime);
661 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
662 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
dde68d36 663 Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
664 THnSparse* curHist=NULL;
665 TString name="";
666
667 name=shortName;
668 name+=event->GetFiredTriggerClasses();
669 name.ToUpper();
670 curHist=(THnSparseF*)fArrayDz->FindObject(name);
74235403 671 if(!curHist){
dde68d36 672 curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
673 fArrayDz->AddLast(curHist);
74235403 674 }
dde68d36 675// curHist=(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
676// if(!curHist){
677// curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
678// fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
679// }
680 curHist->Fill(vecDrift);
74235403 681
dde68d36 682 name=shortName;
683 name+="ALL";
684 name.ToUpper();
685 curHist=(THnSparseF*)fArrayDz->FindObject(name);
74235403 686 if(!curHist){
dde68d36 687 curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
688 fArrayDz->AddLast(curHist);
74235403 689 }
dde68d36 690// curHist=(THnSparseF*)(fMapDz->GetValue("all"));
691// if(!curHist){
692// curHist=new THnSparseF("all","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
693// fMapDz->Add(new TObjString("all"),curHist);
694// }
695 curHist->Fill(vecDrift);
74235403 696 }
697 }
d3ce44cb 698 TTreeSRedirector *cstream = GetDebugStreamer();
699 if (fStreamLevel>0){
700 if (cstream){
701 (*cstream)<<"trackInfo"<<
702 "tr0.="<<track0<<
703 "tr1.="<<track1<<
704 "p0.="<<&param0<<
705 "p1.="<<&param1<<
cc65e4f5 706 "nAC="<<nAC<<
707 "nA0="<<nA0<<
708 "nA1="<<nA1<<
709 "nC0="<<nC0<<
710 "nC1="<<nC1<<
d3ce44cb 711 "isPair="<<isPair<<
712 "isCross="<<isCross<<
dde68d36 713 "isSame="<<isSame<<
d3ce44cb 714 "fDz="<<fDz<<
715 "fRun="<<fRun<<
716 "fTime="<<fTime<<
717 "\n";
718 }
719 }
c74c5f6c 720 } // end 2nd order loop
721 } // end 1st order loop
74235403 722
2be25cc9 723 if (fStreamLevel>0){
724 TTreeSRedirector *cstream = GetDebugStreamer();
725 if (cstream){
2be25cc9 726 (*cstream)<<"timeInfo"<<
727 "run="<<fRun<< // run number
728 "event="<<fEvent<< // event number
729 "time="<<fTime<< // time stamp of event
730 "trigger="<<fTrigger<< // trigger
731 "mag="<<fMagF<< // magnetic field
732 // Environment values
2be25cc9 733 //
734 // accumulated values
735 //
736 "fDz="<<fDz<< //! current delta z
2be25cc9 737 "trigger="<<event->GetFiredTriggerClasses()<<
738 "\n";
739 }
740 }
a390f11f 741 if (GetDebugLevel()>20) printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
c74c5f6c 742}
743
2811495d 744void AliTPCcalibTime::ProcessBeam(const AliESDEvent *const /*event*/){
6ba68367 745 //
746 // Not special treatment yet - the same for cosmic and physic event
747 //
dde68d36 748}
749
6ba68367 750void AliTPCcalibTime::Analyze(){
751 //
752 // Special macro to analyze result of calibration and extract calibration entries
753 // Not yet ported to the Analyze function yet
754 //
755}
74235403 756
2811495d 757THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name) const
758{
6ba68367 759 //
760 // Get histogram for given trigger mask
761 //
dde68d36 762 TIterator* iterator = fArrayDz->MakeIterator();
763 iterator->Reset();
764 TString newName=name;
765 newName.ToUpper();
766 THnSparse* newHist=new THnSparseF(newName,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
767 THnSparse* addHist=NULL;
768 while((addHist=(THnSparseF*)iterator->Next())){
769 if(!addHist) continue;
770 TString histName=addHist->GetName();
771 if(!histName.Contains(newName)) continue;
772 addHist->Print();
773 newHist->Add(addHist);
74235403 774 }
dde68d36 775 return newHist;
74235403 776}
777
2811495d 778TObjArray* AliTPCcalibTime::GetHistoDrift() const
779{
6ba68367 780 //
781 // return array of histograms
782 //
dde68d36 783 return fArrayDz;
74235403 784}
785
786TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){
6ba68367 787 //
788 // Make a drift velocity (delta Z) graph
789 //
dde68d36 790 THnSparse* histoDrift=GetHistoDrift(name);
791 TGraphErrors* graphDrift=NULL;
792 if(histoDrift){
793 graphDrift=FitSlices(histoDrift,2,0,400,100,0.05,0.95, kTRUE);
794 TString end=histoDrift->GetName();
795 Int_t pos=end.Index("_");
796 end=end(pos,end.Capacity()-pos);
797 TString graphName=graphDrift->ClassName();
798 graphName+=end;
799 graphName.ToUpper();
800 graphDrift->SetName(graphName);
74235403 801 }
802 return graphDrift;
803}
804
dde68d36 805TObjArray* AliTPCcalibTime::GetGraphDrift(){
6ba68367 806 //
807 // make a array of drift graphs
808 //
dde68d36 809 TObjArray* arrayGraphDrift=new TObjArray();
810 TIterator* iterator=fArrayDz->MakeIterator();
74235403 811 iterator->Reset();
dde68d36 812 THnSparse* addHist=NULL;
813 while((addHist=(THnSparseF*)iterator->Next())) arrayGraphDrift->AddLast(GetGraphDrift(addHist->GetName()));
814 return arrayGraphDrift;
74235403 815}
816
dde68d36 817AliSplineFit* AliTPCcalibTime::GetFitDrift(const char* name){
6ba68367 818 //
819 // Make a fit AliSplinefit of drift velocity
820 //
d3ce44cb 821 TGraph* graphDrift=GetGraphDrift(name);
dde68d36 822 AliSplineFit* fitDrift=NULL;
74235403 823 if(graphDrift && graphDrift->GetN()){
dde68d36 824 fitDrift=new AliSplineFit();
825 fitDrift->SetGraph(graphDrift);
826 fitDrift->SetMinPoints(graphDrift->GetN()+1);
827 fitDrift->InitKnots(graphDrift,2,0,0.001);
828 fitDrift->SplineFit(0);
829 TString end=graphDrift->GetName();
830 Int_t pos=end.Index("_");
831 end=end(pos,end.Capacity()-pos);
832 TString fitName=fitDrift->ClassName();
833 fitName+=end;
834 fitName.ToUpper();
835 //fitDrift->SetName(fitName);
74235403 836 delete graphDrift;
dde68d36 837 graphDrift=NULL;
74235403 838 }
839 return fitDrift;
840}
841
74235403 842
2811495d 843Long64_t AliTPCcalibTime::Merge(TCollection *const li) {
6ba68367 844 //
845 // Object specific merging procedure
846 //
c74c5f6c 847 TIterator* iter = li->MakeIterator();
848 AliTPCcalibTime* cal = 0;
849
850 while ((cal = (AliTPCcalibTime*)iter->Next())) {
851 if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
852 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
853 return -1;
854 }
2be25cc9 855 for (Int_t imeas=0; imeas<3; imeas++){
856 if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
857 fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
858 fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
859 }
860 }
9112627b 861 //
862 for (Int_t imeas=0; imeas<5; imeas++){
863 if (cal->GetResHistoTPCITS(imeas) && cal->GetResHistoTPCITS(imeas)){
864 fResHistoTPCITS[imeas]->Add(cal->fResHistoTPCITS[imeas]);
865 fResHistoTPCvertex[imeas]->Add(cal->fResHistoTPCvertex[imeas]);
33fc6b9e 866 }
867 if (cal->fResHistoTPCTRD[imeas]){
868 if (fResHistoTPCTRD[imeas])
869 fResHistoTPCTRD[imeas]->Add(cal->fResHistoTPCTRD[imeas]);
870 else
871 fResHistoTPCTRD[imeas]=(THnSparse*)cal->fResHistoTPCTRD[imeas]->Clone();
872 }
873 if (cal->fResHistoTPCTOF[imeas]){
874 if (fResHistoTPCTOF[imeas])
875 fResHistoTPCTOF[imeas]->Add(cal->fResHistoTPCTOF[imeas]);
876 else
877 fResHistoTPCTOF[imeas]=(THnSparse*)cal->fResHistoTPCTOF[imeas]->Clone();
9112627b 878 }
879 }
dde68d36 880 TObjArray* addArray=cal->GetHistoDrift();
881 if(!addArray) return 0;
882 TIterator* iterator = addArray->MakeIterator();
2be25cc9 883 iterator->Reset();
dde68d36 884 THnSparse* addHist=NULL;
885 while((addHist=(THnSparseF*)iterator->Next())){
886 if(!addHist) continue;
2be25cc9 887 addHist->Print();
dde68d36 888 THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName());
2be25cc9 889 if(!localHist){
890 localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
dde68d36 891 fArrayDz->AddLast(localHist);
2be25cc9 892 }
893 localHist->Add(addHist);
894 }
391ffdb2 895
d3ce44cb 896 for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
3e55050f 897 //
898 // Merge alignment
899 //
900 for (Int_t itype=0; itype<3; itype++){
901 //
902 //
903 TObjArray *arr0= 0;
904 TObjArray *arr1= 0;
905 if (itype==0) {arr0=fAlignITSTPC; arr1=cal->fAlignITSTPC;}
906 if (itype==1) {arr0=fAlignTRDTPC; arr1=cal->fAlignTRDTPC;}
907 if (itype==2) {arr0=fAlignTOFTPC; arr1=cal->fAlignTOFTPC;}
908 if (!arr1) continue;
909 if (!arr0) arr0=new TObjArray(arr1->GetEntriesFast());
910 if (arr1->GetEntriesFast()>arr0->GetEntriesFast()){
911 arr0->Expand(arr1->GetEntriesFast());
912 }
913 for (Int_t i=0;i<arr1->GetEntriesFast(); i++){
914 AliRelAlignerKalman *kalman1 = (AliRelAlignerKalman *)arr1->UncheckedAt(i);
915 AliRelAlignerKalman *kalman0 = (AliRelAlignerKalman *)arr0->UncheckedAt(i);
916 if (!kalman1) continue;
917 if (!kalman0) {arr0->AddAt(new AliRelAlignerKalman(*kalman1),i); continue;}
918 kalman0->SetRejectOutliers(kFALSE);
919 kalman0->Merge(kalman1);
920 }
921 }
922
c74c5f6c 923 }
c74c5f6c 924 return 0;
c74c5f6c 925}
926
c74c5f6c 927Bool_t AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
c74c5f6c 928 /*
929 // 0. Same direction - OPOSITE - cutDir +cutT
930 TCut cutDir("cutDir","dir<-0.99")
931 // 1.
932 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
933 //
934 // 2. The same rphi
935 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
936 //
c74c5f6c 937 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
938 // 1/Pt diff cut
939 */
940 const Double_t *p0 = tr0->GetParameter();
941 const Double_t *p1 = tr1->GetParameter();
74235403 942 fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
943 fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
d3ce44cb 944 fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
74235403 945 fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
946 fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
947
c74c5f6c 948 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
949 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
da6c0bc9 950 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
c74c5f6c 951 Double_t d0[3], d1[3];
952 tr0->GetDirection(d0);
953 tr1->GetDirection(d1);
954 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
74235403 955
d3ce44cb 956 fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
957 fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
958 fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
959 fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
960 fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
961
c74c5f6c 962 return kTRUE;
963}
2811495d 964Bool_t AliTPCcalibTime::IsCross(AliESDtrack *const tr0, AliESDtrack *const tr1){
cc65e4f5 965 //
966 // check if the cosmic pair of tracks crossed A/C side
967 //
968 Bool_t result= tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0;
969 if (result==kFALSE) return result;
970 result=kTRUE;
971 return result;
d3ce44cb 972}
973
2811495d 974Bool_t AliTPCcalibTime::IsSame(AliESDtrack *const tr0, AliESDtrack *const tr1){
3e55050f 975 //
976 // track crossing the CE
977 // 0. minimal number of clusters
978 // 1. Same sector +-1
979 // 2. Inner and outer track param on opposite side
980 // 3. Outer and inner track parameter close each to other
981 // 3.
982 Bool_t result=kTRUE;
983 //
984 // inner and outer on opposite sides in z
985 //
986 const Int_t knclCut0 = 30;
987 const Double_t kalphaCut = 0.4;
988 //
989 // 0. minimal number of clusters
990 //
991 if (tr0->GetTPCNcls()<knclCut0) return kFALSE;
992 if (tr1->GetTPCNcls()<knclCut0) return kFALSE;
993 //
994 // 1. alpha cut - sector+-1
995 //
996 if (TMath::Abs(tr0->GetOuterParam()->GetAlpha()-tr1->GetOuterParam()->GetAlpha())>kalphaCut) return kFALSE;
997 //
998 // 2. Z crossing
999 //
1000 if (tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0) result&=kFALSE;
1001 if (tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0) result&=kFALSE;
1002 if (result==kFALSE){
1003 return result;
1004 }
1005 //
1006 //
1007 const Double_t *p0I = tr0->GetInnerParam()->GetParameter();
1008 const Double_t *p1I = tr1->GetInnerParam()->GetParameter();
1009 const Double_t *p0O = tr0->GetOuterParam()->GetParameter();
1010 const Double_t *p1O = tr1->GetOuterParam()->GetParameter();
1011 //
1012 if (TMath::Abs(p0I[0]-p1I[0])>fCutMaxD) result&=kFALSE;
1013 if (TMath::Abs(p0I[1]-p1I[1])>fCutMaxDz) result&=kFALSE;
1014 if (TMath::Abs(p0I[2]-p1I[2])>fCutTheta) result&=kFALSE;
1015 if (TMath::Abs(p0I[3]-p1I[3])>fCutTheta) result&=kFALSE;
1016 if (TMath::Abs(p0O[0]-p1O[0])>fCutMaxD) result&=kFALSE;
1017 if (TMath::Abs(p0O[1]-p1O[1])>fCutMaxDz) result&=kFALSE;
1018 if (TMath::Abs(p0O[2]-p1O[2])>fCutTheta) result&=kFALSE;
1019 if (TMath::Abs(p0O[3]-p1O[3])>fCutTheta) result&=kFALSE;
1020 if (result==kTRUE){
1021 result=kTRUE; // just to put break point here
1022 }
1023 return result;
dde68d36 1024}
1025
31aa7c5c 1026
2811495d 1027void AliTPCcalibTime::ProcessSame(AliESDtrack *const track, AliESDfriendTrack *const friendTrack, const AliESDEvent *const event){
3e55050f 1028 //
1029 // Process TPC tracks crossing CE
1030 //
1031 // 0. Select only track crossing the CE
1032 // 1. Cut on the track length
1033 // 2. Refit the terack on A and C side separatelly
1034 // 3. Fill time histograms
1035 const Int_t kMinNcl=100;
1036 const Int_t kMinNclS=25; // minimul number of clusters on the sides
1037 if (!friendTrack->GetTPCOut()) return;
1038 //
1039 // 0. Select only track crossing the CE
1040 //
1041 if (track->GetInnerParam()->GetZ()*friendTrack->GetTPCOut()->GetZ()>0) return;
1042 //
1043 // 1. cut on track length
1044 //
1045 if (track->GetTPCNcls()<kMinNcl) return;
1046 //
1047 // 2. Refit track sepparatel on A and C side
1048 //
1049 TObject *calibObject;
1050 AliTPCseed *seed = 0;
1051 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
1052 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
1053 }
1054 if (!seed) return;
1055 //
1056 AliExternalTrackParam trackIn(*track->GetInnerParam());
1057 AliExternalTrackParam trackOut(*track->GetOuterParam());
1058 Double_t cov[3]={0.01,0.,0.01}; //use the same errors
1059 Double_t xyz[3]={0,0.,0.0};
1060 Double_t bz =0;
1061 Int_t nclIn=0,nclOut=0;
1062 trackIn.ResetCovariance(30.);
1063 trackOut.ResetCovariance(30.);
1064 //
1065 //2.a Refit inner
1066 //
1067 for (Int_t irow=0;irow<159;irow++) {
1068 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
1069 if (!cl) continue;
1070 if (cl->GetX()<80) continue;
1071 if (track->GetInnerParam()->GetZ()<0 &&(cl->GetDetector()%36)<18) break;
1072 if (track->GetInnerParam()->GetZ()>0 &&(cl->GetDetector()%36)>=18) break;
1073 Int_t sector = cl->GetDetector();
1074 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
1075 if (TMath::Abs(dalpha)>0.01){
1076 if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
1077 }
1078 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
1079 trackIn.GetXYZ(xyz);
1080 bz = AliTracker::GetBz(xyz);
1081 if (!trackIn.PropagateTo(r[0],bz)) break;
1082 nclIn++;
1083 trackIn.Update(&r[1],cov);
1084 }
1085 //
1086 //2.b Refit outer
1087 //
1088 for (Int_t irow=159;irow>0;irow--) {
1089 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
1090 if (!cl) continue;
1091 if (cl->GetX()<80) continue;
1092 if (cl->GetZ()*track->GetOuterParam()->GetZ()<0) break;
1093 if (friendTrack->GetTPCOut()->GetZ()<0 &&(cl->GetDetector()%36)<18) break;
1094 if (friendTrack->GetTPCOut()->GetZ()>0 &&(cl->GetDetector()%36)>=18) break;
1095 Int_t sector = cl->GetDetector();
1096 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
1097 if (TMath::Abs(dalpha)>0.01){
1098 if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
1099 }
1100 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
1101 trackOut.GetXYZ(xyz);
1102 bz = AliTracker::GetBz(xyz);
1103 if (!trackOut.PropagateTo(r[0],bz)) break;
1104 nclOut++;
1105 trackOut.Update(&r[1],cov);
1106 }
1107 trackOut.Rotate(trackIn.GetAlpha());
1108 Double_t meanX = (trackIn.GetX()+trackOut.GetX())*0.5;
1109 trackIn.PropagateTo(meanX,bz);
1110 trackOut.PropagateTo(meanX,bz);
1111 TTreeSRedirector *cstream = GetDebugStreamer();
1112 if (cstream){
1113 TVectorD gxyz(3);
1114 trackIn.GetXYZ(gxyz.GetMatrixArray());
1115 TTimeStamp tstamp(fTime);
3e55050f 1116 (*cstream)<<"tpctpc"<<
1117 "run="<<fRun<< // run number
1118 "event="<<fEvent<< // event number
1119 "time="<<fTime<< // time stamp of event
1120 "trigger="<<fTrigger<< // trigger
1121 "mag="<<fMagF<< // magnetic field
3e55050f 1122 //
1123 "xyz.="<<&gxyz<< // global position
1124 "tIn.="<<&trackIn<< // refitterd track in
1125 "tOut.="<<&trackOut<< // refitter track out
1126 "nclIn="<<nclIn<< //
1127 "nclOut="<<nclOut<< //
1128 "\n";
1129 }
1130 //
1131 // 3. Fill time histograms
1132 // Debug stremaer expression
1133 // chainTPCTPC->Draw("(tIn.fP[1]-tOut.fP[1])*sign(-tIn.fP[3]):tIn.fP[3]","min(nclIn,nclOut)>30","")
1134 if (TMath::Min(nclIn,nclOut)>kMinNclS){
1135 fDz = trackOut.GetZ()-trackIn.GetZ();
1136 if (trackOut.GetTgl()<0) fDz*=-1.;
1137 TTimeStamp tstamp(fTime);
1138 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1139 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1140 Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
1141 //
1142 // fill histograms per trigger class and itegrated
1143 //
1144 THnSparse* curHist=NULL;
1145 for (Int_t itype=0; itype<2; itype++){
1146 TString name="MEAN_VDRIFT_CROSS_";
1147 if (itype==0){
1148 name+=event->GetFiredTriggerClasses();
1149 name.ToUpper();
1150 }else{
1151 name+="ALL";
1152 }
1153 curHist=(THnSparseF*)fArrayDz->FindObject(name);
1154 if(!curHist){
1155 curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
1156 fArrayDz->AddLast(curHist);
1157 }
1158 curHist->Fill(vecDrift);
1159 }
1160 }
1161
1162}
817766d5 1163
2811495d 1164void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, AliESDfriendTrack *const friendTrack, const AliESDEvent *const event, AliESDfriend *const esdFriend){
817766d5 1165 //
3e55050f 1166 // Process track - Update TPC-ITS alignment
1167 // Updates:
1168 // 0. Apply standartd cuts
1169 // 1. Recalucluate the current statistic median/RMS
1170 // 2. Apply median+-rms cut
1171 // 3. Update kalman filter
1172 //
1173 const Int_t kMinTPC = 80; // minimal number of TPC cluster
1174 const Int_t kMinITS = 3; // minimal number of ITS cluster
1175 const Double_t kMinZ = 10; // maximal dz distance
b96c3aef 1176 const Double_t kMaxDy = 2.; // maximal dy distance
1177 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1178 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1179 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1180 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1181 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
b96c3aef 1182 const Double_t kMinPt = 0.3; // minimal pt
be67055b 1183 const Int_t kN=50; // deepnes of history
3e55050f 1184 static Int_t kglast=0;
1185 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
1186 /*
1187 0. Standrd cuts:
1188 TCut cut="abs(pTPC.fP[2]-pITS.fP[2])<0.01&&abs(pTPC.fP[3]-pITS.fP[3])<0.01&&abs(pTPC.fP[2]-pITS.fP[2])<1";
1189 */
817766d5 1190 //
3e55050f 1191 // 0. Apply standard cuts
817766d5 1192 //
1193 Int_t dummycl[1000];
817766d5 1194 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
b842d904 1195 if (track->GetITSclusters(dummycl)<kMinITS) return; // minimal amount of clusters
b96c3aef 1196 if (!track->IsOn(AliESDtrack::kTPCrefit)) return;
3e55050f 1197 if (!friendTrack->GetITSOut()) return;
817766d5 1198 if (!track->GetInnerParam()) return;
1199 if (!track->GetOuterParam()) return;
b842d904 1200 if (track->GetInnerParam()->Pt()<kMinPt) return;
817766d5 1201 // exclude crossing track
1202 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1203 if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ) return;
b96c3aef 1204 if (track->GetInnerParam()->GetX()>90) return;
817766d5 1205 //
1206 AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetInnerParam()));
b842d904 1207 AliExternalTrackParam pITS(*(friendTrack->GetITSOut())); // ITS standalone if possible
1208 AliExternalTrackParam pITS2(*(friendTrack->GetITSOut())); //TPC-ITS track
1209 pITS2.Rotate(pTPC.GetAlpha());
be67055b 1210 // pITS2.PropagateTo(pTPC.GetX(),fMagF);
1211 AliTracker::PropagateTrackToBxByBz(&pITS2,pTPC.GetX(),0.1,0.1,kFALSE);
1212
b842d904 1213 AliESDfriendTrack *itsfriendTrack=0;
1214 //
1215 // try to find standalone ITS track corresponing to the TPC if possible
1216 //
1217 Bool_t hasAlone=kFALSE;
1218 Int_t ntracks=event->GetNumberOfTracks();
1219 for (Int_t i=0; i<ntracks; i++){
1220 AliESDtrack *trackS = event->GetTrack(i);
1221 if (trackS->GetTPCNcls()>0) continue; //continue if has TPC info
6ba68367 1222 itsfriendTrack = esdFriend->GetTrack(i);
b842d904 1223 if (!itsfriendTrack) continue;
1224 if (!itsfriendTrack->GetITSOut()) continue;
1225 if (TMath::Abs(pITS2.GetTgl()-itsfriendTrack->GetITSOut()->GetTgl())> kMaxAngle) continue;
1226 pITS=(*(itsfriendTrack->GetITSOut()));
1227 //
1228 pITS.Rotate(pTPC.GetAlpha());
be67055b 1229 //pITS.PropagateTo(pTPC.GetX(),fMagF);
1230 AliTracker::PropagateTrackToBxByBz(&pITS,pTPC.GetX(),0.1,0.1,kFALSE);
b842d904 1231 if (TMath::Abs(pITS2.GetY()-pITS.GetY())> kMaxDy) continue;
1232 hasAlone=kTRUE;
1233 }
1234 if (!hasAlone) pITS=pITS2;
1235 //
3e55050f 1236 if (TMath::Abs(pITS.GetY()-pTPC.GetY()) >kMaxDy) return;
1237 if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1238 if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1239 //
1240 // 1. Update median and RMS info
1241 //
b322e06a 1242 TVectorD vecDelta(5),vecMedian(5), vecRMS(5);
3e55050f 1243 TVectorD vecDeltaN(5);
1244 Double_t sign=(pITS.GetParameter()[1]>0)? 1.:-1.;
1245 vecDelta[4]=0;
1246 for (Int_t i=0;i<4;i++){
1247 vecDelta[i]=(pITS.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1248 kgdP[i][kglast%kN]=vecDelta[i];
1249 }
1250 kglast=(kglast+1);
1251 Int_t entries=(kglast<kN)?kglast:kN;
1252 for (Int_t i=0;i<4;i++){
1253 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1254 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1255 vecDeltaN[i] = 0;
1256 if (vecRMS[i]>0.){
1257 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1258 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1259 }
1260 }
817766d5 1261 //
3e55050f 1262 // 2. Apply median+-rms cut
817766d5 1263 //
3e55050f 1264 if (kglast<3) return; //median and RMS to be defined
1265 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1266 //
1267 // 3. Update alignment
817766d5 1268 //
1269 Int_t htime = fTime/3600; //time in hours
4486a91f 1270 if (fAlignITSTPC->GetEntriesFast()<htime){
817766d5 1271 fAlignITSTPC->Expand(htime*2+20);
1272 }
1273 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignITSTPC->At(htime);
1274 if (!align){
3e55050f 1275 // make Alignment object if doesn't exist
817766d5 1276 align=new AliRelAlignerKalman();
3e55050f 1277 align->SetRunNumber(fRun);
1278 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1279 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1280 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1281 align->SetRejectOutliers(kFALSE);
1282
1283 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
817766d5 1284 align->SetMagField(fMagF);
1285 fAlignITSTPC->AddAt(align,htime);
1286 }
817766d5 1287 align->AddTrackParams(&pITS,&pTPC);
1288 align->SetTimeStamp(fTime);
3e55050f 1289 align->SetRunNumber(fRun );
391ffdb2 1290 Float_t dca[2],cov[3];
1291 track->GetImpactParameters(dca,cov);
882f5c06 1292 if (TMath::Abs(dca[0])<kMaxDy){
391ffdb2 1293 FillResHistoTPCITS(&pTPC,&pITS);
1294 FillResHistoTPC(track);
1295 }
3e55050f 1296 //
1297 Int_t nupdates=align->GetNUpdates();
1298 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1299 align->SetRejectOutliers(kFALSE);
817766d5 1300 TTreeSRedirector *cstream = GetDebugStreamer();
1301 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
817766d5 1302 TVectorD gpTPC(3), gdTPC(3);
1303 TVectorD gpITS(3), gdITS(3);
1304 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1305 pTPC.GetDirection(gdTPC.GetMatrixArray());
1306 pITS.GetXYZ(gpITS.GetMatrixArray());
1307 pITS.GetDirection(gdITS.GetMatrixArray());
1308 (*cstream)<<"itstpc"<<
1309 "run="<<fRun<< // run number
1310 "event="<<fEvent<< // event number
1311 "time="<<fTime<< // time stamp of event
1312 "trigger="<<fTrigger<< // trigger
1313 "mag="<<fMagF<< // magnetic field
817766d5 1314 //
b842d904 1315 "hasAlone="<<hasAlone<< // has ITS standalone ?
1316 "track.="<<track<< // track info
3e55050f 1317 "nmed="<<kglast<< // number of entries to define median and RMS
1318 "vMed.="<<&vecMedian<< // median of deltas
1319 "vRMS.="<<&vecRMS<< // rms of deltas
1320 "vDelta.="<<&vecDelta<< // delta in respect to median
1321 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1322 "t.="<<track<< // ful track - find proper cuts
1323 "a.="<<align<< // current alignment
b842d904 1324 "pITS.="<<&pITS<< // track param ITS
1325 "pITS2.="<<&pITS2<< // track param ITS+TPC
3e55050f 1326 "pTPC.="<<&pTPC<< // track param TPC
1327 "gpTPC.="<<&gpTPC<< // global position TPC
1328 "gdTPC.="<<&gdTPC<< // global direction TPC
1329 "gpITS.="<<&gpITS<< // global position ITS
1330 "gdITS.="<<&gdITS<< // global position ITS
817766d5 1331 "\n";
1332 }
817766d5 1333}
3e55050f 1334
1335
1336
1337
2811495d 1338void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack *const track, AliESDfriendTrack *const friendTrack){
0dac7b3a 1339 //
3e55050f 1340 // Process track - Update TPC-TRD alignment
1341 // Updates:
1342 // 0. Apply standartd cuts
1343 // 1. Recalucluate the current statistic median/RMS
1344 // 2. Apply median+-rms cut
1345 // 3. Update kalman filter
1346 //
1347 const Int_t kMinTPC = 80; // minimal number of TPC cluster
1348 const Int_t kMinTRD = 50; // minimal number of TRD cluster
1349 const Double_t kMinZ = 20; // maximal dz distance
9112627b 1350 const Double_t kMaxDy = 5.; // maximal dy distance
1351 const Double_t kMaxAngle= 0.1; // maximal angular distance
1352 const Double_t kSigmaCut= 10; // maximal sigma distance to median
3e55050f 1353 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1354 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1355 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
be67055b 1356 const Double_t kRefX = 275; // reference X
1357 const Int_t kN=50; // deepnes of history
3e55050f 1358 static Int_t kglast=0;
1359 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
0dac7b3a 1360 //
3e55050f 1361 // 0. Apply standard cuts
0dac7b3a 1362 //
1363 Int_t dummycl[1000];
1364 if (track->GetTRDclusters(dummycl)<kMinTRD) return; // minimal amount of clusters
1365 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
3e55050f 1366 if (!friendTrack->GetTRDIn()) return;
be67055b 1367 if (!track->IsOn(AliESDtrack::kTRDrefit)) return;
1368 if (!track->IsOn(AliESDtrack::kTRDout)) return;
0dac7b3a 1369 if (!track->GetInnerParam()) return;
be67055b 1370 if (!friendTrack->GetTPCOut()) return;
0dac7b3a 1371 // exclude crossing track
be67055b 1372 if (friendTrack->GetTPCOut()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
3e55050f 1373 if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ) return;
0dac7b3a 1374 //
be67055b 1375 AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(friendTrack->GetTPCOut()));
1376 AliTracker::PropagateTrackToBxByBz(&pTPC,kRefX,0.1,0.1,kFALSE);
0dac7b3a 1377 AliExternalTrackParam pTRD(*(friendTrack->GetTRDIn()));
3e55050f 1378 pTRD.Rotate(pTPC.GetAlpha());
be67055b 1379 // pTRD.PropagateTo(pTPC.GetX(),fMagF);
1380 AliTracker::PropagateTrackToBxByBz(&pTRD,pTPC.GetX(),0.1,0.1,kFALSE);
1381
3e55050f 1382 ((Double_t*)pTRD.GetCovariance())[2]+=3.*3.; // increas sys errors
1383 ((Double_t*)pTRD.GetCovariance())[9]+=0.1*0.1; // increse sys errors
1384
1385 if (TMath::Abs(pTRD.GetY()-pTPC.GetY()) >kMaxDy) return;
1386 if (TMath::Abs(pTRD.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
9112627b 1387 // if (TMath::Abs(pTRD.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
0dac7b3a 1388 //
3e55050f 1389 // 1. Update median and RMS info
0dac7b3a 1390 //
b322e06a 1391 TVectorD vecDelta(5),vecMedian(5), vecRMS(5);
3e55050f 1392 TVectorD vecDeltaN(5);
1393 Double_t sign=(pTRD.GetParameter()[1]>0)? 1.:-1.;
1394 vecDelta[4]=0;
1395 for (Int_t i=0;i<4;i++){
1396 vecDelta[i]=(pTRD.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1397 kgdP[i][kglast%kN]=vecDelta[i];
1398 }
1399 kglast=(kglast+1);
1400 Int_t entries=(kglast<kN)?kglast:kN;
1401 for (Int_t i=0;i<4;i++){
1402 vecMedian[i] = TMath::Median(entries,kgdP[i]);
be67055b 1403
3e55050f 1404 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1405 vecDeltaN[i] = 0;
1406 if (vecRMS[i]>0.){
1407 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1408 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1409 }
1410 }
1411 //
1412 // 2. Apply median+-rms cut
1413 //
1414 if (kglast<3) return; //median and RMS to be defined
1415 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1416 //
1417 // 3. Update alignment
0dac7b3a 1418 //
1419 Int_t htime = fTime/3600; //time in hours
4486a91f 1420 if (fAlignTRDTPC->GetEntriesFast()<htime){
0dac7b3a 1421 fAlignTRDTPC->Expand(htime*2+20);
1422 }
1423 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTRDTPC->At(htime);
1424 if (!align){
3e55050f 1425 // make Alignment object if doesn't exist
0dac7b3a 1426 align=new AliRelAlignerKalman();
3e55050f 1427 align->SetRunNumber(fRun);
1428 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1429 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1430 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1431 align->SetRejectOutliers(kFALSE);
1432 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
0dac7b3a 1433 align->SetMagField(fMagF);
1434 fAlignTRDTPC->AddAt(align,htime);
1435 }
0dac7b3a 1436 align->AddTrackParams(&pTRD,&pTPC);
1437 align->SetTimeStamp(fTime);
3e55050f 1438 align->SetRunNumber(fRun );
0b736a46 1439 Float_t dca[2],cov[3];
1440 track->GetImpactParameters(dca,cov);
1441 if (TMath::Abs(dca[0])<kMaxDy){
1442 FillResHistoTPCTRD(&pTPC,&pTRD); //only primaries
1443 }
3e55050f 1444 //
1445 Int_t nupdates=align->GetNUpdates();
1446 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1447 align->SetRejectOutliers(kFALSE);
0dac7b3a 1448 TTreeSRedirector *cstream = GetDebugStreamer();
1449 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
0dac7b3a 1450 TVectorD gpTPC(3), gdTPC(3);
1451 TVectorD gpTRD(3), gdTRD(3);
1452 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1453 pTPC.GetDirection(gdTPC.GetMatrixArray());
1454 pTRD.GetXYZ(gpTRD.GetMatrixArray());
1455 pTRD.GetDirection(gdTRD.GetMatrixArray());
3e55050f 1456 (*cstream)<<"trdtpc"<<
0dac7b3a 1457 "run="<<fRun<< // run number
1458 "event="<<fEvent<< // event number
1459 "time="<<fTime<< // time stamp of event
1460 "trigger="<<fTrigger<< // trigger
1461 "mag="<<fMagF<< // magnetic field
0dac7b3a 1462 //
3e55050f 1463 "nmed="<<kglast<< // number of entries to define median and RMS
1464 "vMed.="<<&vecMedian<< // median of deltas
1465 "vRMS.="<<&vecRMS<< // rms of deltas
1466 "vDelta.="<<&vecDelta<< // delta in respect to median
1467 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1468 "t.="<<track<< // ful track - find proper cuts
1469 "a.="<<align<< // current alignment
1470 "pTRD.="<<&pTRD<< // track param TRD
1471 "pTPC.="<<&pTPC<< // track param TPC
1472 "gpTPC.="<<&gpTPC<< // global position TPC
1473 "gdTPC.="<<&gdTPC<< // global direction TPC
1474 "gpTRD.="<<&gpTRD<< // global position TRD
1475 "gdTRD.="<<&gdTRD<< // global position TRD
0dac7b3a 1476 "\n";
1477 }
0dac7b3a 1478}
817766d5 1479
1480
2811495d 1481void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack *const track, AliESDfriendTrack *const friendTrack){
817766d5 1482 //
817766d5 1483 //
3e55050f 1484 // Process track - Update TPC-TOF alignment
1485 // Updates:
1486 // -1. Make a TOF "track"
1487 // 0. Apply standartd cuts
1488 // 1. Recalucluate the current statistic median/RMS
1489 // 2. Apply median+-rms cut
1490 // 3. Update kalman filter
1491 //
1492 const Int_t kMinTPC = 80; // minimal number of TPC cluster
b842d904 1493 // const Double_t kMinZ = 10; // maximal dz distance
3e55050f 1494 const Double_t kMaxDy = 5.; // maximal dy distance
b96c3aef 1495 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1496 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1497 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1498 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1499
1500 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
be67055b 1501 const Int_t kN=50; // deepnes of history
3e55050f 1502 static Int_t kglast=0;
1503 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
817766d5 1504 //
3e55050f 1505 // -1. Make a TOF track-
1506 // Clusters are not in friends - use alingment points
1507 //
1508 if (track->GetTOFsignal()<=0) return;
1509 if (!friendTrack->GetTPCOut()) return;
1510 if (!track->GetInnerParam()) return;
be67055b 1511 if (!friendTrack->GetTPCOut()) return;
817766d5 1512 const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
1513 if (!points) return;
be67055b 1514 AliExternalTrackParam pTPC(*(friendTrack->GetTPCOut()));
3e55050f 1515 AliExternalTrackParam pTOF(pTPC);
1516 Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
817766d5 1517 Int_t npoints = points->GetNPoints();
1518 AliTrackPoint point;
3e55050f 1519 Int_t naccept=0;
817766d5 1520 //
1521 for (Int_t ipoint=0;ipoint<npoints;ipoint++){
817766d5 1522 points->GetPoint(point,ipoint);
1523 Float_t xyz[3];
1524 point.GetXYZ(xyz);
1525 Double_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
3e55050f 1526 if (r<350) continue;
817766d5 1527 if (r>400) continue;
3e55050f 1528 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,2.,kTRUE);
1529 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,0.1,kTRUE);
1530 AliTrackPoint lpoint = point.Rotate(pTPC.GetAlpha());
1531 pTPC.PropagateTo(lpoint.GetX(),fMagF);
1532 pTOF=pTPC;
817766d5 1533 ((Double_t*)pTOF.GetParameter())[0] =lpoint.GetY();
1534 ((Double_t*)pTOF.GetParameter())[1] =lpoint.GetZ();
3e55050f 1535 ((Double_t*)pTOF.GetCovariance())[0]+=3.*3./12.;
1536 ((Double_t*)pTOF.GetCovariance())[2]+=3.*3./12.;
0dac7b3a 1537 ((Double_t*)pTOF.GetCovariance())[5]+=0.1*0.1;
1538 ((Double_t*)pTOF.GetCovariance())[9]+=0.1*0.1;
3e55050f 1539 naccept++;
1540 }
1541 if (naccept==0) return; // no tof match clusters
1542 //
1543 // 0. Apply standard cuts
1544 //
1545 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
1546 // exclude crossing track
be67055b 1547 if (friendTrack->GetTPCOut()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
3e55050f 1548 //
1549 if (TMath::Abs(pTOF.GetY()-pTPC.GetY()) >kMaxDy) return;
1550 if (TMath::Abs(pTOF.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1551 if (TMath::Abs(pTOF.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1552 //
1553 // 1. Update median and RMS info
1554 //
b322e06a 1555 TVectorD vecDelta(5),vecMedian(5), vecRMS(5);
3e55050f 1556 TVectorD vecDeltaN(5);
1557 Double_t sign=(pTOF.GetParameter()[1]>0)? 1.:-1.;
1558 vecDelta[4]=0;
1559 for (Int_t i=0;i<4;i++){
1560 vecDelta[i]=(pTOF.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1561 kgdP[i][kglast%kN]=vecDelta[i];
1562 }
1563 kglast=(kglast+1);
1564 Int_t entries=(kglast<kN)?kglast:kN;
1565 Bool_t isOK=kTRUE;
1566 for (Int_t i=0;i<4;i++){
1567 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1568 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1569 vecDeltaN[i] = 0;
1570 if (vecRMS[i]>0.){
1571 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/(vecRMS[i]+1.);
1572 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1573 if (TMath::Abs(vecDeltaN[i])>kSigmaCut) isOK=kFALSE;
817766d5 1574 }
3e55050f 1575 }
1576 //
1577 // 2. Apply median+-rms cut
1578 //
1579 if (kglast<10) return; //median and RMS to be defined
1580 if (!isOK) return;
1581 //
1582 // 3. Update alignment
1583 //
1584 Int_t htime = fTime/3600; //time in hours
4486a91f 1585 if (fAlignTOFTPC->GetEntriesFast()<htime){
3e55050f 1586 fAlignTOFTPC->Expand(htime*2+20);
1587 }
1588 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTOFTPC->At(htime);
1589 if (!align){
1590 // make Alignment object if doesn't exist
1591 align=new AliRelAlignerKalman();
1592 align->SetRunNumber(fRun);
1593 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1594 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1595 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1596 align->SetRejectOutliers(kFALSE);
1597 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
1598 align->SetMagField(fMagF);
1599 fAlignTOFTPC->AddAt(align,htime);
1600 }
1601 align->AddTrackParams(&pTOF,&pTPC);
0b736a46 1602 Float_t dca[2],cov[3];
1603 track->GetImpactParameters(dca,cov);
1604 if (TMath::Abs(dca[0])<kMaxDy){
1605 FillResHistoTPCTOF(&pTPC,&pTOF);
1606 }
3e55050f 1607 align->SetTimeStamp(fTime);
1608 align->SetRunNumber(fRun );
1609 //
1610 Int_t nupdates=align->GetNUpdates();
1611 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1612 align->SetRejectOutliers(kFALSE);
1613 TTreeSRedirector *cstream = GetDebugStreamer();
1614 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
3e55050f 1615 TVectorD gpTPC(3), gdTPC(3);
1616 TVectorD gpTOF(3), gdTOF(3);
1617 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1618 pTPC.GetDirection(gdTPC.GetMatrixArray());
1619 pTOF.GetXYZ(gpTOF.GetMatrixArray());
1620 pTOF.GetDirection(gdTOF.GetMatrixArray());
1621 (*cstream)<<"toftpc"<<
1622 "run="<<fRun<< // run number
1623 "event="<<fEvent<< // event number
1624 "time="<<fTime<< // time stamp of event
1625 "trigger="<<fTrigger<< // trigger
1626 "mag="<<fMagF<< // magnetic field
3e55050f 1627 //
1628 "nmed="<<kglast<< // number of entries to define median and RMS
1629 "vMed.="<<&vecMedian<< // median of deltas
1630 "vRMS.="<<&vecRMS<< // rms of deltas
1631 "vDelta.="<<&vecDelta<< // delta in respect to median
1632 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1633 "t.="<<track<< // ful track - find proper cuts
1634 "a.="<<align<< // current alignment
1635 "pTOF.="<<&pTOF<< // track param TOF
1636 "pTPC.="<<&pTPC<< // track param TPC
1637 "gpTPC.="<<&gpTPC<< // global position TPC
1638 "gdTPC.="<<&gdTPC<< // global direction TPC
1639 "gpTOF.="<<&gpTOF<< // global position TOF
1640 "gdTOF.="<<&gdTOF<< // global position TOF
1641 "\n";
817766d5 1642 }
817766d5 1643}
3e55050f 1644
1645
391ffdb2 1646void AliTPCcalibTime::BookDistortionMaps(){
1647 //
1648 // Book ndimensional histograms of distortions/residuals
1649 // Only primary tracks are selected for analysis
1650 //
1651
1652 Double_t xminTrack[4], xmaxTrack[4];
1653 Int_t binsTrack[4];
1654 TString axisName[4];
1655 TString axisTitle[4];
1656 //
882f5c06 1657 binsTrack[0] =50;
391ffdb2 1658 axisName[0] ="#Delta";
1659 axisTitle[0] ="#Delta";
1660 //
1661 binsTrack[1] =20;
1662 xminTrack[1] =-1.5; xmaxTrack[1]=1.5;
1663 axisName[1] ="tanTheta";
1664 axisTitle[1] ="tan(#Theta)";
1665 //
1666 binsTrack[2] =90;
1667 xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi();
1668 axisName[2] ="phi";
1669 axisTitle[2] ="#phi";
1670 //
1671 binsTrack[3] =20;
1672 xminTrack[3] =-1.; xmaxTrack[3]=1.; // 0.33 GeV cut
1673 axisName[3] ="snp";
1674 //
1675 // delta y
be67055b 1676 xminTrack[0] =-1.5; xmaxTrack[0]=1.5; //
391ffdb2 1677 fResHistoTPCITS[0] = new THnSparseS("TPCITS#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1678 fResHistoTPCvertex[0] = new THnSparseS("TPCVertex#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
be67055b 1679 xminTrack[0] =-1.5; xmaxTrack[0]=1.5; //
391ffdb2 1680 fResHistoTPCTRD[0] = new THnSparseS("TPCTRD#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
0b736a46 1681 xminTrack[0] =-5; xmaxTrack[0]=5; //
1682 fResHistoTPCTOF[0] = new THnSparseS("TPCTOF#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1683 //
1684 // delta z
9112627b 1685 xminTrack[0] =-3.; xmaxTrack[0]=3.; //
391ffdb2 1686 fResHistoTPCITS[1] = new THnSparseS("TPCITS#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1687 fResHistoTPCvertex[1] = new THnSparseS("TPCVertex#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1688 fResHistoTPCTRD[1] = new THnSparseS("TPCTRD#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
0b736a46 1689 xminTrack[0] =-5.; xmaxTrack[0]=5.; //
1690 fResHistoTPCTOF[1] = new THnSparseS("TPCTOF#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1691 //
1692 // delta snp-P2
882f5c06 1693 xminTrack[0] =-0.015; xmaxTrack[0]=0.015; //
391ffdb2 1694 fResHistoTPCITS[2] = new THnSparseS("TPCITS#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
9112627b 1695 fResHistoTPCvertex[2] = new THnSparseS("TPCITSv#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1696 fResHistoTPCTRD[2] = new THnSparseS("TPCTRD#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
0b736a46 1697 fResHistoTPCTOF[2] = new THnSparseS("TPCTOF#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1698 //
1699 // delta theta-P3
9112627b 1700 xminTrack[0] =-0.025; xmaxTrack[0]=0.025; //
391ffdb2 1701 fResHistoTPCITS[3] = new THnSparseS("TPCITS#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
9112627b 1702 fResHistoTPCvertex[3] = new THnSparseS("TPCITSv#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1703 fResHistoTPCTRD[3] = new THnSparseS("TPCTRD#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
0b736a46 1704 fResHistoTPCTOF[3] = new THnSparseS("TPCTOF#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1705 //
1706 // delta theta-P4
9112627b 1707 xminTrack[0] =-0.2; xmaxTrack[0]=0.2; //
391ffdb2 1708 fResHistoTPCITS[4] = new THnSparseS("TPCITS#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
9112627b 1709 fResHistoTPCvertex[4] = new THnSparseS("TPCITSv#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1710 fResHistoTPCTRD[4] = new THnSparseS("TPCTRD#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
0b736a46 1711 fResHistoTPCTOF[4] = new THnSparseS("TPCTOF#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1712 //
1713 for (Int_t ivar=0;ivar<4;ivar++){
1714 for (Int_t ivar2=0;ivar2<4;ivar2++){
1715 fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1716 fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
1717 fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1718 fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
1719 fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1720 fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
1721 }
1722 }
1723}
1724
1725
1726void AliTPCcalibTime::FillResHistoTPCITS(const AliExternalTrackParam * pTPCIn, const AliExternalTrackParam * pITSOut ){
1727 //
1728 // fill residual histograms pTPCIn-pITSOut
1729 // Histogram is filled only for primary tracks
1730 //
1731 Double_t histoX[4];
1732 Double_t xyz[3];
1733 pTPCIn->GetXYZ(xyz);
1734 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1735 histoX[1]= pTPCIn->GetTgl();
1736 histoX[2]= phi;
1737 histoX[3]= pTPCIn->GetSnp();
9112627b 1738 AliExternalTrackParam lits(*pITSOut);
1739 lits.Rotate(pTPCIn->GetAlpha());
1740 lits.PropagateTo(pTPCIn->GetX(),fMagF);
391ffdb2 1741 //
1742 for (Int_t ihisto=0; ihisto<5; ihisto++){
9112627b 1743 histoX[0]=pTPCIn->GetParameter()[ihisto]-lits.GetParameter()[ihisto];
391ffdb2 1744 fResHistoTPCITS[ihisto]->Fill(histoX);
1745 }
1746}
1747
1748
1749void AliTPCcalibTime::FillResHistoTPC(const AliESDtrack * pTrack){
1750 //
1751 // fill residual histograms pTPC - vertex
1752 // Histogram is filled only for primary tracks
1753 //
1754 Double_t histoX[4];
1755 const AliExternalTrackParam * pTPCIn = pTrack->GetInnerParam();
be67055b 1756 AliExternalTrackParam pTPCvertex(*(pTrack->GetInnerParam()));
9112627b 1757 //
1758 AliExternalTrackParam lits(*pTrack);
be67055b 1759 if (TMath::Abs(pTrack->GetY())>3) return; // beam pipe
1760 pTPCvertex.Rotate(lits.GetAlpha());
1761 //pTPCvertex.PropagateTo(pTPCvertex->GetX(),fMagF);
1762 AliTracker::PropagateTrackToBxByBz(&pTPCvertex,lits.GetX(),0.1,2,kFALSE);
1763 AliTracker::PropagateTrackToBxByBz(&pTPCvertex,lits.GetX(),0.1,0.1,kFALSE);
391ffdb2 1764 Double_t xyz[3];
1765 pTPCIn->GetXYZ(xyz);
1766 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1767 histoX[1]= pTPCIn->GetTgl();
1768 histoX[2]= phi;
1769 histoX[3]= pTPCIn->GetSnp();
1770 //
1771 Float_t dca[2], cov[3];
1772 pTrack->GetImpactParametersTPC(dca,cov);
882f5c06 1773 for (Int_t ihisto=0; ihisto<5; ihisto++){
be67055b 1774 histoX[0]=pTPCvertex.GetParameter()[ihisto]-lits.GetParameter()[ihisto];
1775 // if (ihisto<2) histoX[0]=dca[ihisto];
9112627b 1776 fResHistoTPCvertex[ihisto]->Fill(histoX);
391ffdb2 1777 }
1778}
1779
1780
1781void AliTPCcalibTime::FillResHistoTPCTRD(const AliExternalTrackParam * pTPCOut, const AliExternalTrackParam * pTRDIn ){
1782 //
1783 // fill resuidual histogram TPCout-TRDin
1784 //
1785 Double_t histoX[4];
1786 Double_t xyz[3];
1787 pTPCOut->GetXYZ(xyz);
1788 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1789 histoX[1]= pTPCOut->GetTgl();
1790 histoX[2]= phi;
1791 histoX[3]= pTPCOut->GetSnp();
1792 //
882f5c06 1793 AliExternalTrackParam ltrd(*pTRDIn);
1794 ltrd.Rotate(pTPCOut->GetAlpha());
be67055b 1795 // ltrd.PropagateTo(pTPCOut->GetX(),fMagF);
1796 AliTracker::PropagateTrackToBxByBz(&ltrd,pTPCOut->GetX(),0.1,0.1,kFALSE);
882f5c06 1797
391ffdb2 1798 for (Int_t ihisto=0; ihisto<5; ihisto++){
882f5c06 1799 histoX[0]=pTPCOut->GetParameter()[ihisto]-ltrd.GetParameter()[ihisto];
391ffdb2 1800 fResHistoTPCTRD[ihisto]->Fill(histoX);
1801 }
1802
1803}
0b736a46 1804
1805void AliTPCcalibTime::FillResHistoTPCTOF(const AliExternalTrackParam * pTPCOut, const AliExternalTrackParam * pTOFIn ){
1806 //
1807 // fill resuidual histogram TPCout-TOFin
1808 // track propagated to the TOF position
1809 Double_t histoX[4];
1810 Double_t xyz[3];
1811
1812 AliExternalTrackParam ltpc(*pTPCOut);
1813 ltpc.Rotate(pTOFIn->GetAlpha());
1814 AliTracker::PropagateTrackToBxByBz(&ltpc,pTOFIn->GetX(),0.1,0.1,kFALSE);
1815 //
1816 ltpc.GetXYZ(xyz);
1817 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1818 histoX[1]= ltpc.GetTgl();
1819 histoX[2]= phi;
1820 histoX[3]= ltpc.GetSnp();
1821 //
1822 for (Int_t ihisto=0; ihisto<2; ihisto++){
1823 histoX[0]=ltpc.GetParameter()[ihisto]-pTOFIn->GetParameter()[ihisto];
1824 fResHistoTPCTOF[ihisto]->Fill(histoX);
1825 }
1826
1827}