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