]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TPC/AliTPCcalibLaser.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Adding pad type correction (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibLaser.cxx
... / ...
CommitLineData
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/*
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 //
31 // To make laser scan the user interaction neccessary
32 //
33 .x ~/UliStyle.C
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")
41 //
42 // laser track clasification;
43 //
44 TCut cutT("cutT","abs(Tr.fP[3])<0.06");
45 TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
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");
51 //
52 //
53 // Analyze LASER scan
54 //
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();
60 AliTPCcalibLaser::DumpScanInfo(chain)
61 TFile fscan("laserScan.root")
62 TTree * treeT = (TTree*)fscan.Get("Mean")
63 //
64 // Analyze laser
65 //
66 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
67 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
68 AliXRDPROOFtoolkit tool;
69 TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
70 chain->Lookup();
71
72*/
73
74
75
76#include "TLinearFitter.h"
77#include "AliTPCcalibLaser.h"
78#include "AliExternalTrackParam.h"
79#include "AliESDEvent.h"
80#include "AliESDfriend.h"
81#include "AliESDtrack.h"
82#include "AliTPCTracklet.h"
83#include "TH1D.h"
84#include "TH1F.h"
85#include "TProfile.h"
86#include "TVectorD.h"
87#include "TTreeStream.h"
88#include "TFile.h"
89#include "TF1.h"
90#include "TGraphErrors.h"
91#include "AliTPCclusterMI.h"
92#include "AliTPCseed.h"
93#include "AliTracker.h"
94#include "AliLog.h"
95#include "TClonesArray.h"
96#include "TPad.h"
97
98
99#include "TTreeStream.h"
100#include <iostream>
101#include <sstream>
102#include "AliTPCLaserTrack.h"
103
104using namespace std;
105
106ClassImp(AliTPCcalibLaser)
107
108AliTPCcalibLaser::AliTPCcalibLaser():
109 AliTPCcalibBase(),
110 fESD(0),
111 fESDfriend(),
112 fTracksMirror(336),
113 fTracksEsd(336),
114 fTracksEsdParam(336),
115 fTracksTPC(336),
116 fDeltaZ(336),
117 fDeltaPhi(336),
118 fDeltaPhiP(336),
119 fSignals(336),
120 fDeltaYres(336),
121 fDeltaZres(336),
122 fPol2Par2InY(336),
123 fDiffPar1InY(336),
124 fPol2Par2OutY(336),
125 fDiffPar1OutY(336),
126 fPol2Par2InZ(336),
127 fDiffPar1InZ(336),
128 fPol2Par2OutZ(336),
129 fDiffPar1OutZ(336),
130 fFitAside(new TVectorD(3)),
131 fFitCside(new TVectorD(3)),
132 fEdgeXcuts(5),
133 fEdgeYcuts(5),
134 fNClCuts(5),
135 fNcuts(0),
136 fRun(0),
137 fEvent(0)
138{
139 //
140 // Constructor
141 //
142 fTracksEsdParam.SetOwner(kTRUE);
143}
144
145AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
146 AliTPCcalibBase(),
147 fESD(0),
148 fESDfriend(0),
149 fTracksMirror(336),
150 fTracksEsd(336),
151 fTracksEsdParam(336),
152 fTracksTPC(336),
153 fDeltaZ(336), // array of histograms of delta z for each track
154 fDeltaPhi(336), // array of histograms of delta z for each track
155 fDeltaPhiP(336), // array of histograms of delta z for each track
156 fSignals(336), // array of dedx signals
157 fDeltaYres(336),
158 fDeltaZres(336),
159 fPol2Par2InY(336),
160 fDiffPar1InY(336),
161 fPol2Par2OutY(336),
162 fDiffPar1OutY(336),
163 fPol2Par2InZ(336),
164 fDiffPar1InZ(336),
165 fPol2Par2OutZ(336),
166 fDiffPar1OutZ(336),
167 fFitAside(new TVectorD(3)), // drift fit - A side
168 fFitCside(new TVectorD(3)), // drift fit - C- side
169 fEdgeXcuts(5), // cuts in local x direction; used in the refit of the laser tracks
170 fEdgeYcuts(5), // cuts in local y direction; used in the refit of the laser tracks
171 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
172 fNcuts(0), // number of cuts
173 fRun(0),
174 fEvent(0)
175{
176 SetName(name);
177 SetTitle(title);
178 //
179 // Constructor
180 //
181 fTracksEsdParam.SetOwner(kTRUE);
182}
183
184AliTPCcalibLaser::~AliTPCcalibLaser() {
185 //
186 // destructor
187 //
188}
189
190
191
192void AliTPCcalibLaser::Process(AliESDEvent * event) {
193 //
194 //
195 // Loop over tracks and call Process function
196 //
197 fESD = event;
198 if (!fESD) {
199 return;
200 }
201 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
202 if (!fESDfriend) {
203 return;
204 }
205 AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
206 fTracksTPC.Clear();
207 fTracksEsd.Clear();
208 fTracksEsdParam.Delete();
209 //
210 Int_t n=fESD->GetNumberOfTracks();
211 Int_t run = fESD->GetRunNumber();
212 fRun = run;
213 for (Int_t i=0;i<n;++i) {
214 AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
215 AliESDtrack *track=fESD->GetTrack(i);
216 TObject *calibObject=0;
217 AliTPCseed *seed=0;
218 for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
219 if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
220 break;
221 if (track&&seed) FindMirror(track,seed);
222 //
223 }
224
225 FitDriftV();
226 MakeDistHisto();
227 //
228 for (Int_t id=0; id<336; id++){
229 //
230 //
231 if (!fTracksEsdParam.At(id)) continue;
232 DumpLaser(id);
233// RefitLaser(id);
234 RefitLaserJW(id);
235
236 }
237// fEvent++;
238}
239
240void AliTPCcalibLaser::MakeDistHisto(){
241 //
242 //
243 //
244 for (Int_t id=0; id<336; id++){
245 //
246 //
247 if (!fTracksEsdParam.At(id)) continue;
248 if (!AcceptLaser(id)) continue;
249 //
250 //
251 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
252 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
253 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
254 TH1F * hisSignal = (TH1F*)fSignals.At(id);
255
256 if (!hisdz){
257 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
258 hisdz->SetDirectory(0);
259 fDeltaZ.AddAt(hisdz,id);
260 //
261 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
262 hisdphi->SetDirectory(0);
263 fDeltaPhi.AddAt(hisdphi,id);
264 //
265 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
266 hisdphiP->SetDirectory(0);
267 fDeltaPhiP.AddAt(hisdphiP,id);
268 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
269 hisSignal->SetDirectory(0);
270 fSignals.AddAt(hisSignal,id);
271 }
272
273 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
274 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
275 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
276 if (!param) return;
277 if (!ltrp) return;
278 if (!track) return;
279 Double_t xyz[3];
280 Double_t pxyz[3];
281 Double_t lxyz[3];
282 Double_t lpxyz[3];
283 param->GetXYZ(xyz);
284 param->GetPxPyPz(pxyz);
285 ltrp->GetXYZ(lxyz);
286 ltrp->GetPxPyPz(lpxyz);
287 //
288 Float_t dz = param->GetZ()-ltrp->GetZ();
289 Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
290 Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
291 if (hisdz) hisdz->Fill(dz);
292 if (hisdphi) hisdphi->Fill(dphi);
293 if (hisdphiP) hisdphiP->Fill(dphiP);
294 if (hisSignal) hisSignal->Fill(track->GetTPCsignal());
295 }
296}
297
298void AliTPCcalibLaser::FitDriftV(){
299 //
300 // Fit drift velocity - linear approximation in the z and global y
301 //
302 static TLinearFitter fdriftA(3,"hyp2");
303 static TLinearFitter fdriftC(3,"hyp2");
304 fdriftA.ClearPoints();
305 fdriftC.ClearPoints();
306 //
307 for (Int_t id=0; id<336; id++){
308 if (!fTracksEsdParam.At(id)) continue;
309 if (!AcceptLaser(id)) continue;
310 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
311 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
312 Double_t xyz[3];
313 Double_t pxyz[3];
314 Double_t lxyz[3];
315 Double_t lpxyz[3];
316 param->GetXYZ(xyz);
317 param->GetPxPyPz(pxyz);
318 ltrp->GetXYZ(lxyz);
319 ltrp->GetPxPyPz(lpxyz);
320 Double_t xxx[2] = {lxyz[2],lxyz[1]};
321 if (ltrp->GetSide()==0){
322 fdriftA.AddPoint(xxx,xyz[2],1);
323 }else{
324 fdriftC.AddPoint(xxx,xyz[2],1);
325 }
326 }
327 Float_t chi2A = 0;
328 Float_t chi2C = 0;
329 Int_t npointsA=0;
330 Int_t npointsC=0;
331 //
332 if (fdriftA.GetNpoints()>10){
333 fdriftA.Eval();
334 fdriftA.EvalRobust(0.8);
335 fdriftA.GetParameters(*fFitAside);
336 npointsA= fdriftA.GetNpoints();
337 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
338 }
339 if (fdriftC.GetNpoints()>10){
340 fdriftC.Eval();
341 fdriftC.EvalRobust(0.8);
342 fdriftC.GetParameters(*fFitCside);
343 npointsC= fdriftC.GetNpoints();
344 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
345 }
346
347 if (fStreamLevel>0){
348 TTreeSRedirector *cstream = GetDebugStreamer();
349 Int_t time = fESD->GetTimeStamp();
350 if (cstream){
351 (*cstream)<<"driftv"<<
352 "driftA.="<<fFitAside<<
353 "driftC.="<<fFitCside<<
354 "chi2A="<<chi2A<<
355 "chi2C="<<chi2C<<
356 "nA="<<npointsA<<
357 "nC="<<npointsC<<
358 "time="<<time<<
359 "\n";
360 }
361 }
362 //
363}
364
365
366Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
367 //
368 //
369 //
370 /*
371 TCut cutT("cutT","abs(Tr.fP[3])<0.06");
372 TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
373 TCut cutN("cutN","fTPCncls>70");
374 TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
375 TCut cutA = cutT+cutPt+cutP;
376 */
377 AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
378 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
379 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
380
381 if (TMath::Abs(param->GetParameter()[4])>0.03) return kFALSE;
382 if (TMath::Abs(param->GetParameter()[3])>0.06) return kFALSE;
383 if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.06) return kFALSE;
384 if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>10) return kFALSE;
385 //
386 // dedx cut
387 //
388 if (TMath::Abs(track->GetTPCsignal())<20) return kFALSE;
389 if (TMath::Abs(track->GetTPCsignal())>800) return kFALSE;
390 //
391 return kTRUE;
392}
393
394Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
395 //
396 // Find corresponding mirror
397 // add the corresponding tracks
398 //
399 Float_t kRadius0 = 252;
400 Float_t kRadius = 253.4;
401 if (!track->GetOuterParam()) return -1;
402 AliExternalTrackParam param(*(track->GetOuterParam()));
403 AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
404 AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
405 AliTPCLaserTrack ltr;
406 AliTPCLaserTrack *ltrp=0x0;
407 // AliTPCLaserTrack *ltrpjw=0x0;
408 //
409 Int_t id = AliTPCLaserTrack::IdentifyTrack(&param);
410 // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(&param);
411 //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
412
413 if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
414 ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
415 else
416 ltrp=&ltr;
417 /*
418 if (idjw!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(idjw)))
419 ltrpjw=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(idjw);
420 else
421 ltrpjw=&ltr;
422
423
424 if (fStreamLevel>0){
425 TTreeSRedirector *cstream = GetDebugStreamer();
426 if (cstream){
427 (*cstream)<<"idcmp"<<
428 "id=" << id <<
429 "idjw=" << idjw <<
430 "tr.=" << ltrp <<
431 "trjw.="<< ltrpjw <<
432 "seed.="<<seed<<
433 "event="<<fEvent <<
434 "\n";
435 }
436 }
437
438 */
439
440
441 if (id>=0){
442 //
443 //
444 Float_t radius=TMath::Abs(ltrp->GetX());
445 AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
446 //
447 if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
448 fTracksEsdParam.AddAt(param.Clone(),id);
449 fTracksEsd.AddAt(track,id);
450 fTracksTPC.AddAt(seed,id);
451 //
452 }
453
454 return id;
455}
456
457
458
459void AliTPCcalibLaser::DumpLaser(Int_t id) {
460 //
461 // Dump Laser info to the tree
462 //
463 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
464 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
465 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
466 //
467 // Fast laser ID
468 //
469 Double_t xyz[3];
470 Double_t pxyz[3];
471 Double_t lxyz[3];
472 Double_t lpxyz[3];
473 param->GetXYZ(xyz);
474 param->GetPxPyPz(pxyz);
475 ltrp->GetXYZ(lxyz);
476 ltrp->GetPxPyPz(lpxyz);
477
478 if (fStreamLevel>0){
479 TTreeSRedirector *cstream = GetDebugStreamer();
480 Int_t time = fESD->GetTimeStamp();
481 Bool_t accept = AcceptLaser(id);
482 if (cstream){
483 (*cstream)<<"Track"<<
484 "run="<<fRun<<
485 "id="<<id<<
486 "accept="<<accept<<
487 "driftA.="<<fFitAside<<
488 "driftC.="<<fFitCside<<
489 "time="<<time<<
490 //
491 "LTr.="<<ltrp<<
492 "Esd.="<<track<<
493 "Tr.="<<param<<
494 "x0="<<xyz[0]<<
495 "x1="<<xyz[1]<<
496 "x2="<<xyz[2]<<
497 "px0="<<pxyz[0]<<
498 "px1="<<pxyz[1]<<
499 "px2="<<pxyz[2]<<
500 //
501 "lx0="<<lxyz[0]<<
502 "lx1="<<lxyz[1]<<
503 "lx2="<<lxyz[2]<<
504 "lpx0="<<lpxyz[0]<<
505 "lpx1="<<lpxyz[1]<<
506 "lpx2="<<lpxyz[2]<<
507 "\n";
508 }
509 }
510}
511
512void AliTPCcalibLaser::RefitLaserJW(Int_t id){
513 //
514 // Refit the track with different tracklet models:
515 // 1. Per ROC using the kalman filter, different edge cuts
516 // 2. Per ROC linear in y and z
517 // 3. Per ROC quadratic in y and z
518 // 4. Per track offset for each sector, linear for each sector, common quadratic
519 // store x, y, z information for all models and the cluster to calculate the residuals
520 //
521 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
522 AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
523 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
524
525 AliTPCclusterMI dummyCl;
526
527 //two tracklets
528 Int_t kMaxTracklets=2;
529 //=============================================//
530 // Linear Fitters for the Different Approaches //
531 //=============================================//
532 //linear fit model in y and z; inner - outer sector
533 static TLinearFitter fy1I(2,"hyp1");
534 static TLinearFitter fy1O(2,"hyp1");
535 static TLinearFitter fz1I(2,"hyp1");
536 static TLinearFitter fz1O(2,"hyp1");
537 //quadratic fit model in y and z; inner - sector
538 static TLinearFitter fy2I(3,"hyp2");
539 static TLinearFitter fy2O(3,"hyp2");
540 static TLinearFitter fz2I(3,"hyp2");
541 static TLinearFitter fz2O(3,"hyp2");
542 //common quadratic fit for IROC and OROC in y and z
543 static TLinearFitter fy4(5,"hyp4");
544 static TLinearFitter fz4(5,"hyp4");
545
546
547 //set standard cuts
548 if ( fNcuts==0 ){
549 fNcuts=1;
550 fEdgeXcuts[0]=4;
551 fEdgeYcuts[0]=3;
552 fNClCuts[0]=20;
553 }
554 //=============================//
555 // Loop over all Tracklet Cuts //
556 //=============================//
557 for (Int_t icut=0; icut<fNcuts; icut++){
558 AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
559 //cut parameters
560 Double_t edgeCutX = fEdgeXcuts[icut];
561 Double_t edgeCutY = fEdgeYcuts[icut];
562 Int_t nclCut = fNClCuts[icut];
563 //=========================//
564 // Parameters to calculate //
565 //=========================//
566 //fit parameter inner, outer tracklet and combined fit
567 TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
568 TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
569 //
570 TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
571 TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
572 TVectorD vecy4res(5),vecz4res(5);
573 // cluster and track positions for each row - used for residuals
574 TVectorD vecX(159); // x is the same for all (row center)
575 TVectorD vecYkalman(159); // y from kalman fit
576 TVectorD vecZkalman(159); // z from kalman fit
577 TVectorD vecY1(159); // y from pol1 fit per ROC
578 TVectorD vecZ1(159); // z from pol1 fit per ROC
579 TVectorD vecY2(159); // y from pol2 fit per ROC
580 TVectorD vecZ2(159); // z from pol2 fit per ROC
581 TVectorD vecY4(159); // y from sector fit
582 TVectorD vecZ4(159); // z from sector fit
583 TVectorD vecClY(159); // y cluster position
584 TVectorD vecClZ(159); // z cluster position
585 TVectorD vecSec(159); // sector for each row
586 //chi2 of fits
587 Double_t chi2I1z=-1; // chi2 of pol1 fit in z (inner)
588 Double_t chi2I1y=-1; // chi2 of pol1 fit in y (inner)
589 Double_t chi2O1z=-1; // chi2 of pol1 fit in z (outer)
590 Double_t chi2O1y=-1; // chi2 of pol1 fit in y (outer)
591 Double_t chi2I2z=-1; // chi2 of pol2 fit in z (inner)
592 Double_t chi2I2y=-1; // chi2 of pol2 fit in y (inner)
593 Double_t chi2O2z=-1; // chi2 of pol2 fit in z (outer)
594 Double_t chi2O2y=-1; // chi2 of pol2 fit in y (outer)
595 Double_t chi2IOz=-1; // chi2 of hyp4 fit in z (inner+outer)
596 Double_t chi2IOy=-1; // chi2 of hyp4 fit in y (inner+outer)
597 //more
598 Int_t innerSector = -1; // number of inner sector
599 Int_t outerSector = -1; // number of outer sector
600 Int_t nclI=0; // number of clusters (inner)
601 Int_t nclO=0; // number of clusters (outer)
602 //=================================================//
603 // Perform the Kalman Fit using the Tracklet Class //
604 //=================================================//
605 AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
606 TObjArray tracklets=
607 AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
608 kFALSE,nclCut,kMaxTracklets);
609 // tracklet pointers
610 AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
611 AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
612 AliTPCTracklet *tr=0x0;
613 AliTPCTracklet dummy;
614 //continue if we didn't find a tracklet
615 if ( !trInner && !trOuter ) continue;
616 //================================================================================//
617 // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
618 //================================================================================//
619 if ( trInner && trOuter ){
620 if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
621 if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
622 tr = trInner;
623 trInner=trOuter;
624 trOuter=tr;
625 AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
626 }
627 } else {
628 if ( trInner ){
629 if ( !trInner->GetInner() ) continue;
630 trOuter=&dummy;
631 if ( trInner->GetSector()>35 ){
632 trOuter=trInner;
633 trInner=&dummy;
634 }
635 } else { //trOuter
636 if ( !trOuter->GetInner() ) continue;
637 trInner=&dummy;
638 if ( trOuter->GetSector()<36 ){
639 trInner=trOuter;
640 trOuter=&dummy;
641 }
642 }
643 }
644 innerSector = trInner->GetSector();
645 if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
646 outerSector = trOuter->GetSector();
647 if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
648
649 // array of clusters
650 TClonesArray arrCl("AliTPCclusterMI",159);
651 arrCl.ExpandCreateFast(159);
652 //=======================================//
653 // fill fitters with cluster information //
654 //=======================================//
655 AliDebug(3,"Filing Cluster Information");
656 for (Int_t irow=158;irow>-1;--irow) {
657 AliTPCclusterMI *c=track->GetClusterPointer(irow);
658 AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
659 cl=dummyCl;
660 //
661 vecSec[irow]=-1;
662 if (!c) continue;
663 //
664 Int_t roc = static_cast<Int_t>(c->GetDetector());
665 if ( roc!=innerSector && roc!=outerSector ) continue;
666 vecSec[irow]=roc;
667 //store clusters in clones array
668 cl=*c;
669 //cluster position
670 vecX[irow] = c->GetX();
671 vecClY[irow] = c->GetY();
672 vecClZ[irow] = c->GetZ();
673 //
674 Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
675 Double_t y = vecClY[irow];
676 Double_t z = vecClZ[irow];
677 //
678 Double_t x2[2]={x,x*x}; //linear and parabolic parameters
679 Double_t x4[4]={0,0,0,0}; //hyp4 parameters
680 if ( roc == innerSector ) {
681 x4[0]=1; //offset inner - outer sector
682 x4[1]=x; //slope parameter inner sector
683 } else {
684 x4[2]=x; //slope parameter outer sector
685 }
686 x4[3]=x*x; //common parabolic shape
687 //
688 if ( roc==innerSector ){
689 fy1I.AddPoint(x2,y);
690 fz1I.AddPoint(x2,z);
691 fy2I.AddPoint(x2,y);
692 fz2I.AddPoint(x2,z);
693 ++nclI;
694 }
695 if ( roc==outerSector ){
696 fy1O.AddPoint(x2,y);
697 fz1O.AddPoint(x2,z);
698 fy2O.AddPoint(x2,y);
699 fz2O.AddPoint(x2,z);
700 ++nclO;
701 }
702 fy4.AddPoint(x4,y);
703 fz4.AddPoint(x4,z);
704 }
705 //======================================//
706 // Evaluate and retrieve fit parameters //
707 //======================================//
708 AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
709 //inner sector
710 if ( (innerSector>-1) && (fy1I.GetNpoints()>0) ){
711 fy1I.Eval();
712 fz1I.Eval();
713 fy2I.Eval();
714 fz2I.Eval();
715 fy1I.GetParameters(vecy1resInner);
716 fz1I.GetParameters(vecz1resInner);
717 fy2I.GetParameters(vecy2resInner);
718 fz2I.GetParameters(vecz2resInner);
719 chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
720 chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
721 chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
722 chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
723 }
724 //outer sector
725 if ( (outerSector>-1) && (fy1O.GetNpoints()>0) ){
726 fy1O.Eval();
727 fz1O.Eval();
728 fy2O.Eval();
729 fz2O.Eval();
730 fy1O.GetParameters(vecy1resOuter);
731 fz1O.GetParameters(vecz1resOuter);
732 fy2O.GetParameters(vecy2resOuter);
733 fz2O.GetParameters(vecz2resOuter);
734 chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
735 chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
736 chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
737 chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
738 }
739 //combined hyp4 fit
740 if ( innerSector>0 && outerSector>0 ){
741 if (fy4.GetNpoints()>0) {
742 fy4.Eval();
743 fy4.GetParameters(vecy4res);
744 chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
745 }
746 if (fz4.GetNpoints()>0) {
747 fz4.Eval();
748 fz4.GetParameters(vecz4res);
749 chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
750 }
751 }
752 //clear points
753 fy4.ClearPoints(); fz4.ClearPoints();
754 fy1I.ClearPoints(); fy1O.ClearPoints();
755 fz1I.ClearPoints(); fz1O.ClearPoints();
756 fy2I.ClearPoints(); fy2O.ClearPoints();
757 fz2I.ClearPoints(); fz2O.ClearPoints();
758 //==============================//
759 // calculate tracklet positions //
760 //==============================//
761 AliDebug(4,"Calculate tracklet positions");
762 for (Int_t irow=158;irow>-1;--irow) {
763 if ( vecSec[irow]==-1 ) continue; //no cluster info
764 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
765 tr=&dummy;
766 Double_t x=vecX[irow];
767 Double_t xref=x-133.4;
768 //
769 Double_t yoffInner=0;
770 Double_t zoffInner=0;
771 Double_t yslopeInner=0;
772 Double_t yslopeOuter=0;
773 Double_t zslopeInner=0;
774 Double_t zslopeOuter=0;
775 //positions of hyperplane fits
776 if ( vecSec[irow] == outerSector ) {
777 tr=trOuter;
778 vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
779 vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
780 vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
781 vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
782 yslopeOuter=vecy4res[3];
783 zslopeOuter=vecz4res[3];
784 } else {
785 tr=trInner;
786 vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
787 vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
788 vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
789 vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
790 yoffInner=vecy4res[1];
791 zoffInner=vecz4res[1];
792 yslopeInner=vecy4res[2];
793 zslopeInner=vecz4res[2];
794 }
795 vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
796 vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
797 //positions of kalman fits
798 Double_t gxyz[3],xyz[3];
799 AliExternalTrackParam *param = 0x0;
800 //
801 param=tr->GetInner();
802 if (param){
803 param->GetXYZ(gxyz);
804 Float_t bz = AliTracker::GetBz(gxyz);
805 param->GetYAt(x, bz, xyz[1]);
806 param->GetZAt(x, bz, xyz[2]);
807 vecYkalman[irow]=xyz[1];
808 vecZkalman[irow]=xyz[2];
809 }
810 }
811 //=====================================================================//
812 // write results from the different tracklet fits with debug streamers //
813 //=====================================================================//
814 if (fStreamLevel>4){
815 TTreeSRedirector *cstream = GetDebugStreamer();
816 if (cstream){
817 Float_t dedx = track->GetdEdx();
818 (*cstream)<<"FitModels"<<
819 "cutNr=" << icut <<
820 "edgeCutX=" << edgeCutX <<
821 "edgeCutY=" << edgeCutY <<
822 "nclCut=" << nclCut <<
823 "innerSector="<< innerSector <<
824 "outerSector="<< outerSector <<
825 "dEdx=" << dedx <<
826 "LTr.=" << ltrp <<
827 "Tr.=" << extparam <<
828 "yPol1In.=" << &vecy1resInner <<
829 "zPol1In.=" << &vecz1resInner <<
830 "yPol2In.=" << &vecy2resInner <<
831 "zPol2In.=" << &vecz2resInner <<
832 "yPol1Out.=" << &vecy1resOuter <<
833 "zPol1Out.=" << &vecz1resOuter <<
834 "yPol2Out.=" << &vecy2resOuter <<
835 "zPol2Out.=" << &vecz2resOuter <<
836 "yInOut.=" << &vecy4res <<
837 "zInOut.=" << &vecz4res <<
838 "chi2y1In=" << chi2I1y <<
839 "chi2z1In=" << chi2I1z <<
840 "chi2y1Out=" << chi2O1y <<
841 "chi2z1Out=" << chi2O1z <<
842 "chi2y2In=" << chi2I2y <<
843 "chi2z2In=" << chi2I2z <<
844 "chi2y2Out=" << chi2O2y <<
845 "chi2z2Out=" << chi2O2z <<
846 "chi2yInOut=" << chi2IOy <<
847 "chi2zInOut=" << chi2IOz <<
848 "trletIn.=" << trInner <<
849 "trletOut.=" << trOuter <<
850 "nclI=" << nclI <<
851 "nclO=" << nclO <<
852 "\n";
853 }
854 }
855
856 // wirte residuals information
857 if (fStreamLevel>5){
858 TTreeSRedirector *cstream = GetDebugStreamer();
859 if (cstream){
860 Float_t dedx = track->GetdEdx();
861 (*cstream)<<"Residuals"<<
862 "cutNr=" << icut <<
863 "edgeCutX=" << edgeCutX <<
864 "edgeCutY=" << edgeCutY <<
865 "nclCut=" << nclCut <<
866 "LTr.=" << ltrp <<
867 "Tr.=" << extparam<<
868 "dEdx=" << dedx <<
869 "Cl.=" << &arrCl <<
870 "TrX.=" << &vecX <<
871 "TrYpol1.=" << &vecY1 <<
872 "TrZpol1.=" << &vecZ1 <<
873 "TrYpol2.=" << &vecY2 <<
874 "TrZpol2.=" << &vecZ2 <<
875 "TrYInOut.=" << &vecY4 <<
876 "TrZInOut.=" << &vecZ4 <<
877 "ClY.=" << &vecClY <<
878 "ClZ.=" << &vecClZ <<
879 "sec.=" << &vecSec <<
880 "nclI=" << nclI <<
881 "nclO=" << nclO <<
882 "yInOut.=" << &vecy4res <<
883 "zInOut.=" << &vecz4res <<
884 "chi2y1In=" << chi2I1y <<
885 "chi2z1In=" << chi2I1z <<
886 "chi2y1Out=" << chi2O1y <<
887 "chi2z1Out=" << chi2O1z <<
888 "chi2y2In=" << chi2I2y <<
889 "chi2z2In=" << chi2I2z <<
890 "chi2y2Out=" << chi2O2y <<
891 "chi2z2Out=" << chi2O2z <<
892 "chi2yInOut=" << chi2IOy <<
893 "chi2zInOut=" << chi2IOz <<
894 "\n";
895
896 }
897 }
898 //==========================//
899 // Fill Residual Histograms //
900 //==========================//
901 TProfile *profy = (TProfile*)fDeltaYres.UncheckedAt(id);
902 TProfile *profz = (TProfile*)fDeltaZres.UncheckedAt(id);
903 if (!profy){
904 profy=new TProfile(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d",id),115,80,250);
905 profy->SetDirectory(0);
906 fDeltaYres.AddAt(profy,id);
907 }
908 if (!profz){
909 profz=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),115,80,250);
910 profz->SetDirectory(0);
911 fDeltaZres.AddAt(profz,id);
912 }
913 for (Int_t irow=158;irow>-1;--irow) {
914 if (vecSec[irow]==-1)continue; //no cluster info
915 Double_t x = vecX[irow];
916 Double_t ycl = vecClY[irow];
917 Double_t yfit = vecY1[irow];
918 Double_t zcl = vecClZ[irow];
919 Double_t zfit = vecZ1[irow];
920 if (profy)
921 if (profy->GetEntries()<1000000)
922 profy->Fill(x,yfit-ycl);
923 if (profz)
924 if (profz->GetEntries()<1000000)
925 profz->Fill(x,zfit-zcl);
926 }
927 //===============================//
928 // Fill Fit Parameter Histograms //
929 //===============================//
930 TH1F *pol2InnerY =(TH1F*)fPol2Par2InY.UncheckedAt(id);
931 TH1F *diff1InnerY=(TH1F*)fDiffPar1InY.UncheckedAt(id);
932 TH1F *pol2OuterY =(TH1F*)fPol2Par2OutY.UncheckedAt(id);
933 TH1F *diff1OuterY=(TH1F*)fDiffPar1OutY.UncheckedAt(id);
934 TH1F *pol2InnerZ =(TH1F*)fPol2Par2InZ.UncheckedAt(id);
935 TH1F *diff1InnerZ=(TH1F*)fDiffPar1InZ.UncheckedAt(id);
936 TH1F *pol2OuterZ =(TH1F*)fPol2Par2OutZ.UncheckedAt(id);
937 TH1F *diff1OuterZ=(TH1F*)fDiffPar1OutZ.UncheckedAt(id);
938 //create histograms if the do not already exist
939 if (!pol2InnerY){
940 pol2InnerY =new TH1F(Form("pol2par2inY%03d",id),
941 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (inner sector)",id),
942 500,-.005,.005);
943 pol2InnerY->SetDirectory(0);
944 pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
945 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
946 500,0.01,.01);
947 pol2OuterY->SetDirectory(0);
948
949 diff1InnerY=new TH1F(Form("diff1inY%03d",id),
950 Form("diff of 1st derivative from pol1 and pol2 in Y fit for Laser beam %03d (inner sector)",id),
951 500,-.5,.5);
952 diff1InnerY->SetDirectory(0);
953
954 diff1OuterY=new TH1F(Form("diff1outY%03d",id),
955 Form("diff of 1st derivative from pol1 and pol2 in Yfit for Laser beam %03d (outer sector)",id),
956 500,-1,1);
957 diff1OuterY->SetDirectory(0);
958
959 pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
960 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
961 500,-.002,.002);
962 pol2InnerZ->SetDirectory(0);
963
964 pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
965 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
966 500,-.005,.005);
967 pol2OuterZ->SetDirectory(0);
968 diff1InnerZ=new TH1F(Form("diff1inZ%03d",id),
969 Form("diff of 1st derivative from pol1 and pol2 in Z fit for Laser beam %03d (inner sector)",id),
970 500,-.02,.02);
971 diff1InnerZ->SetDirectory(0);
972 diff1OuterZ=new TH1F(Form("diff1outZ%03d",id),
973 Form("diff of 1st derivative from pol1 and pol2 in Zfit for Laser beam %03d (outer sector)",id),
974 500,-.03,.03);
975 diff1OuterZ->SetDirectory(0);
976 //add
977 fPol2Par2InY.AddAt(pol2InnerY,id);
978 fDiffPar1InY.AddAt(diff1InnerY,id);
979 fPol2Par2OutY.AddAt(pol2OuterY,id);
980 fDiffPar1OutY.AddAt(diff1OuterY,id);
981 fPol2Par2InZ.AddAt(pol2InnerZ,id);
982 fDiffPar1InZ.AddAt(diff1InnerZ,id);
983 fPol2Par2OutZ.AddAt(pol2OuterZ,id);
984 fDiffPar1OutZ.AddAt(diff1OuterZ,id);
985 }
986 //fill histograms
987 pol2InnerY ->Fill(vecy2resInner[2]);
988 pol2OuterY ->Fill(vecy2resOuter[2]);
989 diff1InnerY->Fill(vecy2resInner[1]-vecy1resInner[1]);
990 diff1OuterY->Fill(vecy2resOuter[1]-vecy1resOuter[1]);
991 pol2InnerZ ->Fill(vecz2resInner[2]);
992 pol2OuterZ ->Fill(vecz2resOuter[2]);
993 diff1InnerZ->Fill(vecz2resInner[1]-vecz1resInner[1]);
994 diff1OuterZ->Fill(vecz2resOuter[1]-vecz1resOuter[1]);
995
996
997
998 }
999}
1000
1001
1002void AliTPCcalibLaser::RefitLaser(Int_t id){
1003 //
1004 // Refit the track store residuals
1005 //
1006
1007 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
1008 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1009 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1010
1011 //linear fit model in y and z per sector
1012 static TLinearFitter fy1(2,"hyp1");
1013 static TLinearFitter fz1(2,"hyp1");
1014 //quadratic fit model in y and z per sector
1015 static TLinearFitter fy2(3,"hyp2");
1016 static TLinearFitter fz2(3,"hyp2");
1017 static TVectorD vecy2,vecz2,vecy1,vecz1;
1018
1019 const Int_t kMinClusters=20;
1020 Int_t nclusters[72];
1021 //
1022 for (Int_t i=0;i<72;++i) nclusters[i]=0;
1023
1024 for (Int_t i=0;i<160;++i) {
1025 AliTPCclusterMI *c=track->GetClusterPointer(i);
1026 if (c) nclusters[c->GetDetector()]++;
1027 }
1028
1029 for (Int_t isec=0; isec<72;isec++){
1030 if (nclusters[isec]<kMinClusters) continue;
1031 fy2.ClearPoints();
1032 fz2.ClearPoints();
1033 fy1.ClearPoints();
1034 fz1.ClearPoints();
1035 //
1036 for (Int_t irow=0;irow<160;++irow) {
1037 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1038 //if (c && RejectCluster(c)) continue;
1039 Double_t xd = c->GetX()-133.4; // reference x is beteen iroc and oroc
1040 if (c&&c->GetDetector()==isec) {
1041 Double_t x[2]={xd,xd*xd};
1042 fy2.AddPoint(x,c->GetY());
1043 fz2.AddPoint(x,c->GetZ());
1044 //
1045 fy1.AddPoint(x,c->GetY());
1046 fz1.AddPoint(x,c->GetZ());
1047 }
1048 }
1049 fy2.Eval();
1050 fz2.Eval();
1051 fy1.Eval();
1052 fz1.Eval();
1053 fy1.GetParameters(vecy1);
1054 fy2.GetParameters(vecy2);
1055 fz1.GetParameters(vecz1);
1056 fz2.GetParameters(vecz2);
1057
1058 if (fStreamLevel>0){
1059 TTreeSRedirector *cstream = GetDebugStreamer();
1060 if (cstream){
1061 Float_t dedx = track->GetdEdx();
1062 (*cstream)<<"Tracklet"<<
1063 "LTr.="<<ltrp<<
1064 "Tr.="<<param<<
1065 "sec="<<isec<<
1066 "ncl="<<nclusters[isec]<<
1067 "dedx="<<dedx<<
1068 "dedx="<<dedx<<
1069 "vy1.="<<&vecy1<<
1070 "vy2.="<<&vecy2<<
1071 "vz1.="<<&vecz1<<
1072 "vz2.="<<&vecz2<<
1073 "\n";
1074 }
1075 }
1076 }
1077 //
1078 //
1079 //
1080 // for (Int_t irow=0;irow<160;++irow) {
1081 // AliTPCclusterMI *c=track->GetClusterPointer(irow);
1082 // if (c && RejectCluster(c)) continue;
1083 // if (c&&c->GetDetector()==isec) {
1084 // Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
1085 // fy2.AddPoint(&x,c->GetY());
1086 // fz2.AddPoint(&x,c->GetZ());
1087 // //
1088 // fy1.AddPoint(&x,c->GetY());
1089 // fz1.AddPoint(&x,c->GetZ());
1090 // }
1091 // }
1092
1093}
1094
1095
1096void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield,Int_t minEntries){
1097 //
1098 // Dump information about laser beams
1099 // isOK variable indicates usability of the beam
1100 // Beam is not usable if:
1101 // a. No entries in range (krmsCut0)
1102 // b. Big sperad (krmscut1)
1103 // c. RMSto fit sigma bigger then (kmultiCut)
1104 // d. Too big angular spread
1105 //
1106
1107 const Float_t krmsCut0=0.001;
1108 const Float_t krmsCut1=0.16;
1109 const Float_t kmultiCut=2;
1110 const Float_t kcutP0=0.002;
1111 //
1112 AliTPCcalibLaser *laser = this;
1113 TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
1114 TF1 fg("fg","gaus");
1115 //
1116 //
1117 for (Int_t id=0; id<336; id++){
1118 Bool_t isOK=kTRUE;
1119 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
1120 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
1121 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
1122 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
1123 if (!hisphi) continue;;
1124 Double_t entries = hisphi->GetEntries();
1125 if (entries<minEntries) continue;
1126 //
1127 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1128 if (!ltrp) {
1129 AliTPCLaserTrack::LoadTracks();
1130 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1131 }
1132 Float_t meanphi = hisphi->GetMean();
1133 Float_t rmsphi = hisphi->GetRMS();
1134 //
1135 Float_t meanphiP = hisphiP->GetMean();
1136 Float_t rmsphiP = hisphiP->GetRMS();
1137 Float_t meanZ = hisZ->GetMean();
1138 Float_t rmsZ = hisZ->GetRMS();
1139 hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
1140 Double_t gphi1 = fg.GetParameter(1);
1141 Double_t gphi2 = fg.GetParameter(2);
1142 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
1143 Double_t gphiP1 = fg.GetParameter(1);
1144 Double_t gphiP2 = fg.GetParameter(2);
1145 hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
1146 Double_t gz1 = fg.GetParameter(1);
1147 Double_t gz2 = fg.GetParameter(2);
1148 //
1149 Float_t meanS=hisS->GetMean();
1150 //
1151 Double_t lxyz[3];
1152 Double_t lpxyz[3];
1153 ltrp->GetXYZ(lxyz);
1154 ltrp->GetPxPyPz(lpxyz);
1155
1156 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
1157 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
1158 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
1159 if (gphiP2>kcutP0) isOK=kFALSE;
1160 //
1161 (*pcstream)<<"Mean"<<
1162 "isOK="<<isOK<<
1163 "entries="<<entries<< // number of entries
1164 "bz="<<bfield<< // bfield
1165 "LTr.="<<ltrp<< // refernece track
1166 //
1167 "lx0="<<lxyz[0]<< // reference x
1168 "lx1="<<lxyz[1]<< // reference y
1169 "lx2="<<lxyz[2]<< // refernece z
1170 "lpx0="<<lpxyz[0]<< // reference x
1171 "lpx1="<<lpxyz[1]<< // reference y
1172 "lpx2="<<lpxyz[2]<< // refernece z
1173 //
1174 "msig="<<meanS<<
1175 //
1176 "mphi="<<meanphi<< //
1177 "rmsphi="<<rmsphi<< //
1178 "gphi1="<<gphi1<<
1179 "gphi2="<<gphi2<<
1180 //
1181 "mphiP="<<meanphiP<< //
1182 "rmsphiP="<<rmsphiP<< //
1183 "gphiP1="<<gphiP1<<
1184 "gphiP2="<<gphiP2<<
1185 //
1186 "meanZ="<<meanZ<<
1187 "rmsZ="<<rmsZ<<
1188 "gz1="<<gz1<<
1189 "gz2="<<gz2<<
1190
1191 "\n";
1192 }
1193 delete pcstream;
1194}
1195
1196
1197
1198void AliTPCcalibLaser::DumpScanInfo(TTree * chain){
1199 //
1200 //
1201 //
1202 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
1203 TFile * f = pcstream->GetFile();
1204 f->mkdir("dirphi");
1205 f->mkdir("dirphiP");
1206 f->mkdir("dirZ");
1207 TF1 fp("p1","pol1");
1208 //
1209 //
1210 char cut[1000];
1211 char grname[1000];
1212 char grnamefull[1000];
1213
1214 Double_t mphi[100];
1215 Double_t mphiP[100];
1216 Double_t smphi[100];
1217 Double_t smphiP[100];
1218 Double_t mZ[100];
1219 Double_t smZ[100];
1220 Double_t bz[100];
1221 Double_t sbz[100];
1222 // fit parameters
1223 Double_t pphi[3];
1224 Double_t pphiP[3];
1225 Double_t pmZ[3];
1226 //
1227 for (Int_t id=0; id<336; id++){
1228 // id =205;
1229 sprintf(cut,"isOK&&fId==%d",id);
1230 Int_t entries = chain->Draw("bz",cut,"goff");
1231 if (entries<3) continue;
1232 AliTPCLaserTrack *ltrp = 0;;
1233 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
1234 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1235 Double_t lxyz[3];
1236 Double_t lpxyz[3];
1237 ltrp->GetXYZ(lxyz);
1238 ltrp->GetPxPyPz(lpxyz);
1239
1240 chain->Draw("bz",cut,"goff");
1241 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
1242 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
1243 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
1244 //
1245 chain->Draw("gphi1",cut,"goff");
1246 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
1247 chain->Draw("0.05*abs(mphi)+gphi2",cut,"goff");
1248 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
1249 //
1250 chain->Draw("gphiP1",cut,"goff");
1251 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
1252 chain->Draw("0.05*abs(mphiP)+gphiP2",cut,"goff");
1253 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
1254 //
1255 chain->Draw("gz1",cut,"goff");
1256 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
1257 chain->Draw("0.01*abs(meanZ)+gz2",cut,"goff");
1258 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
1259 //
1260 //
1261 sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
1262 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
1263 // store data
1264 // phi
1265 f->cd("dirphi");
1266 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
1267 grphi->Draw("a*");
1268 grphi->Fit(&fp);
1269 pphi[0] = fp.GetParameter(0); // offset
1270 pphi[1] = fp.GetParameter(1); // slope
1271 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1272 sprintf(grname,"phi_id%d",id);
1273 grphi->SetName(grname); grphi->SetTitle(grnamefull);
1274 grphi->GetXaxis()->SetTitle("b_{z} (T)");
1275 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
1276 grphi->SetMaximum(1.2);
1277 grphi->SetMinimum(-1.2);
1278 grphi->Draw("a*");
1279
1280 grphi->Write();
1281 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
1282 // phiP
1283 f->cd("dirphiP)");
1284 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
1285 grphiP->Draw("a*");
1286 grphiP->Fit(&fp);
1287 pphiP[0] = fp.GetParameter(0); // offset
1288 pphiP[1] = fp.GetParameter(1); // slope
1289 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1290 sprintf(grname,"phiP_id%d",id);
1291 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
1292 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
1293 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
1294 grphiP->SetMaximum(pphiP[0]+0.005);
1295 grphiP->SetMinimum(pphiP[0]-0.005);
1296
1297 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
1298 grphiP->Write();
1299 //
1300 //Z
1301 f->cd("dirZ");
1302 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
1303 grmZ->Draw("a*");
1304 grmZ->Fit(&fp);
1305 pmZ[0] = fp.GetParameter(0); // offset
1306 pmZ[1] = fp.GetParameter(1); // slope
1307 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1308 sprintf(grname,"mZ_id%d",id);
1309 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
1310 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
1311 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
1312
1313 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
1314 grmZ->Write();
1315
1316
1317 for (Int_t ientry=0; ientry<entries; ientry++){
1318 (*pcstream)<<"Mean"<<
1319 "id="<<id<<
1320 "LTr.="<<ltrp<<
1321 "entries="<<entries<<
1322 "bz="<<bz[ientry]<<
1323 "lx0="<<lxyz[0]<< // reference x
1324 "lx1="<<lxyz[1]<< // reference y
1325 "lx2="<<lxyz[2]<< // refernece z
1326 "lpx0="<<lpxyz[0]<< // reference x
1327 "lpx1="<<lpxyz[1]<< // reference y
1328 "lpx2="<<lpxyz[2]<< // refernece z
1329 //values
1330 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
1331 "pphi0="<<pphi[0]<< // offset
1332 "pphi1="<<pphi[1]<< // mean
1333 "pphi2="<<pphi[2]<< // norm chi2
1334 //
1335 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
1336 "pphiP0="<<pphiP[0]<< // offset
1337 "pphiP1="<<pphiP[1]<< // mean
1338 "pphiP2="<<pphiP[2]<< // norm chi2
1339 //
1340 "gz1="<<mZ[ientry]<<
1341 "pmZ0="<<pmZ[0]<< // offset
1342 "pmZ1="<<pmZ[1]<< // mean
1343 "pmZ2="<<pmZ[2]<< // norm chi2
1344 "\n";
1345 }
1346 }
1347
1348 delete pcstream;
1349
1350}
1351
1352
1353void AliTPCcalibLaser::Analyze(){
1354 //
1355 //
1356 //
1357}
1358
1359
1360Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
1361
1362 TIterator* iter = li->MakeIterator();
1363 AliTPCcalibLaser* cal = 0;
1364
1365 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
1366 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
1367 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1368 return -1;
1369 }
1370 //
1371 // fHistNTracks->Add(cal->fHistNTracks);
1372// fClusters->Add(cal->fClusters);
1373// fModules->Add(cal->fModules);
1374// fHistPt->Add(cal->fHistPt);
1375// fPtResolution->Add(cal->fPtResolution);
1376// fDeDx->Add(cal->fDeDx);
1377
1378
1379 TH1F *h=0x0;
1380 TH1F *hm=0x0;
1381 TProfile *hp=0x0;
1382 TProfile *hpm=0x0;
1383
1384 for (Int_t id=0; id<336; id++){
1385 // merge fDeltaZ histograms
1386 hm = (TH1F*)cal->fDeltaZ.At(id);
1387 h = (TH1F*)fDeltaZ.At(id);
1388 if (!h) {
1389 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
1390 h->SetDirectory(0);
1391 fDeltaZ.AddAt(h,id);
1392 }
1393 if (hm) h->Add(hm);
1394 // merge fDeltaPhi histograms
1395 hm = (TH1F*)cal->fDeltaPhi.At(id);
1396 h = (TH1F*)fDeltaPhi.At(id);
1397 if (!h) {
1398 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
1399 h->SetDirectory(0);
1400 fDeltaPhi.AddAt(h,id);
1401 }
1402 if (hm) h->Add(hm);
1403 // merge fDeltaPhiP histograms
1404 hm = (TH1F*)cal->fDeltaPhiP.At(id);
1405 h = (TH1F*)fDeltaPhiP.At(id);
1406 if (!h) {
1407 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
1408 h->SetDirectory(0);
1409 fDeltaPhiP.AddAt(h,id);
1410 }
1411 if (hm) h->Add(hm);
1412 // merge fSignals histograms
1413 hm = (TH1F*)cal->fSignals.At(id);
1414 h = (TH1F*)fSignals.At(id);
1415 if (!h) {
1416 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
1417 h->SetDirectory(0);
1418 fSignals.AddAt(h,id);
1419 }
1420 if (hm) h->Add(hm);
1421 //
1422 //
1423 // merge ProfileY histograms
1424 hpm = (TProfile*)cal->fDeltaYres.At(id);
1425 hp = (TProfile*)fDeltaYres.At(id);
1426 if (!hp) {
1427 hp=new TProfile(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d",id),115,80,250);
1428 hp->SetDirectory(0);
1429 fDeltaYres.AddAt(hp,id);
1430 }
1431 if (hpm) hp->Add(hpm);
1432 //
1433 hpm = (TProfile*)cal->fDeltaZres.At(id);
1434 hp = (TProfile*)fDeltaZres.At(id);
1435 if (!hp) {
1436 hp=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),115,80,250);
1437 hp->SetDirectory(0);
1438 fDeltaZres.AddAt(hp,id);
1439 }
1440 if (hpm) hp->Add(hpm);
1441 //
1442 //
1443 //merge fit param histograms
1444 //local hists
1445 TH1F *pol2InnerY =(TH1F*)fPol2Par2InY.UncheckedAt(id);
1446 TH1F *diff1InnerY=(TH1F*)fDiffPar1InY.UncheckedAt(id);
1447 TH1F *pol2OuterY =(TH1F*)fPol2Par2OutY.UncheckedAt(id);
1448 TH1F *diff1OuterY=(TH1F*)fDiffPar1OutY.UncheckedAt(id);
1449 TH1F *pol2InnerZ =(TH1F*)fPol2Par2InZ.UncheckedAt(id);
1450 TH1F *diff1InnerZ=(TH1F*)fDiffPar1InZ.UncheckedAt(id);
1451 TH1F *pol2OuterZ =(TH1F*)fPol2Par2OutZ.UncheckedAt(id);
1452 TH1F *diff1OuterZ=(TH1F*)fDiffPar1OutZ.UncheckedAt(id);
1453 //hists to merge
1454 TH1F *pol2InnerYm =(TH1F*)cal->fPol2Par2InY.UncheckedAt(id);
1455 TH1F *diff1InnerYm=(TH1F*)cal->fDiffPar1InY.UncheckedAt(id);
1456 TH1F *pol2OuterYm =(TH1F*)cal->fPol2Par2OutY.UncheckedAt(id);
1457 TH1F *diff1OuterYm=(TH1F*)cal->fDiffPar1OutY.UncheckedAt(id);
1458 TH1F *pol2InnerZm =(TH1F*)cal->fPol2Par2InZ.UncheckedAt(id);
1459 TH1F *diff1InnerZm=(TH1F*)cal->fDiffPar1InZ.UncheckedAt(id);
1460 TH1F *pol2OuterZm =(TH1F*)cal->fPol2Par2OutZ.UncheckedAt(id);
1461 TH1F *diff1OuterZm=(TH1F*)cal->fDiffPar1OutZ.UncheckedAt(id);
1462 //create histos if they do not exist
1463 if (!pol2InnerY){
1464 pol2InnerY =new TH1F(Form("pol2par2inY%03d",id),
1465 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (inner sector)",id),
1466 500,-.005,.005);
1467 pol2InnerY->SetDirectory(0);
1468
1469 pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
1470 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
1471 500,0.01,.01);
1472 pol2OuterY->SetDirectory(0);
1473 diff1InnerY=new TH1F(Form("diff1inY%03d",id),
1474 Form("diff of 1st derivative from pol1 and pol2 in Y fit for Laser beam %03d (inner sector)",id),
1475 500,-.5,.5);
1476 diff1OuterY=new TH1F(Form("diff1outY%03d",id),
1477 Form("diff of 1st derivative from pol1 and pol2 in Yfit for Laser beam %03d (outer sector)",id),
1478 500,-1,1);
1479 diff1InnerY->SetDirectory(0);
1480
1481
1482 pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
1483 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
1484 500,-.002,.002);
1485 pol2InnerZ->SetDirectory(0);
1486
1487 pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
1488 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
1489 500,-.005,.005);
1490 pol2OuterZ->SetDirectory(0);
1491 diff1InnerZ=new TH1F(Form("diff1inZ%03d",id),
1492 Form("diff of 1st derivative from pol1 and pol2 in Z fit for Laser beam %03d (inner sector)",id),
1493 500,-.02,.02);
1494 diff1InnerZ->SetDirectory(0);
1495 diff1OuterZ=new TH1F(Form("diff1outZ%03d",id),
1496 Form("diff of 1st derivative from pol1 and pol2 in Zfit for Laser beam %03d (outer sector)",id),
1497 500,-.03,.03);
1498 diff1OuterZ->SetDirectory(0);
1499 }
1500 pol2InnerYm =(TH1F*)cal->fPol2Par2InY.UncheckedAt(id);
1501 diff1InnerYm=(TH1F*)cal->fDiffPar1InY.UncheckedAt(id);
1502 pol2OuterYm =(TH1F*)cal->fPol2Par2OutY.UncheckedAt(id);
1503 diff1OuterYm=(TH1F*)cal->fDiffPar1OutY.UncheckedAt(id);
1504 pol2InnerZm =(TH1F*)cal->fPol2Par2InZ.UncheckedAt(id);
1505 diff1InnerZm=(TH1F*)cal->fDiffPar1InZ.UncheckedAt(id);
1506 pol2OuterZm =(TH1F*)cal->fPol2Par2OutZ.UncheckedAt(id);
1507 diff1OuterZm=(TH1F*)cal->fDiffPar1OutZ.UncheckedAt(id);
1508 }
1509 }
1510 return 0;
1511}
1512
1513
1514
1515
1516/*
1517 gSystem->Load("libSTAT.so")
1518 TStatToolkit toolkit;
1519 Double_t chi2;
1520 TVectorD fitParam;
1521 TMatrixD covMatrix;
1522 Int_t npoints;
1523 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
1524
1525
1526TString fstring="";
1527//
1528fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
1529fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
1530fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
1531fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
1532//
1533fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
1534fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
1535fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
1536fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
1537//
1538fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
1539fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
1540fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
1541fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
1542
1543
1544
1545
1546 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
1547
1548 treeT->SetAlias("fit",strq0->Data());
1549
1550
1551 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
1552
1553 treeT->SetAlias("fitP",strqP->Data());
1554
1555
1556 TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
1557 treeT->SetAlias("fitD",strqDrift->Data());
1558
1559
1560treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
1561{
1562for (Int_t i=0; i<6;i++){
1563treeT->SetLineColor(i+2);
1564treeT->SetMarkerSize(1);
1565treeT->SetMarkerStyle(22+i);
1566treeT->SetMarkerColor(i+2);
1567
1568treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
1569}
1570}
1571 */
1572
1573
1574
1575/*
1576 TTree * tree = (TTree*)f.Get("FitModels");
1577
1578 TEventList listLFit0("listLFit0","listLFit0");
1579 TEventList listLFit1("listLFit1","listLFit1");
1580 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
1581 tree->SetEventList(&listLFit0);
1582
1583
1584
1585
1586 gSystem->Load("libSTAT.so")
1587 TStatToolkit toolkit;
1588 Double_t chi2;
1589 TVectorD fitParam;
1590 TMatrixD covMatrix;
1591 Int_t npoints;
1592
1593 chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
1594 chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
1595
1596
1597 TString fstring="";
1598 fstring+="cos(dp)++";
1599 fstring+="sin(dp)++";
1600 fstring+="cos(dt)++";
1601 fstring+="sin(dt)++";
1602
1603 TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
1604
1605
1606*/