new analysis module JCORRAN - Filip Krizek
[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 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
da6c0bc9 410 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
411 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
412 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
413 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
414 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
32100b1d 415 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
31aa7c5c 416 TVectorD vecGoofie(20);
417 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
2be25cc9 418 if (goofieArray){
31aa7c5c 419 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
420 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
421 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
422 }
da6c0bc9 423 }
2be25cc9 424 (*cstream)<<"laserInfo"<<
da6c0bc9 425 "run="<<fRun<< // run number
426 "event="<<fEvent<< // event number
427 "time="<<fTime<< // time stamp of event
428 "trigger="<<fTrigger<< // trigger
429 "mag="<<fMagF<< // magnetic field
430 // Environment values
431 "press0="<<valuePressure0<<
432 "press1="<<valuePressure1<<
433 "pt0="<<ptrelative0<<
434 "pt1="<<ptrelative1<<
435 "temp0="<<temp0<<
436 "temp1="<<temp1<<
817766d5 437 "vecGoofie.="<<&vecGoofie<<
32100b1d 438 "vdcorr="<<vdcorr<<
da6c0bc9 439 //laser
dde68d36 440 "rejectA="<<isReject[0]<<
441 "rejectC="<<isReject[1]<<
2be25cc9 442 "laserA.="<<fLaser->fFitAside<<
443 "laserC.="<<fLaser->fFitCside<<
444 "laserAC.="<<fLaser->fFitACside<<
445 "trigger="<<event->GetFiredTriggerClasses()<<
da6c0bc9 446 "\n";
447 }
448 }
2be25cc9 449 //
2be25cc9 450 // fill histos
451 //
452 TVectorD vdriftA(5), vdriftC(5),vdriftAC(5);
453 vdriftA=*(fLaser->fFitAside);
454 vdriftC=*(fLaser->fFitCside);
455 vdriftAC=*(fLaser->fFitACside);
456 Int_t npointsA=0, npointsC=0;
457 Float_t chi2A=0, chi2C=0;
458 npointsA= TMath::Nint(vdriftA[3]);
459 chi2A= vdriftA[4];
460 npointsC= TMath::Nint(vdriftC[3]);
461 chi2C= vdriftC[4];
2be25cc9 462
2be25cc9 463 TTimeStamp tstamp(fTime);
464 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
465 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
dde68d36 466 Double_t driftA=0, driftC=0;
467 if (vdriftA[1]>1.-kMaxDeltaV) driftA = 1./vdriftA[1]-1.;
468 if (vdriftC[1]>1.-kMaxDeltaV) driftC = 1./vdriftC[1]-1.;
469 //
470 Double_t vecDriftLaserA[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftA,event->GetRunNumber()};
471 Double_t vecDriftLaserC[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftC,event->GetRunNumber()};
472 // Double_t vecDrift[4] ={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitACside))[1])-1,event->GetRunNumber()};
473
474 for (Int_t icalib=0;icalib<3;icalib++){
475 if (icalib==0){ //z0 shift
476 vecDriftLaserA[2]=vdriftA[0]/250.;
477 vecDriftLaserC[2]=vdriftC[0]/250.;
2be25cc9 478 }
dde68d36 479 if (icalib==1){ //vdrel shift
480 vecDriftLaserA[2]=driftA;
481 vecDriftLaserC[2]=driftC;
2be25cc9 482 }
dde68d36 483 if (icalib==2){ //gy shift - full gy - full drift
484 vecDriftLaserA[2]=vdriftA[2]/250.;
485 vecDriftLaserC[2]=vdriftC[2]/250.;
2be25cc9 486 }
391ffdb2 487 //if (isReject[0]==0) fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
488 //if (isReject[1]==0) fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
489 fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
490 fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
2be25cc9 491 }
dde68d36 492
493// THnSparse* curHist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
494// TString shortName=curHist->ClassName();
495// shortName+="_MEAN_DRIFT_LASER_";
496// delete curHist;
497// curHist=NULL;
498// TString name="";
499
500// name=shortName;
501// name+=event->GetFiredTriggerClasses();
502// name.ToUpper();
503// curHist=(THnSparseF*)fArrayDz->FindObject(name);
504// if(!curHist){
505// curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
506// fArrayDz->AddLast(curHist);
507// }
508// curHist->Fill(vecDrift);
509
510// name=shortName;
511// name+="ALL";
512// name.ToUpper();
513// curHist=(THnSparseF*)fArrayDz->FindObject(name);
514// if(!curHist){
515// curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
516// fArrayDz->AddLast(curHist);
517// }
518// curHist->Fill(vecDrift);
2be25cc9 519}
c74c5f6c 520
2811495d 521void AliTPCcalibTime::ProcessCosmic(const AliESDEvent *const event){
6ba68367 522 //
523 // process Cosmic event - track matching A side C side
524 //
c74c5f6c 525 if (!event) {
526 Printf("ERROR: ESD not available");
527 return;
528 }
529 if (event->GetTimeStamp() == 0 ) {
530 Printf("no time stamp!");
531 return;
532 }
74235403 533
2be25cc9 534 //fd
c74c5f6c 535 // Find cosmic pairs
536 //
537 // Track0 is choosen in upper TPC part
538 // Track1 is choosen in lower TPC part
539 //
cc65e4f5 540 const Int_t kMinClustersCross =30;
541 const Int_t kMinClusters =80;
c74c5f6c 542 Int_t ntracks=event->GetNumberOfTracks();
543 if (ntracks==0) return;
74235403 544 if (ntracks > fCutTracks) return;
545
a390f11f 546 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
6ba68367 547 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
c74c5f6c 548
549 TObjArray tpcSeeds(ntracks);
550 Double_t vtxx[3]={0,0,0};
551 Double_t svtxx[3]={0.000001,0.000001,100.};
552 AliESDVertex vtx(vtxx,svtxx);
553 //
554 // track loop
555 //
cc65e4f5 556 TArrayI clusterSideA(ntracks);
557 TArrayI clusterSideC(ntracks);
c74c5f6c 558 for (Int_t i=0;i<ntracks;++i) {
cc65e4f5 559 clusterSideA[i]=0;
560 clusterSideC[i]=0;
74235403 561 AliESDtrack *track = event->GetTrack(i);
562
563 const AliExternalTrackParam * trackIn = track->GetInnerParam();
564 const AliExternalTrackParam * trackOut = track->GetOuterParam();
565 if (!trackIn) continue;
566 if (!trackOut) continue;
567
6ba68367 568 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
b1c27e4f 569 if (!friendTrack) continue;
3e55050f 570 if (friendTrack) ProcessSame(track,friendTrack,event);
6ba68367 571 if (friendTrack) ProcessAlignITS(track,friendTrack,event,esdFriend);
0dac7b3a 572 if (friendTrack) ProcessAlignTRD(track,friendTrack);
817766d5 573 if (friendTrack) ProcessAlignTOF(track,friendTrack);
74235403 574 TObject *calibObject;
575 AliTPCseed *seed = 0;
576 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
cc65e4f5 577 if (seed) {
578 tpcSeeds.AddAt(seed,i);
579 Int_t nA=0, nC=0;
580 for (Int_t irow=159;irow>0;irow--) {
581 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
582 if (!cl) continue;
583 if ((cl->GetDetector()%36)<18) nA++;
584 if ((cl->GetDetector()%36)>=18) nC++;
585 }
586 clusterSideA[i]=nA;
587 clusterSideC[i]=nC;
588 }
c74c5f6c 589 }
c74c5f6c 590 if (ntracks<2) return;
591 //
592 // Find pairs
593 //
d3ce44cb 594
c74c5f6c 595 for (Int_t i=0;i<ntracks;++i) {
74235403 596 AliESDtrack *track0 = event->GetTrack(i);
c74c5f6c 597 // track0 - choosen upper part
598 if (!track0) continue;
599 if (!track0->GetOuterParam()) continue;
600 if (track0->GetOuterParam()->GetAlpha()<0) continue;
601 Double_t d1[3];
602 track0->GetDirection(d1);
603 for (Int_t j=0;j<ntracks;++j) {
604 if (i==j) continue;
74235403 605 AliESDtrack *track1 = event->GetTrack(j);
606 //track 1 lower part
607 if (!track1) continue;
608 if (!track1->GetOuterParam()) continue;
cc65e4f5 609 if (track0->GetTPCNcls()+ track1->GetTPCNcls()< kMinClusters) continue;
610 Int_t nAC = TMath::Max( TMath::Min(clusterSideA[i], clusterSideC[j]),
611 TMath::Min(clusterSideC[i], clusterSideA[j]));
612 if (nAC<kMinClustersCross) continue;
613 Int_t nA0=clusterSideA[i];
614 Int_t nC0=clusterSideC[i];
615 Int_t nA1=clusterSideA[j];
616 Int_t nC1=clusterSideC[j];
3e55050f 617 // if (track1->GetOuterParam()->GetAlpha()>0) continue;
74235403 618 //
619 Double_t d2[3];
620 track1->GetDirection(d2);
621
622 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
623 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
624 if (! seed0) continue;
625 if (! seed1) continue;
626 Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
627 Float_t dist0 = track0->GetLinearD(0,0);
628 Float_t dist1 = track1->GetLinearD(0,0);
629 //
630 // conservative cuts - convergence to be guarantied
631 // applying before track propagation
3e55050f 632 if (TMath::Abs(TMath::Abs(dist0)-TMath::Abs(dist1))>fCutMaxD) continue; // distance to the 0,0
633 if (TMath::Abs(dir)<TMath::Abs(fCutMinDir)) continue; // direction vector product
74235403 634 Float_t bz = AliTracker::GetBz();
635 Float_t dvertex0[2]; //distance to 0,0
636 Float_t dvertex1[2]; //distance to 0,0
637 track0->GetDZ(0,0,0,bz,dvertex0);
638 track1->GetDZ(0,0,0,bz,dvertex1);
639 if (TMath::Abs(dvertex0[1])>250) continue;
640 if (TMath::Abs(dvertex1[1])>250) continue;
641 //
642 //
643 //
644 Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
645 AliExternalTrackParam param0(*track0);
646 AliExternalTrackParam param1(*track1);
647 //
648 // Propagate using Magnetic field and correct fo material budget
649 //
650 AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
651 AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
652 //
653 // Propagate rest to the 0,0 DCA - z should be ignored
654 //
655 //Bool_t b0 = ;
656 param0.PropagateToDCA(&vtx,bz,1000);
657 //Bool_t b1 =
658 param1.PropagateToDCA(&vtx,bz,1000);
74235403 659 param0.GetDZ(0,0,0,bz,dvertex0);
660 param1.GetDZ(0,0,0,bz,dvertex1);
dde68d36 661 Double_t xyz0[3];
662 Double_t xyz1[3];
74235403 663 param0.GetXYZ(xyz0);
664 param1.GetXYZ(xyz1);
665 Bool_t isPair = IsPair(&param0,&param1);
d3ce44cb 666 Bool_t isCross = IsCross(track0, track1);
dde68d36 667 Bool_t isSame = IsSame(track0, track1);
d3ce44cb 668
dde68d36 669 THnSparse* hist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
670 TString shortName=hist->ClassName();
671 shortName+="_MEAN_VDRIFT_COSMICS_";
672 delete hist;
673 hist=NULL;
674
3e55050f 675 if((isSame) || (isCross && isPair)){
676 if (track0->GetTPCNcls()+ track1->GetTPCNcls()> 80) {
74235403 677 fDz = param0.GetZ() - param1.GetZ();
cc65e4f5 678 Double_t sign=(nA0>nA1)? 1:-1;
679 fDz*=sign;
74235403 680 TTimeStamp tstamp(fTime);
681 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
682 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
dde68d36 683 Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
684 THnSparse* curHist=NULL;
685 TString name="";
686
687 name=shortName;
688 name+=event->GetFiredTriggerClasses();
689 name.ToUpper();
690 curHist=(THnSparseF*)fArrayDz->FindObject(name);
74235403 691 if(!curHist){
dde68d36 692 curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
693 fArrayDz->AddLast(curHist);
74235403 694 }
dde68d36 695// curHist=(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
696// if(!curHist){
697// curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
698// fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
699// }
700 curHist->Fill(vecDrift);
74235403 701
dde68d36 702 name=shortName;
703 name+="ALL";
704 name.ToUpper();
705 curHist=(THnSparseF*)fArrayDz->FindObject(name);
74235403 706 if(!curHist){
dde68d36 707 curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
708 fArrayDz->AddLast(curHist);
74235403 709 }
dde68d36 710// curHist=(THnSparseF*)(fMapDz->GetValue("all"));
711// if(!curHist){
712// curHist=new THnSparseF("all","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
713// fMapDz->Add(new TObjString("all"),curHist);
714// }
715 curHist->Fill(vecDrift);
74235403 716 }
717 }
d3ce44cb 718 TTreeSRedirector *cstream = GetDebugStreamer();
719 if (fStreamLevel>0){
720 if (cstream){
721 (*cstream)<<"trackInfo"<<
722 "tr0.="<<track0<<
723 "tr1.="<<track1<<
724 "p0.="<<&param0<<
725 "p1.="<<&param1<<
cc65e4f5 726 "nAC="<<nAC<<
727 "nA0="<<nA0<<
728 "nA1="<<nA1<<
729 "nC0="<<nC0<<
730 "nC1="<<nC1<<
d3ce44cb 731 "isPair="<<isPair<<
732 "isCross="<<isCross<<
dde68d36 733 "isSame="<<isSame<<
d3ce44cb 734 "fDz="<<fDz<<
735 "fRun="<<fRun<<
736 "fTime="<<fTime<<
737 "\n";
738 }
739 }
c74c5f6c 740 } // end 2nd order loop
741 } // end 1st order loop
74235403 742
2be25cc9 743 if (fStreamLevel>0){
744 TTreeSRedirector *cstream = GetDebugStreamer();
745 if (cstream){
746 TTimeStamp tstamp(fTime);
747 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
748 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
749 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
750 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
751 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
752 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
32100b1d 753 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
2be25cc9 754 TVectorD vecGoofie(20);
755 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
756 if (goofieArray){
757 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
758 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
759 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
760 }
761 }
762 (*cstream)<<"timeInfo"<<
763 "run="<<fRun<< // run number
764 "event="<<fEvent<< // event number
765 "time="<<fTime<< // time stamp of event
766 "trigger="<<fTrigger<< // trigger
767 "mag="<<fMagF<< // magnetic field
768 // Environment values
769 "press0="<<valuePressure0<<
770 "press1="<<valuePressure1<<
771 "pt0="<<ptrelative0<<
772 "pt1="<<ptrelative1<<
773 "temp0="<<temp0<<
774 "temp1="<<temp1<<
775 "vecGoofie.=<<"<<&vecGoofie<<
32100b1d 776 "vdcorr="<<vdcorr<<
2be25cc9 777 //
778 // accumulated values
779 //
780 "fDz="<<fDz<< //! current delta z
2be25cc9 781 "trigger="<<event->GetFiredTriggerClasses()<<
782 "\n";
783 }
784 }
a390f11f 785 if (GetDebugLevel()>20) printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
c74c5f6c 786}
787
2811495d 788void AliTPCcalibTime::ProcessBeam(const AliESDEvent *const /*event*/){
6ba68367 789 //
790 // Not special treatment yet - the same for cosmic and physic event
791 //
dde68d36 792}
793
6ba68367 794void AliTPCcalibTime::Analyze(){
795 //
796 // Special macro to analyze result of calibration and extract calibration entries
797 // Not yet ported to the Analyze function yet
798 //
799}
74235403 800
2811495d 801THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name) const
802{
6ba68367 803 //
804 // Get histogram for given trigger mask
805 //
dde68d36 806 TIterator* iterator = fArrayDz->MakeIterator();
807 iterator->Reset();
808 TString newName=name;
809 newName.ToUpper();
810 THnSparse* newHist=new THnSparseF(newName,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
811 THnSparse* addHist=NULL;
812 while((addHist=(THnSparseF*)iterator->Next())){
813 if(!addHist) continue;
814 TString histName=addHist->GetName();
815 if(!histName.Contains(newName)) continue;
816 addHist->Print();
817 newHist->Add(addHist);
74235403 818 }
dde68d36 819 return newHist;
74235403 820}
821
2811495d 822TObjArray* AliTPCcalibTime::GetHistoDrift() const
823{
6ba68367 824 //
825 // return array of histograms
826 //
dde68d36 827 return fArrayDz;
74235403 828}
829
830TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){
6ba68367 831 //
832 // Make a drift velocity (delta Z) graph
833 //
dde68d36 834 THnSparse* histoDrift=GetHistoDrift(name);
835 TGraphErrors* graphDrift=NULL;
836 if(histoDrift){
837 graphDrift=FitSlices(histoDrift,2,0,400,100,0.05,0.95, kTRUE);
838 TString end=histoDrift->GetName();
839 Int_t pos=end.Index("_");
840 end=end(pos,end.Capacity()-pos);
841 TString graphName=graphDrift->ClassName();
842 graphName+=end;
843 graphName.ToUpper();
844 graphDrift->SetName(graphName);
74235403 845 }
846 return graphDrift;
847}
848
dde68d36 849TObjArray* AliTPCcalibTime::GetGraphDrift(){
6ba68367 850 //
851 // make a array of drift graphs
852 //
dde68d36 853 TObjArray* arrayGraphDrift=new TObjArray();
854 TIterator* iterator=fArrayDz->MakeIterator();
74235403 855 iterator->Reset();
dde68d36 856 THnSparse* addHist=NULL;
857 while((addHist=(THnSparseF*)iterator->Next())) arrayGraphDrift->AddLast(GetGraphDrift(addHist->GetName()));
858 return arrayGraphDrift;
74235403 859}
860
dde68d36 861AliSplineFit* AliTPCcalibTime::GetFitDrift(const char* name){
6ba68367 862 //
863 // Make a fit AliSplinefit of drift velocity
864 //
d3ce44cb 865 TGraph* graphDrift=GetGraphDrift(name);
dde68d36 866 AliSplineFit* fitDrift=NULL;
74235403 867 if(graphDrift && graphDrift->GetN()){
dde68d36 868 fitDrift=new AliSplineFit();
869 fitDrift->SetGraph(graphDrift);
870 fitDrift->SetMinPoints(graphDrift->GetN()+1);
871 fitDrift->InitKnots(graphDrift,2,0,0.001);
872 fitDrift->SplineFit(0);
873 TString end=graphDrift->GetName();
874 Int_t pos=end.Index("_");
875 end=end(pos,end.Capacity()-pos);
876 TString fitName=fitDrift->ClassName();
877 fitName+=end;
878 fitName.ToUpper();
879 //fitDrift->SetName(fitName);
74235403 880 delete graphDrift;
dde68d36 881 graphDrift=NULL;
74235403 882 }
883 return fitDrift;
884}
885
74235403 886
2811495d 887Long64_t AliTPCcalibTime::Merge(TCollection *const li) {
6ba68367 888 //
889 // Object specific merging procedure
890 //
c74c5f6c 891 TIterator* iter = li->MakeIterator();
892 AliTPCcalibTime* cal = 0;
893
894 while ((cal = (AliTPCcalibTime*)iter->Next())) {
895 if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
896 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
897 return -1;
898 }
2be25cc9 899 for (Int_t imeas=0; imeas<3; imeas++){
900 if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
901 fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
902 fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
903 }
904 }
dde68d36 905 TObjArray* addArray=cal->GetHistoDrift();
906 if(!addArray) return 0;
907 TIterator* iterator = addArray->MakeIterator();
2be25cc9 908 iterator->Reset();
dde68d36 909 THnSparse* addHist=NULL;
910 while((addHist=(THnSparseF*)iterator->Next())){
911 if(!addHist) continue;
2be25cc9 912 addHist->Print();
dde68d36 913 THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName());
2be25cc9 914 if(!localHist){
915 localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
dde68d36 916 fArrayDz->AddLast(localHist);
2be25cc9 917 }
918 localHist->Add(addHist);
919 }
391ffdb2 920
d3ce44cb 921 for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
3e55050f 922 //
923 // Merge alignment
924 //
925 for (Int_t itype=0; itype<3; itype++){
926 //
927 //
928 TObjArray *arr0= 0;
929 TObjArray *arr1= 0;
930 if (itype==0) {arr0=fAlignITSTPC; arr1=cal->fAlignITSTPC;}
931 if (itype==1) {arr0=fAlignTRDTPC; arr1=cal->fAlignTRDTPC;}
932 if (itype==2) {arr0=fAlignTOFTPC; arr1=cal->fAlignTOFTPC;}
933 if (!arr1) continue;
934 if (!arr0) arr0=new TObjArray(arr1->GetEntriesFast());
935 if (arr1->GetEntriesFast()>arr0->GetEntriesFast()){
936 arr0->Expand(arr1->GetEntriesFast());
937 }
938 for (Int_t i=0;i<arr1->GetEntriesFast(); i++){
939 AliRelAlignerKalman *kalman1 = (AliRelAlignerKalman *)arr1->UncheckedAt(i);
940 AliRelAlignerKalman *kalman0 = (AliRelAlignerKalman *)arr0->UncheckedAt(i);
941 if (!kalman1) continue;
942 if (!kalman0) {arr0->AddAt(new AliRelAlignerKalman(*kalman1),i); continue;}
943 kalman0->SetRejectOutliers(kFALSE);
944 kalman0->Merge(kalman1);
945 }
946 }
947
c74c5f6c 948 }
c74c5f6c 949 return 0;
c74c5f6c 950}
951
c74c5f6c 952Bool_t AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
c74c5f6c 953 /*
954 // 0. Same direction - OPOSITE - cutDir +cutT
955 TCut cutDir("cutDir","dir<-0.99")
956 // 1.
957 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
958 //
959 // 2. The same rphi
960 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
961 //
c74c5f6c 962 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
963 // 1/Pt diff cut
964 */
965 const Double_t *p0 = tr0->GetParameter();
966 const Double_t *p1 = tr1->GetParameter();
74235403 967 fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
968 fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
d3ce44cb 969 fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
74235403 970 fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
971 fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
972
c74c5f6c 973 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
974 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
da6c0bc9 975 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
c74c5f6c 976 Double_t d0[3], d1[3];
977 tr0->GetDirection(d0);
978 tr1->GetDirection(d1);
979 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
74235403 980
d3ce44cb 981 fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
982 fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
983 fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
984 fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
985 fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
986
c74c5f6c 987 return kTRUE;
988}
2811495d 989Bool_t AliTPCcalibTime::IsCross(AliESDtrack *const tr0, AliESDtrack *const tr1){
cc65e4f5 990 //
991 // check if the cosmic pair of tracks crossed A/C side
992 //
993 Bool_t result= tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0;
994 if (result==kFALSE) return result;
995 result=kTRUE;
996 return result;
d3ce44cb 997}
998
2811495d 999Bool_t AliTPCcalibTime::IsSame(AliESDtrack *const tr0, AliESDtrack *const tr1){
3e55050f 1000 //
1001 // track crossing the CE
1002 // 0. minimal number of clusters
1003 // 1. Same sector +-1
1004 // 2. Inner and outer track param on opposite side
1005 // 3. Outer and inner track parameter close each to other
1006 // 3.
1007 Bool_t result=kTRUE;
1008 //
1009 // inner and outer on opposite sides in z
1010 //
1011 const Int_t knclCut0 = 30;
1012 const Double_t kalphaCut = 0.4;
1013 //
1014 // 0. minimal number of clusters
1015 //
1016 if (tr0->GetTPCNcls()<knclCut0) return kFALSE;
1017 if (tr1->GetTPCNcls()<knclCut0) return kFALSE;
1018 //
1019 // 1. alpha cut - sector+-1
1020 //
1021 if (TMath::Abs(tr0->GetOuterParam()->GetAlpha()-tr1->GetOuterParam()->GetAlpha())>kalphaCut) return kFALSE;
1022 //
1023 // 2. Z crossing
1024 //
1025 if (tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0) result&=kFALSE;
1026 if (tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0) result&=kFALSE;
1027 if (result==kFALSE){
1028 return result;
1029 }
1030 //
1031 //
1032 const Double_t *p0I = tr0->GetInnerParam()->GetParameter();
1033 const Double_t *p1I = tr1->GetInnerParam()->GetParameter();
1034 const Double_t *p0O = tr0->GetOuterParam()->GetParameter();
1035 const Double_t *p1O = tr1->GetOuterParam()->GetParameter();
1036 //
1037 if (TMath::Abs(p0I[0]-p1I[0])>fCutMaxD) result&=kFALSE;
1038 if (TMath::Abs(p0I[1]-p1I[1])>fCutMaxDz) result&=kFALSE;
1039 if (TMath::Abs(p0I[2]-p1I[2])>fCutTheta) result&=kFALSE;
1040 if (TMath::Abs(p0I[3]-p1I[3])>fCutTheta) result&=kFALSE;
1041 if (TMath::Abs(p0O[0]-p1O[0])>fCutMaxD) result&=kFALSE;
1042 if (TMath::Abs(p0O[1]-p1O[1])>fCutMaxDz) result&=kFALSE;
1043 if (TMath::Abs(p0O[2]-p1O[2])>fCutTheta) result&=kFALSE;
1044 if (TMath::Abs(p0O[3]-p1O[3])>fCutTheta) result&=kFALSE;
1045 if (result==kTRUE){
1046 result=kTRUE; // just to put break point here
1047 }
1048 return result;
dde68d36 1049}
1050
31aa7c5c 1051
2811495d 1052void AliTPCcalibTime::ProcessSame(AliESDtrack *const track, AliESDfriendTrack *const friendTrack, const AliESDEvent *const event){
3e55050f 1053 //
1054 // Process TPC tracks crossing CE
1055 //
1056 // 0. Select only track crossing the CE
1057 // 1. Cut on the track length
1058 // 2. Refit the terack on A and C side separatelly
1059 // 3. Fill time histograms
1060 const Int_t kMinNcl=100;
1061 const Int_t kMinNclS=25; // minimul number of clusters on the sides
1062 if (!friendTrack->GetTPCOut()) return;
1063 //
1064 // 0. Select only track crossing the CE
1065 //
1066 if (track->GetInnerParam()->GetZ()*friendTrack->GetTPCOut()->GetZ()>0) return;
1067 //
1068 // 1. cut on track length
1069 //
1070 if (track->GetTPCNcls()<kMinNcl) return;
1071 //
1072 // 2. Refit track sepparatel on A and C side
1073 //
1074 TObject *calibObject;
1075 AliTPCseed *seed = 0;
1076 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
1077 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
1078 }
1079 if (!seed) return;
1080 //
1081 AliExternalTrackParam trackIn(*track->GetInnerParam());
1082 AliExternalTrackParam trackOut(*track->GetOuterParam());
1083 Double_t cov[3]={0.01,0.,0.01}; //use the same errors
1084 Double_t xyz[3]={0,0.,0.0};
1085 Double_t bz =0;
1086 Int_t nclIn=0,nclOut=0;
1087 trackIn.ResetCovariance(30.);
1088 trackOut.ResetCovariance(30.);
1089 //
1090 //2.a Refit inner
1091 //
1092 for (Int_t irow=0;irow<159;irow++) {
1093 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
1094 if (!cl) continue;
1095 if (cl->GetX()<80) continue;
1096 if (track->GetInnerParam()->GetZ()<0 &&(cl->GetDetector()%36)<18) break;
1097 if (track->GetInnerParam()->GetZ()>0 &&(cl->GetDetector()%36)>=18) break;
1098 Int_t sector = cl->GetDetector();
1099 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
1100 if (TMath::Abs(dalpha)>0.01){
1101 if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
1102 }
1103 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
1104 trackIn.GetXYZ(xyz);
1105 bz = AliTracker::GetBz(xyz);
1106 if (!trackIn.PropagateTo(r[0],bz)) break;
1107 nclIn++;
1108 trackIn.Update(&r[1],cov);
1109 }
1110 //
1111 //2.b Refit outer
1112 //
1113 for (Int_t irow=159;irow>0;irow--) {
1114 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
1115 if (!cl) continue;
1116 if (cl->GetX()<80) continue;
1117 if (cl->GetZ()*track->GetOuterParam()->GetZ()<0) break;
1118 if (friendTrack->GetTPCOut()->GetZ()<0 &&(cl->GetDetector()%36)<18) break;
1119 if (friendTrack->GetTPCOut()->GetZ()>0 &&(cl->GetDetector()%36)>=18) break;
1120 Int_t sector = cl->GetDetector();
1121 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
1122 if (TMath::Abs(dalpha)>0.01){
1123 if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
1124 }
1125 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
1126 trackOut.GetXYZ(xyz);
1127 bz = AliTracker::GetBz(xyz);
1128 if (!trackOut.PropagateTo(r[0],bz)) break;
1129 nclOut++;
1130 trackOut.Update(&r[1],cov);
1131 }
1132 trackOut.Rotate(trackIn.GetAlpha());
1133 Double_t meanX = (trackIn.GetX()+trackOut.GetX())*0.5;
1134 trackIn.PropagateTo(meanX,bz);
1135 trackOut.PropagateTo(meanX,bz);
1136 TTreeSRedirector *cstream = GetDebugStreamer();
1137 if (cstream){
1138 TVectorD gxyz(3);
1139 trackIn.GetXYZ(gxyz.GetMatrixArray());
1140 TTimeStamp tstamp(fTime);
1141 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1142 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
32100b1d 1143 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1144 (*cstream)<<"tpctpc"<<
1145 "run="<<fRun<< // run number
1146 "event="<<fEvent<< // event number
1147 "time="<<fTime<< // time stamp of event
1148 "trigger="<<fTrigger<< // trigger
1149 "mag="<<fMagF<< // magnetic field
1150 "ptrel0.="<<ptrelative0<<
1151 "ptrel1.="<<ptrelative1<<
32100b1d 1152 "vdcorr="<<vdcorr<< // drift correction applied
3e55050f 1153 //
1154 "xyz.="<<&gxyz<< // global position
1155 "tIn.="<<&trackIn<< // refitterd track in
1156 "tOut.="<<&trackOut<< // refitter track out
1157 "nclIn="<<nclIn<< //
1158 "nclOut="<<nclOut<< //
1159 "\n";
1160 }
1161 //
1162 // 3. Fill time histograms
1163 // Debug stremaer expression
1164 // chainTPCTPC->Draw("(tIn.fP[1]-tOut.fP[1])*sign(-tIn.fP[3]):tIn.fP[3]","min(nclIn,nclOut)>30","")
1165 if (TMath::Min(nclIn,nclOut)>kMinNclS){
1166 fDz = trackOut.GetZ()-trackIn.GetZ();
1167 if (trackOut.GetTgl()<0) fDz*=-1.;
1168 TTimeStamp tstamp(fTime);
1169 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1170 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1171 Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
1172 //
1173 // fill histograms per trigger class and itegrated
1174 //
1175 THnSparse* curHist=NULL;
1176 for (Int_t itype=0; itype<2; itype++){
1177 TString name="MEAN_VDRIFT_CROSS_";
1178 if (itype==0){
1179 name+=event->GetFiredTriggerClasses();
1180 name.ToUpper();
1181 }else{
1182 name+="ALL";
1183 }
1184 curHist=(THnSparseF*)fArrayDz->FindObject(name);
1185 if(!curHist){
1186 curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
1187 fArrayDz->AddLast(curHist);
1188 }
1189 curHist->Fill(vecDrift);
1190 }
1191 }
1192
1193}
817766d5 1194
2811495d 1195void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, AliESDfriendTrack *const friendTrack, const AliESDEvent *const event, AliESDfriend *const esdFriend){
817766d5 1196 //
3e55050f 1197 // Process track - Update TPC-ITS alignment
1198 // Updates:
1199 // 0. Apply standartd cuts
1200 // 1. Recalucluate the current statistic median/RMS
1201 // 2. Apply median+-rms cut
1202 // 3. Update kalman filter
1203 //
1204 const Int_t kMinTPC = 80; // minimal number of TPC cluster
1205 const Int_t kMinITS = 3; // minimal number of ITS cluster
1206 const Double_t kMinZ = 10; // maximal dz distance
b96c3aef 1207 const Double_t kMaxDy = 2.; // maximal dy distance
1208 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1209 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1210 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1211 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1212 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
b96c3aef 1213 const Double_t kMinPt = 0.3; // minimal pt
3e55050f 1214 const Int_t kN=500; // deepnes of history
1215 static Int_t kglast=0;
1216 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
1217 /*
1218 0. Standrd cuts:
1219 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";
1220 */
817766d5 1221 //
3e55050f 1222 // 0. Apply standard cuts
817766d5 1223 //
1224 Int_t dummycl[1000];
817766d5 1225 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
b842d904 1226 if (track->GetITSclusters(dummycl)<kMinITS) return; // minimal amount of clusters
b96c3aef 1227 if (!track->IsOn(AliESDtrack::kTPCrefit)) return;
3e55050f 1228 if (!friendTrack->GetITSOut()) return;
817766d5 1229 if (!track->GetInnerParam()) return;
1230 if (!track->GetOuterParam()) return;
b842d904 1231 if (track->GetInnerParam()->Pt()<kMinPt) return;
817766d5 1232 // exclude crossing track
1233 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1234 if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ) return;
b96c3aef 1235 if (track->GetInnerParam()->GetX()>90) return;
817766d5 1236 //
1237 AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetInnerParam()));
b842d904 1238 AliExternalTrackParam pITS(*(friendTrack->GetITSOut())); // ITS standalone if possible
1239 AliExternalTrackParam pITS2(*(friendTrack->GetITSOut())); //TPC-ITS track
1240 pITS2.Rotate(pTPC.GetAlpha());
1241 pITS2.PropagateTo(pTPC.GetX(),fMagF);
1242 AliESDfriendTrack *itsfriendTrack=0;
1243 //
1244 // try to find standalone ITS track corresponing to the TPC if possible
1245 //
1246 Bool_t hasAlone=kFALSE;
1247 Int_t ntracks=event->GetNumberOfTracks();
1248 for (Int_t i=0; i<ntracks; i++){
1249 AliESDtrack *trackS = event->GetTrack(i);
1250 if (trackS->GetTPCNcls()>0) continue; //continue if has TPC info
6ba68367 1251 itsfriendTrack = esdFriend->GetTrack(i);
b842d904 1252 if (!itsfriendTrack) continue;
1253 if (!itsfriendTrack->GetITSOut()) continue;
1254 if (TMath::Abs(pITS2.GetTgl()-itsfriendTrack->GetITSOut()->GetTgl())> kMaxAngle) continue;
1255 pITS=(*(itsfriendTrack->GetITSOut()));
1256 //
1257 pITS.Rotate(pTPC.GetAlpha());
1258 pITS.PropagateTo(pTPC.GetX(),fMagF);
1259 if (TMath::Abs(pITS2.GetY()-pITS.GetY())> kMaxDy) continue;
1260 hasAlone=kTRUE;
1261 }
1262 if (!hasAlone) pITS=pITS2;
1263 //
3e55050f 1264 if (TMath::Abs(pITS.GetY()-pTPC.GetY()) >kMaxDy) return;
1265 if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1266 if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1267 //
1268 // 1. Update median and RMS info
1269 //
1270 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1271 TVectorD vecDeltaN(5);
1272 Double_t sign=(pITS.GetParameter()[1]>0)? 1.:-1.;
1273 vecDelta[4]=0;
1274 for (Int_t i=0;i<4;i++){
1275 vecDelta[i]=(pITS.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1276 kgdP[i][kglast%kN]=vecDelta[i];
1277 }
1278 kglast=(kglast+1);
1279 Int_t entries=(kglast<kN)?kglast:kN;
1280 for (Int_t i=0;i<4;i++){
1281 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1282 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1283 vecDeltaN[i] = 0;
1284 if (vecRMS[i]>0.){
1285 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1286 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1287 }
1288 }
817766d5 1289 //
3e55050f 1290 // 2. Apply median+-rms cut
817766d5 1291 //
3e55050f 1292 if (kglast<3) return; //median and RMS to be defined
1293 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1294 //
1295 // 3. Update alignment
817766d5 1296 //
1297 Int_t htime = fTime/3600; //time in hours
1298 if (fAlignITSTPC->GetEntries()<htime){
1299 fAlignITSTPC->Expand(htime*2+20);
1300 }
1301 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignITSTPC->At(htime);
1302 if (!align){
3e55050f 1303 // make Alignment object if doesn't exist
817766d5 1304 align=new AliRelAlignerKalman();
3e55050f 1305 align->SetRunNumber(fRun);
1306 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1307 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1308 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1309 align->SetRejectOutliers(kFALSE);
1310
1311 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
817766d5 1312 align->SetMagField(fMagF);
1313 fAlignITSTPC->AddAt(align,htime);
1314 }
817766d5 1315 align->AddTrackParams(&pITS,&pTPC);
1316 align->SetTimeStamp(fTime);
3e55050f 1317 align->SetRunNumber(fRun );
391ffdb2 1318 Float_t dca[2],cov[3];
1319 track->GetImpactParameters(dca,cov);
1320 if (TMath::Abs(dca[0])<kMaxDy&&TMath::Abs(dca[0])<kMaxDy){
1321 FillResHistoTPCITS(&pTPC,&pITS);
1322 FillResHistoTPC(track);
1323 }
3e55050f 1324 //
1325 Int_t nupdates=align->GetNUpdates();
1326 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1327 align->SetRejectOutliers(kFALSE);
817766d5 1328 TTreeSRedirector *cstream = GetDebugStreamer();
1329 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1330 TTimeStamp tstamp(fTime);
1331 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1332 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1333 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1334 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1335 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1336 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1337 TVectorD vecGoofie(20);
1338 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1339 if (goofieArray){
1340 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1341 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1342 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1343 }
1344 }
1345 TVectorD gpTPC(3), gdTPC(3);
1346 TVectorD gpITS(3), gdITS(3);
1347 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1348 pTPC.GetDirection(gdTPC.GetMatrixArray());
1349 pITS.GetXYZ(gpITS.GetMatrixArray());
1350 pITS.GetDirection(gdITS.GetMatrixArray());
32100b1d 1351 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
817766d5 1352 (*cstream)<<"itstpc"<<
1353 "run="<<fRun<< // run number
1354 "event="<<fEvent<< // event number
1355 "time="<<fTime<< // time stamp of event
1356 "trigger="<<fTrigger<< // trigger
1357 "mag="<<fMagF<< // magnetic field
1358 // Environment values
1359 "press0="<<valuePressure0<<
1360 "press1="<<valuePressure1<<
1361 "pt0="<<ptrelative0<<
1362 "pt1="<<ptrelative1<<
1363 "temp0="<<temp0<<
1364 "temp1="<<temp1<<
1365 "vecGoofie.="<<&vecGoofie<<
32100b1d 1366 "vdcorr="<<vdcorr<< // drift correction applied
817766d5 1367 //
b842d904 1368 "hasAlone="<<hasAlone<< // has ITS standalone ?
1369 "track.="<<track<< // track info
3e55050f 1370 "nmed="<<kglast<< // number of entries to define median and RMS
1371 "vMed.="<<&vecMedian<< // median of deltas
1372 "vRMS.="<<&vecRMS<< // rms of deltas
1373 "vDelta.="<<&vecDelta<< // delta in respect to median
1374 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1375 "t.="<<track<< // ful track - find proper cuts
1376 "a.="<<align<< // current alignment
b842d904 1377 "pITS.="<<&pITS<< // track param ITS
1378 "pITS2.="<<&pITS2<< // track param ITS+TPC
3e55050f 1379 "pTPC.="<<&pTPC<< // track param TPC
1380 "gpTPC.="<<&gpTPC<< // global position TPC
1381 "gdTPC.="<<&gdTPC<< // global direction TPC
1382 "gpITS.="<<&gpITS<< // global position ITS
1383 "gdITS.="<<&gdITS<< // global position ITS
817766d5 1384 "\n";
1385 }
817766d5 1386}
3e55050f 1387
1388
1389
1390
2811495d 1391void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack *const track, AliESDfriendTrack *const friendTrack){
0dac7b3a 1392 //
3e55050f 1393 // Process track - Update TPC-TRD alignment
1394 // Updates:
1395 // 0. Apply standartd cuts
1396 // 1. Recalucluate the current statistic median/RMS
1397 // 2. Apply median+-rms cut
1398 // 3. Update kalman filter
1399 //
1400 const Int_t kMinTPC = 80; // minimal number of TPC cluster
1401 const Int_t kMinTRD = 50; // minimal number of TRD cluster
1402 const Double_t kMinZ = 20; // maximal dz distance
b96c3aef 1403 const Double_t kMaxDy = 2.; // maximal dy distance
1404 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1405 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1406 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1407 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1408 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
1409 const Int_t kN=500; // deepnes of history
1410 static Int_t kglast=0;
1411 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
0dac7b3a 1412 //
3e55050f 1413 // 0. Apply standard cuts
0dac7b3a 1414 //
1415 Int_t dummycl[1000];
1416 if (track->GetTRDclusters(dummycl)<kMinTRD) return; // minimal amount of clusters
1417 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
3e55050f 1418 if (!friendTrack->GetTRDIn()) return;
0dac7b3a 1419 if (!track->GetInnerParam()) return;
1420 if (!track->GetOuterParam()) return;
1421 // exclude crossing track
1422 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
3e55050f 1423 if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ) return;
0dac7b3a 1424 //
1425 AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetOuterParam()));
1426 AliExternalTrackParam pTRD(*(friendTrack->GetTRDIn()));
3e55050f 1427 pTRD.Rotate(pTPC.GetAlpha());
1428 pTRD.PropagateTo(pTPC.GetX(),fMagF);
1429 ((Double_t*)pTRD.GetCovariance())[2]+=3.*3.; // increas sys errors
1430 ((Double_t*)pTRD.GetCovariance())[9]+=0.1*0.1; // increse sys errors
1431
1432 if (TMath::Abs(pTRD.GetY()-pTPC.GetY()) >kMaxDy) return;
1433 if (TMath::Abs(pTRD.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1434 if (TMath::Abs(pTRD.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
0dac7b3a 1435 //
3e55050f 1436 // 1. Update median and RMS info
0dac7b3a 1437 //
3e55050f 1438 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1439 TVectorD vecDeltaN(5);
1440 Double_t sign=(pTRD.GetParameter()[1]>0)? 1.:-1.;
1441 vecDelta[4]=0;
1442 for (Int_t i=0;i<4;i++){
1443 vecDelta[i]=(pTRD.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1444 kgdP[i][kglast%kN]=vecDelta[i];
1445 }
1446 kglast=(kglast+1);
1447 Int_t entries=(kglast<kN)?kglast:kN;
1448 for (Int_t i=0;i<4;i++){
1449 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1450 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1451 vecDeltaN[i] = 0;
1452 if (vecRMS[i]>0.){
1453 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1454 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1455 }
1456 }
1457 //
1458 // 2. Apply median+-rms cut
1459 //
1460 if (kglast<3) return; //median and RMS to be defined
1461 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1462 //
1463 // 3. Update alignment
0dac7b3a 1464 //
1465 Int_t htime = fTime/3600; //time in hours
1466 if (fAlignTRDTPC->GetEntries()<htime){
1467 fAlignTRDTPC->Expand(htime*2+20);
1468 }
1469 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTRDTPC->At(htime);
1470 if (!align){
3e55050f 1471 // make Alignment object if doesn't exist
0dac7b3a 1472 align=new AliRelAlignerKalman();
3e55050f 1473 align->SetRunNumber(fRun);
1474 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1475 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1476 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1477 align->SetRejectOutliers(kFALSE);
1478 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
0dac7b3a 1479 align->SetMagField(fMagF);
1480 fAlignTRDTPC->AddAt(align,htime);
1481 }
0dac7b3a 1482 align->AddTrackParams(&pTRD,&pTPC);
1483 align->SetTimeStamp(fTime);
3e55050f 1484 align->SetRunNumber(fRun );
391ffdb2 1485 FillResHistoTPCITS(&pTPC,&pTRD);
3e55050f 1486 //
1487 Int_t nupdates=align->GetNUpdates();
1488 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1489 align->SetRejectOutliers(kFALSE);
0dac7b3a 1490 TTreeSRedirector *cstream = GetDebugStreamer();
1491 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1492 TTimeStamp tstamp(fTime);
1493 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1494 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1495 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1496 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1497 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1498 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1499 TVectorD vecGoofie(20);
1500 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1501 if (goofieArray){
1502 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1503 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1504 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1505 }
1506 }
1507 TVectorD gpTPC(3), gdTPC(3);
1508 TVectorD gpTRD(3), gdTRD(3);
1509 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1510 pTPC.GetDirection(gdTPC.GetMatrixArray());
1511 pTRD.GetXYZ(gpTRD.GetMatrixArray());
1512 pTRD.GetDirection(gdTRD.GetMatrixArray());
32100b1d 1513 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1514 (*cstream)<<"trdtpc"<<
0dac7b3a 1515 "run="<<fRun<< // run number
1516 "event="<<fEvent<< // event number
1517 "time="<<fTime<< // time stamp of event
1518 "trigger="<<fTrigger<< // trigger
1519 "mag="<<fMagF<< // magnetic field
1520 // Environment values
1521 "press0="<<valuePressure0<<
1522 "press1="<<valuePressure1<<
1523 "pt0="<<ptrelative0<<
1524 "pt1="<<ptrelative1<<
1525 "temp0="<<temp0<<
1526 "temp1="<<temp1<<
1527 "vecGoofie.="<<&vecGoofie<<
32100b1d 1528 "vdcorr="<<vdcorr<< // drift correction applied
0dac7b3a 1529 //
3e55050f 1530 "nmed="<<kglast<< // number of entries to define median and RMS
1531 "vMed.="<<&vecMedian<< // median of deltas
1532 "vRMS.="<<&vecRMS<< // rms of deltas
1533 "vDelta.="<<&vecDelta<< // delta in respect to median
1534 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1535 "t.="<<track<< // ful track - find proper cuts
1536 "a.="<<align<< // current alignment
1537 "pTRD.="<<&pTRD<< // track param TRD
1538 "pTPC.="<<&pTPC<< // track param TPC
1539 "gpTPC.="<<&gpTPC<< // global position TPC
1540 "gdTPC.="<<&gdTPC<< // global direction TPC
1541 "gpTRD.="<<&gpTRD<< // global position TRD
1542 "gdTRD.="<<&gdTRD<< // global position TRD
0dac7b3a 1543 "\n";
1544 }
0dac7b3a 1545}
817766d5 1546
1547
2811495d 1548void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack *const track, AliESDfriendTrack *const friendTrack){
817766d5 1549 //
817766d5 1550 //
3e55050f 1551 // Process track - Update TPC-TOF alignment
1552 // Updates:
1553 // -1. Make a TOF "track"
1554 // 0. Apply standartd cuts
1555 // 1. Recalucluate the current statistic median/RMS
1556 // 2. Apply median+-rms cut
1557 // 3. Update kalman filter
1558 //
1559 const Int_t kMinTPC = 80; // minimal number of TPC cluster
b842d904 1560 // const Double_t kMinZ = 10; // maximal dz distance
3e55050f 1561 const Double_t kMaxDy = 5.; // maximal dy distance
b96c3aef 1562 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1563 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1564 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1565 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1566
1567 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
1568 const Int_t kN=1000; // deepnes of history
1569 static Int_t kglast=0;
1570 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
817766d5 1571 //
3e55050f 1572 // -1. Make a TOF track-
1573 // Clusters are not in friends - use alingment points
1574 //
1575 if (track->GetTOFsignal()<=0) return;
1576 if (!friendTrack->GetTPCOut()) return;
1577 if (!track->GetInnerParam()) return;
1578 if (!track->GetOuterParam()) return;
817766d5 1579 const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
1580 if (!points) return;
3e55050f 1581 AliExternalTrackParam pTPC(*(track->GetOuterParam()));
1582 AliExternalTrackParam pTOF(pTPC);
1583 Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
817766d5 1584 Int_t npoints = points->GetNPoints();
1585 AliTrackPoint point;
3e55050f 1586 Int_t naccept=0;
817766d5 1587 //
1588 for (Int_t ipoint=0;ipoint<npoints;ipoint++){
817766d5 1589 points->GetPoint(point,ipoint);
1590 Float_t xyz[3];
1591 point.GetXYZ(xyz);
1592 Double_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
3e55050f 1593 if (r<350) continue;
817766d5 1594 if (r>400) continue;
3e55050f 1595 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,2.,kTRUE);
1596 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,0.1,kTRUE);
1597 AliTrackPoint lpoint = point.Rotate(pTPC.GetAlpha());
1598 pTPC.PropagateTo(lpoint.GetX(),fMagF);
1599 pTOF=pTPC;
817766d5 1600 ((Double_t*)pTOF.GetParameter())[0] =lpoint.GetY();
1601 ((Double_t*)pTOF.GetParameter())[1] =lpoint.GetZ();
3e55050f 1602 ((Double_t*)pTOF.GetCovariance())[0]+=3.*3./12.;
1603 ((Double_t*)pTOF.GetCovariance())[2]+=3.*3./12.;
0dac7b3a 1604 ((Double_t*)pTOF.GetCovariance())[5]+=0.1*0.1;
1605 ((Double_t*)pTOF.GetCovariance())[9]+=0.1*0.1;
3e55050f 1606 naccept++;
1607 }
1608 if (naccept==0) return; // no tof match clusters
1609 //
1610 // 0. Apply standard cuts
1611 //
1612 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
1613 // exclude crossing track
1614 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1615 //
1616 if (TMath::Abs(pTOF.GetY()-pTPC.GetY()) >kMaxDy) return;
1617 if (TMath::Abs(pTOF.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1618 if (TMath::Abs(pTOF.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1619 //
1620 // 1. Update median and RMS info
1621 //
1622 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1623 TVectorD vecDeltaN(5);
1624 Double_t sign=(pTOF.GetParameter()[1]>0)? 1.:-1.;
1625 vecDelta[4]=0;
1626 for (Int_t i=0;i<4;i++){
1627 vecDelta[i]=(pTOF.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1628 kgdP[i][kglast%kN]=vecDelta[i];
1629 }
1630 kglast=(kglast+1);
1631 Int_t entries=(kglast<kN)?kglast:kN;
1632 Bool_t isOK=kTRUE;
1633 for (Int_t i=0;i<4;i++){
1634 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1635 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1636 vecDeltaN[i] = 0;
1637 if (vecRMS[i]>0.){
1638 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/(vecRMS[i]+1.);
1639 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1640 if (TMath::Abs(vecDeltaN[i])>kSigmaCut) isOK=kFALSE;
817766d5 1641 }
3e55050f 1642 }
1643 //
1644 // 2. Apply median+-rms cut
1645 //
1646 if (kglast<10) return; //median and RMS to be defined
1647 if (!isOK) return;
1648 //
1649 // 3. Update alignment
1650 //
1651 Int_t htime = fTime/3600; //time in hours
1652 if (fAlignTOFTPC->GetEntries()<htime){
1653 fAlignTOFTPC->Expand(htime*2+20);
1654 }
1655 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTOFTPC->At(htime);
1656 if (!align){
1657 // make Alignment object if doesn't exist
1658 align=new AliRelAlignerKalman();
1659 align->SetRunNumber(fRun);
1660 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1661 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1662 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1663 align->SetRejectOutliers(kFALSE);
1664 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
1665 align->SetMagField(fMagF);
1666 fAlignTOFTPC->AddAt(align,htime);
1667 }
1668 align->AddTrackParams(&pTOF,&pTPC);
1669 align->SetTimeStamp(fTime);
1670 align->SetRunNumber(fRun );
1671 //
1672 Int_t nupdates=align->GetNUpdates();
1673 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1674 align->SetRejectOutliers(kFALSE);
1675 TTreeSRedirector *cstream = GetDebugStreamer();
1676 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1677 TTimeStamp tstamp(fTime);
1678 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1679 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1680 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1681 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1682 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1683 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1684 TVectorD vecGoofie(20);
1685 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1686 if (goofieArray){
1687 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1688 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1689 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1690 }
817766d5 1691 }
3e55050f 1692 TVectorD gpTPC(3), gdTPC(3);
1693 TVectorD gpTOF(3), gdTOF(3);
1694 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1695 pTPC.GetDirection(gdTPC.GetMatrixArray());
1696 pTOF.GetXYZ(gpTOF.GetMatrixArray());
1697 pTOF.GetDirection(gdTOF.GetMatrixArray());
32100b1d 1698 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1699 (*cstream)<<"toftpc"<<
1700 "run="<<fRun<< // run number
1701 "event="<<fEvent<< // event number
1702 "time="<<fTime<< // time stamp of event
1703 "trigger="<<fTrigger<< // trigger
1704 "mag="<<fMagF<< // magnetic field
1705 // Environment values
1706 "press0="<<valuePressure0<<
1707 "press1="<<valuePressure1<<
1708 "pt0="<<ptrelative0<<
1709 "pt1="<<ptrelative1<<
1710 "temp0="<<temp0<<
1711 "temp1="<<temp1<<
1712 "vecGoofie.="<<&vecGoofie<<
32100b1d 1713 "vdcorr="<<vdcorr<< // drift correction applied
3e55050f 1714 //
1715 "nmed="<<kglast<< // number of entries to define median and RMS
1716 "vMed.="<<&vecMedian<< // median of deltas
1717 "vRMS.="<<&vecRMS<< // rms of deltas
1718 "vDelta.="<<&vecDelta<< // delta in respect to median
1719 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1720 "t.="<<track<< // ful track - find proper cuts
1721 "a.="<<align<< // current alignment
1722 "pTOF.="<<&pTOF<< // track param TOF
1723 "pTPC.="<<&pTPC<< // track param TPC
1724 "gpTPC.="<<&gpTPC<< // global position TPC
1725 "gdTPC.="<<&gdTPC<< // global direction TPC
1726 "gpTOF.="<<&gpTOF<< // global position TOF
1727 "gdTOF.="<<&gdTOF<< // global position TOF
1728 "\n";
817766d5 1729 }
817766d5 1730}
3e55050f 1731
1732
391ffdb2 1733void AliTPCcalibTime::BookDistortionMaps(){
1734 //
1735 // Book ndimensional histograms of distortions/residuals
1736 // Only primary tracks are selected for analysis
1737 //
1738
1739 Double_t xminTrack[4], xmaxTrack[4];
1740 Int_t binsTrack[4];
1741 TString axisName[4];
1742 TString axisTitle[4];
1743 //
1744 binsTrack[0] =100;
1745 axisName[0] ="#Delta";
1746 axisTitle[0] ="#Delta";
1747 //
1748 binsTrack[1] =20;
1749 xminTrack[1] =-1.5; xmaxTrack[1]=1.5;
1750 axisName[1] ="tanTheta";
1751 axisTitle[1] ="tan(#Theta)";
1752 //
1753 binsTrack[2] =90;
1754 xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi();
1755 axisName[2] ="phi";
1756 axisTitle[2] ="#phi";
1757 //
1758 binsTrack[3] =20;
1759 xminTrack[3] =-1.; xmaxTrack[3]=1.; // 0.33 GeV cut
1760 axisName[3] ="snp";
1761 //
1762 // delta y
1763 xminTrack[0] =-1.5; xmaxTrack[0]=1.5; //
1764 fResHistoTPCITS[0] = new THnSparseS("TPCITS#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1765 fResHistoTPCvertex[0] = new THnSparseS("TPCVertex#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1766 fResHistoTPCTRD[0] = new THnSparseS("TPCTRD#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1767 //
1768 // delta z
1769 xminTrack[0] =-1.5; xmaxTrack[0]=1.5; //
1770 fResHistoTPCITS[1] = new THnSparseS("TPCITS#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1771 fResHistoTPCvertex[1] = new THnSparseS("TPCVertex#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1772 fResHistoTPCTRD[1] = new THnSparseS("TPCTRD#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1773 //
1774 // delta snp-P2
1775 xminTrack[0] =-0.01; xmaxTrack[0]=0.01; //
1776 fResHistoTPCITS[2] = new THnSparseS("TPCITS#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
1777 fResHistoTPCvertex[2] = new THnSparseS("TPCITS#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
1778 fResHistoTPCTRD[2] = new THnSparseS("TPCTRD#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
1779 //
1780 // delta theta-P3
1781 xminTrack[0] =-0.01; xmaxTrack[0]=0.01; //
1782 fResHistoTPCITS[3] = new THnSparseS("TPCITS#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
1783 fResHistoTPCvertex[3] = new THnSparseS("TPCITS#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
1784 fResHistoTPCTRD[3] = new THnSparseS("TPCTRD#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
1785 //
1786 // delta theta-P4
1787 xminTrack[0] =-0.05; xmaxTrack[0]=0.05; //
1788 fResHistoTPCITS[4] = new THnSparseS("TPCITS#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
1789 fResHistoTPCvertex[4] = new THnSparseS("TPCITS#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
1790 fResHistoTPCTRD[4] = new THnSparseS("TPCTRD#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
1791 //
1792 for (Int_t ivar=0;ivar<4;ivar++){
1793 for (Int_t ivar2=0;ivar2<4;ivar2++){
1794 fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1795 fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
1796 fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1797 fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
1798 fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1799 fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
1800 }
1801 }
1802}
1803
1804
1805void AliTPCcalibTime::FillResHistoTPCITS(const AliExternalTrackParam * pTPCIn, const AliExternalTrackParam * pITSOut ){
1806 //
1807 // fill residual histograms pTPCIn-pITSOut
1808 // Histogram is filled only for primary tracks
1809 //
1810 Double_t histoX[4];
1811 Double_t xyz[3];
1812 pTPCIn->GetXYZ(xyz);
1813 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1814 histoX[1]= pTPCIn->GetTgl();
1815 histoX[2]= phi;
1816 histoX[3]= pTPCIn->GetSnp();
1817 //
1818 for (Int_t ihisto=0; ihisto<5; ihisto++){
1819 histoX[0]=pTPCIn->GetParameter()[ihisto]-pITSOut->GetParameter()[ihisto];
1820 fResHistoTPCITS[ihisto]->Fill(histoX);
1821 }
1822}
1823
1824
1825void AliTPCcalibTime::FillResHistoTPC(const AliESDtrack * pTrack){
1826 //
1827 // fill residual histograms pTPC - vertex
1828 // Histogram is filled only for primary tracks
1829 //
1830 Double_t histoX[4];
1831 const AliExternalTrackParam * pTPCIn = pTrack->GetInnerParam();
1832 const AliExternalTrackParam * pTPCvertex = pTrack->GetTPCInnerParam();
1833 Double_t xyz[3];
1834 pTPCIn->GetXYZ(xyz);
1835 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1836 histoX[1]= pTPCIn->GetTgl();
1837 histoX[2]= phi;
1838 histoX[3]= pTPCIn->GetSnp();
1839 //
1840 Float_t dca[2], cov[3];
1841 pTrack->GetImpactParametersTPC(dca,cov);
1842 for (Int_t ihisto=0; ihisto<2; ihisto++){
1843 histoX[0]=pTPCvertex->GetParameter()[ihisto]-pTrack->GetParameter()[ihisto];
1844 fResHistoTPCITS[ihisto]->Fill(histoX);
1845 }
1846}
1847
1848
1849void AliTPCcalibTime::FillResHistoTPCTRD(const AliExternalTrackParam * pTPCOut, const AliExternalTrackParam * pTRDIn ){
1850 //
1851 // fill resuidual histogram TPCout-TRDin
1852 //
1853 Double_t histoX[4];
1854 Double_t xyz[3];
1855 pTPCOut->GetXYZ(xyz);
1856 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1857 histoX[1]= pTPCOut->GetTgl();
1858 histoX[2]= phi;
1859 histoX[3]= pTPCOut->GetSnp();
1860 //
1861 for (Int_t ihisto=0; ihisto<5; ihisto++){
1862 histoX[0]=pTPCOut->GetParameter()[ihisto]-pTRDIn->GetParameter()[ihisto];
1863 fResHistoTPCTRD[ihisto]->Fill(histoX);
1864 }
1865
1866}