]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibTime.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Stronger cuts before filling histogram
[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
dde68d36 45 TFile f("CalibObjectsTrain1.root");
46 AliTPCcalibTime *calib = (AliTPCcalibTime *)f->Get("calibTime");
47 calib->GetHistoDrift("all")->Projection(2,0)->Draw()
48 calib->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();
abae1b29 57 TChain * chainLaser = tool.MakeChain("time.txt.Good","laserInfo",0,10200);
dde68d36 58 chainLaser->Lookup();
c74c5f6c 59*/
60
c74c5f6c 61#include "Riostream.h"
62#include "TChain.h"
63#include "TTree.h"
64#include "TH1F.h"
65#include "TH2F.h"
66#include "TH3F.h"
67#include "THnSparse.h"
68#include "TList.h"
69#include "TMath.h"
70#include "TCanvas.h"
71#include "TFile.h"
72#include "TF1.h"
73#include "TVectorD.h"
74#include "TProfile.h"
75#include "TGraphErrors.h"
76#include "TCanvas.h"
77
78#include "AliTPCclusterMI.h"
79#include "AliTPCseed.h"
80#include "AliESDVertex.h"
81#include "AliESDEvent.h"
82#include "AliESDfriend.h"
83#include "AliESDInputHandler.h"
84#include "AliAnalysisManager.h"
85
86#include "AliTracker.h"
f7a1cc68 87#include "AliMagF.h"
c74c5f6c 88#include "AliTPCCalROC.h"
89
90#include "AliLog.h"
91
92#include "AliTPCcalibTime.h"
93
94#include "TTreeStream.h"
95#include "AliTPCTracklet.h"
da6c0bc9 96#include "TTimeStamp.h"
97#include "AliTPCcalibDB.h"
98#include "AliTPCcalibLaser.h"
31aa7c5c 99#include "AliDCSSensorArray.h"
100#include "AliDCSSensor.h"
c74c5f6c 101
102ClassImp(AliTPCcalibTime)
103
104
105AliTPCcalibTime::AliTPCcalibTime()
da6c0bc9 106 :AliTPCcalibBase(),
da6c0bc9 107 fLaser(0), // pointer to laser calibration
108 fDz(0), // current delta z
d3ce44cb 109 fCutMaxD(3), // maximal distance in rfi ditection
110 fCutMaxDz(25), // maximal distance in rfi ditection
c74c5f6c 111 fCutTheta(0.03), // maximal distan theta
2be25cc9 112 fCutMinDir(-0.99), // direction vector products
74235403 113 fCutTracks(10),
dde68d36 114 fArrayDz(0), //NEW! Tmap of V drifts for different triggers
2be25cc9 115 fTimeBins(0),
116 fTimeStart(0),
117 fTimeEnd(0),
118 fPtBins(0),
119 fPtStart(0),
120 fPtEnd(0),
121 fVdriftBins(0),
122 fVdriftStart(0),
123 fVdriftEnd(0),
124 fRunBins(0),
125 fRunStart(0),
126 fRunEnd(0)
127// fBinsVdrift(fTimeBins,fPtBins,fVdriftBins),
128// fXminVdrift(fTimeStart,fPtStart,fVdriftStart),
129// fXmaxVdrift(fTimeEnd,fPtEnd,fVdriftEnd)
c74c5f6c 130{
131 AliInfo("Default Constructor");
2be25cc9 132 for (Int_t i=0;i<3;i++) {
133 fHistVdriftLaserA[i]=0;
134 fHistVdriftLaserC[i]=0;
135 }
d3ce44cb 136 for (Int_t i=0;i<10;i++) {
74235403 137 fCosmiMatchingHisto[i]=0;
138 }
c74c5f6c 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
dde68d36 146 fCutMaxDz(40), // maximal distance in rfi ditection
d3ce44cb 147 fCutTheta(5*0.004644),// maximal distan theta
148 fCutMinDir(-0.99), // direction vector products
74235403 149 fCutTracks(10),
dde68d36 150 fArrayDz(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{
c74c5f6c 164 SetName(name);
165 SetTitle(title);
2be25cc9 166 for (Int_t i=0;i<3;i++) {
167 fHistVdriftLaserA[i]=0;
168 fHistVdriftLaserC[i]=0;
169 }
c74c5f6c 170
171 AliInfo("Non Default Constructor");
74235403 172 fTimeBins =(EndTime-StartTime)/deltaIntegrationTimeVdrift;
2be25cc9 173 fTimeStart =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi();
174 fTimeEnd =EndTime; //(((TObjString*)(mapGRP->GetValue("fAliceStopTime")))->GetString()).Atoi();
dde68d36 175 fPtBins = 400;
2be25cc9 176 fPtStart = -0.04;
177 fPtEnd = 0.04;
dde68d36 178 fVdriftBins = 500;
179 fVdriftStart= -0.1;
180 fVdriftEnd = 0.1;
181 fRunBins = 100001;
182 fRunStart = -1.5;
d3ce44cb 183 fRunEnd = 99999.5;
2be25cc9 184
185 Int_t binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins };
186 Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart , fRunStart};
187 Double_t xmaxVdriftLaser[4] = {fTimeEnd , fPtEnd , fVdriftEnd , fRunEnd };
dde68d36 188 TString axisTitle[4]={
189 "T",
190 "#delta_{P/T}",
191 "value",
192 "run"
193 };
194 TString histoName[3]={
195 "Loffset",
196 "Lcorr",
197 "Lgy"
198 };
2be25cc9 199
dde68d36 200
2be25cc9 201 for (Int_t i=0;i<3;i++) {
202 fHistVdriftLaserA[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
203 fHistVdriftLaserC[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
dde68d36 204 fHistVdriftLaserA[i]->SetName(histoName[i]);
205 fHistVdriftLaserC[i]->SetName(histoName[i]);
206 for (Int_t iaxis=0; iaxis<4;iaxis++){
207 fHistVdriftLaserA[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
208 fHistVdriftLaserC[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
209 }
2be25cc9 210 }
211 fBinsVdrift[0] = fTimeBins;
212 fBinsVdrift[1] = fPtBins;
213 fBinsVdrift[2] = fVdriftBins;
214 fBinsVdrift[3] = fRunBins;
215 fXminVdrift[0] = fTimeStart;
216 fXminVdrift[1] = fPtStart;
217 fXminVdrift[2] = fVdriftStart;
218 fXminVdrift[3] = fRunStart;
219 fXmaxVdrift[0] = fTimeEnd;
220 fXmaxVdrift[1] = fPtEnd;
221 fXmaxVdrift[2] = fVdriftEnd;
222 fXmaxVdrift[3] = fRunEnd;
223
dde68d36 224 fArrayDz=new TObjArray();
abae1b29 225 // fArrayDz->AddLast(fHistVdriftLaserA[0]);
226// fArrayDz->AddLast(fHistVdriftLaserA[1]);
227// fArrayDz->AddLast(fHistVdriftLaserA[2]);
228// fArrayDz->AddLast(fHistVdriftLaserC[0]);
229// fArrayDz->AddLast(fHistVdriftLaserC[1]);
230// fArrayDz->AddLast(fHistVdriftLaserC[2]);
c74c5f6c 231
d3ce44cb 232 fCosmiMatchingHisto[0]=new TH1F("Cosmics matching","p0-all" ,100,-10*0.5356 ,10*0.5356 );
233 fCosmiMatchingHisto[1]=new TH1F("Cosmics matching","p1-all" ,100,-10*4.541 ,10*4.541 );
234 fCosmiMatchingHisto[2]=new TH1F("Cosmics matching","p2-all" ,100,-10*0.01134 ,10*0.01134 );
235 fCosmiMatchingHisto[3]=new TH1F("Cosmics matching","p3-all" ,100,-10*0.004644,10*0.004644);
236 fCosmiMatchingHisto[4]=new TH1F("Cosmics matching","p4-all" ,100,-10*0.03773 ,10*0.03773 );
237 fCosmiMatchingHisto[5]=new TH1F("Cosmics matching","p0-isPair",100,-10*0.5356 ,10*0.5356 );
238 fCosmiMatchingHisto[6]=new TH1F("Cosmics matching","p1-isPair",100,-10*4.541 ,10*4.541 );
239 fCosmiMatchingHisto[7]=new TH1F("Cosmics matching","p2-isPair",100,-10*0.01134 ,10*0.01134 );
240 fCosmiMatchingHisto[8]=new TH1F("Cosmics matching","p3-isPair",100,-10*0.004644,10*0.004644);
241 fCosmiMatchingHisto[9]=new TH1F("Cosmics matching","p4-isPair",100,-10*0.03773 ,10*0.03773 );
242// Char_t nameHisto[3]={'p','0','\n'};
243// for (Int_t i=0;i<10;i++){
244// fCosmiMatchingHisto[i]=new TH1F("Cosmics matching",nameHisto,8192,0,0);
245// nameHisto[1]++;
246// if(i==4) nameHisto[1]='0';
247// }
74235403 248}
c74c5f6c 249
250AliTPCcalibTime::~AliTPCcalibTime(){
251 //
2be25cc9 252 // Destructor
c74c5f6c 253 //
d3ce44cb 254 for(Int_t i=0;i<3;i++){
255 if(fHistVdriftLaserA[i]){
256 delete fHistVdriftLaserA[i];
257 fHistVdriftLaserA[i]=NULL;
258 }
259 if(fHistVdriftLaserC[i]){
260 delete fHistVdriftLaserC[i];
261 fHistVdriftLaserC[i]=NULL;
262 }
263 }
dde68d36 264 if(fArrayDz){
265 fArrayDz->SetOwner();
266 fArrayDz->Delete();
267 delete fArrayDz;
268 fArrayDz=NULL;
2be25cc9 269 }
d3ce44cb 270 for(Int_t i=0;i<5;i++){
271 if(fCosmiMatchingHisto[i]){
272 delete fCosmiMatchingHisto[i];
273 fCosmiMatchingHisto[i]=NULL;
274 }
74235403 275 }
c74c5f6c 276}
2be25cc9 277
dde68d36 278Bool_t AliTPCcalibTime::IsLaser(AliESDEvent */*event*/){
279 return kTRUE; //More accurate creteria to be added
da6c0bc9 280}
dde68d36 281Bool_t AliTPCcalibTime::IsCosmics(AliESDEvent */*event*/){
282 return kTRUE; //More accurate creteria to be added
283}
284Bool_t AliTPCcalibTime::IsBeam(AliESDEvent */*event*/){
285 return kTRUE; //More accurate creteria to be added
286}
287void AliTPCcalibTime::ResetCurrent(){
288 fDz=0; //Reset current dz
2be25cc9 289}
2be25cc9 290void AliTPCcalibTime::Process(AliESDEvent *event){
291 if(!event) return;
292 if (event->GetNumberOfTracks()<2) return;
293 ResetCurrent();
dde68d36 294 if(IsLaser (event)) ProcessLaser (event);
295 if(IsCosmics(event)) ProcessCosmic(event);
296 if(IsBeam (event)) ProcessBeam (event);
2be25cc9 297}
298
299void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){
c74c5f6c 300 //
2be25cc9 301 // Fit drift velocity using laser
302 //
303 // 0. cuts
dde68d36 304 const Int_t kMinTracks = 40; // minimal number of laser tracks
305 const Int_t kMinTracksSide = 20; // minimal number of tracks per side
306 const Float_t kMaxDeltaZ = 30.; // maximal trigger delay
307 const Float_t kMaxDeltaV = 0.05; // maximal deltaV
2be25cc9 308 const Float_t kMaxRMS = 0.1; // maximal RMS of tracks
dde68d36 309 //
2be25cc9 310 /*
311 TCut cutRMS("sqrt(laserA.fElements[4])<0.1&&sqrt(laserC.fElements[4])<0.1");
312 TCut cutZ("abs(laserA.fElements[0]-laserC.fElements[0])<3");
313 TCut cutV("abs(laserA.fElements[1]-laserC.fElements[1])<0.01");
314 TCut cutY("abs(laserA.fElements[2]-laserC.fElements[2])<2");
315 TCut cutAll = cutRMS+cutZ+cutV+cutY;
316 */
317 if (event->GetNumberOfTracks()<kMinTracks) return;
c74c5f6c 318 //
2be25cc9 319 if(!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
320 fLaser->Process(event);
321 if (fLaser->GetNtracks()<kMinTracks) return; // small amount of tracks cut
dde68d36 322 if (fLaser->fFitAside->GetNrows()==0 && fLaser->fFitCside->GetNrows()==0) return; // no fit neither a or C side
c74c5f6c 323 //
2be25cc9 324 // debug streamer - activate stream level
325 // Use it for tuning of the cuts
da6c0bc9 326 //
dde68d36 327 // cuts to be applied
328 //
329 Int_t isReject[2]={0,0};
330 //
331 // not enough tracks
332 if (TMath::Abs((*fLaser->fFitAside)[3]) < kMinTracksSide) isReject[0]|=1;
333 if (TMath::Abs((*fLaser->fFitCside)[3]) < kMinTracksSide) isReject[1]|=1;
334 // unreasonable z offset
335 if (TMath::Abs((*fLaser->fFitAside)[0])>kMaxDeltaZ) isReject[0]|=2;
336 if (TMath::Abs((*fLaser->fFitCside)[0])>kMaxDeltaZ) isReject[1]|=2;
337 // unreasonable drift velocity
338 if (TMath::Abs((*fLaser->fFitAside)[1]-1)>kMaxDeltaV) isReject[0]|=4;
339 if (TMath::Abs((*fLaser->fFitCside)[1]-1)>kMaxDeltaV) isReject[1]|=4;
340 // big chi2
341 if (TMath::Sqrt(TMath::Abs((*fLaser->fFitAside)[4]))>kMaxRMS ) isReject[0]|=8;
342 if (TMath::Sqrt(TMath::Abs((*fLaser->fFitCside)[4]))>kMaxRMS ) isReject[1]|=8;
343
344
345
346
2be25cc9 347 if (fStreamLevel>0){
348 printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
349
da6c0bc9 350 TTreeSRedirector *cstream = GetDebugStreamer();
351 if (cstream){
352 TTimeStamp tstamp(fTime);
2be25cc9 353 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
da6c0bc9 354 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
355 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
356 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
357 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
358 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
31aa7c5c 359 TVectorD vecGoofie(20);
360 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
2be25cc9 361 if (goofieArray){
31aa7c5c 362 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
363 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
364 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
365 }
da6c0bc9 366 }
2be25cc9 367 (*cstream)<<"laserInfo"<<
da6c0bc9 368 "run="<<fRun<< // run number
369 "event="<<fEvent<< // event number
370 "time="<<fTime<< // time stamp of event
371 "trigger="<<fTrigger<< // trigger
372 "mag="<<fMagF<< // magnetic field
373 // Environment values
374 "press0="<<valuePressure0<<
375 "press1="<<valuePressure1<<
376 "pt0="<<ptrelative0<<
377 "pt1="<<ptrelative1<<
378 "temp0="<<temp0<<
379 "temp1="<<temp1<<
31aa7c5c 380 "vecGoofie.=<<"<<&vecGoofie<<
da6c0bc9 381 //laser
dde68d36 382 "rejectA="<<isReject[0]<<
383 "rejectC="<<isReject[1]<<
2be25cc9 384 "laserA.="<<fLaser->fFitAside<<
385 "laserC.="<<fLaser->fFitCside<<
386 "laserAC.="<<fLaser->fFitACside<<
387 "trigger="<<event->GetFiredTriggerClasses()<<
da6c0bc9 388 "\n";
389 }
390 }
2be25cc9 391 //
2be25cc9 392 // fill histos
393 //
394 TVectorD vdriftA(5), vdriftC(5),vdriftAC(5);
395 vdriftA=*(fLaser->fFitAside);
396 vdriftC=*(fLaser->fFitCside);
397 vdriftAC=*(fLaser->fFitACside);
398 Int_t npointsA=0, npointsC=0;
399 Float_t chi2A=0, chi2C=0;
400 npointsA= TMath::Nint(vdriftA[3]);
401 chi2A= vdriftA[4];
402 npointsC= TMath::Nint(vdriftC[3]);
403 chi2C= vdriftC[4];
2be25cc9 404
2be25cc9 405 TTimeStamp tstamp(fTime);
406 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
407 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
dde68d36 408 Double_t driftA=0, driftC=0;
409 if (vdriftA[1]>1.-kMaxDeltaV) driftA = 1./vdriftA[1]-1.;
410 if (vdriftC[1]>1.-kMaxDeltaV) driftC = 1./vdriftC[1]-1.;
411 //
412 Double_t vecDriftLaserA[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftA,event->GetRunNumber()};
413 Double_t vecDriftLaserC[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftC,event->GetRunNumber()};
414 // Double_t vecDrift[4] ={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitACside))[1])-1,event->GetRunNumber()};
415
416 for (Int_t icalib=0;icalib<3;icalib++){
417 if (icalib==0){ //z0 shift
418 vecDriftLaserA[2]=vdriftA[0]/250.;
419 vecDriftLaserC[2]=vdriftC[0]/250.;
2be25cc9 420 }
dde68d36 421 if (icalib==1){ //vdrel shift
422 vecDriftLaserA[2]=driftA;
423 vecDriftLaserC[2]=driftC;
2be25cc9 424 }
dde68d36 425 if (icalib==2){ //gy shift - full gy - full drift
426 vecDriftLaserA[2]=vdriftA[2]/250.;
427 vecDriftLaserC[2]=vdriftC[2]/250.;
2be25cc9 428 }
abae1b29 429 if (isReject[0]==0) fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
430 if (isReject[1]==0) fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
2be25cc9 431 }
dde68d36 432
433// THnSparse* curHist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
434// TString shortName=curHist->ClassName();
435// shortName+="_MEAN_DRIFT_LASER_";
436// delete curHist;
437// curHist=NULL;
438// TString name="";
439
440// name=shortName;
441// name+=event->GetFiredTriggerClasses();
442// name.ToUpper();
443// curHist=(THnSparseF*)fArrayDz->FindObject(name);
444// if(!curHist){
445// curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
446// fArrayDz->AddLast(curHist);
447// }
448// curHist->Fill(vecDrift);
449
450// name=shortName;
451// name+="ALL";
452// name.ToUpper();
453// curHist=(THnSparseF*)fArrayDz->FindObject(name);
454// if(!curHist){
455// curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
456// fArrayDz->AddLast(curHist);
457// }
458// curHist->Fill(vecDrift);
2be25cc9 459}
c74c5f6c 460
2be25cc9 461void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){
c74c5f6c 462 if (!event) {
463 Printf("ERROR: ESD not available");
464 return;
465 }
466 if (event->GetTimeStamp() == 0 ) {
467 Printf("no time stamp!");
468 return;
469 }
74235403 470
2be25cc9 471 //fd
c74c5f6c 472 // Find cosmic pairs
473 //
474 // Track0 is choosen in upper TPC part
475 // Track1 is choosen in lower TPC part
476 //
477 Int_t ntracks=event->GetNumberOfTracks();
478 if (ntracks==0) return;
74235403 479 if (ntracks > fCutTracks) return;
480
c74c5f6c 481 if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
482 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
483
484 TObjArray tpcSeeds(ntracks);
485 Double_t vtxx[3]={0,0,0};
486 Double_t svtxx[3]={0.000001,0.000001,100.};
487 AliESDVertex vtx(vtxx,svtxx);
488 //
489 // track loop
490 //
491 for (Int_t i=0;i<ntracks;++i) {
74235403 492 AliESDtrack *track = event->GetTrack(i);
493
494 const AliExternalTrackParam * trackIn = track->GetInnerParam();
495 const AliExternalTrackParam * trackOut = track->GetOuterParam();
496 if (!trackIn) continue;
497 if (!trackOut) continue;
498
499 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
500 TObject *calibObject;
501 AliTPCseed *seed = 0;
502 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
503 if (seed) tpcSeeds.AddAt(seed,i);
c74c5f6c 504 }
c74c5f6c 505 if (ntracks<2) return;
506 //
507 // Find pairs
508 //
d3ce44cb 509
c74c5f6c 510 for (Int_t i=0;i<ntracks;++i) {
74235403 511 AliESDtrack *track0 = event->GetTrack(i);
c74c5f6c 512 // track0 - choosen upper part
513 if (!track0) continue;
514 if (!track0->GetOuterParam()) continue;
515 if (track0->GetOuterParam()->GetAlpha()<0) continue;
516 Double_t d1[3];
517 track0->GetDirection(d1);
518 for (Int_t j=0;j<ntracks;++j) {
519 if (i==j) continue;
74235403 520 AliESDtrack *track1 = event->GetTrack(j);
521 //track 1 lower part
522 if (!track1) continue;
523 if (!track1->GetOuterParam()) continue;
524 if (track1->GetOuterParam()->GetAlpha()>0) continue;
525 //
526 Double_t d2[3];
527 track1->GetDirection(d2);
528
529 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
530 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
531 if (! seed0) continue;
532 if (! seed1) continue;
533 Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
534 Float_t dist0 = track0->GetLinearD(0,0);
535 Float_t dist1 = track1->GetLinearD(0,0);
536 //
537 // conservative cuts - convergence to be guarantied
538 // applying before track propagation
539 if (TMath::Abs(dist0+dist1)>fCutMaxD) continue; // distance to the 0,0
540 if (dir>fCutMinDir) continue; // direction vector product
541 Float_t bz = AliTracker::GetBz();
542 Float_t dvertex0[2]; //distance to 0,0
543 Float_t dvertex1[2]; //distance to 0,0
544 track0->GetDZ(0,0,0,bz,dvertex0);
545 track1->GetDZ(0,0,0,bz,dvertex1);
546 if (TMath::Abs(dvertex0[1])>250) continue;
547 if (TMath::Abs(dvertex1[1])>250) continue;
548 //
549 //
550 //
551 Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
552 AliExternalTrackParam param0(*track0);
553 AliExternalTrackParam param1(*track1);
554 //
555 // Propagate using Magnetic field and correct fo material budget
556 //
557 AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
558 AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
559 //
560 // Propagate rest to the 0,0 DCA - z should be ignored
561 //
562 //Bool_t b0 = ;
563 param0.PropagateToDCA(&vtx,bz,1000);
564 //Bool_t b1 =
565 param1.PropagateToDCA(&vtx,bz,1000);
74235403 566 param0.GetDZ(0,0,0,bz,dvertex0);
567 param1.GetDZ(0,0,0,bz,dvertex1);
dde68d36 568 Double_t xyz0[3];
569 Double_t xyz1[3];
74235403 570 param0.GetXYZ(xyz0);
571 param1.GetXYZ(xyz1);
572 Bool_t isPair = IsPair(&param0,&param1);
d3ce44cb 573 Bool_t isCross = IsCross(track0, track1);
dde68d36 574 Bool_t isSame = IsSame(track0, track1);
d3ce44cb 575
dde68d36 576 THnSparse* hist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
577 TString shortName=hist->ClassName();
578 shortName+="_MEAN_VDRIFT_COSMICS_";
579 delete hist;
580 hist=NULL;
581
582 if(isSame || (isCross && isPair)){
74235403 583 if (track0->GetTPCNcls() > 80) {
584 fDz = param0.GetZ() - param1.GetZ();
d3ce44cb 585 if(track0->GetOuterParam()->GetZ()<0) fDz=-fDz;
74235403 586 TTimeStamp tstamp(fTime);
587 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
588 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
dde68d36 589 Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
590 THnSparse* curHist=NULL;
591 TString name="";
592
593 name=shortName;
594 name+=event->GetFiredTriggerClasses();
595 name.ToUpper();
596 curHist=(THnSparseF*)fArrayDz->FindObject(name);
74235403 597 if(!curHist){
dde68d36 598 curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
599 fArrayDz->AddLast(curHist);
74235403 600 }
dde68d36 601// curHist=(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
602// if(!curHist){
603// curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
604// fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
605// }
606 curHist->Fill(vecDrift);
74235403 607
dde68d36 608 name=shortName;
609 name+="ALL";
610 name.ToUpper();
611 curHist=(THnSparseF*)fArrayDz->FindObject(name);
74235403 612 if(!curHist){
dde68d36 613 curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
614 fArrayDz->AddLast(curHist);
74235403 615 }
dde68d36 616// curHist=(THnSparseF*)(fMapDz->GetValue("all"));
617// if(!curHist){
618// curHist=new THnSparseF("all","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
619// fMapDz->Add(new TObjString("all"),curHist);
620// }
621 curHist->Fill(vecDrift);
74235403 622 }
623 }
d3ce44cb 624 TTreeSRedirector *cstream = GetDebugStreamer();
625 if (fStreamLevel>0){
626 if (cstream){
627 (*cstream)<<"trackInfo"<<
628 "tr0.="<<track0<<
629 "tr1.="<<track1<<
630 "p0.="<<&param0<<
631 "p1.="<<&param1<<
632 "isPair="<<isPair<<
633 "isCross="<<isCross<<
dde68d36 634 "isSame="<<isSame<<
d3ce44cb 635 "fDz="<<fDz<<
636 "fRun="<<fRun<<
637 "fTime="<<fTime<<
638 "\n";
639 }
640 }
c74c5f6c 641 } // end 2nd order loop
642 } // end 1st order loop
74235403 643
2be25cc9 644 if (fStreamLevel>0){
645 TTreeSRedirector *cstream = GetDebugStreamer();
646 if (cstream){
647 TTimeStamp tstamp(fTime);
648 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
649 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
650 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
651 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
652 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
653 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
654 TVectorD vecGoofie(20);
655 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
656 if (goofieArray){
657 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
658 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
659 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
660 }
661 }
662 (*cstream)<<"timeInfo"<<
663 "run="<<fRun<< // run number
664 "event="<<fEvent<< // event number
665 "time="<<fTime<< // time stamp of event
666 "trigger="<<fTrigger<< // trigger
667 "mag="<<fMagF<< // magnetic field
668 // Environment values
669 "press0="<<valuePressure0<<
670 "press1="<<valuePressure1<<
671 "pt0="<<ptrelative0<<
672 "pt1="<<ptrelative1<<
673 "temp0="<<temp0<<
674 "temp1="<<temp1<<
675 "vecGoofie.=<<"<<&vecGoofie<<
676 //
677 // accumulated values
678 //
679 "fDz="<<fDz<< //! current delta z
2be25cc9 680 "trigger="<<event->GetFiredTriggerClasses()<<
681 "\n";
682 }
683 }
684 printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
c74c5f6c 685}
686
dde68d36 687void AliTPCcalibTime::ProcessBeam(AliESDEvent */*event*/){
688}
689
74235403 690void AliTPCcalibTime::Analyze(){}
691
74235403 692THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name){
dde68d36 693 TIterator* iterator = fArrayDz->MakeIterator();
694 iterator->Reset();
695 TString newName=name;
696 newName.ToUpper();
697 THnSparse* newHist=new THnSparseF(newName,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
698 THnSparse* addHist=NULL;
699 while((addHist=(THnSparseF*)iterator->Next())){
700 if(!addHist) continue;
701 TString histName=addHist->GetName();
702 if(!histName.Contains(newName)) continue;
703 addHist->Print();
704 newHist->Add(addHist);
74235403 705 }
dde68d36 706 return newHist;
74235403 707}
708
dde68d36 709TObjArray* AliTPCcalibTime::GetHistoDrift(){
710 return fArrayDz;
74235403 711}
712
713TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){
dde68d36 714 THnSparse* histoDrift=GetHistoDrift(name);
715 TGraphErrors* graphDrift=NULL;
716 if(histoDrift){
717 graphDrift=FitSlices(histoDrift,2,0,400,100,0.05,0.95, kTRUE);
718 TString end=histoDrift->GetName();
719 Int_t pos=end.Index("_");
720 end=end(pos,end.Capacity()-pos);
721 TString graphName=graphDrift->ClassName();
722 graphName+=end;
723 graphName.ToUpper();
724 graphDrift->SetName(graphName);
74235403 725 }
726 return graphDrift;
727}
728
dde68d36 729TObjArray* AliTPCcalibTime::GetGraphDrift(){
730 TObjArray* arrayGraphDrift=new TObjArray();
731 TIterator* iterator=fArrayDz->MakeIterator();
74235403 732 iterator->Reset();
dde68d36 733 THnSparse* addHist=NULL;
734 while((addHist=(THnSparseF*)iterator->Next())) arrayGraphDrift->AddLast(GetGraphDrift(addHist->GetName()));
735 return arrayGraphDrift;
74235403 736}
737
dde68d36 738AliSplineFit* AliTPCcalibTime::GetFitDrift(const char* name){
d3ce44cb 739 TGraph* graphDrift=GetGraphDrift(name);
dde68d36 740 AliSplineFit* fitDrift=NULL;
74235403 741 if(graphDrift && graphDrift->GetN()){
dde68d36 742 fitDrift=new AliSplineFit();
743 fitDrift->SetGraph(graphDrift);
744 fitDrift->SetMinPoints(graphDrift->GetN()+1);
745 fitDrift->InitKnots(graphDrift,2,0,0.001);
746 fitDrift->SplineFit(0);
747 TString end=graphDrift->GetName();
748 Int_t pos=end.Index("_");
749 end=end(pos,end.Capacity()-pos);
750 TString fitName=fitDrift->ClassName();
751 fitName+=end;
752 fitName.ToUpper();
753 //fitDrift->SetName(fitName);
74235403 754 delete graphDrift;
dde68d36 755 graphDrift=NULL;
74235403 756 }
757 return fitDrift;
758}
759
dde68d36 760//TObjArray* AliTPCcalibTime::GetFitDrift(){
761// TObjArray* arrayFitDrift=new TObjArray();
762// TIterator* iterator = fArrayDz->MakeIterator();
763// iterator->Reset();
764// THnSparse* addHist=NULL;
765// while((addHist=(THnSparseF*)iterator->Next())) arrayFitDrift->AddLast(GetFitDrift(addHist->GetName()));
766// return arrayFitDrift;
767//}
74235403 768
c74c5f6c 769Long64_t AliTPCcalibTime::Merge(TCollection *li) {
c74c5f6c 770 TIterator* iter = li->MakeIterator();
771 AliTPCcalibTime* cal = 0;
772
773 while ((cal = (AliTPCcalibTime*)iter->Next())) {
774 if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
775 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
776 return -1;
777 }
2be25cc9 778 for (Int_t imeas=0; imeas<3; imeas++){
779 if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
780 fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
781 fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
782 }
783 }
dde68d36 784 TObjArray* addArray=cal->GetHistoDrift();
785 if(!addArray) return 0;
786 TIterator* iterator = addArray->MakeIterator();
2be25cc9 787 iterator->Reset();
dde68d36 788 THnSparse* addHist=NULL;
789 while((addHist=(THnSparseF*)iterator->Next())){
790 if(!addHist) continue;
2be25cc9 791 addHist->Print();
dde68d36 792 THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName());
2be25cc9 793 if(!localHist){
794 localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
dde68d36 795 fArrayDz->AddLast(localHist);
2be25cc9 796 }
797 localHist->Add(addHist);
798 }
dde68d36 799// TMap * addMap=cal->GetHistoDrift();
800// if(!addMap) return 0;
801// TIterator* iterator = addMap->MakeIterator();
802// iterator->Reset();
803// TPair* addPair=0;
804// while((addPair=(TPair *)(addMap->FindObject(iterator->Next())))){
805// THnSparse* addHist=dynamic_cast<THnSparseF*>(addPair->Value());
806// if (!addHist) continue;
807// addHist->Print();
808// THnSparse* localHist=dynamic_cast<THnSparseF*>(fMapDz->GetValue(addHist->GetName()));
809// if(!localHist){
810// localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
811// fMapDz->Add(new TObjString(addHist->GetName()),localHist);
812// }
813// localHist->Add(addHist);
814// }
d3ce44cb 815 for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
c74c5f6c 816 }
c74c5f6c 817 return 0;
c74c5f6c 818}
819
c74c5f6c 820Bool_t AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
c74c5f6c 821 /*
822 // 0. Same direction - OPOSITE - cutDir +cutT
823 TCut cutDir("cutDir","dir<-0.99")
824 // 1.
825 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
826 //
827 // 2. The same rphi
828 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
829 //
c74c5f6c 830 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
831 // 1/Pt diff cut
832 */
833 const Double_t *p0 = tr0->GetParameter();
834 const Double_t *p1 = tr1->GetParameter();
74235403 835 fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
836 fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
d3ce44cb 837 fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
74235403 838 fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
839 fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
840
c74c5f6c 841 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
842 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
da6c0bc9 843 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
c74c5f6c 844 Double_t d0[3], d1[3];
845 tr0->GetDirection(d0);
846 tr1->GetDirection(d1);
847 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
74235403 848
d3ce44cb 849 fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
850 fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
851 fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
852 fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
853 fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
854
c74c5f6c 855 return kTRUE;
856}
d3ce44cb 857Bool_t AliTPCcalibTime::IsCross(AliESDtrack *tr0, AliESDtrack *tr1){
858 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;
859}
860
dde68d36 861Bool_t AliTPCcalibTime::IsSame(AliESDtrack */*tr0*/, AliESDtrack */*tr1*/){
862 // To be implemented
863 return kFALSE;
864}
865
d3ce44cb 866/*
867chainDrift->Draw("p0.fP[0]+p1.fP[0]","isPair");
868 mean ~-0.02 ~-0.03913
869 RMS ~ 0.5 ~ 0.5356 --> 3 (fCutMaxD)
870
871chainDrift->Draw("p0.fP[1]-p1.fP[1]","isPair");
872 mean ~ 0.1855
873 RMS ~ 4.541 -->25 (fCutMaxDz)
874
875chainDrift->Draw("p1.fAlpha-p0.fAlpha+pi","isPair");
876//chainDrift->Draw("p1.fAlpha+p0.fAlpha","isPair");
877//chainDrift->Draw("p1.fP[2]-p0.fP[2]+pi","isPair");
878//chainDrift->Draw("p1.fP[2]+p0.fP[2]","isPair");
879 mean ~ 0 ~ 0.001898
880 RMS ~ 0.009 ~ 0.01134 --> 0.06
881
882chainDrift->Draw("p0.fP[3]+p1.fP[3]","isPair");
883 mean ~ 0.0013 ~ 0.001539
884 RMS ~ 0.003 ~ 0.004644 --> 0.03 (fCutTheta)
885
886chainDrift->Draw("p1.fP[4]+p0.fP[4]>>his(100,-0.2,0.2)","isPair")
887 mean ~ 0.012 ~-0.0009729
888 RMS ~ 0.036 ~ 0.03773 --> 0.2
889*/
31aa7c5c 890