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