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