Completely misleading method name changed.
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTime.cxx
CommitLineData
d3ce44cb 1
c74c5f6c 2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
17/*
74235403 18Comments to be written here:
c74c5f6c 19
74235403 201. What do we calibrate.
21
22 Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
c74c5f6c 23
24 AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
25
74235403 262. How to interpret results
27
283. Simple example
c74c5f6c 29
74235403 30 a) determine the required time range:
c74c5f6c 31
74235403 32 AliXRDPROOFtoolkit tool;
33 TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
34 chain->Draw("GetTimeStamp()")
c74c5f6c 35
74235403 36 b) analyse calibration object on Proof in calibration train
c74c5f6c 37
74235403 38 AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift);
c74c5f6c 39
74235403 40 c) plot results
41 .x ~/NimStyle.C
2be25cc9 42 gSystem->Load("libANALYSIS");
43 gSystem->Load("libTPCcalib");
c74c5f6c 44
74235403 45 TFile f("CalibObjects.root");
46 AliTPCcalibTime *cal = (AliTPCcalibTime *)f->Get("TPCCalib")->FindObject("calibTime");
47 cal->GetHistoDrift("all")->Projection(2,0)->Draw()
48 cal->GetFitDrift("all")->Draw("lp")
da6c0bc9 49
74235403 504. Analysis using debug streamers.
da6c0bc9 51
74235403 52 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
53 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
54 AliXRDPROOFtoolkit tool;
55 TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
56 chainTime->Lookup();
c74c5f6c 57*/
58
c74c5f6c 59#include "Riostream.h"
60#include "TChain.h"
61#include "TTree.h"
62#include "TH1F.h"
63#include "TH2F.h"
64#include "TH3F.h"
65#include "THnSparse.h"
66#include "TList.h"
67#include "TMath.h"
68#include "TCanvas.h"
69#include "TFile.h"
70#include "TF1.h"
71#include "TVectorD.h"
72#include "TProfile.h"
73#include "TGraphErrors.h"
74#include "TCanvas.h"
75
76#include "AliTPCclusterMI.h"
77#include "AliTPCseed.h"
78#include "AliESDVertex.h"
79#include "AliESDEvent.h"
80#include "AliESDfriend.h"
81#include "AliESDInputHandler.h"
82#include "AliAnalysisManager.h"
83
84#include "AliTracker.h"
f7a1cc68 85#include "AliMagF.h"
c74c5f6c 86#include "AliTPCCalROC.h"
87
88#include "AliLog.h"
89
90#include "AliTPCcalibTime.h"
91
92#include "TTreeStream.h"
93#include "AliTPCTracklet.h"
da6c0bc9 94#include "TTimeStamp.h"
95#include "AliTPCcalibDB.h"
96#include "AliTPCcalibLaser.h"
31aa7c5c 97#include "AliDCSSensorArray.h"
98#include "AliDCSSensor.h"
c74c5f6c 99
100ClassImp(AliTPCcalibTime)
101
102
103AliTPCcalibTime::AliTPCcalibTime()
da6c0bc9 104 :AliTPCcalibBase(),
da6c0bc9 105 fLaser(0), // pointer to laser calibration
106 fDz(0), // current delta z
d3ce44cb 107 fCutMaxD(3), // maximal distance in rfi ditection
108 fCutMaxDz(25), // maximal distance in rfi ditection
c74c5f6c 109 fCutTheta(0.03), // maximal distan theta
2be25cc9 110 fCutMinDir(-0.99), // direction vector products
74235403 111 fCutTracks(10),
2be25cc9 112 fMapDz(0), //NEW! Tmap of V drifts for different triggers
113 fTimeBins(0),
114 fTimeStart(0),
115 fTimeEnd(0),
116 fPtBins(0),
117 fPtStart(0),
118 fPtEnd(0),
119 fVdriftBins(0),
120 fVdriftStart(0),
121 fVdriftEnd(0),
122 fRunBins(0),
123 fRunStart(0),
124 fRunEnd(0)
125// fBinsVdrift(fTimeBins,fPtBins,fVdriftBins),
126// fXminVdrift(fTimeStart,fPtStart,fVdriftStart),
127// fXmaxVdrift(fTimeEnd,fPtEnd,fVdriftEnd)
c74c5f6c 128
129{
130 AliInfo("Default Constructor");
2be25cc9 131 for (Int_t i=0;i<3;i++) {
132 fHistVdriftLaserA[i]=0;
133 fHistVdriftLaserC[i]=0;
134 }
d3ce44cb 135 for (Int_t i=0;i<10;i++) {
74235403 136 fCosmiMatchingHisto[i]=0;
137 }
c74c5f6c 138}
139
140
74235403 141AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeVdrift)
c74c5f6c 142 :AliTPCcalibBase(),
d3ce44cb 143 fLaser(0), // pointer to laser calibration
144 fDz(0), // current delta z
145 fCutMaxD(5*0.5356), // maximal distance in rfi ditection
146 fCutMaxDz(5*4.541), // maximal distance in rfi ditection
147 fCutTheta(5*0.004644),// maximal distan theta
148 fCutMinDir(-0.99), // direction vector products
74235403 149 fCutTracks(10),
d3ce44cb 150 fMapDz(0), //Tmap of V drifts for different triggers
2be25cc9 151 fTimeBins(0),
152 fTimeStart(0),
153 fTimeEnd(0),
154 fPtBins(0),
155 fPtStart(0),
156 fPtEnd(0),
157 fVdriftBins(0),
158 fVdriftStart(0),
159 fVdriftEnd(0),
160 fRunBins(0),
161 fRunStart(0),
162 fRunEnd(0)
c74c5f6c 163{
164
165 SetName(name);
166 SetTitle(title);
2be25cc9 167 for (Int_t i=0;i<3;i++) {
168 fHistVdriftLaserA[i]=0;
169 fHistVdriftLaserC[i]=0;
170 }
c74c5f6c 171
172 AliInfo("Non Default Constructor");
173
74235403 174 fTimeBins =(EndTime-StartTime)/deltaIntegrationTimeVdrift;
2be25cc9 175 fTimeStart =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi();
176 fTimeEnd =EndTime; //(((TObjString*)(mapGRP->GetValue("fAliceStopTime")))->GetString()).Atoi();
177 fPtBins = 200;
178 fPtStart = -0.04;
179 fPtEnd = 0.04;
180 fVdriftBins = 200;
181 fVdriftStart= -20.0/500.0;
182 fVdriftEnd = 20.0/500.0;
74235403 183 fRunBins = 100000;
2be25cc9 184 fRunStart = -0.5;
d3ce44cb 185 fRunEnd = 99999.5;
2be25cc9 186
187 Int_t binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins };
188 Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart , fRunStart};
189 Double_t xmaxVdriftLaser[4] = {fTimeEnd , fPtEnd , fVdriftEnd , fRunEnd };
190
191 for (Int_t i=0;i<3;i++) {
192 fHistVdriftLaserA[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
193 fHistVdriftLaserC[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
194 }
195 fBinsVdrift[0] = fTimeBins;
196 fBinsVdrift[1] = fPtBins;
197 fBinsVdrift[2] = fVdriftBins;
198 fBinsVdrift[3] = fRunBins;
199 fXminVdrift[0] = fTimeStart;
200 fXminVdrift[1] = fPtStart;
201 fXminVdrift[2] = fVdriftStart;
202 fXminVdrift[3] = fRunStart;
203 fXmaxVdrift[0] = fTimeEnd;
204 fXmaxVdrift[1] = fPtEnd;
205 fXmaxVdrift[2] = fVdriftEnd;
206 fXmaxVdrift[3] = fRunEnd;
207
208 fMapDz=new TMap();
c74c5f6c 209
d3ce44cb 210 fCosmiMatchingHisto[0]=new TH1F("Cosmics matching","p0-all" ,100,-10*0.5356 ,10*0.5356 );
211 fCosmiMatchingHisto[1]=new TH1F("Cosmics matching","p1-all" ,100,-10*4.541 ,10*4.541 );
212 fCosmiMatchingHisto[2]=new TH1F("Cosmics matching","p2-all" ,100,-10*0.01134 ,10*0.01134 );
213 fCosmiMatchingHisto[3]=new TH1F("Cosmics matching","p3-all" ,100,-10*0.004644,10*0.004644);
214 fCosmiMatchingHisto[4]=new TH1F("Cosmics matching","p4-all" ,100,-10*0.03773 ,10*0.03773 );
215 fCosmiMatchingHisto[5]=new TH1F("Cosmics matching","p0-isPair",100,-10*0.5356 ,10*0.5356 );
216 fCosmiMatchingHisto[6]=new TH1F("Cosmics matching","p1-isPair",100,-10*4.541 ,10*4.541 );
217 fCosmiMatchingHisto[7]=new TH1F("Cosmics matching","p2-isPair",100,-10*0.01134 ,10*0.01134 );
218 fCosmiMatchingHisto[8]=new TH1F("Cosmics matching","p3-isPair",100,-10*0.004644,10*0.004644);
219 fCosmiMatchingHisto[9]=new TH1F("Cosmics matching","p4-isPair",100,-10*0.03773 ,10*0.03773 );
220// Char_t nameHisto[3]={'p','0','\n'};
221// for (Int_t i=0;i<10;i++){
222// fCosmiMatchingHisto[i]=new TH1F("Cosmics matching",nameHisto,8192,0,0);
223// nameHisto[1]++;
224// if(i==4) nameHisto[1]='0';
225// }
74235403 226}
c74c5f6c 227
228AliTPCcalibTime::~AliTPCcalibTime(){
229 //
2be25cc9 230 // Destructor
c74c5f6c 231 //
d3ce44cb 232 for(Int_t i=0;i<3;i++){
233 if(fHistVdriftLaserA[i]){
234 delete fHistVdriftLaserA[i];
235 fHistVdriftLaserA[i]=NULL;
236 }
237 if(fHistVdriftLaserC[i]){
238 delete fHistVdriftLaserC[i];
239 fHistVdriftLaserC[i]=NULL;
240 }
241 }
242 if(fMapDz){
243 fMapDz->SetOwner();
244 fMapDz->Delete();
245 delete fMapDz;
246 fMapDz=NULL;
2be25cc9 247 }
d3ce44cb 248 for(Int_t i=0;i<5;i++){
249 if(fCosmiMatchingHisto[i]){
250 delete fCosmiMatchingHisto[i];
251 fCosmiMatchingHisto[i]=NULL;
252 }
74235403 253 }
c74c5f6c 254}
2be25cc9 255
da6c0bc9 256void AliTPCcalibTime::ResetCurrent(){
257 //
258 // reset current values
259 //
260 fDz=0; // current delta z
da6c0bc9 261}
262
2be25cc9 263Bool_t AliTPCcalibTime::IsLaser(AliESDEvent *event){
264 return ((event->GetFiredTriggerClasses()).Contains("0LSR")==1);
265}
c74c5f6c 266
2be25cc9 267void AliTPCcalibTime::Process(AliESDEvent *event){
268 if(!event) return;
269 if (event->GetNumberOfTracks()<2) return;
270 ResetCurrent();
271
272// if(IsLaser(event))
273 ProcessLaser (event);
274// else
275 ProcessCosmic(event);
276}
277
278void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){
c74c5f6c 279 //
2be25cc9 280 // Fit drift velocity using laser
281 //
282 // 0. cuts
283 const Int_t kMinTracks = 20; // minimal number of laser tracks
284 const Int_t kMinTracksSide = 10; // minimal number of tracks per side
285 const Float_t kMaxRMS = 0.1; // maximal RMS of tracks
286 const Float_t kMaxDeltaZ = 3.; // maximal deltaZ A-C side
287 const Float_t kMaxDeltaV = 0.01; // maximal deltaV A-C side
288 const Float_t kMaxDeltaY = 2.; // maximal deltaY A-C side
289 /*
290 TCut cutRMS("sqrt(laserA.fElements[4])<0.1&&sqrt(laserC.fElements[4])<0.1");
291 TCut cutZ("abs(laserA.fElements[0]-laserC.fElements[0])<3");
292 TCut cutV("abs(laserA.fElements[1]-laserC.fElements[1])<0.01");
293 TCut cutY("abs(laserA.fElements[2]-laserC.fElements[2])<2");
294 TCut cutAll = cutRMS+cutZ+cutV+cutY;
295 */
296 if (event->GetNumberOfTracks()<kMinTracks) return;
c74c5f6c 297 //
2be25cc9 298 if(!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
299 fLaser->Process(event);
300 if (fLaser->GetNtracks()<kMinTracks) return; // small amount of tracks cut
301 if (fLaser->fFitAside->GetNrows()==0) return; // no fit A side
302 if (fLaser->fFitCside->GetNrows()==0) return; // no fit C side
c74c5f6c 303 //
2be25cc9 304 // debug streamer - activate stream level
305 // Use it for tuning of the cuts
da6c0bc9 306 //
2be25cc9 307 if (fStreamLevel>0){
308 printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
309
da6c0bc9 310 TTreeSRedirector *cstream = GetDebugStreamer();
311 if (cstream){
312 TTimeStamp tstamp(fTime);
2be25cc9 313 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
da6c0bc9 314 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
315 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
316 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
317 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
318 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
31aa7c5c 319 TVectorD vecGoofie(20);
320 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
2be25cc9 321 if (goofieArray){
31aa7c5c 322 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
323 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
324 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
325 }
da6c0bc9 326 }
2be25cc9 327 (*cstream)<<"laserInfo"<<
da6c0bc9 328 "run="<<fRun<< // run number
329 "event="<<fEvent<< // event number
330 "time="<<fTime<< // time stamp of event
331 "trigger="<<fTrigger<< // trigger
332 "mag="<<fMagF<< // magnetic field
333 // Environment values
334 "press0="<<valuePressure0<<
335 "press1="<<valuePressure1<<
336 "pt0="<<ptrelative0<<
337 "pt1="<<ptrelative1<<
338 "temp0="<<temp0<<
339 "temp1="<<temp1<<
31aa7c5c 340 "vecGoofie.=<<"<<&vecGoofie<<
da6c0bc9 341 //laser
2be25cc9 342 "laserA.="<<fLaser->fFitAside<<
343 "laserC.="<<fLaser->fFitCside<<
344 "laserAC.="<<fLaser->fFitACside<<
345 "trigger="<<event->GetFiredTriggerClasses()<<
da6c0bc9 346 "\n";
347 }
348 }
2be25cc9 349 //
350 // Apply custs
351 //
352 if ((*fLaser->fFitAside)[3] <kMinTracksSide) return; // enough tracks A side
353 if ((*fLaser->fFitCside)[3]<kMinTracksSide) return; // enough tracks C side
354 //
355 if (TMath::Abs((*fLaser->fFitAside)[0]-(*fLaser->fFitCside)[0])>kMaxDeltaZ) return;
356 if (TMath::Abs((*fLaser->fFitAside)[2]-(*fLaser->fFitCside)[2])>kMaxDeltaY) return;
357 if (TMath::Abs((*fLaser->fFitAside)[1]-(*fLaser->fFitCside)[1])>kMaxDeltaV) return;
358 if (TMath::Sqrt(TMath::Abs((*fLaser->fFitAside)[4]))>kMaxRMS) return;
359 if (TMath::Sqrt(TMath::Abs((*fLaser->fFitCside)[4]))>kMaxRMS) return;
360 //
361 // fill histos
362 //
363 TVectorD vdriftA(5), vdriftC(5),vdriftAC(5);
364 vdriftA=*(fLaser->fFitAside);
365 vdriftC=*(fLaser->fFitCside);
366 vdriftAC=*(fLaser->fFitACside);
367 Int_t npointsA=0, npointsC=0;
368 Float_t chi2A=0, chi2C=0;
369 npointsA= TMath::Nint(vdriftA[3]);
370 chi2A= vdriftA[4];
371 npointsC= TMath::Nint(vdriftC[3]);
372 chi2C= vdriftC[4];
2be25cc9 373
2be25cc9 374 TTimeStamp tstamp(fTime);
375 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
376 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
74235403 377
2be25cc9 378 Double_t vecVdriftLaserA[4]={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitAside))[1])-1,event->GetRunNumber()};
379 Double_t vecVdriftLaserC[4]={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitCside))[1])-1,event->GetRunNumber()};
380
381 for (Int_t i=0;i<3;i++){
382 if (i==0){ //z0 shift
383 vecVdriftLaserA[3]=(*(fLaser->fFitAside))[0]/250.;
384 vecVdriftLaserA[3]=(*(fLaser->fFitCside))[0]/250.;
385 }
386 if (i==1){ //vdrel shift
387 vecVdriftLaserA[3]=1./(*(fLaser->fFitAside))[1]-1.;
388 vecVdriftLaserA[3]=1./(*(fLaser->fFitCside))[1]-1.;
389 }
390 if (i==2){ //gy shift - full gy - full drift
391 vecVdriftLaserA[3]=(*(fLaser->fFitAside))[2]/250.;
392 vecVdriftLaserA[3]=(*(fLaser->fFitCside))[2]/250.;
393 }
394 fHistVdriftLaserA[i]->Fill(vecVdriftLaserA);
395 fHistVdriftLaserC[i]->Fill(vecVdriftLaserC);
396 }
2be25cc9 397}
c74c5f6c 398
2be25cc9 399void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){
c74c5f6c 400 if (!event) {
401 Printf("ERROR: ESD not available");
402 return;
403 }
404 if (event->GetTimeStamp() == 0 ) {
405 Printf("no time stamp!");
406 return;
407 }
74235403 408
2be25cc9 409 //fd
c74c5f6c 410 // Find cosmic pairs
411 //
412 // Track0 is choosen in upper TPC part
413 // Track1 is choosen in lower TPC part
414 //
415 Int_t ntracks=event->GetNumberOfTracks();
416 if (ntracks==0) return;
74235403 417 if (ntracks > fCutTracks) return;
418
c74c5f6c 419 if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
420 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
421
422 TObjArray tpcSeeds(ntracks);
423 Double_t vtxx[3]={0,0,0};
424 Double_t svtxx[3]={0.000001,0.000001,100.};
425 AliESDVertex vtx(vtxx,svtxx);
426 //
427 // track loop
428 //
429 for (Int_t i=0;i<ntracks;++i) {
74235403 430 AliESDtrack *track = event->GetTrack(i);
431
432 const AliExternalTrackParam * trackIn = track->GetInnerParam();
433 const AliExternalTrackParam * trackOut = track->GetOuterParam();
434 if (!trackIn) continue;
435 if (!trackOut) continue;
436
437 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
438 TObject *calibObject;
439 AliTPCseed *seed = 0;
440 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
441 if (seed) tpcSeeds.AddAt(seed,i);
c74c5f6c 442 }
c74c5f6c 443 if (ntracks<2) return;
444 //
445 // Find pairs
446 //
d3ce44cb 447
c74c5f6c 448 for (Int_t i=0;i<ntracks;++i) {
74235403 449 AliESDtrack *track0 = event->GetTrack(i);
c74c5f6c 450 // track0 - choosen upper part
451 if (!track0) continue;
452 if (!track0->GetOuterParam()) continue;
453 if (track0->GetOuterParam()->GetAlpha()<0) continue;
454 Double_t d1[3];
455 track0->GetDirection(d1);
456 for (Int_t j=0;j<ntracks;++j) {
457 if (i==j) continue;
74235403 458 AliESDtrack *track1 = event->GetTrack(j);
459 //track 1 lower part
460 if (!track1) continue;
461 if (!track1->GetOuterParam()) continue;
462 if (track1->GetOuterParam()->GetAlpha()>0) continue;
463 //
464 Double_t d2[3];
465 track1->GetDirection(d2);
466
467 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
468 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
469 if (! seed0) continue;
470 if (! seed1) continue;
471 Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
472 Float_t dist0 = track0->GetLinearD(0,0);
473 Float_t dist1 = track1->GetLinearD(0,0);
474 //
475 // conservative cuts - convergence to be guarantied
476 // applying before track propagation
477 if (TMath::Abs(dist0+dist1)>fCutMaxD) continue; // distance to the 0,0
478 if (dir>fCutMinDir) continue; // direction vector product
479 Float_t bz = AliTracker::GetBz();
480 Float_t dvertex0[2]; //distance to 0,0
481 Float_t dvertex1[2]; //distance to 0,0
482 track0->GetDZ(0,0,0,bz,dvertex0);
483 track1->GetDZ(0,0,0,bz,dvertex1);
484 if (TMath::Abs(dvertex0[1])>250) continue;
485 if (TMath::Abs(dvertex1[1])>250) continue;
486 //
487 //
488 //
489 Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
490 AliExternalTrackParam param0(*track0);
491 AliExternalTrackParam param1(*track1);
492 //
493 // Propagate using Magnetic field and correct fo material budget
494 //
495 AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
496 AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
497 //
498 // Propagate rest to the 0,0 DCA - z should be ignored
499 //
500 //Bool_t b0 = ;
501 param0.PropagateToDCA(&vtx,bz,1000);
502 //Bool_t b1 =
503 param1.PropagateToDCA(&vtx,bz,1000);
504 //
505 param0.GetDZ(0,0,0,bz,dvertex0);
506 param1.GetDZ(0,0,0,bz,dvertex1);
507 //
508 Double_t xyz0[3];//,pxyz0[3];
509 Double_t xyz1[3];//,pxyz1[3];
510 param0.GetXYZ(xyz0);
511 param1.GetXYZ(xyz1);
512 Bool_t isPair = IsPair(&param0,&param1);
d3ce44cb 513 Bool_t isCross = IsCross(track0, track1);
514
515// Double_t z0 = track0->GetOuterParam()->GetZ();
516// Double_t z1 = track1->GetOuterParam()->GetZ();
74235403 517
d3ce44cb 518// Double_t z0inner = track0->GetInnerParam()->GetZ();
519// Double_t z1inner = track1->GetInnerParam()->GetZ();
74235403 520
d3ce44cb 521 if (isPair && isCross){
74235403 522 if (track0->GetTPCNcls() > 80) {
523 fDz = param0.GetZ() - param1.GetZ();
d3ce44cb 524 if(track0->GetOuterParam()->GetZ()<0) fDz=-fDz;
74235403 525 TTimeStamp tstamp(fTime);
526 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
527 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
528 Double_t vecVdrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
529 THnSparse* curHist=0;
530
531 curHist=(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
532 if(!curHist){
533 curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
534 fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
535 }
536 curHist->Fill(vecVdrift);
537
538 curHist=(THnSparseF*)(fMapDz->GetValue("all"));
539 if(!curHist){
540 curHist=new THnSparseF("all","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
541 fMapDz->Add(new TObjString("all"),curHist);
542 }
543 curHist->Fill(vecVdrift);
544 }
545 }
d3ce44cb 546 TTreeSRedirector *cstream = GetDebugStreamer();
547 if (fStreamLevel>0){
548 if (cstream){
549 (*cstream)<<"trackInfo"<<
550 "tr0.="<<track0<<
551 "tr1.="<<track1<<
552 "p0.="<<&param0<<
553 "p1.="<<&param1<<
554 "isPair="<<isPair<<
555 "isCross="<<isCross<<
556 "fDz="<<fDz<<
557 "fRun="<<fRun<<
558 "fTime="<<fTime<<
559 "\n";
560 }
561 }
c74c5f6c 562 } // end 2nd order loop
563 } // end 1st order loop
74235403 564
2be25cc9 565 if (fStreamLevel>0){
566 TTreeSRedirector *cstream = GetDebugStreamer();
567 if (cstream){
568 TTimeStamp tstamp(fTime);
569 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
570 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
571 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
572 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
573 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
574 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
575 TVectorD vecGoofie(20);
576 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
577 if (goofieArray){
578 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
579 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
580 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
581 }
582 }
583 (*cstream)<<"timeInfo"<<
584 "run="<<fRun<< // run number
585 "event="<<fEvent<< // event number
586 "time="<<fTime<< // time stamp of event
587 "trigger="<<fTrigger<< // trigger
588 "mag="<<fMagF<< // magnetic field
589 // Environment values
590 "press0="<<valuePressure0<<
591 "press1="<<valuePressure1<<
592 "pt0="<<ptrelative0<<
593 "pt1="<<ptrelative1<<
594 "temp0="<<temp0<<
595 "temp1="<<temp1<<
596 "vecGoofie.=<<"<<&vecGoofie<<
597 //
598 // accumulated values
599 //
600 "fDz="<<fDz<< //! current delta z
2be25cc9 601 "trigger="<<event->GetFiredTriggerClasses()<<
602 "\n";
603 }
604 }
605 printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
c74c5f6c 606}
607
74235403 608void AliTPCcalibTime::Analyze(){}
609
610THnSparse* AliTPCcalibTime::GetHistoDrift(TObjString* name){
611 return (THnSparseF*)(fMapDz->GetValue(name));
612}
74235403 613THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name){
614 TObjString* objName=new TObjString(name);
615 THnSparse* histoDrift=0;
616 if(objName){
617 histoDrift=GetHistoDrift(objName);
618 delete objName;
619 objName=0;
620 }
621 return histoDrift;
622}
623
624TMap* AliTPCcalibTime::GetHistoDrift(){
625 return fMapDz;
626}
627
628TGraphErrors* AliTPCcalibTime::GetGraphDrift(TObjString* name){
629 THnSparse* histoDrift=GetHistoDrift(name);
630 TGraphErrors* graphDrift=0;
631 if(histoDrift) graphDrift=FitSlices(histoDrift,2,0,400,100,0.05,0.95, kTRUE);
632 return graphDrift;
633}
634
635TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){
636 TObjString* objName=new TObjString(name);
637 TGraphErrors* graphDrift=0;
638 if(objName){
639 graphDrift=GetGraphDrift(objName);
640 delete objName;
641 objName=0;
642 }
643 return graphDrift;
644}
645
646TMap* AliTPCcalibTime::GetGraphDrift(){
647 TMap* mapGraphDrift=new TMap();
648 TIterator* iterator=fMapDz->MakeIterator();
649 iterator->Reset();
650 TPair* addPair=0;
651 while((addPair=(TPair*)(fMapDz->FindObject(iterator->Next())))) mapGraphDrift->Add((TObjString*)addPair->Key(), GetGraphDrift((TObjString*)addPair->Key()));
652 return mapGraphDrift;
653}
654
655TGraph* AliTPCcalibTime::GetFitDrift(TObjString* name){
d3ce44cb 656 TGraph* graphDrift=GetGraphDrift(name);
74235403 657 TGraph* fitDrift=0;
658 if(graphDrift && graphDrift->GetN()){
659 AliSplineFit fit;
660 fit.SetGraph(graphDrift);
661 fit.SetMinPoints(graphDrift->GetN()+1);
662 fit.InitKnots(graphDrift,2,0,0.001);
663 fit.SplineFit(0);
664 fitDrift=fit.MakeGraph(graphDrift->GetX()[0],graphDrift->GetX()[graphDrift->GetN()-1],50000,0);
665 delete graphDrift;
666 graphDrift=0;
667 }
668 return fitDrift;
669}
670
671TGraph* AliTPCcalibTime::GetFitDrift(const char* name){
672 TObjString* objName=new TObjString(name);
673 TGraph* fitDrift=0;
674 if(objName){
675 fitDrift=GetFitDrift(objName);
676 delete objName;
677 objName=0;
678 }
679 return fitDrift;
680}
681
682TMap* AliTPCcalibTime::GetFitDrift(){
683 TMap* mapFitDrift=new TMap();
684 TIterator* iterator = fMapDz->MakeIterator();
685 iterator->Reset();
686 TPair* addPair=0;
687 while((addPair=(TPair*)(fMapDz->FindObject(iterator->Next())))) mapFitDrift->Add((TObjString*)addPair->Key(), GetFitDrift((TObjString*)addPair->Key()));
688 return mapFitDrift;
c74c5f6c 689}
c74c5f6c 690Long64_t AliTPCcalibTime::Merge(TCollection *li) {
691
692 TIterator* iter = li->MakeIterator();
693 AliTPCcalibTime* cal = 0;
694
695 while ((cal = (AliTPCcalibTime*)iter->Next())) {
696 if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
697 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
698 return -1;
699 }
2be25cc9 700 for (Int_t imeas=0; imeas<3; imeas++){
701 if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
702 fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
703 fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
704 }
705 }
74235403 706 TMap * addMap=cal->GetHistoDrift();
2be25cc9 707 if(!addMap) return 0;
708 TIterator* iterator = addMap->MakeIterator();
709 iterator->Reset();
710 TPair* addPair=0;
711 while((addPair=(TPair *)(addMap->FindObject(iterator->Next())))){
712 THnSparse* addHist=dynamic_cast<THnSparseF*>(addPair->Value());
713 if (!addHist) continue;
714 addHist->Print();
715 THnSparse* localHist=dynamic_cast<THnSparseF*>(fMapDz->GetValue(addHist->GetName()));
716 if(!localHist){
717 localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
718 fMapDz->Add(new TObjString(addHist->GetName()),localHist);
719 }
720 localHist->Add(addHist);
721 }
d3ce44cb 722 for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
c74c5f6c 723 }
c74c5f6c 724 return 0;
c74c5f6c 725}
726
c74c5f6c 727Bool_t AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
c74c5f6c 728 /*
729 // 0. Same direction - OPOSITE - cutDir +cutT
730 TCut cutDir("cutDir","dir<-0.99")
731 // 1.
732 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
733 //
734 // 2. The same rphi
735 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
736 //
c74c5f6c 737 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
738 // 1/Pt diff cut
739 */
740 const Double_t *p0 = tr0->GetParameter();
741 const Double_t *p1 = tr1->GetParameter();
74235403 742 fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
743 fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
d3ce44cb 744 fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
74235403 745 fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
746 fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
747
c74c5f6c 748 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
749 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
da6c0bc9 750 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
c74c5f6c 751 Double_t d0[3], d1[3];
752 tr0->GetDirection(d0);
753 tr1->GetDirection(d1);
754 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
74235403 755
d3ce44cb 756 fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
757 fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
758 fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
759 fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
760 fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
761
c74c5f6c 762 return kTRUE;
763}
d3ce44cb 764Bool_t AliTPCcalibTime::IsCross(AliESDtrack *tr0, AliESDtrack *tr1){
765 return tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0 && tr0->GetInnerParam()->GetZ()*tr1->GetInnerParam()->GetZ()<0 && tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0 && tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0;
766}
767
768/*
769chainDrift->Draw("p0.fP[0]+p1.fP[0]","isPair");
770 mean ~-0.02 ~-0.03913
771 RMS ~ 0.5 ~ 0.5356 --> 3 (fCutMaxD)
772
773chainDrift->Draw("p0.fP[1]-p1.fP[1]","isPair");
774 mean ~ 0.1855
775 RMS ~ 4.541 -->25 (fCutMaxDz)
776
777chainDrift->Draw("p1.fAlpha-p0.fAlpha+pi","isPair");
778//chainDrift->Draw("p1.fAlpha+p0.fAlpha","isPair");
779//chainDrift->Draw("p1.fP[2]-p0.fP[2]+pi","isPair");
780//chainDrift->Draw("p1.fP[2]+p0.fP[2]","isPair");
781 mean ~ 0 ~ 0.001898
782 RMS ~ 0.009 ~ 0.01134 --> 0.06
783
784chainDrift->Draw("p0.fP[3]+p1.fP[3]","isPair");
785 mean ~ 0.0013 ~ 0.001539
786 RMS ~ 0.003 ~ 0.004644 --> 0.03 (fCutTheta)
787
788chainDrift->Draw("p1.fP[4]+p0.fP[4]>>his(100,-0.2,0.2)","isPair")
789 mean ~ 0.012 ~-0.0009729
790 RMS ~ 0.036 ~ 0.03773 --> 0.2
791*/
31aa7c5c 792