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