Adding the 3 different fit function of laser tracks (Jens)
[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)");
ae69f16f 715 grphi->SetMaximum(1.2);
716 grphi->SetMinimum(-1.2);
1fd56785 717 grphi->Draw("a*");
718
719 grphi->Write();
720 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
721 // phiP
722 f->cd("dirphiP)");
723 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
724 grphiP->Draw("a*");
725 grphiP->Fit(&fp);
726 pphiP[0] = fp.GetParameter(0); // offset
727 pphiP[1] = fp.GetParameter(1); // slope
728 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
729 sprintf(grname,"phiP_id%d",id);
730 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
731 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
732 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
ae69f16f 733 grphiP->SetMaximum(pphiP[0]+0.005);
734 grphiP->SetMinimum(pphiP[0]-0.005);
1fd56785 735
736 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
737 grphiP->Write();
738 //
739 //Z
740 f->cd("dirZ");
741 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
742 grmZ->Draw("a*");
743 grmZ->Fit(&fp);
744 pmZ[0] = fp.GetParameter(0); // offset
745 pmZ[1] = fp.GetParameter(1); // slope
746 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
747 sprintf(grname,"mZ_id%d",id);
748 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
749 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
750 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
751
752 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
753 grmZ->Write();
754
755
756 for (Int_t ientry=0; ientry<entries; ientry++){
757 (*pcstream)<<"Mean"<<
758 "id="<<id<<
759 "LTr.="<<ltrp<<
760 "entries="<<entries<<
761 "bz="<<bz[ientry]<<
762 "lx0="<<lxyz[0]<< // reference x
763 "lx1="<<lxyz[1]<< // reference y
764 "lx2="<<lxyz[2]<< // refernece z
765 "lpx0="<<lpxyz[0]<< // reference x
766 "lpx1="<<lpxyz[1]<< // reference y
767 "lpx2="<<lpxyz[2]<< // refernece z
768 //values
769 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
770 "pphi0="<<pphi[0]<< // offset
771 "pphi1="<<pphi[1]<< // mean
772 "pphi2="<<pphi[2]<< // norm chi2
773 //
774 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
775 "pphiP0="<<pphiP[0]<< // offset
776 "pphiP1="<<pphiP[1]<< // mean
777 "pphiP2="<<pphiP[2]<< // norm chi2
778 //
779 "gz1="<<mZ[ientry]<<
780 "pmZ0="<<pmZ[0]<< // offset
781 "pmZ1="<<pmZ[1]<< // mean
782 "pmZ2="<<pmZ[2]<< // norm chi2
783 "\n";
784 }
785 }
786
787 delete pcstream;
788
789}
790
791
c6914c83 792void AliTPCcalibLaser::Analyze(){
793 //
794 //
795 //
796}
797
798
c03e3250 799Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
c6914c83 800
c03e3250 801 TIterator* iter = li->MakeIterator();
802 AliTPCcalibLaser* cal = 0;
c6914c83 803
c03e3250 804 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
805 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
806 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
807 return -1;
808 }
809 //
810 // fHistNTracks->Add(cal->fHistNTracks);
811// fClusters->Add(cal->fClusters);
812// fModules->Add(cal->fModules);
813// fHistPt->Add(cal->fHistPt);
814// fPtResolution->Add(cal->fPtResolution);
815// fDeDx->Add(cal->fDeDx);
816
817
818 TH1F *h=0x0;
819 TH1F *hm=0x0;
820
821 for (Int_t id=0; id<336; id++){
822 // merge fDeltaZ histograms
823 hm = (TH1F*)cal->fDeltaZ.At(id);
824 h = (TH1F*)fDeltaZ.At(id);
825 if (!h) {
826 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
827 fDeltaZ.AddAt(h,id);
828 }
829 if (hm) h->Add(hm);
830 // merge fDeltaPhi histograms
831 hm = (TH1F*)cal->fDeltaPhi.At(id);
832 h = (TH1F*)fDeltaPhi.At(id);
833 if (!h) {
834 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
835 fDeltaPhi.AddAt(h,id);
836 }
837 if (hm) h->Add(hm);
838 // merge fDeltaPhiP histograms
839 hm = (TH1F*)cal->fDeltaPhiP.At(id);
840 h = (TH1F*)fDeltaPhiP.At(id);
841 if (!h) {
842 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
843 fDeltaPhiP.AddAt(h,id);
844 }
845 if (hm) h->Add(hm);
846 // merge fSignals histograms
847 hm = (TH1F*)cal->fSignals.At(id);
848 h = (TH1F*)fSignals.At(id);
849 if (!h) {
850 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
851 fSignals.AddAt(h,id);
852 }
853 if (hm) h->Add(hm);
854 }
855 }
856 return 0;
857}
858
859
860
861/*
862 gSystem->Load("libSTAT.so")
863 TStatToolkit toolkit;
864 Double_t chi2;
865 TVectorD fitParam;
866 TMatrixD covMatrix;
867 Int_t npoints;
ae69f16f 868 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
c03e3250 869
870
871TString fstring="";
1fd56785 872//
873fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
874fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
875fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
876fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
877//
878fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
879fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
880fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
881fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
882//
883fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
884fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
885fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
886fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
c5595838 887
c5595838 888
c5595838 889
c5595838 890
1fd56785 891 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c03e3250 892
1fd56785 893 treeT->SetAlias("fit",strq0->Data());
894
c03e3250 895
1fd56785 896 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c5595838 897
1fd56785 898 treeT->SetAlias("fitP",strqP->Data());
c5595838 899
900
ae69f16f 901 TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
902 treeT->SetAlias("fitD",strqDrift->Data());
c5595838 903
c03e3250 904
ae69f16f 905treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
906{
907for (Int_t i=0; i<6;i++){
908treeT->SetLineColor(i+2);
909treeT->SetMarkerSize(1);
910treeT->SetMarkerStyle(22+i);
911treeT->SetMarkerColor(i+2);
c03e3250 912
ae69f16f 913treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
914}
915}
c03e3250 916 */
68ff0583 917
918
919
920/*
921 TTree * tree = (TTree*)f.Get("FitModels");
922
923 TEventList listLFit0("listLFit0","listLFit0");
924 TEventList listLFit1("listLFit1","listLFit1");
925
926 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
927 tree->SetEventList(&listLFit0);
928
929
930*/