]> 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 pad type correction (Marian)
[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;
5732d662 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);
5732d662 905 profy->SetDirectory(0);
95a0e09b 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);
5732d662 910 profz->SetDirectory(0);
95a0e09b 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];
5732d662 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);
95a0e09b 926 }
ecb7e705 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),
674ee45c 942 500,-.005,.005);
5732d662 943 pol2InnerY->SetDirectory(0);
ecb7e705 944 pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
945 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
674ee45c 946 500,0.01,.01);
5732d662 947 pol2OuterY->SetDirectory(0);
948
ecb7e705 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),
674ee45c 951 500,-.5,.5);
5732d662 952 diff1InnerY->SetDirectory(0);
953
ecb7e705 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),
674ee45c 956 500,-1,1);
5732d662 957 diff1OuterY->SetDirectory(0);
958
ecb7e705 959 pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
960 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
674ee45c 961 500,-.002,.002);
5732d662 962 pol2InnerZ->SetDirectory(0);
963
ecb7e705 964 pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
965 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
674ee45c 966 500,-.005,.005);
5732d662 967 pol2OuterZ->SetDirectory(0);
ecb7e705 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),
674ee45c 970 500,-.02,.02);
5732d662 971 diff1InnerZ->SetDirectory(0);
ecb7e705 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),
674ee45c 974 500,-.03,.03);
5732d662 975 diff1OuterZ->SetDirectory(0);
ecb7e705 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]);
95a0e09b 995
263d466a 996
ecb7e705 997
998 }
263d466a 999}
1000
e9f38a4b 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);
263d466a 1010
1011 //linear fit model in y and z per sector
e9f38a4b 1012 static TLinearFitter fy1(2,"hyp1");
1013 static TLinearFitter fz1(2,"hyp1");
263d466a 1014 //quadratic fit model in y and z per sector
1015 static TLinearFitter fy2(3,"hyp2");
1016 static TLinearFitter fz2(3,"hyp2");
e9f38a4b 1017 static TVectorD vecy2,vecz2,vecy1,vecz1;
1018
1019 const Int_t kMinClusters=20;
263d466a 1020 Int_t nclusters[72];
e9f38a4b 1021 //
1022 for (Int_t i=0;i<72;++i) nclusters[i]=0;
1023
263d466a 1024 for (Int_t i=0;i<160;++i) {
e9f38a4b 1025 AliTPCclusterMI *c=track->GetClusterPointer(i);
1026 if (c) nclusters[c->GetDetector()]++;
1027 }
263d466a 1028
e9f38a4b 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 //
263d466a 1036 for (Int_t irow=0;irow<160;++irow) {
e9f38a4b 1037 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1038 //if (c && RejectCluster(c)) continue;
263d466a 1039 Double_t xd = c->GetX()-133.4; // reference x is beteen iroc and oroc
e9f38a4b 1040 if (c&&c->GetDetector()==isec) {
e9f38a4b 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);
263d466a 1057
e9f38a4b 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 //
263d466a 1080 // for (Int_t irow=0;irow<160;++irow) {
e9f38a4b 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 // }
263d466a 1091 // }
1092
e9f38a4b 1093}
1094
1095
1fd56785 1096void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield,Int_t minEntries){
2b35e8e6 1097 //
1fd56785 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;
2b35e8e6 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++){
1fd56785 1118 Bool_t isOK=kTRUE;
2b35e8e6 1119 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
1120 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
1121 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
1fd56785 1122 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
2b35e8e6 1123 if (!hisphi) continue;;
1124 Double_t entries = hisphi->GetEntries();
1fd56785 1125 if (entries<minEntries) continue;
2b35e8e6 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();
1fd56785 1134 //
2b35e8e6 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());
263d466a 1140 Double_t gphi1 = fg.GetParameter(1);
1141 Double_t gphi2 = fg.GetParameter(2);
2b35e8e6 1142 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
263d466a 1143 Double_t gphiP1 = fg.GetParameter(1);
1144 Double_t gphiP2 = fg.GetParameter(2);
2b35e8e6 1145 hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
263d466a 1146 Double_t gz1 = fg.GetParameter(1);
1147 Double_t gz2 = fg.GetParameter(2);
1fd56785 1148 //
1149 Float_t meanS=hisS->GetMean();
2b35e8e6 1150 //
1151 Double_t lxyz[3];
1152 Double_t lpxyz[3];
1153 ltrp->GetXYZ(lxyz);
1154 ltrp->GetPxPyPz(lpxyz);
1fd56785 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;
2b35e8e6 1160 //
1161 (*pcstream)<<"Mean"<<
1fd56785 1162 "isOK="<<isOK<<
2b35e8e6 1163 "entries="<<entries<< // number of entries
1164 "bz="<<bfield<< // bfield
1165 "LTr.="<<ltrp<< // refernece track
1fd56785 1166 //
2b35e8e6 1167 "lx0="<<lxyz[0]<< // reference x
1168 "lx1="<<lxyz[1]<< // reference y
263d466a 1169 "lx2="<<lxyz[2]<< // refernece z
2b35e8e6 1170 "lpx0="<<lpxyz[0]<< // reference x
1171 "lpx1="<<lpxyz[1]<< // reference y
263d466a 1172 "lpx2="<<lpxyz[2]<< // refernece z
1fd56785 1173 //
1174 "msig="<<meanS<<
1175 //
2b35e8e6 1176 "mphi="<<meanphi<< //
1177 "rmsphi="<<rmsphi<< //
1178 "gphi1="<<gphi1<<
1179 "gphi2="<<gphi2<<
1fd56785 1180 //
2b35e8e6 1181 "mphiP="<<meanphiP<< //
1182 "rmsphiP="<<rmsphiP<< //
1183 "gphiP1="<<gphiP1<<
1184 "gphiP2="<<gphiP2<<
1fd56785 1185 //
2b35e8e6 1186 "meanZ="<<meanZ<<
1187 "rmsZ="<<rmsZ<<
1188 "gz1="<<gz1<<
1189 "gz2="<<gz2<<
1190
1191 "\n";
1192 }
1193 delete pcstream;
1194}
1195
1fd56785 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)");
ae69f16f 1276 grphi->SetMaximum(1.2);
1277 grphi->SetMinimum(-1.2);
1fd56785 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)");
ae69f16f 1294 grphiP->SetMaximum(pphiP[0]+0.005);
1295 grphiP->SetMinimum(pphiP[0]-0.005);
1fd56785 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
c6914c83 1353void AliTPCcalibLaser::Analyze(){
1354 //
1355 //
1356 //
1357}
1358
1359
c03e3250 1360Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
c6914c83 1361
c03e3250 1362 TIterator* iter = li->MakeIterator();
1363 AliTPCcalibLaser* cal = 0;
c6914c83 1364
c03e3250 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);
263d466a 1377
c03e3250 1378
1379 TH1F *h=0x0;
1380 TH1F *hm=0x0;
95a0e09b 1381 TProfile *hp=0x0;
1382 TProfile *hpm=0x0;
c03e3250 1383
1384 for (Int_t id=0; id<336; id++){
1385 // merge fDeltaZ histograms
263d466a 1386 hm = (TH1F*)cal->fDeltaZ.At(id);
c03e3250 1387 h = (TH1F*)fDeltaZ.At(id);
1388 if (!h) {
1389 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
5732d662 1390 h->SetDirectory(0);
c03e3250 1391 fDeltaZ.AddAt(h,id);
1392 }
1393 if (hm) h->Add(hm);
1394 // merge fDeltaPhi histograms
263d466a 1395 hm = (TH1F*)cal->fDeltaPhi.At(id);
c03e3250 1396 h = (TH1F*)fDeltaPhi.At(id);
1397 if (!h) {
1398 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
5732d662 1399 h->SetDirectory(0);
c03e3250 1400 fDeltaPhi.AddAt(h,id);
1401 }
1402 if (hm) h->Add(hm);
1403 // merge fDeltaPhiP histograms
263d466a 1404 hm = (TH1F*)cal->fDeltaPhiP.At(id);
c03e3250 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);
5732d662 1408 h->SetDirectory(0);
c03e3250 1409 fDeltaPhiP.AddAt(h,id);
1410 }
1411 if (hm) h->Add(hm);
1412 // merge fSignals histograms
263d466a 1413 hm = (TH1F*)cal->fSignals.At(id);
c03e3250 1414 h = (TH1F*)fSignals.At(id);
1415 if (!h) {
1416 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
5732d662 1417 h->SetDirectory(0);
c03e3250 1418 fSignals.AddAt(h,id);
1419 }
263d466a 1420 if (hm) h->Add(hm);
95a0e09b 1421 //
1422 //
1423 // merge ProfileY histograms
1424 hpm = (TProfile*)cal->fDeltaYres.At(id);
1425 hp = (TProfile*)fDeltaYres.At(id);
1426 if (!hp) {
5732d662 1427 hp=new TProfile(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d",id),115,80,250);
1428 hp->SetDirectory(0);
95a0e09b 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) {
5732d662 1436 hp=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),115,80,250);
1437 hp->SetDirectory(0);
95a0e09b 1438 fDeltaZres.AddAt(hp,id);
1439 }
1440 if (hpm) hp->Add(hpm);
1441 //
1442 //
ecb7e705 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),
674ee45c 1466 500,-.005,.005);
5732d662 1467 pol2InnerY->SetDirectory(0);
1468
ecb7e705 1469 pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
1470 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
674ee45c 1471 500,0.01,.01);
5732d662 1472 pol2OuterY->SetDirectory(0);
ecb7e705 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),
674ee45c 1475 500,-.5,.5);
ecb7e705 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),
674ee45c 1478 500,-1,1);
5732d662 1479 diff1InnerY->SetDirectory(0);
1480
1481
ecb7e705 1482 pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
1483 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
674ee45c 1484 500,-.002,.002);
5732d662 1485 pol2InnerZ->SetDirectory(0);
1486
ecb7e705 1487 pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
1488 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
674ee45c 1489 500,-.005,.005);
5732d662 1490 pol2OuterZ->SetDirectory(0);
ecb7e705 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),
674ee45c 1493 500,-.02,.02);
5732d662 1494 diff1InnerZ->SetDirectory(0);
ecb7e705 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),
674ee45c 1497 500,-.03,.03);
5732d662 1498 diff1OuterZ->SetDirectory(0);
ecb7e705 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);
c03e3250 1508 }
1509 }
1510 return 0;
1511}
1512
1513
1514
263d466a 1515
c03e3250 1516/*
1517 gSystem->Load("libSTAT.so")
1518 TStatToolkit toolkit;
1519 Double_t chi2;
1520 TVectorD fitParam;
1521 TMatrixD covMatrix;
1522 Int_t npoints;
ae69f16f 1523 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
c03e3250 1524
1525
1526TString fstring="";
1fd56785 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
c5595838 1542
c5595838 1543
c5595838 1544
c5595838 1545
1fd56785 1546 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c03e3250 1547
1fd56785 1548 treeT->SetAlias("fit",strq0->Data());
1549
c03e3250 1550
1fd56785 1551 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c5595838 1552
1fd56785 1553 treeT->SetAlias("fitP",strqP->Data());
c5595838 1554
1555
ae69f16f 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());
c5595838 1558
c03e3250 1559
ae69f16f 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);
c03e3250 1567
ae69f16f 1568treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
1569}
1570}
c03e3250 1571 */
68ff0583 1572
1573
1574
1575/*
1576 TTree * tree = (TTree*)f.Get("FitModels");
1577
1578 TEventList listLFit0("listLFit0","listLFit0");
1579 TEventList listLFit1("listLFit1","listLFit1");
68ff0583 1580 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
1581 tree->SetEventList(&listLFit0);
1582
1583
95a0e09b 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
68ff0583 1606*/