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