]> 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
Adding the normalized dE/dx vs P plot
[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
b96c3aef 1174 const Double_t kMaxDy = 2.; // maximal dy distance
1175 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 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
b96c3aef 1180 const Double_t kMinPt = 0.3; // minimal pt
3e55050f 1181 const Int_t kN=500; // deepnes of history
1182 static Int_t kglast=0;
1183 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
1184 /*
1185 0. Standrd cuts:
1186 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";
1187 */
817766d5 1188 //
3e55050f 1189 // 0. Apply standard cuts
817766d5 1190 //
1191 Int_t dummycl[1000];
1192 if (track->GetITSclusters(dummycl)<kMinITS) return; // minimal amount of clusters
1193 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
b96c3aef 1194 if (!track->IsOn(AliESDtrack::kTPCrefit)) return;
3e55050f 1195 if (!friendTrack->GetITSOut()) return;
817766d5 1196 if (!track->GetInnerParam()) return;
1197 if (!track->GetOuterParam()) return;
b96c3aef 1198 if (track->Pt()<kMinPt) return;
817766d5 1199 // exclude crossing track
1200 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1201 if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ) return;
b96c3aef 1202 if (track->GetInnerParam()->GetX()>90) return;
817766d5 1203 //
1204 AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetInnerParam()));
1205 AliExternalTrackParam pITS(*(friendTrack->GetITSOut()));
3e55050f 1206 pITS.Rotate(pTPC.GetAlpha());
1207 pITS.PropagateTo(pTPC.GetX(),fMagF);
1208 if (TMath::Abs(pITS.GetY()-pTPC.GetY()) >kMaxDy) return;
1209 if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1210 if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1211 //
1212 // 1. Update median and RMS info
1213 //
1214 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1215 TVectorD vecDeltaN(5);
1216 Double_t sign=(pITS.GetParameter()[1]>0)? 1.:-1.;
1217 vecDelta[4]=0;
1218 for (Int_t i=0;i<4;i++){
1219 vecDelta[i]=(pITS.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1220 kgdP[i][kglast%kN]=vecDelta[i];
1221 }
1222 kglast=(kglast+1);
1223 Int_t entries=(kglast<kN)?kglast:kN;
1224 for (Int_t i=0;i<4;i++){
1225 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1226 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1227 vecDeltaN[i] = 0;
1228 if (vecRMS[i]>0.){
1229 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1230 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1231 }
1232 }
817766d5 1233 //
3e55050f 1234 // 2. Apply median+-rms cut
817766d5 1235 //
3e55050f 1236 if (kglast<3) return; //median and RMS to be defined
1237 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1238 //
1239 // 3. Update alignment
817766d5 1240 //
1241 Int_t htime = fTime/3600; //time in hours
1242 if (fAlignITSTPC->GetEntries()<htime){
1243 fAlignITSTPC->Expand(htime*2+20);
1244 }
1245 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignITSTPC->At(htime);
1246 if (!align){
3e55050f 1247 // make Alignment object if doesn't exist
817766d5 1248 align=new AliRelAlignerKalman();
3e55050f 1249 align->SetRunNumber(fRun);
1250 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1251 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1252 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1253 align->SetRejectOutliers(kFALSE);
1254
1255 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
817766d5 1256 align->SetMagField(fMagF);
1257 fAlignITSTPC->AddAt(align,htime);
1258 }
817766d5 1259 align->AddTrackParams(&pITS,&pTPC);
1260 align->SetTimeStamp(fTime);
3e55050f 1261 align->SetRunNumber(fRun );
1262 //
1263 Int_t nupdates=align->GetNUpdates();
1264 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1265 align->SetRejectOutliers(kFALSE);
817766d5 1266 TTreeSRedirector *cstream = GetDebugStreamer();
1267 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1268 TTimeStamp tstamp(fTime);
1269 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1270 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1271 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1272 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1273 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1274 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1275 TVectorD vecGoofie(20);
1276 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1277 if (goofieArray){
1278 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1279 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1280 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1281 }
1282 }
1283 TVectorD gpTPC(3), gdTPC(3);
1284 TVectorD gpITS(3), gdITS(3);
1285 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1286 pTPC.GetDirection(gdTPC.GetMatrixArray());
1287 pITS.GetXYZ(gpITS.GetMatrixArray());
1288 pITS.GetDirection(gdITS.GetMatrixArray());
32100b1d 1289 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
817766d5 1290 (*cstream)<<"itstpc"<<
1291 "run="<<fRun<< // run number
1292 "event="<<fEvent<< // event number
1293 "time="<<fTime<< // time stamp of event
1294 "trigger="<<fTrigger<< // trigger
1295 "mag="<<fMagF<< // magnetic field
1296 // Environment values
1297 "press0="<<valuePressure0<<
1298 "press1="<<valuePressure1<<
1299 "pt0="<<ptrelative0<<
1300 "pt1="<<ptrelative1<<
1301 "temp0="<<temp0<<
1302 "temp1="<<temp1<<
1303 "vecGoofie.="<<&vecGoofie<<
32100b1d 1304 "vdcorr="<<vdcorr<< // drift correction applied
817766d5 1305 //
3e55050f 1306 "nmed="<<kglast<< // number of entries to define median and RMS
1307 "vMed.="<<&vecMedian<< // median of deltas
1308 "vRMS.="<<&vecRMS<< // rms of deltas
1309 "vDelta.="<<&vecDelta<< // delta in respect to median
1310 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1311 "t.="<<track<< // ful track - find proper cuts
1312 "a.="<<align<< // current alignment
1313 "pITS.="<<&pITS<< // track param ITS
1314 "pTPC.="<<&pTPC<< // track param TPC
1315 "gpTPC.="<<&gpTPC<< // global position TPC
1316 "gdTPC.="<<&gdTPC<< // global direction TPC
1317 "gpITS.="<<&gpITS<< // global position ITS
1318 "gdITS.="<<&gdITS<< // global position ITS
817766d5 1319 "\n";
1320 }
817766d5 1321}
3e55050f 1322
1323
1324
1325
0dac7b3a 1326void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack* track, AliESDfriendTrack *friendTrack){
1327 //
3e55050f 1328 // Process track - Update TPC-TRD alignment
1329 // Updates:
1330 // 0. Apply standartd cuts
1331 // 1. Recalucluate the current statistic median/RMS
1332 // 2. Apply median+-rms cut
1333 // 3. Update kalman filter
1334 //
1335 const Int_t kMinTPC = 80; // minimal number of TPC cluster
1336 const Int_t kMinTRD = 50; // minimal number of TRD cluster
1337 const Double_t kMinZ = 20; // maximal dz distance
b96c3aef 1338 const Double_t kMaxDy = 2.; // maximal dy distance
1339 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1340 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1341 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1342 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1343 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
1344 const Int_t kN=500; // deepnes of history
1345 static Int_t kglast=0;
1346 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
0dac7b3a 1347 //
3e55050f 1348 // 0. Apply standard cuts
0dac7b3a 1349 //
1350 Int_t dummycl[1000];
1351 if (track->GetTRDclusters(dummycl)<kMinTRD) return; // minimal amount of clusters
1352 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
3e55050f 1353 if (!friendTrack->GetTRDIn()) return;
0dac7b3a 1354 if (!track->GetInnerParam()) return;
1355 if (!track->GetOuterParam()) return;
1356 // exclude crossing track
1357 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
3e55050f 1358 if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ) return;
0dac7b3a 1359 //
1360 AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetOuterParam()));
1361 AliExternalTrackParam pTRD(*(friendTrack->GetTRDIn()));
3e55050f 1362 pTRD.Rotate(pTPC.GetAlpha());
1363 pTRD.PropagateTo(pTPC.GetX(),fMagF);
1364 ((Double_t*)pTRD.GetCovariance())[2]+=3.*3.; // increas sys errors
1365 ((Double_t*)pTRD.GetCovariance())[9]+=0.1*0.1; // increse sys errors
1366
1367 if (TMath::Abs(pTRD.GetY()-pTPC.GetY()) >kMaxDy) return;
1368 if (TMath::Abs(pTRD.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1369 if (TMath::Abs(pTRD.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
0dac7b3a 1370 //
3e55050f 1371 // 1. Update median and RMS info
0dac7b3a 1372 //
3e55050f 1373 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1374 TVectorD vecDeltaN(5);
1375 Double_t sign=(pTRD.GetParameter()[1]>0)? 1.:-1.;
1376 vecDelta[4]=0;
1377 for (Int_t i=0;i<4;i++){
1378 vecDelta[i]=(pTRD.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1379 kgdP[i][kglast%kN]=vecDelta[i];
1380 }
1381 kglast=(kglast+1);
1382 Int_t entries=(kglast<kN)?kglast:kN;
1383 for (Int_t i=0;i<4;i++){
1384 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1385 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1386 vecDeltaN[i] = 0;
1387 if (vecRMS[i]>0.){
1388 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1389 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1390 }
1391 }
1392 //
1393 // 2. Apply median+-rms cut
1394 //
1395 if (kglast<3) return; //median and RMS to be defined
1396 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1397 //
1398 // 3. Update alignment
0dac7b3a 1399 //
1400 Int_t htime = fTime/3600; //time in hours
1401 if (fAlignTRDTPC->GetEntries()<htime){
1402 fAlignTRDTPC->Expand(htime*2+20);
1403 }
1404 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTRDTPC->At(htime);
1405 if (!align){
3e55050f 1406 // make Alignment object if doesn't exist
0dac7b3a 1407 align=new AliRelAlignerKalman();
3e55050f 1408 align->SetRunNumber(fRun);
1409 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1410 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1411 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1412 align->SetRejectOutliers(kFALSE);
1413 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
0dac7b3a 1414 align->SetMagField(fMagF);
1415 fAlignTRDTPC->AddAt(align,htime);
1416 }
0dac7b3a 1417 align->AddTrackParams(&pTRD,&pTPC);
1418 align->SetTimeStamp(fTime);
3e55050f 1419 align->SetRunNumber(fRun );
1420 //
1421 Int_t nupdates=align->GetNUpdates();
1422 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1423 align->SetRejectOutliers(kFALSE);
0dac7b3a 1424 TTreeSRedirector *cstream = GetDebugStreamer();
1425 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1426 TTimeStamp tstamp(fTime);
1427 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1428 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1429 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1430 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1431 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1432 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1433 TVectorD vecGoofie(20);
1434 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1435 if (goofieArray){
1436 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1437 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1438 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1439 }
1440 }
1441 TVectorD gpTPC(3), gdTPC(3);
1442 TVectorD gpTRD(3), gdTRD(3);
1443 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1444 pTPC.GetDirection(gdTPC.GetMatrixArray());
1445 pTRD.GetXYZ(gpTRD.GetMatrixArray());
1446 pTRD.GetDirection(gdTRD.GetMatrixArray());
32100b1d 1447 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1448 (*cstream)<<"trdtpc"<<
0dac7b3a 1449 "run="<<fRun<< // run number
1450 "event="<<fEvent<< // event number
1451 "time="<<fTime<< // time stamp of event
1452 "trigger="<<fTrigger<< // trigger
1453 "mag="<<fMagF<< // magnetic field
1454 // Environment values
1455 "press0="<<valuePressure0<<
1456 "press1="<<valuePressure1<<
1457 "pt0="<<ptrelative0<<
1458 "pt1="<<ptrelative1<<
1459 "temp0="<<temp0<<
1460 "temp1="<<temp1<<
1461 "vecGoofie.="<<&vecGoofie<<
32100b1d 1462 "vdcorr="<<vdcorr<< // drift correction applied
0dac7b3a 1463 //
3e55050f 1464 "nmed="<<kglast<< // number of entries to define median and RMS
1465 "vMed.="<<&vecMedian<< // median of deltas
1466 "vRMS.="<<&vecRMS<< // rms of deltas
1467 "vDelta.="<<&vecDelta<< // delta in respect to median
1468 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1469 "t.="<<track<< // ful track - find proper cuts
1470 "a.="<<align<< // current alignment
1471 "pTRD.="<<&pTRD<< // track param TRD
1472 "pTPC.="<<&pTPC<< // track param TPC
1473 "gpTPC.="<<&gpTPC<< // global position TPC
1474 "gdTPC.="<<&gdTPC<< // global direction TPC
1475 "gpTRD.="<<&gpTRD<< // global position TRD
1476 "gdTRD.="<<&gdTRD<< // global position TRD
0dac7b3a 1477 "\n";
1478 }
0dac7b3a 1479}
817766d5 1480
1481
1482void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack* track, AliESDfriendTrack *friendTrack){
1483 //
817766d5 1484 //
3e55050f 1485 // Process track - Update TPC-TOF alignment
1486 // Updates:
1487 // -1. Make a TOF "track"
1488 // 0. Apply standartd cuts
1489 // 1. Recalucluate the current statistic median/RMS
1490 // 2. Apply median+-rms cut
1491 // 3. Update kalman filter
1492 //
1493 const Int_t kMinTPC = 80; // minimal number of TPC cluster
1494 const Double_t kMinZ = 10; // maximal dz distance
1495 const Double_t kMaxDy = 5.; // maximal dy distance
b96c3aef 1496 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1497 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1498 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1499 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1500
1501 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
1502 const Int_t kN=1000; // deepnes of history
1503 static Int_t kglast=0;
1504 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
817766d5 1505 //
3e55050f 1506 // -1. Make a TOF track-
1507 // Clusters are not in friends - use alingment points
1508 //
1509 if (track->GetTOFsignal()<=0) return;
1510 if (!friendTrack->GetTPCOut()) return;
1511 if (!track->GetInnerParam()) return;
1512 if (!track->GetOuterParam()) return;
817766d5 1513 const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
1514 if (!points) return;
3e55050f 1515 AliExternalTrackParam pTPC(*(track->GetOuterParam()));
1516 AliExternalTrackParam pTOF(pTPC);
1517 Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
817766d5 1518 Int_t npoints = points->GetNPoints();
1519 AliTrackPoint point;
3e55050f 1520 Int_t naccept=0;
817766d5 1521 //
1522 for (Int_t ipoint=0;ipoint<npoints;ipoint++){
817766d5 1523 points->GetPoint(point,ipoint);
1524 Float_t xyz[3];
1525 point.GetXYZ(xyz);
1526 Double_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
3e55050f 1527 if (r<350) continue;
817766d5 1528 if (r>400) continue;
3e55050f 1529 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,2.,kTRUE);
1530 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,0.1,kTRUE);
1531 AliTrackPoint lpoint = point.Rotate(pTPC.GetAlpha());
1532 pTPC.PropagateTo(lpoint.GetX(),fMagF);
1533 pTOF=pTPC;
817766d5 1534 ((Double_t*)pTOF.GetParameter())[0] =lpoint.GetY();
1535 ((Double_t*)pTOF.GetParameter())[1] =lpoint.GetZ();
3e55050f 1536 ((Double_t*)pTOF.GetCovariance())[0]+=3.*3./12.;
1537 ((Double_t*)pTOF.GetCovariance())[2]+=3.*3./12.;
0dac7b3a 1538 ((Double_t*)pTOF.GetCovariance())[5]+=0.1*0.1;
1539 ((Double_t*)pTOF.GetCovariance())[9]+=0.1*0.1;
3e55050f 1540 naccept++;
1541 }
1542 if (naccept==0) return; // no tof match clusters
1543 //
1544 // 0. Apply standard cuts
1545 //
1546 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
1547 // exclude crossing track
1548 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1549 //
1550 if (TMath::Abs(pTOF.GetY()-pTPC.GetY()) >kMaxDy) return;
1551 if (TMath::Abs(pTOF.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1552 if (TMath::Abs(pTOF.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1553 //
1554 // 1. Update median and RMS info
1555 //
1556 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1557 TVectorD vecDeltaN(5);
1558 Double_t sign=(pTOF.GetParameter()[1]>0)? 1.:-1.;
1559 vecDelta[4]=0;
1560 for (Int_t i=0;i<4;i++){
1561 vecDelta[i]=(pTOF.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1562 kgdP[i][kglast%kN]=vecDelta[i];
1563 }
1564 kglast=(kglast+1);
1565 Int_t entries=(kglast<kN)?kglast:kN;
1566 Bool_t isOK=kTRUE;
1567 for (Int_t i=0;i<4;i++){
1568 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1569 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1570 vecDeltaN[i] = 0;
1571 if (vecRMS[i]>0.){
1572 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/(vecRMS[i]+1.);
1573 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1574 if (TMath::Abs(vecDeltaN[i])>kSigmaCut) isOK=kFALSE;
817766d5 1575 }
3e55050f 1576 }
1577 //
1578 // 2. Apply median+-rms cut
1579 //
1580 if (kglast<10) return; //median and RMS to be defined
1581 if (!isOK) return;
1582 //
1583 // 3. Update alignment
1584 //
1585 Int_t htime = fTime/3600; //time in hours
1586 if (fAlignTOFTPC->GetEntries()<htime){
1587 fAlignTOFTPC->Expand(htime*2+20);
1588 }
1589 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTOFTPC->At(htime);
1590 if (!align){
1591 // make Alignment object if doesn't exist
1592 align=new AliRelAlignerKalman();
1593 align->SetRunNumber(fRun);
1594 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1595 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1596 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1597 align->SetRejectOutliers(kFALSE);
1598 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
1599 align->SetMagField(fMagF);
1600 fAlignTOFTPC->AddAt(align,htime);
1601 }
1602 align->AddTrackParams(&pTOF,&pTPC);
1603 align->SetTimeStamp(fTime);
1604 align->SetRunNumber(fRun );
1605 //
1606 Int_t nupdates=align->GetNUpdates();
1607 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1608 align->SetRejectOutliers(kFALSE);
1609 TTreeSRedirector *cstream = GetDebugStreamer();
1610 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1611 TTimeStamp tstamp(fTime);
1612 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1613 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1614 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1615 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1616 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1617 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1618 TVectorD vecGoofie(20);
1619 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1620 if (goofieArray){
1621 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1622 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1623 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1624 }
817766d5 1625 }
3e55050f 1626 TVectorD gpTPC(3), gdTPC(3);
1627 TVectorD gpTOF(3), gdTOF(3);
1628 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1629 pTPC.GetDirection(gdTPC.GetMatrixArray());
1630 pTOF.GetXYZ(gpTOF.GetMatrixArray());
1631 pTOF.GetDirection(gdTOF.GetMatrixArray());
32100b1d 1632 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1633 (*cstream)<<"toftpc"<<
1634 "run="<<fRun<< // run number
1635 "event="<<fEvent<< // event number
1636 "time="<<fTime<< // time stamp of event
1637 "trigger="<<fTrigger<< // trigger
1638 "mag="<<fMagF<< // magnetic field
1639 // Environment values
1640 "press0="<<valuePressure0<<
1641 "press1="<<valuePressure1<<
1642 "pt0="<<ptrelative0<<
1643 "pt1="<<ptrelative1<<
1644 "temp0="<<temp0<<
1645 "temp1="<<temp1<<
1646 "vecGoofie.="<<&vecGoofie<<
32100b1d 1647 "vdcorr="<<vdcorr<< // drift correction applied
3e55050f 1648 //
1649 "nmed="<<kglast<< // number of entries to define median and RMS
1650 "vMed.="<<&vecMedian<< // median of deltas
1651 "vRMS.="<<&vecRMS<< // rms of deltas
1652 "vDelta.="<<&vecDelta<< // delta in respect to median
1653 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1654 "t.="<<track<< // ful track - find proper cuts
1655 "a.="<<align<< // current alignment
1656 "pTOF.="<<&pTOF<< // track param TOF
1657 "pTPC.="<<&pTPC<< // track param TPC
1658 "gpTPC.="<<&gpTPC<< // global position TPC
1659 "gdTPC.="<<&gdTPC<< // global direction TPC
1660 "gpTOF.="<<&gpTOF<< // global position TOF
1661 "gdTOF.="<<&gdTOF<< // global position TOF
1662 "\n";
817766d5 1663 }
817766d5 1664}
3e55050f 1665
1666