]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibLaser.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
ProcessDCS() { return kFALSE; } is removed following the discussion at
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibLaser.cxx
CommitLineData
c6914c83 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/*
2b35e8e6 17 //
18 // FUNCTIONALITY:
19 //
20 // 1. The laser track is associated with the mirror
21 // see function FindMirror
22 //
23 // 2. The laser track is accepted for the analysis under certain condition
24 // (see function Accpet laser)
25 //
26 // 3. The drift velocity and jitter is calculated event by event
27 // (see function drift velocity)
28 //
29 //
30 //
1fd56785 31 // To make laser scan the user interaction neccessary
32 //
33 .x ~/UliStyle.C
2b35e8e6 34 gSystem->Load("libANALYSIS");
35 gSystem->Load("libTPCcalib");
36 TFile fcalib("CalibObjects.root");
37 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
38 AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)array->FindObject("laserTPC");
39 laser->DumpMeanInfo(-0.4)
40 TFile fmean("laserMean.root")
1fd56785 41 //
2b35e8e6 42 // laser track clasification;
1fd56785 43 //
7b18d067 44 TCut cutT("cutT","abs(Tr.fP[3])<0.06");
45 TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
c18f4385 46 TCut cutN("cutN","fTPCncls>70");
47 TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
48 TCut cutA = cutT+cutPt+cutP;
49 TFile f("laserTPCDebug.root");
50 TTree * treeT = (TTree*)f.Get("Track");
1fd56785 51 //
52 //
53 // Analyze LASER scan
54 //
2b35e8e6 55 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
56 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
57 AliXRDPROOFtoolkit tool;
58 TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
59 chain->Lookup();
1fd56785 60 AliTPCcalibLaser::DumpScanInfo(chain)
61 TFile fscan("laserScan.root")
62 TTree * treeT = (TTree*)fscan.Get("Mean")
63
c6914c83 64*/
65
66
67
68#include "TLinearFitter.h"
69#include "AliTPCcalibLaser.h"
70#include "AliExternalTrackParam.h"
e9f38a4b 71#include "AliESDEvent.h"
72#include "AliESDfriend.h"
c6914c83 73#include "AliESDtrack.h"
74#include "AliTPCTracklet.h"
75#include "TH1D.h"
76#include "TVectorD.h"
77#include "TTreeStream.h"
78#include "TFile.h"
79#include "TF1.h"
80#include "TGraphErrors.h"
81#include "AliTPCclusterMI.h"
82#include "AliTPCseed.h"
83#include "AliTracker.h"
84#include "TClonesArray.h"
1fd56785 85#include "TPad.h"
c6914c83 86
87
88#include "TTreeStream.h"
89#include <iostream>
90#include <sstream>
7b18d067 91#include "AliTPCLaserTrack.h"
92
c6914c83 93using namespace std;
94
95ClassImp(AliTPCcalibLaser)
96
97AliTPCcalibLaser::AliTPCcalibLaser():
e9f38a4b 98 AliTPCcalibBase(),
99 fESD(0),
100 fESDfriend(),
c18f4385 101 fTracksMirror(336),
102 fTracksEsd(336),
103 fTracksEsdParam(336),
104 fTracksTPC(336),
105 fDeltaZ(336), // array of histograms of delta z for each track
106 fDeltaPhi(336), // array of histograms of delta z for each track
107 fDeltaPhiP(336), // array of histograms of delta z for each track
2b35e8e6 108 fSignals(336), // array of dedx signals
c18f4385 109 fFitAside(new TVectorD(3)), // drift fit - A side
110 fFitCside(new TVectorD(3)), // drift fit - C- side
e9f38a4b 111 fRun(0)
c6914c83 112{
113 //
114 // Constructor
115 //
e9f38a4b 116 fTracksEsdParam.SetOwner(kTRUE);
c6914c83 117}
118
119AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
e9f38a4b 120 AliTPCcalibBase(),
121 fESD(0),
122 fESDfriend(0),
c18f4385 123 fTracksMirror(336),
124 fTracksEsd(336),
125 fTracksEsdParam(336),
126 fTracksTPC(336),
127 fDeltaZ(336), // array of histograms of delta z for each track
128 fDeltaPhi(336), // array of histograms of delta z for each track
129 fDeltaPhiP(336), // array of histograms of delta z for each track
2b35e8e6 130 fSignals(336), // array of dedx signals
c18f4385 131 fFitAside(new TVectorD(3)), // drift fit - A side
132 fFitCside(new TVectorD(3)), // drift fit - C- side
e9f38a4b 133 fRun(0)
c6914c83 134{
135 SetName(name);
136 SetTitle(title);
137 //
138 // Constructor
139 //
e9f38a4b 140 fTracksEsdParam.SetOwner(kTRUE);
c6914c83 141}
142
143AliTPCcalibLaser::~AliTPCcalibLaser() {
144 //
145 // destructor
146 //
147}
148
e9f38a4b 149
150
151void AliTPCcalibLaser::Process(AliESDEvent * event) {
c6914c83 152 //
153 //
e9f38a4b 154 // Loop over tracks and call Process function
155 //
156 fESD = event;
157 if (!fESD) {
158 return;
159 }
160 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
161 if (!fESDfriend) {
162 return;
163 }
164 fTracksTPC.Clear();
165 fTracksEsd.Clear();
166 fTracksEsdParam.Delete();
167 //
168 Int_t n=fESD->GetNumberOfTracks();
169 Int_t run = fESD->GetRunNumber();
170 fRun = run;
171 for (Int_t i=0;i<n;++i) {
172 AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
173 AliESDtrack *track=fESD->GetTrack(i);
174 TObject *calibObject=0;
175 AliTPCseed *seed=0;
176 for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
177 if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
178 break;
179 if (track&&seed) FindMirror(track,seed);
180 //
181 }
182
183 FitDriftV();
c18f4385 184 MakeDistHisto();
e9f38a4b 185 //
c18f4385 186 for (Int_t id=0; id<336; id++){
e9f38a4b 187 //
188 //
189 if (!fTracksEsdParam.At(id)) continue;
190 DumpLaser(id);
191 RefitLaser(id);
c18f4385 192
193 }
194}
195
196void AliTPCcalibLaser::MakeDistHisto(){
197 //
198 //
199 //
200 for (Int_t id=0; id<336; id++){
201 //
202 //
203 if (!fTracksEsdParam.At(id)) continue;
204 if (!AcceptLaser(id)) continue;
205 //
206 //
207 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
208 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
209 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
2b35e8e6 210 TH1F * hisSignal = (TH1F*)fSignals.At(id);
211
c18f4385 212 if (!hisdz){
c03e3250 213 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
214 hisdz->SetDirectory(0);
c18f4385 215 fDeltaZ.AddAt(hisdz,id);
216 //
c03e3250 217 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
218 hisdphi->SetDirectory(0);
c18f4385 219 fDeltaPhi.AddAt(hisdphi,id);
220 //
c03e3250 221 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
222 hisdphiP->SetDirectory(0);
c18f4385 223 fDeltaPhiP.AddAt(hisdphiP,id);
c03e3250 224 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
225 hisSignal->SetDirectory(0);
2b35e8e6 226 fSignals.AddAt(hisSignal,id);
c18f4385 227 }
228
229 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
230 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
2b35e8e6 231 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
c03e3250 232 if (!param) return;
233 if (!ltrp) return;
234 if (!track) return;
c18f4385 235 Double_t xyz[3];
236 Double_t pxyz[3];
237 Double_t lxyz[3];
238 Double_t lpxyz[3];
239 param->GetXYZ(xyz);
240 param->GetPxPyPz(pxyz);
241 ltrp->GetXYZ(lxyz);
242 ltrp->GetPxPyPz(lpxyz);
243 //
244 Float_t dz = param->GetZ()-ltrp->GetZ();
245 Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
246 Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
c03e3250 247 if (hisdz) hisdz->Fill(dz);
248 if (hisdphi) hisdphi->Fill(dphi);
249 if (hisdphiP) hisdphiP->Fill(dphiP);
250 if (hisSignal) hisSignal->Fill(track->GetTPCsignal());
c18f4385 251 }
252}
253
254void AliTPCcalibLaser::FitDriftV(){
255 //
256 // Fit drift velocity - linear approximation in the z and global y
257 //
258 static TLinearFitter fdriftA(3,"hyp2");
259 static TLinearFitter fdriftC(3,"hyp2");
260 fdriftA.ClearPoints();
261 fdriftC.ClearPoints();
262 //
263 for (Int_t id=0; id<336; id++){
264 if (!fTracksEsdParam.At(id)) continue;
265 if (!AcceptLaser(id)) continue;
266 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
267 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
268 Double_t xyz[3];
269 Double_t pxyz[3];
270 Double_t lxyz[3];
271 Double_t lpxyz[3];
272 param->GetXYZ(xyz);
273 param->GetPxPyPz(pxyz);
274 ltrp->GetXYZ(lxyz);
275 ltrp->GetPxPyPz(lpxyz);
276 Double_t xxx[2] = {lxyz[2],lxyz[1]};
277 if (ltrp->GetSide()==0){
278 fdriftA.AddPoint(xxx,xyz[2],1);
279 }else{
280 fdriftC.AddPoint(xxx,xyz[2],1);
281 }
282 }
283 Float_t chi2A = 0;
284 Float_t chi2C = 0;
285 Int_t npointsA=0;
286 Int_t npointsC=0;
287 //
288 if (fdriftA.GetNpoints()>10){
289 fdriftA.Eval();
290 fdriftA.EvalRobust(0.8);
291 fdriftA.GetParameters(*fFitAside);
292 npointsA= fdriftA.GetNpoints();
293 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
294 }
295 if (fdriftC.GetNpoints()>10){
296 fdriftC.Eval();
297 fdriftC.EvalRobust(0.8);
298 fdriftC.GetParameters(*fFitCside);
299 npointsC= fdriftC.GetNpoints();
300 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
e9f38a4b 301 }
c18f4385 302
303 if (fStreamLevel>0){
304 TTreeSRedirector *cstream = GetDebugStreamer();
305 Int_t time = fESD->GetTimeStamp();
306 if (cstream){
307 (*cstream)<<"driftv"<<
308 "driftA.="<<fFitAside<<
309 "driftC.="<<fFitCside<<
310 "chi2A="<<chi2A<<
311 "chi2C="<<chi2C<<
312 "nA="<<npointsA<<
313 "nC="<<npointsC<<
314 "time="<<time<<
315 "\n";
316 }
317 }
318 //
319}
320
321
322Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
323 //
324 //
325 //
326 /*
327 TCut cutT("cutT","abs(Tr.fP[3])<0.06");
328 TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
329 TCut cutN("cutN","fTPCncls>70");
330 TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
331 TCut cutA = cutT+cutPt+cutP;
332 */
333 AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
334 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
335 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
336
337 if (TMath::Abs(param->GetParameter()[4])>0.03) return kFALSE;
338 if (TMath::Abs(param->GetParameter()[3])>0.06) return kFALSE;
339 if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.06) return kFALSE;
340 if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>10) return kFALSE;
341 //
342 // dedx cut
343 //
344 if (TMath::Abs(track->GetTPCsignal())<20) return kFALSE;
345 if (TMath::Abs(track->GetTPCsignal())>800) return kFALSE;
346 //
347 return kTRUE;
e9f38a4b 348}
349
350Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
351 //
352 // Find corresponding mirror
353 // add the corresponding tracks
c6914c83 354 //
1bd37cc0 355 Float_t kRadius0 = 252;
c18f4385 356 Float_t kRadius = 253.4;
e9f38a4b 357 if (!track->GetOuterParam()) return -1;
c6914c83 358 AliExternalTrackParam param(*(track->GetOuterParam()));
1bd37cc0 359 AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
360 AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
7b18d067 361 AliTPCLaserTrack ltr;
362 AliTPCLaserTrack *ltrp=0x0;
1bd37cc0 363 //
7b18d067 364 Int_t id = AliTPCLaserTrack::IdentifyTrack(&param);
1bd37cc0 365 if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
366 ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
367 else
368 ltrp=&ltr;
e9f38a4b 369
370 if (id>=0){
c18f4385 371 //
372 //
373 Float_t radius=TMath::Abs(ltrp->GetX());
374 AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
375 //
e9f38a4b 376 if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
377 fTracksEsdParam.AddAt(param.Clone(),id);
378 fTracksEsd.AddAt(track,id);
379 fTracksTPC.AddAt(seed,id);
c18f4385 380 //
e9f38a4b 381 }
c18f4385 382
e9f38a4b 383 return id;
384}
385
386
387
388void AliTPCcalibLaser::DumpLaser(Int_t id) {
389 //
390 // Dump Laser info to the tree
391 //
392 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
393 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
394 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1bd37cc0 395 //
396 // Fast laser ID
397 //
398 Double_t xyz[3];
399 Double_t pxyz[3];
c18f4385 400 Double_t lxyz[3];
401 Double_t lpxyz[3];
e9f38a4b 402 param->GetXYZ(xyz);
403 param->GetPxPyPz(pxyz);
c18f4385 404 ltrp->GetXYZ(lxyz);
405 ltrp->GetPxPyPz(lpxyz);
406
c6914c83 407 if (fStreamLevel>0){
408 TTreeSRedirector *cstream = GetDebugStreamer();
c18f4385 409 Int_t time = fESD->GetTimeStamp();
410 Bool_t accept = AcceptLaser(id);
c6914c83 411 if (cstream){
412 (*cstream)<<"Track"<<
e9f38a4b 413 "run="<<fRun<<
7b18d067 414 "id="<<id<<
c18f4385 415 "accept="<<accept<<
416 "driftA.="<<fFitAside<<
417 "driftC.="<<fFitCside<<
418 "time="<<time<<
7b18d067 419 //
420 "LTr.="<<ltrp<<
421 "Esd.="<<track<<
e9f38a4b 422 "Tr.="<<param<<
c6914c83 423 "x0="<<xyz[0]<<
424 "x1="<<xyz[1]<<
425 "x2="<<xyz[2]<<
426 "px0="<<pxyz[0]<<
427 "px1="<<pxyz[1]<<
428 "px2="<<pxyz[2]<<
c18f4385 429 //
430 "lx0="<<lxyz[0]<<
431 "lx1="<<lxyz[1]<<
432 "lx2="<<lxyz[2]<<
433 "lpx0="<<lpxyz[0]<<
434 "lpx1="<<lpxyz[1]<<
435 "lpx2="<<lpxyz[2]<<
c6914c83 436 "\n";
437 }
438 }
439}
440
e9f38a4b 441
442void AliTPCcalibLaser::RefitLaser(Int_t id){
443 //
444 // Refit the track store residuals
445 //
446
447 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
448 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
449 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
450
451 //
452 static TLinearFitter fy2(3,"hyp2");
453 static TLinearFitter fz2(3,"hyp2");
454 static TLinearFitter fy1(2,"hyp1");
455 static TLinearFitter fz1(2,"hyp1");
456 static TVectorD vecy2,vecz2,vecy1,vecz1;
457
458 const Int_t kMinClusters=20;
459 Int_t nclusters[72];
460 //
461 for (Int_t i=0;i<72;++i) nclusters[i]=0;
462
463 for (Int_t i=0;i<160;++i) {
464 AliTPCclusterMI *c=track->GetClusterPointer(i);
465 if (c) nclusters[c->GetDetector()]++;
466 }
467
468 for (Int_t isec=0; isec<72;isec++){
469 if (nclusters[isec]<kMinClusters) continue;
470 fy2.ClearPoints();
471 fz2.ClearPoints();
472 fy1.ClearPoints();
473 fz1.ClearPoints();
474 //
475 for (Int_t irow=0;irow<160;++irow) {
476 AliTPCclusterMI *c=track->GetClusterPointer(irow);
477 //if (c && RejectCluster(c)) continue;
478 if (c&&c->GetDetector()==isec) {
479 Double_t xd = c->GetX()-120;;
480 Double_t x[2]={xd,xd*xd};
481 fy2.AddPoint(x,c->GetY());
482 fz2.AddPoint(x,c->GetZ());
483 //
484 fy1.AddPoint(x,c->GetY());
485 fz1.AddPoint(x,c->GetZ());
486 }
487 }
488 fy2.Eval();
489 fz2.Eval();
490 fy1.Eval();
491 fz1.Eval();
492 fy1.GetParameters(vecy1);
493 fy2.GetParameters(vecy2);
494 fz1.GetParameters(vecz1);
495 fz2.GetParameters(vecz2);
496
497 if (fStreamLevel>0){
498 TTreeSRedirector *cstream = GetDebugStreamer();
499 if (cstream){
500 Float_t dedx = track->GetdEdx();
501 (*cstream)<<"Tracklet"<<
502 "LTr.="<<ltrp<<
503 "Tr.="<<param<<
504 "sec="<<isec<<
505 "ncl="<<nclusters[isec]<<
506 "dedx="<<dedx<<
507 "dedx="<<dedx<<
508 "vy1.="<<&vecy1<<
509 "vy2.="<<&vecy2<<
510 "vz1.="<<&vecz1<<
511 "vz2.="<<&vecz2<<
512 "\n";
513 }
514 }
515 }
516 //
517 //
518 //
519 // for (Int_t irow=0;irow<160;++irow) {
520 // AliTPCclusterMI *c=track->GetClusterPointer(irow);
521 // if (c && RejectCluster(c)) continue;
522 // if (c&&c->GetDetector()==isec) {
523 // Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
524 // fy2.AddPoint(&x,c->GetY());
525 // fz2.AddPoint(&x,c->GetZ());
526 // //
527 // fy1.AddPoint(&x,c->GetY());
528 // fz1.AddPoint(&x,c->GetZ());
529 // }
530 // }
531
532}
533
534
1fd56785 535void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield,Int_t minEntries){
2b35e8e6 536 //
1fd56785 537 // Dump information about laser beams
538 // isOK variable indicates usability of the beam
539 // Beam is not usable if:
540 // a. No entries in range (krmsCut0)
541 // b. Big sperad (krmscut1)
542 // c. RMSto fit sigma bigger then (kmultiCut)
543 // d. Too big angular spread
544 //
545
546 const Float_t krmsCut0=0.001;
547 const Float_t krmsCut1=0.16;
548 const Float_t kmultiCut=2;
549 const Float_t kcutP0=0.002;
2b35e8e6 550 //
551 AliTPCcalibLaser *laser = this;
552 TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
553 TF1 fg("fg","gaus");
554 //
555 //
556 for (Int_t id=0; id<336; id++){
1fd56785 557 Bool_t isOK=kTRUE;
2b35e8e6 558 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
559 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
560 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
1fd56785 561 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
2b35e8e6 562 if (!hisphi) continue;;
563 Double_t entries = hisphi->GetEntries();
1fd56785 564 if (entries<minEntries) continue;
2b35e8e6 565 //
566 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
567 if (!ltrp) {
568 AliTPCLaserTrack::LoadTracks();
569 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
570 }
571 Float_t meanphi = hisphi->GetMean();
572 Float_t rmsphi = hisphi->GetRMS();
1fd56785 573 //
2b35e8e6 574 Float_t meanphiP = hisphiP->GetMean();
575 Float_t rmsphiP = hisphiP->GetRMS();
576 Float_t meanZ = hisZ->GetMean();
577 Float_t rmsZ = hisZ->GetRMS();
578 hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
579 Double_t gphi1 = fg.GetParameter(1);
580 Double_t gphi2 = fg.GetParameter(2);
581 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
582 Double_t gphiP1 = fg.GetParameter(1);
583 Double_t gphiP2 = fg.GetParameter(2);
584 hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
585 Double_t gz1 = fg.GetParameter(1);
586 Double_t gz2 = fg.GetParameter(2);
1fd56785 587 //
588 Float_t meanS=hisS->GetMean();
2b35e8e6 589 //
590 Double_t lxyz[3];
591 Double_t lpxyz[3];
592 ltrp->GetXYZ(lxyz);
593 ltrp->GetPxPyPz(lpxyz);
1fd56785 594
595 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
596 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
597 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
598 if (gphiP2>kcutP0) isOK=kFALSE;
2b35e8e6 599 //
600 (*pcstream)<<"Mean"<<
1fd56785 601 "isOK="<<isOK<<
2b35e8e6 602 "entries="<<entries<< // number of entries
603 "bz="<<bfield<< // bfield
604 "LTr.="<<ltrp<< // refernece track
1fd56785 605 //
2b35e8e6 606 "lx0="<<lxyz[0]<< // reference x
607 "lx1="<<lxyz[1]<< // reference y
608 "lx2="<<lxyz[2]<< // refernece z
609 "lpx0="<<lpxyz[0]<< // reference x
610 "lpx1="<<lpxyz[1]<< // reference y
611 "lpx2="<<lpxyz[2]<< // refernece z
1fd56785 612 //
613 "msig="<<meanS<<
614 //
2b35e8e6 615 "mphi="<<meanphi<< //
616 "rmsphi="<<rmsphi<< //
617 "gphi1="<<gphi1<<
618 "gphi2="<<gphi2<<
1fd56785 619 //
2b35e8e6 620 "mphiP="<<meanphiP<< //
621 "rmsphiP="<<rmsphiP<< //
622 "gphiP1="<<gphiP1<<
623 "gphiP2="<<gphiP2<<
1fd56785 624 //
2b35e8e6 625 "meanZ="<<meanZ<<
626 "rmsZ="<<rmsZ<<
627 "gz1="<<gz1<<
628 "gz2="<<gz2<<
629
630 "\n";
631 }
632 delete pcstream;
633}
634
1fd56785 635
636
637void AliTPCcalibLaser::DumpScanInfo(TTree * chain){
638 //
639 //
640 //
641 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
642 TFile * f = pcstream->GetFile();
643 f->mkdir("dirphi");
644 f->mkdir("dirphiP");
645 f->mkdir("dirZ");
646 TF1 fp("p1","pol1");
647 //
648 //
649 char cut[1000];
650 char grname[1000];
651 char grnamefull[1000];
652
653 Double_t mphi[100];
654 Double_t mphiP[100];
655 Double_t smphi[100];
656 Double_t smphiP[100];
657 Double_t mZ[100];
658 Double_t smZ[100];
659 Double_t bz[100];
660 Double_t sbz[100];
661 // fit parameters
662 Double_t pphi[3];
663 Double_t pphiP[3];
664 Double_t pmZ[3];
665 //
666 for (Int_t id=0; id<336; id++){
667 // id =205;
668 sprintf(cut,"isOK&&fId==%d",id);
669 Int_t entries = chain->Draw("bz",cut,"goff");
670 if (entries<3) continue;
671 AliTPCLaserTrack *ltrp = 0;;
672 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
673 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
674 Double_t lxyz[3];
675 Double_t lpxyz[3];
676 ltrp->GetXYZ(lxyz);
677 ltrp->GetPxPyPz(lpxyz);
678
679 chain->Draw("bz",cut,"goff");
680 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
681 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
682 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
683 //
684 chain->Draw("gphi1",cut,"goff");
685 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
686 chain->Draw("0.05*abs(mphi)+gphi2",cut,"goff");
687 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
688 //
689 chain->Draw("gphiP1",cut,"goff");
690 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
691 chain->Draw("0.05*abs(mphiP)+gphiP2",cut,"goff");
692 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
693 //
694 chain->Draw("gz1",cut,"goff");
695 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
696 chain->Draw("0.01*abs(meanZ)+gz2",cut,"goff");
697 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
698 //
699 //
700 sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
701 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
702 // store data
703 // phi
704 f->cd("dirphi");
705 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
706 grphi->Draw("a*");
707 grphi->Fit(&fp);
708 pphi[0] = fp.GetParameter(0); // offset
709 pphi[1] = fp.GetParameter(1); // slope
710 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
711 sprintf(grname,"phi_id%d",id);
712 grphi->SetName(grname); grphi->SetTitle(grnamefull);
713 grphi->GetXaxis()->SetTitle("b_{z} (T)");
714 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
715 grphi->Draw("a*");
716
717 grphi->Write();
718 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
719 // phiP
720 f->cd("dirphiP)");
721 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
722 grphiP->Draw("a*");
723 grphiP->Fit(&fp);
724 pphiP[0] = fp.GetParameter(0); // offset
725 pphiP[1] = fp.GetParameter(1); // slope
726 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
727 sprintf(grname,"phiP_id%d",id);
728 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
729 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
730 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
731
732 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
733 grphiP->Write();
734 //
735 //Z
736 f->cd("dirZ");
737 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
738 grmZ->Draw("a*");
739 grmZ->Fit(&fp);
740 pmZ[0] = fp.GetParameter(0); // offset
741 pmZ[1] = fp.GetParameter(1); // slope
742 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
743 sprintf(grname,"mZ_id%d",id);
744 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
745 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
746 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
747
748 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
749 grmZ->Write();
750
751
752 for (Int_t ientry=0; ientry<entries; ientry++){
753 (*pcstream)<<"Mean"<<
754 "id="<<id<<
755 "LTr.="<<ltrp<<
756 "entries="<<entries<<
757 "bz="<<bz[ientry]<<
758 "lx0="<<lxyz[0]<< // reference x
759 "lx1="<<lxyz[1]<< // reference y
760 "lx2="<<lxyz[2]<< // refernece z
761 "lpx0="<<lpxyz[0]<< // reference x
762 "lpx1="<<lpxyz[1]<< // reference y
763 "lpx2="<<lpxyz[2]<< // refernece z
764 //values
765 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
766 "pphi0="<<pphi[0]<< // offset
767 "pphi1="<<pphi[1]<< // mean
768 "pphi2="<<pphi[2]<< // norm chi2
769 //
770 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
771 "pphiP0="<<pphiP[0]<< // offset
772 "pphiP1="<<pphiP[1]<< // mean
773 "pphiP2="<<pphiP[2]<< // norm chi2
774 //
775 "gz1="<<mZ[ientry]<<
776 "pmZ0="<<pmZ[0]<< // offset
777 "pmZ1="<<pmZ[1]<< // mean
778 "pmZ2="<<pmZ[2]<< // norm chi2
779 "\n";
780 }
781 }
782
783 delete pcstream;
784
785}
786
787
c6914c83 788void AliTPCcalibLaser::Analyze(){
789 //
790 //
791 //
792}
793
794
c03e3250 795Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
c6914c83 796
c03e3250 797 TIterator* iter = li->MakeIterator();
798 AliTPCcalibLaser* cal = 0;
c6914c83 799
c03e3250 800 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
801 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
802 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
803 return -1;
804 }
805 //
806 // fHistNTracks->Add(cal->fHistNTracks);
807// fClusters->Add(cal->fClusters);
808// fModules->Add(cal->fModules);
809// fHistPt->Add(cal->fHistPt);
810// fPtResolution->Add(cal->fPtResolution);
811// fDeDx->Add(cal->fDeDx);
812
813
814 TH1F *h=0x0;
815 TH1F *hm=0x0;
816
817 for (Int_t id=0; id<336; id++){
818 // merge fDeltaZ histograms
819 hm = (TH1F*)cal->fDeltaZ.At(id);
820 h = (TH1F*)fDeltaZ.At(id);
821 if (!h) {
822 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
823 fDeltaZ.AddAt(h,id);
824 }
825 if (hm) h->Add(hm);
826 // merge fDeltaPhi histograms
827 hm = (TH1F*)cal->fDeltaPhi.At(id);
828 h = (TH1F*)fDeltaPhi.At(id);
829 if (!h) {
830 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
831 fDeltaPhi.AddAt(h,id);
832 }
833 if (hm) h->Add(hm);
834 // merge fDeltaPhiP histograms
835 hm = (TH1F*)cal->fDeltaPhiP.At(id);
836 h = (TH1F*)fDeltaPhiP.At(id);
837 if (!h) {
838 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
839 fDeltaPhiP.AddAt(h,id);
840 }
841 if (hm) h->Add(hm);
842 // merge fSignals histograms
843 hm = (TH1F*)cal->fSignals.At(id);
844 h = (TH1F*)fSignals.At(id);
845 if (!h) {
846 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
847 fSignals.AddAt(h,id);
848 }
849 if (hm) h->Add(hm);
850 }
851 }
852 return 0;
853}
854
855
856
857/*
858 gSystem->Load("libSTAT.so")
859 TStatToolkit toolkit;
860 Double_t chi2;
861 TVectorD fitParam;
862 TMatrixD covMatrix;
863 Int_t npoints;
1fd56785 864 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003");
c03e3250 865
866
867TString fstring="";
1fd56785 868//
869fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
870fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
871fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
872fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
873//
874fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
875fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
876fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
877fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
878//
879fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
880fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
881fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
882fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
c5595838 883
c5595838 884
c5595838 885
c5595838 886
1fd56785 887 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c03e3250 888
1fd56785 889 treeT->SetAlias("fit",strq0->Data());
890
c03e3250 891
1fd56785 892 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c5595838 893
1fd56785 894 treeT->SetAlias("fitP",strqP->Data());
c5595838 895
896
c5595838 897
c03e3250 898
899
900
901
1fd56785 902
c03e3250 903 */