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