]> 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
Adding histogram
[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)
263d466a 25 //
2b35e8e6 26 // 3. The drift velocity and jitter is calculated event by event
263d466a 27 // (see function drift velocity)
2b35e8e6 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+")
0728a4f6 57 AliXRDPROOFtoolkit tool;
2b35e8e6 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")
95a0e09b 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
c6914c83 72*/
73
74
75
76#include "TLinearFitter.h"
77#include "AliTPCcalibLaser.h"
78#include "AliExternalTrackParam.h"
e9f38a4b 79#include "AliESDEvent.h"
80#include "AliESDfriend.h"
c6914c83 81#include "AliESDtrack.h"
82#include "AliTPCTracklet.h"
83#include "TH1D.h"
ecb7e705 84#include "TH1F.h"
95a0e09b 85#include "TProfile.h"
c6914c83 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"
263d466a 94#include "AliLog.h"
c6914c83 95#include "TClonesArray.h"
1fd56785 96#include "TPad.h"
c6914c83 97
98
99#include "TTreeStream.h"
100#include <iostream>
101#include <sstream>
7b18d067 102#include "AliTPCLaserTrack.h"
103
c6914c83 104using namespace std;
105
106ClassImp(AliTPCcalibLaser)
107
108AliTPCcalibLaser::AliTPCcalibLaser():
e9f38a4b 109 AliTPCcalibBase(),
110 fESD(0),
111 fESDfriend(),
c18f4385 112 fTracksMirror(336),
113 fTracksEsd(336),
114 fTracksEsdParam(336),
115 fTracksTPC(336),
263d466a 116 fDeltaZ(336),
117 fDeltaPhi(336),
118 fDeltaPhiP(336),
95a0e09b 119 fSignals(336),
120 fDeltaYres(336),
121 fDeltaZres(336),
ecb7e705 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)),
263d466a 131 fFitCside(new TVectorD(3)),
132 fEdgeXcuts(5),
133 fEdgeYcuts(5),
134 fNClCuts(5),
135 fNcuts(0),
136 fRun(0),
137 fEvent(0)
c6914c83 138{
139 //
140 // Constructor
141 //
e9f38a4b 142 fTracksEsdParam.SetOwner(kTRUE);
c6914c83 143}
144
145AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
e9f38a4b 146 AliTPCcalibBase(),
147 fESD(0),
148 fESDfriend(0),
c18f4385 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
2b35e8e6 156 fSignals(336), // array of dedx signals
95a0e09b 157 fDeltaYres(336),
158 fDeltaZres(336),
ecb7e705 159 fPol2Par2InY(336),
160 fDiffPar1InY(336),
161 fPol2Par2OutY(336),
162 fDiffPar1OutY(336),
163 fPol2Par2InZ(336),
164 fDiffPar1InZ(336),
165 fPol2Par2OutZ(336),
166 fDiffPar1OutZ(336),
c18f4385 167 fFitAside(new TVectorD(3)), // drift fit - A side
168 fFitCside(new TVectorD(3)), // drift fit - C- side
263d466a 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)
c6914c83 175{
176 SetName(name);
177 SetTitle(title);
178 //
179 // Constructor
180 //
263d466a 181 fTracksEsdParam.SetOwner(kTRUE);
c6914c83 182}
183
184AliTPCcalibLaser::~AliTPCcalibLaser() {
185 //
186 // destructor
187 //
188}
189
e9f38a4b 190
191
192void AliTPCcalibLaser::Process(AliESDEvent * event) {
c6914c83 193 //
263d466a 194 //
e9f38a4b 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 }
263d466a 205 AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
e9f38a4b 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 }
263d466a 224
e9f38a4b 225 FitDriftV();
c18f4385 226 MakeDistHisto();
e9f38a4b 227 //
c18f4385 228 for (Int_t id=0; id<336; id++){
e9f38a4b 229 //
230 //
231 if (!fTracksEsdParam.At(id)) continue;
232 DumpLaser(id);
263d466a 233// RefitLaser(id);
234 RefitLaserJW(id);
235
c18f4385 236 }
263d466a 237// fEvent++;
c18f4385 238}
239
240void AliTPCcalibLaser::MakeDistHisto(){
241 //
242 //
243 //
244 for (Int_t id=0; id<336; id++){
245 //
246 //
263d466a 247 if (!fTracksEsdParam.At(id)) continue;
c18f4385 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);
2b35e8e6 254 TH1F * hisSignal = (TH1F*)fSignals.At(id);
255
263d466a 256 if (!hisdz){
c03e3250 257 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
258 hisdz->SetDirectory(0);
c18f4385 259 fDeltaZ.AddAt(hisdz,id);
260 //
c03e3250 261 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
262 hisdphi->SetDirectory(0);
c18f4385 263 fDeltaPhi.AddAt(hisdphi,id);
264 //
c03e3250 265 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
266 hisdphiP->SetDirectory(0);
c18f4385 267 fDeltaPhiP.AddAt(hisdphiP,id);
c03e3250 268 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
269 hisSignal->SetDirectory(0);
2b35e8e6 270 fSignals.AddAt(hisSignal,id);
c18f4385 271 }
272
273 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
274 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
2b35e8e6 275 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
c03e3250 276 if (!param) return;
277 if (!ltrp) return;
278 if (!track) return;
c18f4385 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];
c03e3250 291 if (hisdz) hisdz->Fill(dz);
292 if (hisdphi) hisdphi->Fill(dphi);
263d466a 293 if (hisdphiP) hisdphiP->Fill(dphiP);
c03e3250 294 if (hisSignal) hisSignal->Fill(track->GetTPCsignal());
c18f4385 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++){
263d466a 308 if (!fTracksEsdParam.At(id)) continue;
c18f4385 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();
e9f38a4b 345 }
263d466a 346
c18f4385 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;
263d466a 382 if (TMath::Abs(param->GetParameter()[3])>0.06) return kFALSE;
c18f4385 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 //
263d466a 391 return kTRUE;
e9f38a4b 392}
393
394Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
395 //
396 // Find corresponding mirror
397 // add the corresponding tracks
c6914c83 398 //
1bd37cc0 399 Float_t kRadius0 = 252;
c18f4385 400 Float_t kRadius = 253.4;
e9f38a4b 401 if (!track->GetOuterParam()) return -1;
c6914c83 402 AliExternalTrackParam param(*(track->GetOuterParam()));
1bd37cc0 403 AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
404 AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
7b18d067 405 AliTPCLaserTrack ltr;
406 AliTPCLaserTrack *ltrp=0x0;
263d466a 407 AliTPCLaserTrack *ltrpjw=0x0;
1bd37cc0 408 //
263d466a 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)))
1bd37cc0 414 ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
263d466a 415 else
1bd37cc0 416 ltrp=&ltr;
263d466a 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
e9f38a4b 441 if (id>=0){
c18f4385 442 //
443 //
444 Float_t radius=TMath::Abs(ltrp->GetX());
445 AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
446 //
e9f38a4b 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);
c18f4385 451 //
e9f38a4b 452 }
c18f4385 453
e9f38a4b 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);
1bd37cc0 466 //
467 // Fast laser ID
468 //
469 Double_t xyz[3];
470 Double_t pxyz[3];
c18f4385 471 Double_t lxyz[3];
472 Double_t lpxyz[3];
e9f38a4b 473 param->GetXYZ(xyz);
474 param->GetPxPyPz(pxyz);
c18f4385 475 ltrp->GetXYZ(lxyz);
476 ltrp->GetPxPyPz(lpxyz);
477
c6914c83 478 if (fStreamLevel>0){
479 TTreeSRedirector *cstream = GetDebugStreamer();
c18f4385 480 Int_t time = fESD->GetTimeStamp();
481 Bool_t accept = AcceptLaser(id);
c6914c83 482 if (cstream){
483 (*cstream)<<"Track"<<
e9f38a4b 484 "run="<<fRun<<
7b18d067 485 "id="<<id<<
c18f4385 486 "accept="<<accept<<
487 "driftA.="<<fFitAside<<
488 "driftC.="<<fFitCside<<
489 "time="<<time<<
7b18d067 490 //
491 "LTr.="<<ltrp<<
492 "Esd.="<<track<<
e9f38a4b 493 "Tr.="<<param<<
c6914c83 494 "x0="<<xyz[0]<<
495 "x1="<<xyz[1]<<
496 "x2="<<xyz[2]<<
497 "px0="<<pxyz[0]<<
498 "px1="<<pxyz[1]<<
499 "px2="<<pxyz[2]<<
c18f4385 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]<<
c6914c83 507 "\n";
508 }
509 }
510}
511
263d466a 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;
0728a4f6 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");
263d466a 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 }
0728a4f6 554 //=============================//
555 // Loop over all Tracklet Cuts //
556 //=============================//
263d466a 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];
0728a4f6 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
263d466a 569 //
0728a4f6 570 TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
571 TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
263d466a 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
0728a4f6 586 //chi2 of fits
263d466a 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)
0728a4f6 597 //more
598 Int_t innerSector = -1; // number of inner sector
599 Int_t outerSector = -1; // number of outer sector
263d466a 600 Int_t nclI=0; // number of clusters (inner)
601 Int_t nclO=0; // number of clusters (outer)
0728a4f6 602 //=================================================//
603 // Perform the Kalman Fit using the Tracklet Class //
604 //=================================================//
263d466a 605 AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
606 TObjArray tracklets=
607 AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
608 kFALSE,nclCut,kMaxTracklets);
263d466a 609 // tracklet pointers
263d466a 610 AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
611 AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
0728a4f6 612 AliTPCTracklet *tr=0x0;
613 AliTPCTracklet dummy;
614 //continue if we didn't find a tracklet
263d466a 615 if ( !trInner && !trOuter ) continue;
0728a4f6 616 //================================================================================//
617 // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
618 //================================================================================//
263d466a 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
263d466a 649 // array of clusters
650 TClonesArray arrCl("AliTPCclusterMI",159);
651 arrCl.ExpandCreateFast(159);
0728a4f6 652 //=======================================//
653 // fill fitters with cluster information //
654 //=======================================//
655 AliDebug(3,"Filing Cluster Information");
263d466a 656 for (Int_t irow=158;irow>-1;--irow) {
657 AliTPCclusterMI *c=track->GetClusterPointer(irow);
0728a4f6 658 AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
659 cl=dummyCl;
660 //
263d466a 661 vecSec[irow]=-1;
662 if (!c) continue;
0728a4f6 663 //
664 Int_t roc = static_cast<Int_t>(c->GetDetector());
665 if ( roc!=innerSector && roc!=outerSector ) continue;
666 vecSec[irow]=roc;
263d466a 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();
263d466a 673 //
0728a4f6 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
263d466a 680 if ( roc == innerSector ) {
0728a4f6 681 x4[0]=1; //offset inner - outer sector
682 x4[1]=x; //slope parameter inner sector
263d466a 683 } else {
0728a4f6 684 x4[2]=x; //slope parameter outer sector
263d466a 685 }
0728a4f6 686 x4[3]=x*x; //common parabolic shape
263d466a 687 //
0728a4f6 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 }
263d466a 702 fy4.AddPoint(x4,y);
703 fz4.AddPoint(x4,z);
704 }
0728a4f6 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
263d466a 740 if ( innerSector>0 && outerSector>0 ){
0728a4f6 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 }
263d466a 751 }
0728a4f6 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 //==============================//
263d466a 761 AliDebug(4,"Calculate tracklet positions");
762 for (Int_t irow=158;irow>-1;--irow) {
0728a4f6 763 if ( vecSec[irow]==-1 ) continue; //no cluster info
263d466a 764 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
765 tr=&dummy;
766 Double_t x=vecX[irow];
0728a4f6 767 Double_t xref=x-133.4;
768 //
263d466a 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;
0728a4f6 774 Double_t zslopeOuter=0;
775 //positions of hyperplane fits
263d466a 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];
263d466a 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];
263d466a 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;
0728a4f6 797 //positions of kalman fits
263d466a 798 Double_t gxyz[3],xyz[3];
799 AliExternalTrackParam *param = 0x0;
0728a4f6 800 //
263d466a 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 }
263d466a 810 }
0728a4f6 811 //=====================================================================//
812 // write results from the different tracklet fits with debug streamers //
813 //=====================================================================//
263d466a 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 <<
0728a4f6 879 "sec.=" << &vecSec <<
263d466a 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 }
95a0e09b 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 fDeltaYres.AddAt(profy,id);
906 }
907 if (!profz){
908 profz=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),115,80,250);
909 fDeltaZres.AddAt(profz,id);
910 }
911 for (Int_t irow=158;irow>-1;--irow) {
912 if (vecSec[irow]==-1)continue; //no cluster info
913 Double_t x = vecX[irow];
914 Double_t ycl = vecClY[irow];
915 Double_t yfit = vecY1[irow];
916 Double_t zcl = vecClZ[irow];
917 Double_t zfit = vecZ1[irow];
918 profy->Fill(x,yfit-ycl);
919 profz->Fill(x,zfit-zcl);
920 }
ecb7e705 921 //===============================//
922 // Fill Fit Parameter Histograms //
923 //===============================//
924 TH1F *pol2InnerY =(TH1F*)fPol2Par2InY.UncheckedAt(id);
925 TH1F *diff1InnerY=(TH1F*)fDiffPar1InY.UncheckedAt(id);
926 TH1F *pol2OuterY =(TH1F*)fPol2Par2OutY.UncheckedAt(id);
927 TH1F *diff1OuterY=(TH1F*)fDiffPar1OutY.UncheckedAt(id);
928 TH1F *pol2InnerZ =(TH1F*)fPol2Par2InZ.UncheckedAt(id);
929 TH1F *diff1InnerZ=(TH1F*)fDiffPar1InZ.UncheckedAt(id);
930 TH1F *pol2OuterZ =(TH1F*)fPol2Par2OutZ.UncheckedAt(id);
931 TH1F *diff1OuterZ=(TH1F*)fDiffPar1OutZ.UncheckedAt(id);
932 //create histograms if the do not already exist
933 if (!pol2InnerY){
934 pol2InnerY =new TH1F(Form("pol2par2inY%03d",id),
935 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (inner sector)",id),
936 100,-.005,.005);
937 pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
938 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
939 100,0.01,.01);
940 diff1InnerY=new TH1F(Form("diff1inY%03d",id),
941 Form("diff of 1st derivative from pol1 and pol2 in Y fit for Laser beam %03d (inner sector)",id),
942 100,-.5,.5);
943 diff1OuterY=new TH1F(Form("diff1outY%03d",id),
944 Form("diff of 1st derivative from pol1 and pol2 in Yfit for Laser beam %03d (outer sector)",id),
945 100,-1,1);
946 pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
947 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
948 100,-.002,.002);
949 pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
950 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
951 100,-.005,.005);
952 diff1InnerZ=new TH1F(Form("diff1inZ%03d",id),
953 Form("diff of 1st derivative from pol1 and pol2 in Z fit for Laser beam %03d (inner sector)",id),
954 100,-.02,.02);
955 diff1OuterZ=new TH1F(Form("diff1outZ%03d",id),
956 Form("diff of 1st derivative from pol1 and pol2 in Zfit for Laser beam %03d (outer sector)",id),
957 100,-.03,.03);
958 //add
959 fPol2Par2InY.AddAt(pol2InnerY,id);
960 fDiffPar1InY.AddAt(diff1InnerY,id);
961 fPol2Par2OutY.AddAt(pol2OuterY,id);
962 fDiffPar1OutY.AddAt(diff1OuterY,id);
963 fPol2Par2InZ.AddAt(pol2InnerZ,id);
964 fDiffPar1InZ.AddAt(diff1InnerZ,id);
965 fPol2Par2OutZ.AddAt(pol2OuterZ,id);
966 fDiffPar1OutZ.AddAt(diff1OuterZ,id);
967 }
968 //fill histograms
969 pol2InnerY ->Fill(vecy2resInner[2]);
970 pol2OuterY ->Fill(vecy2resOuter[2]);
971 diff1InnerY->Fill(vecy2resInner[1]-vecy1resInner[1]);
972 diff1OuterY->Fill(vecy2resOuter[1]-vecy1resOuter[1]);
973 pol2InnerZ ->Fill(vecz2resInner[2]);
974 pol2OuterZ ->Fill(vecz2resOuter[2]);
975 diff1InnerZ->Fill(vecz2resInner[1]-vecz1resInner[1]);
976 diff1OuterZ->Fill(vecz2resOuter[1]-vecz1resOuter[1]);
95a0e09b 977
263d466a 978
ecb7e705 979
980 }
263d466a 981}
982
e9f38a4b 983
984void AliTPCcalibLaser::RefitLaser(Int_t id){
985 //
986 // Refit the track store residuals
987 //
988
989 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
990 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
991 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
263d466a 992
993 //linear fit model in y and z per sector
e9f38a4b 994 static TLinearFitter fy1(2,"hyp1");
995 static TLinearFitter fz1(2,"hyp1");
263d466a 996 //quadratic fit model in y and z per sector
997 static TLinearFitter fy2(3,"hyp2");
998 static TLinearFitter fz2(3,"hyp2");
e9f38a4b 999 static TVectorD vecy2,vecz2,vecy1,vecz1;
1000
1001 const Int_t kMinClusters=20;
263d466a 1002 Int_t nclusters[72];
e9f38a4b 1003 //
1004 for (Int_t i=0;i<72;++i) nclusters[i]=0;
1005
263d466a 1006 for (Int_t i=0;i<160;++i) {
e9f38a4b 1007 AliTPCclusterMI *c=track->GetClusterPointer(i);
1008 if (c) nclusters[c->GetDetector()]++;
1009 }
263d466a 1010
e9f38a4b 1011 for (Int_t isec=0; isec<72;isec++){
1012 if (nclusters[isec]<kMinClusters) continue;
1013 fy2.ClearPoints();
1014 fz2.ClearPoints();
1015 fy1.ClearPoints();
1016 fz1.ClearPoints();
1017 //
263d466a 1018 for (Int_t irow=0;irow<160;++irow) {
e9f38a4b 1019 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1020 //if (c && RejectCluster(c)) continue;
263d466a 1021 Double_t xd = c->GetX()-133.4; // reference x is beteen iroc and oroc
e9f38a4b 1022 if (c&&c->GetDetector()==isec) {
e9f38a4b 1023 Double_t x[2]={xd,xd*xd};
1024 fy2.AddPoint(x,c->GetY());
1025 fz2.AddPoint(x,c->GetZ());
1026 //
1027 fy1.AddPoint(x,c->GetY());
1028 fz1.AddPoint(x,c->GetZ());
1029 }
1030 }
1031 fy2.Eval();
1032 fz2.Eval();
1033 fy1.Eval();
1034 fz1.Eval();
1035 fy1.GetParameters(vecy1);
1036 fy2.GetParameters(vecy2);
1037 fz1.GetParameters(vecz1);
1038 fz2.GetParameters(vecz2);
263d466a 1039
e9f38a4b 1040 if (fStreamLevel>0){
1041 TTreeSRedirector *cstream = GetDebugStreamer();
1042 if (cstream){
1043 Float_t dedx = track->GetdEdx();
1044 (*cstream)<<"Tracklet"<<
1045 "LTr.="<<ltrp<<
1046 "Tr.="<<param<<
1047 "sec="<<isec<<
1048 "ncl="<<nclusters[isec]<<
1049 "dedx="<<dedx<<
1050 "dedx="<<dedx<<
1051 "vy1.="<<&vecy1<<
1052 "vy2.="<<&vecy2<<
1053 "vz1.="<<&vecz1<<
1054 "vz2.="<<&vecz2<<
1055 "\n";
1056 }
1057 }
1058 }
1059 //
1060 //
1061 //
263d466a 1062 // for (Int_t irow=0;irow<160;++irow) {
e9f38a4b 1063 // AliTPCclusterMI *c=track->GetClusterPointer(irow);
1064 // if (c && RejectCluster(c)) continue;
1065 // if (c&&c->GetDetector()==isec) {
1066 // Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
1067 // fy2.AddPoint(&x,c->GetY());
1068 // fz2.AddPoint(&x,c->GetZ());
1069 // //
1070 // fy1.AddPoint(&x,c->GetY());
1071 // fz1.AddPoint(&x,c->GetZ());
1072 // }
263d466a 1073 // }
1074
e9f38a4b 1075}
1076
1077
1fd56785 1078void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield,Int_t minEntries){
2b35e8e6 1079 //
1fd56785 1080 // Dump information about laser beams
1081 // isOK variable indicates usability of the beam
1082 // Beam is not usable if:
1083 // a. No entries in range (krmsCut0)
1084 // b. Big sperad (krmscut1)
1085 // c. RMSto fit sigma bigger then (kmultiCut)
1086 // d. Too big angular spread
1087 //
1088
1089 const Float_t krmsCut0=0.001;
1090 const Float_t krmsCut1=0.16;
1091 const Float_t kmultiCut=2;
1092 const Float_t kcutP0=0.002;
2b35e8e6 1093 //
1094 AliTPCcalibLaser *laser = this;
1095 TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
1096 TF1 fg("fg","gaus");
1097 //
1098 //
1099 for (Int_t id=0; id<336; id++){
1fd56785 1100 Bool_t isOK=kTRUE;
2b35e8e6 1101 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
1102 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
1103 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
1fd56785 1104 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
2b35e8e6 1105 if (!hisphi) continue;;
1106 Double_t entries = hisphi->GetEntries();
1fd56785 1107 if (entries<minEntries) continue;
2b35e8e6 1108 //
1109 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1110 if (!ltrp) {
1111 AliTPCLaserTrack::LoadTracks();
1112 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1113 }
1114 Float_t meanphi = hisphi->GetMean();
1115 Float_t rmsphi = hisphi->GetRMS();
1fd56785 1116 //
2b35e8e6 1117 Float_t meanphiP = hisphiP->GetMean();
1118 Float_t rmsphiP = hisphiP->GetRMS();
1119 Float_t meanZ = hisZ->GetMean();
1120 Float_t rmsZ = hisZ->GetRMS();
1121 hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
263d466a 1122 Double_t gphi1 = fg.GetParameter(1);
1123 Double_t gphi2 = fg.GetParameter(2);
2b35e8e6 1124 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
263d466a 1125 Double_t gphiP1 = fg.GetParameter(1);
1126 Double_t gphiP2 = fg.GetParameter(2);
2b35e8e6 1127 hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
263d466a 1128 Double_t gz1 = fg.GetParameter(1);
1129 Double_t gz2 = fg.GetParameter(2);
1fd56785 1130 //
1131 Float_t meanS=hisS->GetMean();
2b35e8e6 1132 //
1133 Double_t lxyz[3];
1134 Double_t lpxyz[3];
1135 ltrp->GetXYZ(lxyz);
1136 ltrp->GetPxPyPz(lpxyz);
1fd56785 1137
1138 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
1139 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
1140 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
1141 if (gphiP2>kcutP0) isOK=kFALSE;
2b35e8e6 1142 //
1143 (*pcstream)<<"Mean"<<
1fd56785 1144 "isOK="<<isOK<<
2b35e8e6 1145 "entries="<<entries<< // number of entries
1146 "bz="<<bfield<< // bfield
1147 "LTr.="<<ltrp<< // refernece track
1fd56785 1148 //
2b35e8e6 1149 "lx0="<<lxyz[0]<< // reference x
1150 "lx1="<<lxyz[1]<< // reference y
263d466a 1151 "lx2="<<lxyz[2]<< // refernece z
2b35e8e6 1152 "lpx0="<<lpxyz[0]<< // reference x
1153 "lpx1="<<lpxyz[1]<< // reference y
263d466a 1154 "lpx2="<<lpxyz[2]<< // refernece z
1fd56785 1155 //
1156 "msig="<<meanS<<
1157 //
2b35e8e6 1158 "mphi="<<meanphi<< //
1159 "rmsphi="<<rmsphi<< //
1160 "gphi1="<<gphi1<<
1161 "gphi2="<<gphi2<<
1fd56785 1162 //
2b35e8e6 1163 "mphiP="<<meanphiP<< //
1164 "rmsphiP="<<rmsphiP<< //
1165 "gphiP1="<<gphiP1<<
1166 "gphiP2="<<gphiP2<<
1fd56785 1167 //
2b35e8e6 1168 "meanZ="<<meanZ<<
1169 "rmsZ="<<rmsZ<<
1170 "gz1="<<gz1<<
1171 "gz2="<<gz2<<
1172
1173 "\n";
1174 }
1175 delete pcstream;
1176}
1177
1fd56785 1178
1179
1180void AliTPCcalibLaser::DumpScanInfo(TTree * chain){
1181 //
1182 //
1183 //
1184 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
1185 TFile * f = pcstream->GetFile();
1186 f->mkdir("dirphi");
1187 f->mkdir("dirphiP");
1188 f->mkdir("dirZ");
1189 TF1 fp("p1","pol1");
1190 //
1191 //
1192 char cut[1000];
1193 char grname[1000];
1194 char grnamefull[1000];
1195
1196 Double_t mphi[100];
1197 Double_t mphiP[100];
1198 Double_t smphi[100];
1199 Double_t smphiP[100];
1200 Double_t mZ[100];
1201 Double_t smZ[100];
1202 Double_t bz[100];
1203 Double_t sbz[100];
1204 // fit parameters
1205 Double_t pphi[3];
1206 Double_t pphiP[3];
1207 Double_t pmZ[3];
1208 //
1209 for (Int_t id=0; id<336; id++){
1210 // id =205;
1211 sprintf(cut,"isOK&&fId==%d",id);
1212 Int_t entries = chain->Draw("bz",cut,"goff");
1213 if (entries<3) continue;
1214 AliTPCLaserTrack *ltrp = 0;;
1215 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
1216 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1217 Double_t lxyz[3];
1218 Double_t lpxyz[3];
1219 ltrp->GetXYZ(lxyz);
1220 ltrp->GetPxPyPz(lpxyz);
1221
1222 chain->Draw("bz",cut,"goff");
1223 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
1224 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
1225 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
1226 //
1227 chain->Draw("gphi1",cut,"goff");
1228 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
1229 chain->Draw("0.05*abs(mphi)+gphi2",cut,"goff");
1230 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
1231 //
1232 chain->Draw("gphiP1",cut,"goff");
1233 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
1234 chain->Draw("0.05*abs(mphiP)+gphiP2",cut,"goff");
1235 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
1236 //
1237 chain->Draw("gz1",cut,"goff");
1238 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
1239 chain->Draw("0.01*abs(meanZ)+gz2",cut,"goff");
1240 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
1241 //
1242 //
1243 sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
1244 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
1245 // store data
1246 // phi
1247 f->cd("dirphi");
1248 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
1249 grphi->Draw("a*");
1250 grphi->Fit(&fp);
1251 pphi[0] = fp.GetParameter(0); // offset
1252 pphi[1] = fp.GetParameter(1); // slope
1253 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1254 sprintf(grname,"phi_id%d",id);
1255 grphi->SetName(grname); grphi->SetTitle(grnamefull);
1256 grphi->GetXaxis()->SetTitle("b_{z} (T)");
1257 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
ae69f16f 1258 grphi->SetMaximum(1.2);
1259 grphi->SetMinimum(-1.2);
1fd56785 1260 grphi->Draw("a*");
1261
1262 grphi->Write();
1263 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
1264 // phiP
1265 f->cd("dirphiP)");
1266 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
1267 grphiP->Draw("a*");
1268 grphiP->Fit(&fp);
1269 pphiP[0] = fp.GetParameter(0); // offset
1270 pphiP[1] = fp.GetParameter(1); // slope
1271 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1272 sprintf(grname,"phiP_id%d",id);
1273 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
1274 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
1275 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
ae69f16f 1276 grphiP->SetMaximum(pphiP[0]+0.005);
1277 grphiP->SetMinimum(pphiP[0]-0.005);
1fd56785 1278
1279 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
1280 grphiP->Write();
1281 //
1282 //Z
1283 f->cd("dirZ");
1284 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
1285 grmZ->Draw("a*");
1286 grmZ->Fit(&fp);
1287 pmZ[0] = fp.GetParameter(0); // offset
1288 pmZ[1] = fp.GetParameter(1); // slope
1289 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1290 sprintf(grname,"mZ_id%d",id);
1291 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
1292 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
1293 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
1294
1295 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
1296 grmZ->Write();
1297
1298
1299 for (Int_t ientry=0; ientry<entries; ientry++){
1300 (*pcstream)<<"Mean"<<
1301 "id="<<id<<
1302 "LTr.="<<ltrp<<
1303 "entries="<<entries<<
1304 "bz="<<bz[ientry]<<
1305 "lx0="<<lxyz[0]<< // reference x
1306 "lx1="<<lxyz[1]<< // reference y
1307 "lx2="<<lxyz[2]<< // refernece z
1308 "lpx0="<<lpxyz[0]<< // reference x
1309 "lpx1="<<lpxyz[1]<< // reference y
1310 "lpx2="<<lpxyz[2]<< // refernece z
1311 //values
1312 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
1313 "pphi0="<<pphi[0]<< // offset
1314 "pphi1="<<pphi[1]<< // mean
1315 "pphi2="<<pphi[2]<< // norm chi2
1316 //
1317 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
1318 "pphiP0="<<pphiP[0]<< // offset
1319 "pphiP1="<<pphiP[1]<< // mean
1320 "pphiP2="<<pphiP[2]<< // norm chi2
1321 //
1322 "gz1="<<mZ[ientry]<<
1323 "pmZ0="<<pmZ[0]<< // offset
1324 "pmZ1="<<pmZ[1]<< // mean
1325 "pmZ2="<<pmZ[2]<< // norm chi2
1326 "\n";
1327 }
1328 }
1329
1330 delete pcstream;
1331
1332}
1333
1334
c6914c83 1335void AliTPCcalibLaser::Analyze(){
1336 //
1337 //
1338 //
1339}
1340
1341
c03e3250 1342Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
c6914c83 1343
c03e3250 1344 TIterator* iter = li->MakeIterator();
1345 AliTPCcalibLaser* cal = 0;
c6914c83 1346
c03e3250 1347 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
1348 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
1349 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1350 return -1;
1351 }
1352 //
1353 // fHistNTracks->Add(cal->fHistNTracks);
1354// fClusters->Add(cal->fClusters);
1355// fModules->Add(cal->fModules);
1356// fHistPt->Add(cal->fHistPt);
1357// fPtResolution->Add(cal->fPtResolution);
1358// fDeDx->Add(cal->fDeDx);
263d466a 1359
c03e3250 1360
1361 TH1F *h=0x0;
1362 TH1F *hm=0x0;
95a0e09b 1363 TProfile *hp=0x0;
1364 TProfile *hpm=0x0;
c03e3250 1365
1366 for (Int_t id=0; id<336; id++){
1367 // merge fDeltaZ histograms
263d466a 1368 hm = (TH1F*)cal->fDeltaZ.At(id);
c03e3250 1369 h = (TH1F*)fDeltaZ.At(id);
1370 if (!h) {
1371 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
1372 fDeltaZ.AddAt(h,id);
1373 }
1374 if (hm) h->Add(hm);
1375 // merge fDeltaPhi histograms
263d466a 1376 hm = (TH1F*)cal->fDeltaPhi.At(id);
c03e3250 1377 h = (TH1F*)fDeltaPhi.At(id);
1378 if (!h) {
1379 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
1380 fDeltaPhi.AddAt(h,id);
1381 }
1382 if (hm) h->Add(hm);
1383 // merge fDeltaPhiP histograms
263d466a 1384 hm = (TH1F*)cal->fDeltaPhiP.At(id);
c03e3250 1385 h = (TH1F*)fDeltaPhiP.At(id);
1386 if (!h) {
1387 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
1388 fDeltaPhiP.AddAt(h,id);
1389 }
1390 if (hm) h->Add(hm);
1391 // merge fSignals histograms
263d466a 1392 hm = (TH1F*)cal->fSignals.At(id);
c03e3250 1393 h = (TH1F*)fSignals.At(id);
1394 if (!h) {
1395 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
1396 fSignals.AddAt(h,id);
1397 }
263d466a 1398 if (hm) h->Add(hm);
95a0e09b 1399 //
1400 //
1401 // merge ProfileY histograms
1402 hpm = (TProfile*)cal->fDeltaYres.At(id);
1403 hp = (TProfile*)fDeltaYres.At(id);
1404 if (!hp) {
1405 hp=new TProfile(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d",id),115,80,250);;
1406 fDeltaYres.AddAt(hp,id);
1407 }
1408 if (hpm) hp->Add(hpm);
1409 //
1410 hpm = (TProfile*)cal->fDeltaZres.At(id);
1411 hp = (TProfile*)fDeltaZres.At(id);
1412 if (!hp) {
1413 hp=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),115,80,250);;
1414 fDeltaZres.AddAt(hp,id);
1415 }
1416 if (hpm) hp->Add(hpm);
1417 //
1418 //
ecb7e705 1419 //merge fit param histograms
1420 //local hists
1421 TH1F *pol2InnerY =(TH1F*)fPol2Par2InY.UncheckedAt(id);
1422 TH1F *diff1InnerY=(TH1F*)fDiffPar1InY.UncheckedAt(id);
1423 TH1F *pol2OuterY =(TH1F*)fPol2Par2OutY.UncheckedAt(id);
1424 TH1F *diff1OuterY=(TH1F*)fDiffPar1OutY.UncheckedAt(id);
1425 TH1F *pol2InnerZ =(TH1F*)fPol2Par2InZ.UncheckedAt(id);
1426 TH1F *diff1InnerZ=(TH1F*)fDiffPar1InZ.UncheckedAt(id);
1427 TH1F *pol2OuterZ =(TH1F*)fPol2Par2OutZ.UncheckedAt(id);
1428 TH1F *diff1OuterZ=(TH1F*)fDiffPar1OutZ.UncheckedAt(id);
1429 //hists to merge
1430 TH1F *pol2InnerYm =(TH1F*)cal->fPol2Par2InY.UncheckedAt(id);
1431 TH1F *diff1InnerYm=(TH1F*)cal->fDiffPar1InY.UncheckedAt(id);
1432 TH1F *pol2OuterYm =(TH1F*)cal->fPol2Par2OutY.UncheckedAt(id);
1433 TH1F *diff1OuterYm=(TH1F*)cal->fDiffPar1OutY.UncheckedAt(id);
1434 TH1F *pol2InnerZm =(TH1F*)cal->fPol2Par2InZ.UncheckedAt(id);
1435 TH1F *diff1InnerZm=(TH1F*)cal->fDiffPar1InZ.UncheckedAt(id);
1436 TH1F *pol2OuterZm =(TH1F*)cal->fPol2Par2OutZ.UncheckedAt(id);
1437 TH1F *diff1OuterZm=(TH1F*)cal->fDiffPar1OutZ.UncheckedAt(id);
1438 //create histos if they do not exist
1439 if (!pol2InnerY){
1440 pol2InnerY =new TH1F(Form("pol2par2inY%03d",id),
1441 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (inner sector)",id),
1442 100,-.005,.005);
1443 pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
1444 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
1445 100,0.01,.01);
1446 diff1InnerY=new TH1F(Form("diff1inY%03d",id),
1447 Form("diff of 1st derivative from pol1 and pol2 in Y fit for Laser beam %03d (inner sector)",id),
1448 100,-.5,.5);
1449 diff1OuterY=new TH1F(Form("diff1outY%03d",id),
1450 Form("diff of 1st derivative from pol1 and pol2 in Yfit for Laser beam %03d (outer sector)",id),
1451 100,-1,1);
1452 pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
1453 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
1454 100,-.002,.002);
1455 pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
1456 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
1457 100,-.005,.005);
1458 diff1InnerZ=new TH1F(Form("diff1inZ%03d",id),
1459 Form("diff of 1st derivative from pol1 and pol2 in Z fit for Laser beam %03d (inner sector)",id),
1460 100,-.02,.02);
1461 diff1OuterZ=new TH1F(Form("diff1outZ%03d",id),
1462 Form("diff of 1st derivative from pol1 and pol2 in Zfit for Laser beam %03d (outer sector)",id),
1463 100,-.03,.03);
1464 }
1465 pol2InnerYm =(TH1F*)cal->fPol2Par2InY.UncheckedAt(id);
1466 diff1InnerYm=(TH1F*)cal->fDiffPar1InY.UncheckedAt(id);
1467 pol2OuterYm =(TH1F*)cal->fPol2Par2OutY.UncheckedAt(id);
1468 diff1OuterYm=(TH1F*)cal->fDiffPar1OutY.UncheckedAt(id);
1469 pol2InnerZm =(TH1F*)cal->fPol2Par2InZ.UncheckedAt(id);
1470 diff1InnerZm=(TH1F*)cal->fDiffPar1InZ.UncheckedAt(id);
1471 pol2OuterZm =(TH1F*)cal->fPol2Par2OutZ.UncheckedAt(id);
1472 diff1OuterZm=(TH1F*)cal->fDiffPar1OutZ.UncheckedAt(id);
c03e3250 1473 }
1474 }
1475 return 0;
1476}
1477
1478
1479
263d466a 1480
c03e3250 1481/*
1482 gSystem->Load("libSTAT.so")
1483 TStatToolkit toolkit;
1484 Double_t chi2;
1485 TVectorD fitParam;
1486 TMatrixD covMatrix;
1487 Int_t npoints;
ae69f16f 1488 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
c03e3250 1489
1490
1491TString fstring="";
1fd56785 1492//
1493fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
1494fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
1495fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
1496fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
1497//
1498fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
1499fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
1500fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
1501fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
1502//
1503fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
1504fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
1505fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
1506fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
c5595838 1507
c5595838 1508
c5595838 1509
c5595838 1510
1fd56785 1511 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c03e3250 1512
1fd56785 1513 treeT->SetAlias("fit",strq0->Data());
1514
c03e3250 1515
1fd56785 1516 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c5595838 1517
1fd56785 1518 treeT->SetAlias("fitP",strqP->Data());
c5595838 1519
1520
ae69f16f 1521 TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
1522 treeT->SetAlias("fitD",strqDrift->Data());
c5595838 1523
c03e3250 1524
ae69f16f 1525treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
1526{
1527for (Int_t i=0; i<6;i++){
1528treeT->SetLineColor(i+2);
1529treeT->SetMarkerSize(1);
1530treeT->SetMarkerStyle(22+i);
1531treeT->SetMarkerColor(i+2);
c03e3250 1532
ae69f16f 1533treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
1534}
1535}
c03e3250 1536 */
68ff0583 1537
1538
1539
1540/*
1541 TTree * tree = (TTree*)f.Get("FitModels");
1542
1543 TEventList listLFit0("listLFit0","listLFit0");
1544 TEventList listLFit1("listLFit1","listLFit1");
68ff0583 1545 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
1546 tree->SetEventList(&listLFit0);
1547
1548
95a0e09b 1549
1550
1551 gSystem->Load("libSTAT.so")
1552 TStatToolkit toolkit;
1553 Double_t chi2;
1554 TVectorD fitParam;
1555 TMatrixD covMatrix;
1556 Int_t npoints;
1557
1558 chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
1559 chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
1560
1561
1562 TString fstring="";
1563 fstring+="cos(dp)++";
1564 fstring+="sin(dp)++";
1565 fstring+="cos(dt)++";
1566 fstring+="sin(dt)++";
1567
1568 TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
1569
1570
68ff0583 1571*/