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