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