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