leave the QA.root file open after the first time it has been opened
[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 }
9112627b 905 //
906 for (Int_t imeas=0; imeas<5; imeas++){
907 if (cal->GetResHistoTPCITS(imeas) && cal->GetResHistoTPCITS(imeas)){
908 fResHistoTPCITS[imeas]->Add(cal->fResHistoTPCITS[imeas]);
909 fResHistoTPCvertex[imeas]->Add(cal->fResHistoTPCvertex[imeas]);
910 fResHistoTPCTRD[imeas]->Add(cal->fResHistoTPCTRD[imeas]);
911 }
912 }
dde68d36 913 TObjArray* addArray=cal->GetHistoDrift();
914 if(!addArray) return 0;
915 TIterator* iterator = addArray->MakeIterator();
2be25cc9 916 iterator->Reset();
dde68d36 917 THnSparse* addHist=NULL;
918 while((addHist=(THnSparseF*)iterator->Next())){
919 if(!addHist) continue;
2be25cc9 920 addHist->Print();
dde68d36 921 THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName());
2be25cc9 922 if(!localHist){
923 localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
dde68d36 924 fArrayDz->AddLast(localHist);
2be25cc9 925 }
926 localHist->Add(addHist);
927 }
391ffdb2 928
d3ce44cb 929 for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
3e55050f 930 //
931 // Merge alignment
932 //
933 for (Int_t itype=0; itype<3; itype++){
934 //
935 //
936 TObjArray *arr0= 0;
937 TObjArray *arr1= 0;
938 if (itype==0) {arr0=fAlignITSTPC; arr1=cal->fAlignITSTPC;}
939 if (itype==1) {arr0=fAlignTRDTPC; arr1=cal->fAlignTRDTPC;}
940 if (itype==2) {arr0=fAlignTOFTPC; arr1=cal->fAlignTOFTPC;}
941 if (!arr1) continue;
942 if (!arr0) arr0=new TObjArray(arr1->GetEntriesFast());
943 if (arr1->GetEntriesFast()>arr0->GetEntriesFast()){
944 arr0->Expand(arr1->GetEntriesFast());
945 }
946 for (Int_t i=0;i<arr1->GetEntriesFast(); i++){
947 AliRelAlignerKalman *kalman1 = (AliRelAlignerKalman *)arr1->UncheckedAt(i);
948 AliRelAlignerKalman *kalman0 = (AliRelAlignerKalman *)arr0->UncheckedAt(i);
949 if (!kalman1) continue;
950 if (!kalman0) {arr0->AddAt(new AliRelAlignerKalman(*kalman1),i); continue;}
951 kalman0->SetRejectOutliers(kFALSE);
952 kalman0->Merge(kalman1);
953 }
954 }
955
c74c5f6c 956 }
c74c5f6c 957 return 0;
c74c5f6c 958}
959
c74c5f6c 960Bool_t AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
c74c5f6c 961 /*
962 // 0. Same direction - OPOSITE - cutDir +cutT
963 TCut cutDir("cutDir","dir<-0.99")
964 // 1.
965 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
966 //
967 // 2. The same rphi
968 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
969 //
c74c5f6c 970 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
971 // 1/Pt diff cut
972 */
973 const Double_t *p0 = tr0->GetParameter();
974 const Double_t *p1 = tr1->GetParameter();
74235403 975 fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
976 fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
d3ce44cb 977 fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
74235403 978 fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
979 fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
980
c74c5f6c 981 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
982 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
da6c0bc9 983 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
c74c5f6c 984 Double_t d0[3], d1[3];
985 tr0->GetDirection(d0);
986 tr1->GetDirection(d1);
987 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
74235403 988
d3ce44cb 989 fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
990 fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
991 fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
992 fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
993 fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
994
c74c5f6c 995 return kTRUE;
996}
2811495d 997Bool_t AliTPCcalibTime::IsCross(AliESDtrack *const tr0, AliESDtrack *const tr1){
cc65e4f5 998 //
999 // check if the cosmic pair of tracks crossed A/C side
1000 //
1001 Bool_t result= tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0;
1002 if (result==kFALSE) return result;
1003 result=kTRUE;
1004 return result;
d3ce44cb 1005}
1006
2811495d 1007Bool_t AliTPCcalibTime::IsSame(AliESDtrack *const tr0, AliESDtrack *const tr1){
3e55050f 1008 //
1009 // track crossing the CE
1010 // 0. minimal number of clusters
1011 // 1. Same sector +-1
1012 // 2. Inner and outer track param on opposite side
1013 // 3. Outer and inner track parameter close each to other
1014 // 3.
1015 Bool_t result=kTRUE;
1016 //
1017 // inner and outer on opposite sides in z
1018 //
1019 const Int_t knclCut0 = 30;
1020 const Double_t kalphaCut = 0.4;
1021 //
1022 // 0. minimal number of clusters
1023 //
1024 if (tr0->GetTPCNcls()<knclCut0) return kFALSE;
1025 if (tr1->GetTPCNcls()<knclCut0) return kFALSE;
1026 //
1027 // 1. alpha cut - sector+-1
1028 //
1029 if (TMath::Abs(tr0->GetOuterParam()->GetAlpha()-tr1->GetOuterParam()->GetAlpha())>kalphaCut) return kFALSE;
1030 //
1031 // 2. Z crossing
1032 //
1033 if (tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0) result&=kFALSE;
1034 if (tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0) result&=kFALSE;
1035 if (result==kFALSE){
1036 return result;
1037 }
1038 //
1039 //
1040 const Double_t *p0I = tr0->GetInnerParam()->GetParameter();
1041 const Double_t *p1I = tr1->GetInnerParam()->GetParameter();
1042 const Double_t *p0O = tr0->GetOuterParam()->GetParameter();
1043 const Double_t *p1O = tr1->GetOuterParam()->GetParameter();
1044 //
1045 if (TMath::Abs(p0I[0]-p1I[0])>fCutMaxD) result&=kFALSE;
1046 if (TMath::Abs(p0I[1]-p1I[1])>fCutMaxDz) result&=kFALSE;
1047 if (TMath::Abs(p0I[2]-p1I[2])>fCutTheta) result&=kFALSE;
1048 if (TMath::Abs(p0I[3]-p1I[3])>fCutTheta) result&=kFALSE;
1049 if (TMath::Abs(p0O[0]-p1O[0])>fCutMaxD) result&=kFALSE;
1050 if (TMath::Abs(p0O[1]-p1O[1])>fCutMaxDz) result&=kFALSE;
1051 if (TMath::Abs(p0O[2]-p1O[2])>fCutTheta) result&=kFALSE;
1052 if (TMath::Abs(p0O[3]-p1O[3])>fCutTheta) result&=kFALSE;
1053 if (result==kTRUE){
1054 result=kTRUE; // just to put break point here
1055 }
1056 return result;
dde68d36 1057}
1058
31aa7c5c 1059
2811495d 1060void AliTPCcalibTime::ProcessSame(AliESDtrack *const track, AliESDfriendTrack *const friendTrack, const AliESDEvent *const event){
3e55050f 1061 //
1062 // Process TPC tracks crossing CE
1063 //
1064 // 0. Select only track crossing the CE
1065 // 1. Cut on the track length
1066 // 2. Refit the terack on A and C side separatelly
1067 // 3. Fill time histograms
1068 const Int_t kMinNcl=100;
1069 const Int_t kMinNclS=25; // minimul number of clusters on the sides
1070 if (!friendTrack->GetTPCOut()) return;
1071 //
1072 // 0. Select only track crossing the CE
1073 //
1074 if (track->GetInnerParam()->GetZ()*friendTrack->GetTPCOut()->GetZ()>0) return;
1075 //
1076 // 1. cut on track length
1077 //
1078 if (track->GetTPCNcls()<kMinNcl) return;
1079 //
1080 // 2. Refit track sepparatel on A and C side
1081 //
1082 TObject *calibObject;
1083 AliTPCseed *seed = 0;
1084 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
1085 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
1086 }
1087 if (!seed) return;
1088 //
1089 AliExternalTrackParam trackIn(*track->GetInnerParam());
1090 AliExternalTrackParam trackOut(*track->GetOuterParam());
1091 Double_t cov[3]={0.01,0.,0.01}; //use the same errors
1092 Double_t xyz[3]={0,0.,0.0};
1093 Double_t bz =0;
1094 Int_t nclIn=0,nclOut=0;
1095 trackIn.ResetCovariance(30.);
1096 trackOut.ResetCovariance(30.);
1097 //
1098 //2.a Refit inner
1099 //
1100 for (Int_t irow=0;irow<159;irow++) {
1101 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
1102 if (!cl) continue;
1103 if (cl->GetX()<80) continue;
1104 if (track->GetInnerParam()->GetZ()<0 &&(cl->GetDetector()%36)<18) break;
1105 if (track->GetInnerParam()->GetZ()>0 &&(cl->GetDetector()%36)>=18) break;
1106 Int_t sector = cl->GetDetector();
1107 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
1108 if (TMath::Abs(dalpha)>0.01){
1109 if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
1110 }
1111 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
1112 trackIn.GetXYZ(xyz);
1113 bz = AliTracker::GetBz(xyz);
1114 if (!trackIn.PropagateTo(r[0],bz)) break;
1115 nclIn++;
1116 trackIn.Update(&r[1],cov);
1117 }
1118 //
1119 //2.b Refit outer
1120 //
1121 for (Int_t irow=159;irow>0;irow--) {
1122 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
1123 if (!cl) continue;
1124 if (cl->GetX()<80) continue;
1125 if (cl->GetZ()*track->GetOuterParam()->GetZ()<0) break;
1126 if (friendTrack->GetTPCOut()->GetZ()<0 &&(cl->GetDetector()%36)<18) break;
1127 if (friendTrack->GetTPCOut()->GetZ()>0 &&(cl->GetDetector()%36)>=18) break;
1128 Int_t sector = cl->GetDetector();
1129 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
1130 if (TMath::Abs(dalpha)>0.01){
1131 if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
1132 }
1133 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
1134 trackOut.GetXYZ(xyz);
1135 bz = AliTracker::GetBz(xyz);
1136 if (!trackOut.PropagateTo(r[0],bz)) break;
1137 nclOut++;
1138 trackOut.Update(&r[1],cov);
1139 }
1140 trackOut.Rotate(trackIn.GetAlpha());
1141 Double_t meanX = (trackIn.GetX()+trackOut.GetX())*0.5;
1142 trackIn.PropagateTo(meanX,bz);
1143 trackOut.PropagateTo(meanX,bz);
1144 TTreeSRedirector *cstream = GetDebugStreamer();
1145 if (cstream){
1146 TVectorD gxyz(3);
1147 trackIn.GetXYZ(gxyz.GetMatrixArray());
1148 TTimeStamp tstamp(fTime);
1149 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1150 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
32100b1d 1151 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1152 (*cstream)<<"tpctpc"<<
1153 "run="<<fRun<< // run number
1154 "event="<<fEvent<< // event number
1155 "time="<<fTime<< // time stamp of event
1156 "trigger="<<fTrigger<< // trigger
1157 "mag="<<fMagF<< // magnetic field
1158 "ptrel0.="<<ptrelative0<<
1159 "ptrel1.="<<ptrelative1<<
32100b1d 1160 "vdcorr="<<vdcorr<< // drift correction applied
3e55050f 1161 //
1162 "xyz.="<<&gxyz<< // global position
1163 "tIn.="<<&trackIn<< // refitterd track in
1164 "tOut.="<<&trackOut<< // refitter track out
1165 "nclIn="<<nclIn<< //
1166 "nclOut="<<nclOut<< //
1167 "\n";
1168 }
1169 //
1170 // 3. Fill time histograms
1171 // Debug stremaer expression
1172 // chainTPCTPC->Draw("(tIn.fP[1]-tOut.fP[1])*sign(-tIn.fP[3]):tIn.fP[3]","min(nclIn,nclOut)>30","")
1173 if (TMath::Min(nclIn,nclOut)>kMinNclS){
1174 fDz = trackOut.GetZ()-trackIn.GetZ();
1175 if (trackOut.GetTgl()<0) fDz*=-1.;
1176 TTimeStamp tstamp(fTime);
1177 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1178 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1179 Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
1180 //
1181 // fill histograms per trigger class and itegrated
1182 //
1183 THnSparse* curHist=NULL;
1184 for (Int_t itype=0; itype<2; itype++){
1185 TString name="MEAN_VDRIFT_CROSS_";
1186 if (itype==0){
1187 name+=event->GetFiredTriggerClasses();
1188 name.ToUpper();
1189 }else{
1190 name+="ALL";
1191 }
1192 curHist=(THnSparseF*)fArrayDz->FindObject(name);
1193 if(!curHist){
1194 curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
1195 fArrayDz->AddLast(curHist);
1196 }
1197 curHist->Fill(vecDrift);
1198 }
1199 }
1200
1201}
817766d5 1202
2811495d 1203void AliTPCcalibTime::ProcessAlignITS(AliESDtrack *const track, AliESDfriendTrack *const friendTrack, const AliESDEvent *const event, AliESDfriend *const esdFriend){
817766d5 1204 //
3e55050f 1205 // Process track - Update TPC-ITS alignment
1206 // Updates:
1207 // 0. Apply standartd cuts
1208 // 1. Recalucluate the current statistic median/RMS
1209 // 2. Apply median+-rms cut
1210 // 3. Update kalman filter
1211 //
1212 const Int_t kMinTPC = 80; // minimal number of TPC cluster
1213 const Int_t kMinITS = 3; // minimal number of ITS cluster
1214 const Double_t kMinZ = 10; // maximal dz distance
b96c3aef 1215 const Double_t kMaxDy = 2.; // maximal dy distance
1216 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1217 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1218 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1219 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1220 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
b96c3aef 1221 const Double_t kMinPt = 0.3; // minimal pt
3e55050f 1222 const Int_t kN=500; // deepnes of history
1223 static Int_t kglast=0;
1224 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
1225 /*
1226 0. Standrd cuts:
1227 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";
1228 */
817766d5 1229 //
3e55050f 1230 // 0. Apply standard cuts
817766d5 1231 //
1232 Int_t dummycl[1000];
817766d5 1233 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
b842d904 1234 if (track->GetITSclusters(dummycl)<kMinITS) return; // minimal amount of clusters
b96c3aef 1235 if (!track->IsOn(AliESDtrack::kTPCrefit)) return;
3e55050f 1236 if (!friendTrack->GetITSOut()) return;
817766d5 1237 if (!track->GetInnerParam()) return;
1238 if (!track->GetOuterParam()) return;
b842d904 1239 if (track->GetInnerParam()->Pt()<kMinPt) return;
817766d5 1240 // exclude crossing track
1241 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1242 if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ) return;
b96c3aef 1243 if (track->GetInnerParam()->GetX()>90) return;
817766d5 1244 //
1245 AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetInnerParam()));
b842d904 1246 AliExternalTrackParam pITS(*(friendTrack->GetITSOut())); // ITS standalone if possible
1247 AliExternalTrackParam pITS2(*(friendTrack->GetITSOut())); //TPC-ITS track
1248 pITS2.Rotate(pTPC.GetAlpha());
1249 pITS2.PropagateTo(pTPC.GetX(),fMagF);
1250 AliESDfriendTrack *itsfriendTrack=0;
1251 //
1252 // try to find standalone ITS track corresponing to the TPC if possible
1253 //
1254 Bool_t hasAlone=kFALSE;
1255 Int_t ntracks=event->GetNumberOfTracks();
1256 for (Int_t i=0; i<ntracks; i++){
1257 AliESDtrack *trackS = event->GetTrack(i);
1258 if (trackS->GetTPCNcls()>0) continue; //continue if has TPC info
6ba68367 1259 itsfriendTrack = esdFriend->GetTrack(i);
b842d904 1260 if (!itsfriendTrack) continue;
1261 if (!itsfriendTrack->GetITSOut()) continue;
1262 if (TMath::Abs(pITS2.GetTgl()-itsfriendTrack->GetITSOut()->GetTgl())> kMaxAngle) continue;
1263 pITS=(*(itsfriendTrack->GetITSOut()));
1264 //
1265 pITS.Rotate(pTPC.GetAlpha());
1266 pITS.PropagateTo(pTPC.GetX(),fMagF);
1267 if (TMath::Abs(pITS2.GetY()-pITS.GetY())> kMaxDy) continue;
1268 hasAlone=kTRUE;
1269 }
1270 if (!hasAlone) pITS=pITS2;
1271 //
3e55050f 1272 if (TMath::Abs(pITS.GetY()-pTPC.GetY()) >kMaxDy) return;
1273 if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1274 if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1275 //
1276 // 1. Update median and RMS info
1277 //
1278 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1279 TVectorD vecDeltaN(5);
1280 Double_t sign=(pITS.GetParameter()[1]>0)? 1.:-1.;
1281 vecDelta[4]=0;
1282 for (Int_t i=0;i<4;i++){
1283 vecDelta[i]=(pITS.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1284 kgdP[i][kglast%kN]=vecDelta[i];
1285 }
1286 kglast=(kglast+1);
1287 Int_t entries=(kglast<kN)?kglast:kN;
1288 for (Int_t i=0;i<4;i++){
1289 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1290 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1291 vecDeltaN[i] = 0;
1292 if (vecRMS[i]>0.){
1293 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1294 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1295 }
1296 }
817766d5 1297 //
3e55050f 1298 // 2. Apply median+-rms cut
817766d5 1299 //
3e55050f 1300 if (kglast<3) return; //median and RMS to be defined
1301 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1302 //
1303 // 3. Update alignment
817766d5 1304 //
1305 Int_t htime = fTime/3600; //time in hours
1306 if (fAlignITSTPC->GetEntries()<htime){
1307 fAlignITSTPC->Expand(htime*2+20);
1308 }
1309 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignITSTPC->At(htime);
1310 if (!align){
3e55050f 1311 // make Alignment object if doesn't exist
817766d5 1312 align=new AliRelAlignerKalman();
3e55050f 1313 align->SetRunNumber(fRun);
1314 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1315 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1316 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1317 align->SetRejectOutliers(kFALSE);
1318
1319 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
817766d5 1320 align->SetMagField(fMagF);
1321 fAlignITSTPC->AddAt(align,htime);
1322 }
817766d5 1323 align->AddTrackParams(&pITS,&pTPC);
1324 align->SetTimeStamp(fTime);
3e55050f 1325 align->SetRunNumber(fRun );
391ffdb2 1326 Float_t dca[2],cov[3];
1327 track->GetImpactParameters(dca,cov);
1328 if (TMath::Abs(dca[0])<kMaxDy&&TMath::Abs(dca[0])<kMaxDy){
1329 FillResHistoTPCITS(&pTPC,&pITS);
1330 FillResHistoTPC(track);
1331 }
3e55050f 1332 //
1333 Int_t nupdates=align->GetNUpdates();
1334 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1335 align->SetRejectOutliers(kFALSE);
817766d5 1336 TTreeSRedirector *cstream = GetDebugStreamer();
1337 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1338 TTimeStamp tstamp(fTime);
1339 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1340 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1341 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1342 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1343 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1344 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1345 TVectorD vecGoofie(20);
1346 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1347 if (goofieArray){
1348 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1349 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1350 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1351 }
1352 }
1353 TVectorD gpTPC(3), gdTPC(3);
1354 TVectorD gpITS(3), gdITS(3);
1355 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1356 pTPC.GetDirection(gdTPC.GetMatrixArray());
1357 pITS.GetXYZ(gpITS.GetMatrixArray());
1358 pITS.GetDirection(gdITS.GetMatrixArray());
32100b1d 1359 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
817766d5 1360 (*cstream)<<"itstpc"<<
1361 "run="<<fRun<< // run number
1362 "event="<<fEvent<< // event number
1363 "time="<<fTime<< // time stamp of event
1364 "trigger="<<fTrigger<< // trigger
1365 "mag="<<fMagF<< // magnetic field
1366 // Environment values
1367 "press0="<<valuePressure0<<
1368 "press1="<<valuePressure1<<
1369 "pt0="<<ptrelative0<<
1370 "pt1="<<ptrelative1<<
1371 "temp0="<<temp0<<
1372 "temp1="<<temp1<<
1373 "vecGoofie.="<<&vecGoofie<<
32100b1d 1374 "vdcorr="<<vdcorr<< // drift correction applied
817766d5 1375 //
b842d904 1376 "hasAlone="<<hasAlone<< // has ITS standalone ?
1377 "track.="<<track<< // track info
3e55050f 1378 "nmed="<<kglast<< // number of entries to define median and RMS
1379 "vMed.="<<&vecMedian<< // median of deltas
1380 "vRMS.="<<&vecRMS<< // rms of deltas
1381 "vDelta.="<<&vecDelta<< // delta in respect to median
1382 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1383 "t.="<<track<< // ful track - find proper cuts
1384 "a.="<<align<< // current alignment
b842d904 1385 "pITS.="<<&pITS<< // track param ITS
1386 "pITS2.="<<&pITS2<< // track param ITS+TPC
3e55050f 1387 "pTPC.="<<&pTPC<< // track param TPC
1388 "gpTPC.="<<&gpTPC<< // global position TPC
1389 "gdTPC.="<<&gdTPC<< // global direction TPC
1390 "gpITS.="<<&gpITS<< // global position ITS
1391 "gdITS.="<<&gdITS<< // global position ITS
817766d5 1392 "\n";
1393 }
817766d5 1394}
3e55050f 1395
1396
1397
1398
2811495d 1399void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack *const track, AliESDfriendTrack *const friendTrack){
0dac7b3a 1400 //
3e55050f 1401 // Process track - Update TPC-TRD alignment
1402 // Updates:
1403 // 0. Apply standartd cuts
1404 // 1. Recalucluate the current statistic median/RMS
1405 // 2. Apply median+-rms cut
1406 // 3. Update kalman filter
1407 //
1408 const Int_t kMinTPC = 80; // minimal number of TPC cluster
1409 const Int_t kMinTRD = 50; // minimal number of TRD cluster
1410 const Double_t kMinZ = 20; // maximal dz distance
9112627b 1411 const Double_t kMaxDy = 5.; // maximal dy distance
1412 const Double_t kMaxAngle= 0.1; // maximal angular distance
1413 const Double_t kSigmaCut= 10; // maximal sigma distance to median
3e55050f 1414 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1415 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1416 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
1417 const Int_t kN=500; // deepnes of history
1418 static Int_t kglast=0;
1419 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
0dac7b3a 1420 //
3e55050f 1421 // 0. Apply standard cuts
0dac7b3a 1422 //
1423 Int_t dummycl[1000];
1424 if (track->GetTRDclusters(dummycl)<kMinTRD) return; // minimal amount of clusters
1425 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
3e55050f 1426 if (!friendTrack->GetTRDIn()) return;
0dac7b3a 1427 if (!track->GetInnerParam()) return;
1428 if (!track->GetOuterParam()) return;
1429 // exclude crossing track
1430 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
3e55050f 1431 if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ) return;
0dac7b3a 1432 //
1433 AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetOuterParam()));
1434 AliExternalTrackParam pTRD(*(friendTrack->GetTRDIn()));
3e55050f 1435 pTRD.Rotate(pTPC.GetAlpha());
1436 pTRD.PropagateTo(pTPC.GetX(),fMagF);
1437 ((Double_t*)pTRD.GetCovariance())[2]+=3.*3.; // increas sys errors
1438 ((Double_t*)pTRD.GetCovariance())[9]+=0.1*0.1; // increse sys errors
1439
1440 if (TMath::Abs(pTRD.GetY()-pTPC.GetY()) >kMaxDy) return;
1441 if (TMath::Abs(pTRD.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
9112627b 1442 // if (TMath::Abs(pTRD.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
0dac7b3a 1443 //
3e55050f 1444 // 1. Update median and RMS info
0dac7b3a 1445 //
3e55050f 1446 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1447 TVectorD vecDeltaN(5);
1448 Double_t sign=(pTRD.GetParameter()[1]>0)? 1.:-1.;
1449 vecDelta[4]=0;
1450 for (Int_t i=0;i<4;i++){
1451 vecDelta[i]=(pTRD.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1452 kgdP[i][kglast%kN]=vecDelta[i];
1453 }
1454 kglast=(kglast+1);
1455 Int_t entries=(kglast<kN)?kglast:kN;
1456 for (Int_t i=0;i<4;i++){
1457 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1458 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1459 vecDeltaN[i] = 0;
1460 if (vecRMS[i]>0.){
1461 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1462 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1463 }
1464 }
1465 //
1466 // 2. Apply median+-rms cut
1467 //
1468 if (kglast<3) return; //median and RMS to be defined
1469 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1470 //
1471 // 3. Update alignment
0dac7b3a 1472 //
1473 Int_t htime = fTime/3600; //time in hours
1474 if (fAlignTRDTPC->GetEntries()<htime){
1475 fAlignTRDTPC->Expand(htime*2+20);
1476 }
1477 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTRDTPC->At(htime);
1478 if (!align){
3e55050f 1479 // make Alignment object if doesn't exist
0dac7b3a 1480 align=new AliRelAlignerKalman();
3e55050f 1481 align->SetRunNumber(fRun);
1482 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1483 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1484 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1485 align->SetRejectOutliers(kFALSE);
1486 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
0dac7b3a 1487 align->SetMagField(fMagF);
1488 fAlignTRDTPC->AddAt(align,htime);
1489 }
0dac7b3a 1490 align->AddTrackParams(&pTRD,&pTPC);
1491 align->SetTimeStamp(fTime);
3e55050f 1492 align->SetRunNumber(fRun );
9112627b 1493 FillResHistoTPCTRD(&pTPC,&pTRD);
3e55050f 1494 //
1495 Int_t nupdates=align->GetNUpdates();
1496 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1497 align->SetRejectOutliers(kFALSE);
0dac7b3a 1498 TTreeSRedirector *cstream = GetDebugStreamer();
1499 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1500 TTimeStamp tstamp(fTime);
1501 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1502 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1503 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1504 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1505 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1506 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1507 TVectorD vecGoofie(20);
1508 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1509 if (goofieArray){
1510 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1511 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1512 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1513 }
1514 }
1515 TVectorD gpTPC(3), gdTPC(3);
1516 TVectorD gpTRD(3), gdTRD(3);
1517 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1518 pTPC.GetDirection(gdTPC.GetMatrixArray());
1519 pTRD.GetXYZ(gpTRD.GetMatrixArray());
1520 pTRD.GetDirection(gdTRD.GetMatrixArray());
32100b1d 1521 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1522 (*cstream)<<"trdtpc"<<
0dac7b3a 1523 "run="<<fRun<< // run number
1524 "event="<<fEvent<< // event number
1525 "time="<<fTime<< // time stamp of event
1526 "trigger="<<fTrigger<< // trigger
1527 "mag="<<fMagF<< // magnetic field
1528 // Environment values
1529 "press0="<<valuePressure0<<
1530 "press1="<<valuePressure1<<
1531 "pt0="<<ptrelative0<<
1532 "pt1="<<ptrelative1<<
1533 "temp0="<<temp0<<
1534 "temp1="<<temp1<<
1535 "vecGoofie.="<<&vecGoofie<<
32100b1d 1536 "vdcorr="<<vdcorr<< // drift correction applied
0dac7b3a 1537 //
3e55050f 1538 "nmed="<<kglast<< // number of entries to define median and RMS
1539 "vMed.="<<&vecMedian<< // median of deltas
1540 "vRMS.="<<&vecRMS<< // rms of deltas
1541 "vDelta.="<<&vecDelta<< // delta in respect to median
1542 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1543 "t.="<<track<< // ful track - find proper cuts
1544 "a.="<<align<< // current alignment
1545 "pTRD.="<<&pTRD<< // track param TRD
1546 "pTPC.="<<&pTPC<< // track param TPC
1547 "gpTPC.="<<&gpTPC<< // global position TPC
1548 "gdTPC.="<<&gdTPC<< // global direction TPC
1549 "gpTRD.="<<&gpTRD<< // global position TRD
1550 "gdTRD.="<<&gdTRD<< // global position TRD
0dac7b3a 1551 "\n";
1552 }
0dac7b3a 1553}
817766d5 1554
1555
2811495d 1556void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack *const track, AliESDfriendTrack *const friendTrack){
817766d5 1557 //
817766d5 1558 //
3e55050f 1559 // Process track - Update TPC-TOF alignment
1560 // Updates:
1561 // -1. Make a TOF "track"
1562 // 0. Apply standartd cuts
1563 // 1. Recalucluate the current statistic median/RMS
1564 // 2. Apply median+-rms cut
1565 // 3. Update kalman filter
1566 //
1567 const Int_t kMinTPC = 80; // minimal number of TPC cluster
b842d904 1568 // const Double_t kMinZ = 10; // maximal dz distance
3e55050f 1569 const Double_t kMaxDy = 5.; // maximal dy distance
b96c3aef 1570 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1571 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1572 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1573 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1574
1575 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
1576 const Int_t kN=1000; // deepnes of history
1577 static Int_t kglast=0;
1578 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
817766d5 1579 //
3e55050f 1580 // -1. Make a TOF track-
1581 // Clusters are not in friends - use alingment points
1582 //
1583 if (track->GetTOFsignal()<=0) return;
1584 if (!friendTrack->GetTPCOut()) return;
1585 if (!track->GetInnerParam()) return;
1586 if (!track->GetOuterParam()) return;
817766d5 1587 const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
1588 if (!points) return;
3e55050f 1589 AliExternalTrackParam pTPC(*(track->GetOuterParam()));
1590 AliExternalTrackParam pTOF(pTPC);
1591 Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
817766d5 1592 Int_t npoints = points->GetNPoints();
1593 AliTrackPoint point;
3e55050f 1594 Int_t naccept=0;
817766d5 1595 //
1596 for (Int_t ipoint=0;ipoint<npoints;ipoint++){
817766d5 1597 points->GetPoint(point,ipoint);
1598 Float_t xyz[3];
1599 point.GetXYZ(xyz);
1600 Double_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
3e55050f 1601 if (r<350) continue;
817766d5 1602 if (r>400) continue;
3e55050f 1603 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,2.,kTRUE);
1604 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,0.1,kTRUE);
1605 AliTrackPoint lpoint = point.Rotate(pTPC.GetAlpha());
1606 pTPC.PropagateTo(lpoint.GetX(),fMagF);
1607 pTOF=pTPC;
817766d5 1608 ((Double_t*)pTOF.GetParameter())[0] =lpoint.GetY();
1609 ((Double_t*)pTOF.GetParameter())[1] =lpoint.GetZ();
3e55050f 1610 ((Double_t*)pTOF.GetCovariance())[0]+=3.*3./12.;
1611 ((Double_t*)pTOF.GetCovariance())[2]+=3.*3./12.;
0dac7b3a 1612 ((Double_t*)pTOF.GetCovariance())[5]+=0.1*0.1;
1613 ((Double_t*)pTOF.GetCovariance())[9]+=0.1*0.1;
3e55050f 1614 naccept++;
1615 }
1616 if (naccept==0) return; // no tof match clusters
1617 //
1618 // 0. Apply standard cuts
1619 //
1620 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
1621 // exclude crossing track
1622 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1623 //
1624 if (TMath::Abs(pTOF.GetY()-pTPC.GetY()) >kMaxDy) return;
1625 if (TMath::Abs(pTOF.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1626 if (TMath::Abs(pTOF.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1627 //
1628 // 1. Update median and RMS info
1629 //
1630 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1631 TVectorD vecDeltaN(5);
1632 Double_t sign=(pTOF.GetParameter()[1]>0)? 1.:-1.;
1633 vecDelta[4]=0;
1634 for (Int_t i=0;i<4;i++){
1635 vecDelta[i]=(pTOF.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1636 kgdP[i][kglast%kN]=vecDelta[i];
1637 }
1638 kglast=(kglast+1);
1639 Int_t entries=(kglast<kN)?kglast:kN;
1640 Bool_t isOK=kTRUE;
1641 for (Int_t i=0;i<4;i++){
1642 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1643 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1644 vecDeltaN[i] = 0;
1645 if (vecRMS[i]>0.){
1646 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/(vecRMS[i]+1.);
1647 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1648 if (TMath::Abs(vecDeltaN[i])>kSigmaCut) isOK=kFALSE;
817766d5 1649 }
3e55050f 1650 }
1651 //
1652 // 2. Apply median+-rms cut
1653 //
1654 if (kglast<10) return; //median and RMS to be defined
1655 if (!isOK) return;
1656 //
1657 // 3. Update alignment
1658 //
1659 Int_t htime = fTime/3600; //time in hours
1660 if (fAlignTOFTPC->GetEntries()<htime){
1661 fAlignTOFTPC->Expand(htime*2+20);
1662 }
1663 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTOFTPC->At(htime);
1664 if (!align){
1665 // make Alignment object if doesn't exist
1666 align=new AliRelAlignerKalman();
1667 align->SetRunNumber(fRun);
1668 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1669 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1670 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1671 align->SetRejectOutliers(kFALSE);
1672 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
1673 align->SetMagField(fMagF);
1674 fAlignTOFTPC->AddAt(align,htime);
1675 }
1676 align->AddTrackParams(&pTOF,&pTPC);
1677 align->SetTimeStamp(fTime);
1678 align->SetRunNumber(fRun );
1679 //
1680 Int_t nupdates=align->GetNUpdates();
1681 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1682 align->SetRejectOutliers(kFALSE);
1683 TTreeSRedirector *cstream = GetDebugStreamer();
1684 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1685 TTimeStamp tstamp(fTime);
1686 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1687 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1688 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1689 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1690 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1691 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1692 TVectorD vecGoofie(20);
1693 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1694 if (goofieArray){
1695 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1696 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1697 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1698 }
817766d5 1699 }
3e55050f 1700 TVectorD gpTPC(3), gdTPC(3);
1701 TVectorD gpTOF(3), gdTOF(3);
1702 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1703 pTPC.GetDirection(gdTPC.GetMatrixArray());
1704 pTOF.GetXYZ(gpTOF.GetMatrixArray());
1705 pTOF.GetDirection(gdTOF.GetMatrixArray());
32100b1d 1706 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1707 (*cstream)<<"toftpc"<<
1708 "run="<<fRun<< // run number
1709 "event="<<fEvent<< // event number
1710 "time="<<fTime<< // time stamp of event
1711 "trigger="<<fTrigger<< // trigger
1712 "mag="<<fMagF<< // magnetic field
1713 // Environment values
1714 "press0="<<valuePressure0<<
1715 "press1="<<valuePressure1<<
1716 "pt0="<<ptrelative0<<
1717 "pt1="<<ptrelative1<<
1718 "temp0="<<temp0<<
1719 "temp1="<<temp1<<
1720 "vecGoofie.="<<&vecGoofie<<
32100b1d 1721 "vdcorr="<<vdcorr<< // drift correction applied
3e55050f 1722 //
1723 "nmed="<<kglast<< // number of entries to define median and RMS
1724 "vMed.="<<&vecMedian<< // median of deltas
1725 "vRMS.="<<&vecRMS<< // rms of deltas
1726 "vDelta.="<<&vecDelta<< // delta in respect to median
1727 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1728 "t.="<<track<< // ful track - find proper cuts
1729 "a.="<<align<< // current alignment
1730 "pTOF.="<<&pTOF<< // track param TOF
1731 "pTPC.="<<&pTPC<< // track param TPC
1732 "gpTPC.="<<&gpTPC<< // global position TPC
1733 "gdTPC.="<<&gdTPC<< // global direction TPC
1734 "gpTOF.="<<&gpTOF<< // global position TOF
1735 "gdTOF.="<<&gdTOF<< // global position TOF
1736 "\n";
817766d5 1737 }
817766d5 1738}
3e55050f 1739
1740
391ffdb2 1741void AliTPCcalibTime::BookDistortionMaps(){
1742 //
1743 // Book ndimensional histograms of distortions/residuals
1744 // Only primary tracks are selected for analysis
1745 //
1746
1747 Double_t xminTrack[4], xmaxTrack[4];
1748 Int_t binsTrack[4];
1749 TString axisName[4];
1750 TString axisTitle[4];
1751 //
1752 binsTrack[0] =100;
1753 axisName[0] ="#Delta";
1754 axisTitle[0] ="#Delta";
1755 //
1756 binsTrack[1] =20;
1757 xminTrack[1] =-1.5; xmaxTrack[1]=1.5;
1758 axisName[1] ="tanTheta";
1759 axisTitle[1] ="tan(#Theta)";
1760 //
1761 binsTrack[2] =90;
1762 xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi();
1763 axisName[2] ="phi";
1764 axisTitle[2] ="#phi";
1765 //
1766 binsTrack[3] =20;
1767 xminTrack[3] =-1.; xmaxTrack[3]=1.; // 0.33 GeV cut
1768 axisName[3] ="snp";
1769 //
1770 // delta y
9112627b 1771 xminTrack[0] =-1.0; xmaxTrack[0]=1.0; //
391ffdb2 1772 fResHistoTPCITS[0] = new THnSparseS("TPCITS#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1773 fResHistoTPCvertex[0] = new THnSparseS("TPCVertex#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1774 fResHistoTPCTRD[0] = new THnSparseS("TPCTRD#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1775 //
1776 // delta z
9112627b 1777 xminTrack[0] =-3.; xmaxTrack[0]=3.; //
391ffdb2 1778 fResHistoTPCITS[1] = new THnSparseS("TPCITS#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1779 fResHistoTPCvertex[1] = new THnSparseS("TPCVertex#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1780 fResHistoTPCTRD[1] = new THnSparseS("TPCTRD#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
1781 //
1782 // delta snp-P2
9112627b 1783 xminTrack[0] =-0.025; xmaxTrack[0]=0.025; //
391ffdb2 1784 fResHistoTPCITS[2] = new THnSparseS("TPCITS#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
9112627b 1785 fResHistoTPCvertex[2] = new THnSparseS("TPCITSv#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1786 fResHistoTPCTRD[2] = new THnSparseS("TPCTRD#Delta_{#phi}","#Delta_{#phi}", 4, binsTrack,xminTrack, xmaxTrack);
1787 //
1788 // delta theta-P3
9112627b 1789 xminTrack[0] =-0.025; xmaxTrack[0]=0.025; //
391ffdb2 1790 fResHistoTPCITS[3] = new THnSparseS("TPCITS#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
9112627b 1791 fResHistoTPCvertex[3] = new THnSparseS("TPCITSv#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1792 fResHistoTPCTRD[3] = new THnSparseS("TPCTRD#Delta_{#theta}","#Delta_{#theta}", 4, binsTrack,xminTrack, xmaxTrack);
1793 //
1794 // delta theta-P4
9112627b 1795 xminTrack[0] =-0.2; xmaxTrack[0]=0.2; //
391ffdb2 1796 fResHistoTPCITS[4] = new THnSparseS("TPCITS#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
9112627b 1797 fResHistoTPCvertex[4] = new THnSparseS("TPCITSv#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
391ffdb2 1798 fResHistoTPCTRD[4] = new THnSparseS("TPCTRD#Delta_{1/pt}","#Delta_{1/pt}", 4, binsTrack,xminTrack, xmaxTrack);
1799 //
1800 for (Int_t ivar=0;ivar<4;ivar++){
1801 for (Int_t ivar2=0;ivar2<4;ivar2++){
1802 fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1803 fResHistoTPCITS[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
1804 fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1805 fResHistoTPCTRD[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
1806 fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
1807 fResHistoTPCvertex[ivar]->GetAxis(ivar2)->SetTitle(axisTitle[ivar2].Data());
1808 }
1809 }
1810}
1811
1812
1813void AliTPCcalibTime::FillResHistoTPCITS(const AliExternalTrackParam * pTPCIn, const AliExternalTrackParam * pITSOut ){
1814 //
1815 // fill residual histograms pTPCIn-pITSOut
1816 // Histogram is filled only for primary tracks
1817 //
1818 Double_t histoX[4];
1819 Double_t xyz[3];
1820 pTPCIn->GetXYZ(xyz);
1821 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1822 histoX[1]= pTPCIn->GetTgl();
1823 histoX[2]= phi;
1824 histoX[3]= pTPCIn->GetSnp();
9112627b 1825 AliExternalTrackParam lits(*pITSOut);
1826 lits.Rotate(pTPCIn->GetAlpha());
1827 lits.PropagateTo(pTPCIn->GetX(),fMagF);
391ffdb2 1828 //
1829 for (Int_t ihisto=0; ihisto<5; ihisto++){
9112627b 1830 histoX[0]=pTPCIn->GetParameter()[ihisto]-lits.GetParameter()[ihisto];
391ffdb2 1831 fResHistoTPCITS[ihisto]->Fill(histoX);
1832 }
1833}
1834
1835
1836void AliTPCcalibTime::FillResHistoTPC(const AliESDtrack * pTrack){
1837 //
1838 // fill residual histograms pTPC - vertex
1839 // Histogram is filled only for primary tracks
1840 //
1841 Double_t histoX[4];
1842 const AliExternalTrackParam * pTPCIn = pTrack->GetInnerParam();
1843 const AliExternalTrackParam * pTPCvertex = pTrack->GetTPCInnerParam();
9112627b 1844 //
1845 AliExternalTrackParam lits(*pTrack);
1846 lits.Rotate(pTPCvertex->GetAlpha());
1847 lits.PropagateTo(pTPCvertex->GetX(),fMagF);
391ffdb2 1848 Double_t xyz[3];
1849 pTPCIn->GetXYZ(xyz);
1850 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1851 histoX[1]= pTPCIn->GetTgl();
1852 histoX[2]= phi;
1853 histoX[3]= pTPCIn->GetSnp();
1854 //
1855 Float_t dca[2], cov[3];
1856 pTrack->GetImpactParametersTPC(dca,cov);
1857 for (Int_t ihisto=0; ihisto<2; ihisto++){
9112627b 1858 histoX[0]=pTPCvertex->GetParameter()[ihisto]-lits.GetParameter()[ihisto];
1859 fResHistoTPCvertex[ihisto]->Fill(histoX);
391ffdb2 1860 }
1861}
1862
1863
1864void AliTPCcalibTime::FillResHistoTPCTRD(const AliExternalTrackParam * pTPCOut, const AliExternalTrackParam * pTRDIn ){
1865 //
1866 // fill resuidual histogram TPCout-TRDin
1867 //
1868 Double_t histoX[4];
1869 Double_t xyz[3];
1870 pTPCOut->GetXYZ(xyz);
1871 Double_t phi= TMath::ATan2(xyz[1],xyz[0]);
1872 histoX[1]= pTPCOut->GetTgl();
1873 histoX[2]= phi;
1874 histoX[3]= pTPCOut->GetSnp();
1875 //
1876 for (Int_t ihisto=0; ihisto<5; ihisto++){
1877 histoX[0]=pTPCOut->GetParameter()[ihisto]-pTRDIn->GetParameter()[ihisto];
1878 fResHistoTPCTRD[ihisto]->Fill(histoX);
1879 }
1880
1881}