Warning removal.
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTime.cxx
CommitLineData
c74c5f6c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
74235403 17Comments to be written here:
c74c5f6c 18
74235403 191. What do we calibrate.
20
21 Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
c74c5f6c 22
23 AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
24
74235403 252. How to interpret results
26
273. Simple example
c74c5f6c 28
74235403 29 a) determine the required time range:
c74c5f6c 30
74235403 31 AliXRDPROOFtoolkit tool;
32 TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
33 chain->Draw("GetTimeStamp()")
c74c5f6c 34
74235403 35 b) analyse calibration object on Proof in calibration train
c74c5f6c 36
74235403 37 AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift);
c74c5f6c 38
74235403 39 c) plot results
40 .x ~/NimStyle.C
2be25cc9 41 gSystem->Load("libANALYSIS");
42 gSystem->Load("libTPCcalib");
c74c5f6c 43
dde68d36 44 TFile f("CalibObjectsTrain1.root");
45 AliTPCcalibTime *calib = (AliTPCcalibTime *)f->Get("calibTime");
46 calib->GetHistoDrift("all")->Projection(2,0)->Draw()
47 calib->GetFitDrift("all")->Draw("lp")
da6c0bc9 48
74235403 494. Analysis using debug streamers.
da6c0bc9 50
74235403 51 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
52 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
53 AliXRDPROOFtoolkit tool;
3e55050f 54 TChain * chainTime = tool.MakeChainRandom("time.txt","trackInfo",0,10000);
817766d5 55
3e55050f 56 AliXRDPROOFtoolkit::FilterList("timetpctpc.txt","* tpctpc",1)
57 AliXRDPROOFtoolkit::FilterList("timetoftpc.txt","* toftpc",1)
58 AliXRDPROOFtoolkit::FilterList("timeitstpc.txt","* itstpc",1)
59 AliXRDPROOFtoolkit::FilterList("timelaser.txt","* laserInfo",1)
60 TChain * chainTPCTPC = tool.MakeChainRandom("timetpctpc.txt.Good","tpctpc",0,10000);
817766d5 61 TChain * chainTPCITS = tool.MakeChainRandom("timeitstpc.txt.Good","itstpc",0,10000);
3e55050f 62 TChain * chainTPCTOF = tool.MakeChainRandom("timetoftpc.txt.Good","toftpc",0,10000);
817766d5 63 TChain * chainLaser = tool.MakeChainRandom("timelaser.txt.Good","laserInfo",0,10000);
74235403 64 chainTime->Lookup();
dde68d36 65 chainLaser->Lookup();
c74c5f6c 66*/
67
c74c5f6c 68#include "Riostream.h"
69#include "TChain.h"
70#include "TTree.h"
71#include "TH1F.h"
72#include "TH2F.h"
73#include "TH3F.h"
74#include "THnSparse.h"
75#include "TList.h"
76#include "TMath.h"
77#include "TCanvas.h"
78#include "TFile.h"
79#include "TF1.h"
80#include "TVectorD.h"
81#include "TProfile.h"
82#include "TGraphErrors.h"
83#include "TCanvas.h"
c74c5f6c 84#include "AliTPCclusterMI.h"
85#include "AliTPCseed.h"
86#include "AliESDVertex.h"
87#include "AliESDEvent.h"
88#include "AliESDfriend.h"
89#include "AliESDInputHandler.h"
90#include "AliAnalysisManager.h"
91
92#include "AliTracker.h"
f7a1cc68 93#include "AliMagF.h"
c74c5f6c 94#include "AliTPCCalROC.h"
3e55050f 95#include "AliTPCParam.h"
c74c5f6c 96
97#include "AliLog.h"
98
99#include "AliTPCcalibTime.h"
817766d5 100#include "AliRelAlignerKalman.h"
c74c5f6c 101
102#include "TTreeStream.h"
103#include "AliTPCTracklet.h"
da6c0bc9 104#include "TTimeStamp.h"
105#include "AliTPCcalibDB.h"
106#include "AliTPCcalibLaser.h"
31aa7c5c 107#include "AliDCSSensorArray.h"
108#include "AliDCSSensor.h"
c74c5f6c 109
817766d5 110#include "TDatabasePDG.h"
111#include "AliTrackPointArray.h"
112
c74c5f6c 113ClassImp(AliTPCcalibTime)
114
115
116AliTPCcalibTime::AliTPCcalibTime()
da6c0bc9 117 :AliTPCcalibBase(),
da6c0bc9 118 fLaser(0), // pointer to laser calibration
119 fDz(0), // current delta z
d3ce44cb 120 fCutMaxD(3), // maximal distance in rfi ditection
121 fCutMaxDz(25), // maximal distance in rfi ditection
c74c5f6c 122 fCutTheta(0.03), // maximal distan theta
2be25cc9 123 fCutMinDir(-0.99), // direction vector products
b842d904 124 fCutTracks(100),
dde68d36 125 fArrayDz(0), //NEW! Tmap of V drifts for different triggers
817766d5 126 fAlignITSTPC(0), //alignemnt array ITS TPC match
127 fAlignTRDTPC(0), //alignemnt array TRD TPC match
128 fAlignTOFTPC(0), //alignemnt array TOF TPC match
2be25cc9 129 fTimeBins(0),
130 fTimeStart(0),
131 fTimeEnd(0),
132 fPtBins(0),
133 fPtStart(0),
134 fPtEnd(0),
135 fVdriftBins(0),
136 fVdriftStart(0),
137 fVdriftEnd(0),
138 fRunBins(0),
139 fRunStart(0),
140 fRunEnd(0)
141// fBinsVdrift(fTimeBins,fPtBins,fVdriftBins),
142// fXminVdrift(fTimeStart,fPtStart,fVdriftStart),
143// fXmaxVdrift(fTimeEnd,fPtEnd,fVdriftEnd)
c74c5f6c 144{
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
b842d904 163 fCutTracks(100),
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);
b842d904 543 if (friendTrack) ProcessAlignITS(track,friendTrack,event,ESDfriend);
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
b842d904 1162void AliTPCcalibTime::ProcessAlignITS(AliESDtrack* track, AliESDfriendTrack *friendTrack, AliESDEvent *event,AliESDfriend *ESDfriend){
817766d5 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];
817766d5 1192 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
b842d904 1193 if (track->GetITSclusters(dummycl)<kMinITS) return; // minimal amount of clusters
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;
b842d904 1198 if (track->GetInnerParam()->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()));
b842d904 1205 AliExternalTrackParam pITS(*(friendTrack->GetITSOut())); // ITS standalone if possible
1206 AliExternalTrackParam pITS2(*(friendTrack->GetITSOut())); //TPC-ITS track
1207 pITS2.Rotate(pTPC.GetAlpha());
1208 pITS2.PropagateTo(pTPC.GetX(),fMagF);
1209 AliESDfriendTrack *itsfriendTrack=0;
1210 //
1211 // try to find standalone ITS track corresponing to the TPC if possible
1212 //
1213 Bool_t hasAlone=kFALSE;
1214 Int_t ntracks=event->GetNumberOfTracks();
1215 for (Int_t i=0; i<ntracks; i++){
1216 AliESDtrack *trackS = event->GetTrack(i);
1217 if (trackS->GetTPCNcls()>0) continue; //continue if has TPC info
1218 itsfriendTrack = ESDfriend->GetTrack(i);
1219 if (!itsfriendTrack) continue;
1220 if (!itsfriendTrack->GetITSOut()) continue;
1221 if (TMath::Abs(pITS2.GetTgl()-itsfriendTrack->GetITSOut()->GetTgl())> kMaxAngle) continue;
1222 pITS=(*(itsfriendTrack->GetITSOut()));
1223 //
1224 pITS.Rotate(pTPC.GetAlpha());
1225 pITS.PropagateTo(pTPC.GetX(),fMagF);
1226 if (TMath::Abs(pITS2.GetY()-pITS.GetY())> kMaxDy) continue;
1227 hasAlone=kTRUE;
1228 }
1229 if (!hasAlone) pITS=pITS2;
1230 //
3e55050f 1231 if (TMath::Abs(pITS.GetY()-pTPC.GetY()) >kMaxDy) return;
1232 if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1233 if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1234 //
1235 // 1. Update median and RMS info
1236 //
1237 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1238 TVectorD vecDeltaN(5);
1239 Double_t sign=(pITS.GetParameter()[1]>0)? 1.:-1.;
1240 vecDelta[4]=0;
1241 for (Int_t i=0;i<4;i++){
1242 vecDelta[i]=(pITS.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1243 kgdP[i][kglast%kN]=vecDelta[i];
1244 }
1245 kglast=(kglast+1);
1246 Int_t entries=(kglast<kN)?kglast:kN;
1247 for (Int_t i=0;i<4;i++){
1248 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1249 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1250 vecDeltaN[i] = 0;
1251 if (vecRMS[i]>0.){
1252 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1253 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1254 }
1255 }
817766d5 1256 //
3e55050f 1257 // 2. Apply median+-rms cut
817766d5 1258 //
3e55050f 1259 if (kglast<3) return; //median and RMS to be defined
1260 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1261 //
1262 // 3. Update alignment
817766d5 1263 //
1264 Int_t htime = fTime/3600; //time in hours
1265 if (fAlignITSTPC->GetEntries()<htime){
1266 fAlignITSTPC->Expand(htime*2+20);
1267 }
1268 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignITSTPC->At(htime);
1269 if (!align){
3e55050f 1270 // make Alignment object if doesn't exist
817766d5 1271 align=new AliRelAlignerKalman();
3e55050f 1272 align->SetRunNumber(fRun);
1273 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1274 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1275 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1276 align->SetRejectOutliers(kFALSE);
1277
1278 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
817766d5 1279 align->SetMagField(fMagF);
1280 fAlignITSTPC->AddAt(align,htime);
1281 }
817766d5 1282 align->AddTrackParams(&pITS,&pTPC);
1283 align->SetTimeStamp(fTime);
3e55050f 1284 align->SetRunNumber(fRun );
1285 //
1286 Int_t nupdates=align->GetNUpdates();
1287 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1288 align->SetRejectOutliers(kFALSE);
817766d5 1289 TTreeSRedirector *cstream = GetDebugStreamer();
1290 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1291 TTimeStamp tstamp(fTime);
1292 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1293 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1294 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1295 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1296 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1297 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1298 TVectorD vecGoofie(20);
1299 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1300 if (goofieArray){
1301 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1302 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1303 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1304 }
1305 }
1306 TVectorD gpTPC(3), gdTPC(3);
1307 TVectorD gpITS(3), gdITS(3);
1308 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1309 pTPC.GetDirection(gdTPC.GetMatrixArray());
1310 pITS.GetXYZ(gpITS.GetMatrixArray());
1311 pITS.GetDirection(gdITS.GetMatrixArray());
32100b1d 1312 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
817766d5 1313 (*cstream)<<"itstpc"<<
1314 "run="<<fRun<< // run number
1315 "event="<<fEvent<< // event number
1316 "time="<<fTime<< // time stamp of event
1317 "trigger="<<fTrigger<< // trigger
1318 "mag="<<fMagF<< // magnetic field
1319 // Environment values
1320 "press0="<<valuePressure0<<
1321 "press1="<<valuePressure1<<
1322 "pt0="<<ptrelative0<<
1323 "pt1="<<ptrelative1<<
1324 "temp0="<<temp0<<
1325 "temp1="<<temp1<<
1326 "vecGoofie.="<<&vecGoofie<<
32100b1d 1327 "vdcorr="<<vdcorr<< // drift correction applied
817766d5 1328 //
b842d904 1329 "hasAlone="<<hasAlone<< // has ITS standalone ?
1330 "track.="<<track<< // track info
3e55050f 1331 "nmed="<<kglast<< // number of entries to define median and RMS
1332 "vMed.="<<&vecMedian<< // median of deltas
1333 "vRMS.="<<&vecRMS<< // rms of deltas
1334 "vDelta.="<<&vecDelta<< // delta in respect to median
1335 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1336 "t.="<<track<< // ful track - find proper cuts
1337 "a.="<<align<< // current alignment
b842d904 1338 "pITS.="<<&pITS<< // track param ITS
1339 "pITS2.="<<&pITS2<< // track param ITS+TPC
3e55050f 1340 "pTPC.="<<&pTPC<< // track param TPC
1341 "gpTPC.="<<&gpTPC<< // global position TPC
1342 "gdTPC.="<<&gdTPC<< // global direction TPC
1343 "gpITS.="<<&gpITS<< // global position ITS
1344 "gdITS.="<<&gdITS<< // global position ITS
817766d5 1345 "\n";
1346 }
817766d5 1347}
3e55050f 1348
1349
1350
1351
0dac7b3a 1352void AliTPCcalibTime::ProcessAlignTRD(AliESDtrack* track, AliESDfriendTrack *friendTrack){
1353 //
3e55050f 1354 // Process track - Update TPC-TRD alignment
1355 // Updates:
1356 // 0. Apply standartd cuts
1357 // 1. Recalucluate the current statistic median/RMS
1358 // 2. Apply median+-rms cut
1359 // 3. Update kalman filter
1360 //
1361 const Int_t kMinTPC = 80; // minimal number of TPC cluster
1362 const Int_t kMinTRD = 50; // minimal number of TRD cluster
1363 const Double_t kMinZ = 20; // maximal dz distance
b96c3aef 1364 const Double_t kMaxDy = 2.; // maximal dy distance
1365 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1366 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1367 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1368 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1369 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
1370 const Int_t kN=500; // deepnes of history
1371 static Int_t kglast=0;
1372 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
0dac7b3a 1373 //
3e55050f 1374 // 0. Apply standard cuts
0dac7b3a 1375 //
1376 Int_t dummycl[1000];
1377 if (track->GetTRDclusters(dummycl)<kMinTRD) return; // minimal amount of clusters
1378 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
3e55050f 1379 if (!friendTrack->GetTRDIn()) return;
0dac7b3a 1380 if (!track->GetInnerParam()) return;
1381 if (!track->GetOuterParam()) return;
1382 // exclude crossing track
1383 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
3e55050f 1384 if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ) return;
0dac7b3a 1385 //
1386 AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetOuterParam()));
1387 AliExternalTrackParam pTRD(*(friendTrack->GetTRDIn()));
3e55050f 1388 pTRD.Rotate(pTPC.GetAlpha());
1389 pTRD.PropagateTo(pTPC.GetX(),fMagF);
1390 ((Double_t*)pTRD.GetCovariance())[2]+=3.*3.; // increas sys errors
1391 ((Double_t*)pTRD.GetCovariance())[9]+=0.1*0.1; // increse sys errors
1392
1393 if (TMath::Abs(pTRD.GetY()-pTPC.GetY()) >kMaxDy) return;
1394 if (TMath::Abs(pTRD.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1395 if (TMath::Abs(pTRD.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
0dac7b3a 1396 //
3e55050f 1397 // 1. Update median and RMS info
0dac7b3a 1398 //
3e55050f 1399 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1400 TVectorD vecDeltaN(5);
1401 Double_t sign=(pTRD.GetParameter()[1]>0)? 1.:-1.;
1402 vecDelta[4]=0;
1403 for (Int_t i=0;i<4;i++){
1404 vecDelta[i]=(pTRD.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1405 kgdP[i][kglast%kN]=vecDelta[i];
1406 }
1407 kglast=(kglast+1);
1408 Int_t entries=(kglast<kN)?kglast:kN;
1409 for (Int_t i=0;i<4;i++){
1410 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1411 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1412 vecDeltaN[i] = 0;
1413 if (vecRMS[i]>0.){
1414 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
1415 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1416 }
1417 }
1418 //
1419 // 2. Apply median+-rms cut
1420 //
1421 if (kglast<3) return; //median and RMS to be defined
1422 if ( vecDeltaN[4]/4.>kSigmaCut) return;
1423 //
1424 // 3. Update alignment
0dac7b3a 1425 //
1426 Int_t htime = fTime/3600; //time in hours
1427 if (fAlignTRDTPC->GetEntries()<htime){
1428 fAlignTRDTPC->Expand(htime*2+20);
1429 }
1430 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTRDTPC->At(htime);
1431 if (!align){
3e55050f 1432 // make Alignment object if doesn't exist
0dac7b3a 1433 align=new AliRelAlignerKalman();
3e55050f 1434 align->SetRunNumber(fRun);
1435 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1436 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1437 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1438 align->SetRejectOutliers(kFALSE);
1439 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
0dac7b3a 1440 align->SetMagField(fMagF);
1441 fAlignTRDTPC->AddAt(align,htime);
1442 }
0dac7b3a 1443 align->AddTrackParams(&pTRD,&pTPC);
1444 align->SetTimeStamp(fTime);
3e55050f 1445 align->SetRunNumber(fRun );
1446 //
1447 Int_t nupdates=align->GetNUpdates();
1448 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1449 align->SetRejectOutliers(kFALSE);
0dac7b3a 1450 TTreeSRedirector *cstream = GetDebugStreamer();
1451 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1452 TTimeStamp tstamp(fTime);
1453 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1454 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1455 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1456 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1457 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1458 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1459 TVectorD vecGoofie(20);
1460 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1461 if (goofieArray){
1462 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1463 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1464 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1465 }
1466 }
1467 TVectorD gpTPC(3), gdTPC(3);
1468 TVectorD gpTRD(3), gdTRD(3);
1469 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1470 pTPC.GetDirection(gdTPC.GetMatrixArray());
1471 pTRD.GetXYZ(gpTRD.GetMatrixArray());
1472 pTRD.GetDirection(gdTRD.GetMatrixArray());
32100b1d 1473 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1474 (*cstream)<<"trdtpc"<<
0dac7b3a 1475 "run="<<fRun<< // run number
1476 "event="<<fEvent<< // event number
1477 "time="<<fTime<< // time stamp of event
1478 "trigger="<<fTrigger<< // trigger
1479 "mag="<<fMagF<< // magnetic field
1480 // Environment values
1481 "press0="<<valuePressure0<<
1482 "press1="<<valuePressure1<<
1483 "pt0="<<ptrelative0<<
1484 "pt1="<<ptrelative1<<
1485 "temp0="<<temp0<<
1486 "temp1="<<temp1<<
1487 "vecGoofie.="<<&vecGoofie<<
32100b1d 1488 "vdcorr="<<vdcorr<< // drift correction applied
0dac7b3a 1489 //
3e55050f 1490 "nmed="<<kglast<< // number of entries to define median and RMS
1491 "vMed.="<<&vecMedian<< // median of deltas
1492 "vRMS.="<<&vecRMS<< // rms of deltas
1493 "vDelta.="<<&vecDelta<< // delta in respect to median
1494 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1495 "t.="<<track<< // ful track - find proper cuts
1496 "a.="<<align<< // current alignment
1497 "pTRD.="<<&pTRD<< // track param TRD
1498 "pTPC.="<<&pTPC<< // track param TPC
1499 "gpTPC.="<<&gpTPC<< // global position TPC
1500 "gdTPC.="<<&gdTPC<< // global direction TPC
1501 "gpTRD.="<<&gpTRD<< // global position TRD
1502 "gdTRD.="<<&gdTRD<< // global position TRD
0dac7b3a 1503 "\n";
1504 }
0dac7b3a 1505}
817766d5 1506
1507
1508void AliTPCcalibTime::ProcessAlignTOF(AliESDtrack* track, AliESDfriendTrack *friendTrack){
1509 //
817766d5 1510 //
3e55050f 1511 // Process track - Update TPC-TOF alignment
1512 // Updates:
1513 // -1. Make a TOF "track"
1514 // 0. Apply standartd cuts
1515 // 1. Recalucluate the current statistic median/RMS
1516 // 2. Apply median+-rms cut
1517 // 3. Update kalman filter
1518 //
1519 const Int_t kMinTPC = 80; // minimal number of TPC cluster
b842d904 1520 // const Double_t kMinZ = 10; // maximal dz distance
3e55050f 1521 const Double_t kMaxDy = 5.; // maximal dy distance
b96c3aef 1522 const Double_t kMaxAngle= 0.015; // maximal angular distance
3e55050f 1523 const Double_t kSigmaCut= 5; // maximal sigma distance to median
1524 const Double_t kVdErr = 0.1; // initial uncertainty of the vd correction
1525 const Double_t kVdYErr = 0.05; // initial uncertainty of the vd correction
1526
1527 const Double_t kOutCut = 1.0; // outlyer cut in AliRelAlgnmentKalman
1528 const Int_t kN=1000; // deepnes of history
1529 static Int_t kglast=0;
1530 static Double_t* kgdP[4]={new Double_t[kN], new Double_t[kN], new Double_t[kN], new Double_t[kN]};
817766d5 1531 //
3e55050f 1532 // -1. Make a TOF track-
1533 // Clusters are not in friends - use alingment points
1534 //
1535 if (track->GetTOFsignal()<=0) return;
1536 if (!friendTrack->GetTPCOut()) return;
1537 if (!track->GetInnerParam()) return;
1538 if (!track->GetOuterParam()) return;
817766d5 1539 const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
1540 if (!points) return;
3e55050f 1541 AliExternalTrackParam pTPC(*(track->GetOuterParam()));
1542 AliExternalTrackParam pTOF(pTPC);
1543 Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
817766d5 1544 Int_t npoints = points->GetNPoints();
1545 AliTrackPoint point;
3e55050f 1546 Int_t naccept=0;
817766d5 1547 //
1548 for (Int_t ipoint=0;ipoint<npoints;ipoint++){
817766d5 1549 points->GetPoint(point,ipoint);
1550 Float_t xyz[3];
1551 point.GetXYZ(xyz);
1552 Double_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
3e55050f 1553 if (r<350) continue;
817766d5 1554 if (r>400) continue;
3e55050f 1555 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,2.,kTRUE);
1556 AliTracker::PropagateTrackToBxByBz(&pTPC,r,mass,0.1,kTRUE);
1557 AliTrackPoint lpoint = point.Rotate(pTPC.GetAlpha());
1558 pTPC.PropagateTo(lpoint.GetX(),fMagF);
1559 pTOF=pTPC;
817766d5 1560 ((Double_t*)pTOF.GetParameter())[0] =lpoint.GetY();
1561 ((Double_t*)pTOF.GetParameter())[1] =lpoint.GetZ();
3e55050f 1562 ((Double_t*)pTOF.GetCovariance())[0]+=3.*3./12.;
1563 ((Double_t*)pTOF.GetCovariance())[2]+=3.*3./12.;
0dac7b3a 1564 ((Double_t*)pTOF.GetCovariance())[5]+=0.1*0.1;
1565 ((Double_t*)pTOF.GetCovariance())[9]+=0.1*0.1;
3e55050f 1566 naccept++;
1567 }
1568 if (naccept==0) return; // no tof match clusters
1569 //
1570 // 0. Apply standard cuts
1571 //
1572 if (track->GetTPCNcls()<kMinTPC) return; // minimal amount of clusters cut
1573 // exclude crossing track
1574 if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0) return;
1575 //
1576 if (TMath::Abs(pTOF.GetY()-pTPC.GetY()) >kMaxDy) return;
1577 if (TMath::Abs(pTOF.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1578 if (TMath::Abs(pTOF.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1579 //
1580 // 1. Update median and RMS info
1581 //
1582 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
1583 TVectorD vecDeltaN(5);
1584 Double_t sign=(pTOF.GetParameter()[1]>0)? 1.:-1.;
1585 vecDelta[4]=0;
1586 for (Int_t i=0;i<4;i++){
1587 vecDelta[i]=(pTOF.GetParameter()[i]-pTPC.GetParameter()[i])*sign;
1588 kgdP[i][kglast%kN]=vecDelta[i];
1589 }
1590 kglast=(kglast+1);
1591 Int_t entries=(kglast<kN)?kglast:kN;
1592 Bool_t isOK=kTRUE;
1593 for (Int_t i=0;i<4;i++){
1594 vecMedian[i] = TMath::Median(entries,kgdP[i]);
1595 vecRMS[i] = TMath::RMS(entries,kgdP[i]);
1596 vecDeltaN[i] = 0;
1597 if (vecRMS[i]>0.){
1598 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/(vecRMS[i]+1.);
1599 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
1600 if (TMath::Abs(vecDeltaN[i])>kSigmaCut) isOK=kFALSE;
817766d5 1601 }
3e55050f 1602 }
1603 //
1604 // 2. Apply median+-rms cut
1605 //
1606 if (kglast<10) return; //median and RMS to be defined
1607 if (!isOK) return;
1608 //
1609 // 3. Update alignment
1610 //
1611 Int_t htime = fTime/3600; //time in hours
1612 if (fAlignTOFTPC->GetEntries()<htime){
1613 fAlignTOFTPC->Expand(htime*2+20);
1614 }
1615 AliRelAlignerKalman* align = (AliRelAlignerKalman*)fAlignTOFTPC->At(htime);
1616 if (!align){
1617 // make Alignment object if doesn't exist
1618 align=new AliRelAlignerKalman();
1619 align->SetRunNumber(fRun);
1620 (*align->GetStateCov())(6,6)=kVdErr*kVdErr;
1621 (*align->GetStateCov())(8,8)=kVdYErr*kVdYErr;
1622 align->SetOutRejSigma(kOutCut+kOutCut*kN);
1623 align->SetRejectOutliers(kFALSE);
1624 align->SetTPCvd(AliTPCcalibDB::Instance()->GetParameters()->GetDriftV()/1000000.);
1625 align->SetMagField(fMagF);
1626 fAlignTOFTPC->AddAt(align,htime);
1627 }
1628 align->AddTrackParams(&pTOF,&pTPC);
1629 align->SetTimeStamp(fTime);
1630 align->SetRunNumber(fRun );
1631 //
1632 Int_t nupdates=align->GetNUpdates();
1633 align->SetOutRejSigma(kOutCut+kOutCut*kN/Double_t(nupdates));
1634 align->SetRejectOutliers(kFALSE);
1635 TTreeSRedirector *cstream = GetDebugStreamer();
1636 if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1637 TTimeStamp tstamp(fTime);
1638 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1639 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1640 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1641 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1642 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1643 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1644 TVectorD vecGoofie(20);
1645 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1646 if (goofieArray){
1647 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1648 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1649 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1650 }
817766d5 1651 }
3e55050f 1652 TVectorD gpTPC(3), gdTPC(3);
1653 TVectorD gpTOF(3), gdTOF(3);
1654 pTPC.GetXYZ(gpTPC.GetMatrixArray());
1655 pTPC.GetDirection(gdTPC.GetMatrixArray());
1656 pTOF.GetXYZ(gpTOF.GetMatrixArray());
1657 pTOF.GetDirection(gdTOF.GetMatrixArray());
32100b1d 1658 Double_t vdcorr = AliTPCcalibDB::Instance()->GetVDriftCorrectionTime(tstamp,fRun,0,1);
3e55050f 1659 (*cstream)<<"toftpc"<<
1660 "run="<<fRun<< // run number
1661 "event="<<fEvent<< // event number
1662 "time="<<fTime<< // time stamp of event
1663 "trigger="<<fTrigger<< // trigger
1664 "mag="<<fMagF<< // magnetic field
1665 // Environment values
1666 "press0="<<valuePressure0<<
1667 "press1="<<valuePressure1<<
1668 "pt0="<<ptrelative0<<
1669 "pt1="<<ptrelative1<<
1670 "temp0="<<temp0<<
1671 "temp1="<<temp1<<
1672 "vecGoofie.="<<&vecGoofie<<
32100b1d 1673 "vdcorr="<<vdcorr<< // drift correction applied
3e55050f 1674 //
1675 "nmed="<<kglast<< // number of entries to define median and RMS
1676 "vMed.="<<&vecMedian<< // median of deltas
1677 "vRMS.="<<&vecRMS<< // rms of deltas
1678 "vDelta.="<<&vecDelta<< // delta in respect to median
1679 "vDeltaN.="<<&vecDeltaN<< // normalized delta in respect to median
1680 "t.="<<track<< // ful track - find proper cuts
1681 "a.="<<align<< // current alignment
1682 "pTOF.="<<&pTOF<< // track param TOF
1683 "pTPC.="<<&pTPC<< // track param TPC
1684 "gpTPC.="<<&gpTPC<< // global position TPC
1685 "gdTPC.="<<&gdTPC<< // global direction TPC
1686 "gpTOF.="<<&gpTOF<< // global position TOF
1687 "gdTOF.="<<&gdTOF<< // global position TOF
1688 "\n";
817766d5 1689 }
817766d5 1690}
3e55050f 1691
1692