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