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