]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibLaser.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removed obsolete classes
[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");
e5c59218 39 laser->DumpMeanInfo(-0,0,10)
2b35e8e6 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;
e5c59218 49 TChain * chainTrL = tool.MakeChain("laser.txt","Track",0,10200);
50
1fd56785 51 //
52 //
53 // Analyze LASER scan
54 //
e5c59218 55
2b35e8e6 56 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
57 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
0728a4f6 58 AliXRDPROOFtoolkit tool;
2b35e8e6 59 TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
60 chain->Lookup();
1fd56785 61 AliTPCcalibLaser::DumpScanInfo(chain)
62 TFile fscan("laserScan.root")
63 TTree * treeT = (TTree*)fscan.Get("Mean")
95a0e09b 64 //
65 // Analyze laser
66 //
67 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
68 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
69 AliXRDPROOFtoolkit tool;
70 TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
71 chain->Lookup();
e5c59218 72 TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
73 chainFit->Lookup();
95a0e09b 74
c6914c83 75*/
76
77
78
79#include "TLinearFitter.h"
80#include "AliTPCcalibLaser.h"
81#include "AliExternalTrackParam.h"
e9f38a4b 82#include "AliESDEvent.h"
83#include "AliESDfriend.h"
c6914c83 84#include "AliESDtrack.h"
85#include "AliTPCTracklet.h"
86#include "TH1D.h"
ecb7e705 87#include "TH1F.h"
95a0e09b 88#include "TProfile.h"
c6914c83 89#include "TVectorD.h"
90#include "TTreeStream.h"
91#include "TFile.h"
92#include "TF1.h"
93#include "TGraphErrors.h"
94#include "AliTPCclusterMI.h"
95#include "AliTPCseed.h"
96#include "AliTracker.h"
263d466a 97#include "AliLog.h"
c6914c83 98#include "TClonesArray.h"
1fd56785 99#include "TPad.h"
e5c59218 100#include "TSystem.h"
101#include "TCut.h"
102#include "TChain.h"
103#include "TH2F.h"
104#include "TStatToolkit.h"
105#include "TROOT.h"
c6914c83 106
107
108#include "TTreeStream.h"
109#include <iostream>
110#include <sstream>
7b18d067 111#include "AliTPCLaserTrack.h"
112
c6914c83 113using namespace std;
114
115ClassImp(AliTPCcalibLaser)
116
117AliTPCcalibLaser::AliTPCcalibLaser():
e9f38a4b 118 AliTPCcalibBase(),
119 fESD(0),
120 fESDfriend(),
c18f4385 121 fTracksMirror(336),
122 fTracksEsd(336),
123 fTracksEsdParam(336),
124 fTracksTPC(336),
263d466a 125 fDeltaZ(336),
e5c59218 126 fDeltaP3(336),
127 fDeltaP4(336),
263d466a 128 fDeltaPhi(336),
129 fDeltaPhiP(336),
95a0e09b 130 fSignals(336),
131 fDeltaYres(336),
132 fDeltaZres(336),
ecb7e705 133 fPol2Par2InY(336),
134 fDiffPar1InY(336),
135 fPol2Par2OutY(336),
136 fDiffPar1OutY(336),
137 fPol2Par2InZ(336),
138 fDiffPar1InZ(336),
139 fPol2Par2OutZ(336),
140 fDiffPar1OutZ(336),
141 fFitAside(new TVectorD(3)),
263d466a 142 fFitCside(new TVectorD(3)),
143 fEdgeXcuts(5),
144 fEdgeYcuts(5),
145 fNClCuts(5),
146 fNcuts(0),
147 fRun(0),
148 fEvent(0)
c6914c83 149{
150 //
151 // Constructor
152 //
e9f38a4b 153 fTracksEsdParam.SetOwner(kTRUE);
c6914c83 154}
155
156AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
e9f38a4b 157 AliTPCcalibBase(),
158 fESD(0),
159 fESDfriend(0),
c18f4385 160 fTracksMirror(336),
161 fTracksEsd(336),
162 fTracksEsdParam(336),
163 fTracksTPC(336),
164 fDeltaZ(336), // array of histograms of delta z for each track
e5c59218 165 fDeltaP3(336), // array of histograms of delta z for each track
166 fDeltaP4(336), // array of histograms of P3 for each track
167 fDeltaPhi(336), // array of histograms of P4 for each track
c18f4385 168 fDeltaPhiP(336), // array of histograms of delta z for each track
2b35e8e6 169 fSignals(336), // array of dedx signals
95a0e09b 170 fDeltaYres(336),
171 fDeltaZres(336),
ecb7e705 172 fPol2Par2InY(336),
173 fDiffPar1InY(336),
174 fPol2Par2OutY(336),
175 fDiffPar1OutY(336),
176 fPol2Par2InZ(336),
177 fDiffPar1InZ(336),
178 fPol2Par2OutZ(336),
179 fDiffPar1OutZ(336),
c18f4385 180 fFitAside(new TVectorD(3)), // drift fit - A side
181 fFitCside(new TVectorD(3)), // drift fit - C- side
263d466a 182 fEdgeXcuts(5), // cuts in local x direction; used in the refit of the laser tracks
183 fEdgeYcuts(5), // cuts in local y direction; used in the refit of the laser tracks
184 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
185 fNcuts(0), // number of cuts
186 fRun(0),
187 fEvent(0)
c6914c83 188{
189 SetName(name);
190 SetTitle(title);
191 //
192 // Constructor
193 //
263d466a 194 fTracksEsdParam.SetOwner(kTRUE);
c6914c83 195}
196
197AliTPCcalibLaser::~AliTPCcalibLaser() {
198 //
199 // destructor
200 //
201}
202
e9f38a4b 203
204
205void AliTPCcalibLaser::Process(AliESDEvent * event) {
c6914c83 206 //
263d466a 207 //
e9f38a4b 208 // Loop over tracks and call Process function
209 //
210 fESD = event;
211 if (!fESD) {
212 return;
213 }
214 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
215 if (!fESDfriend) {
216 return;
217 }
263d466a 218 AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
e9f38a4b 219 fTracksTPC.Clear();
220 fTracksEsd.Clear();
221 fTracksEsdParam.Delete();
222 //
223 Int_t n=fESD->GetNumberOfTracks();
224 Int_t run = fESD->GetRunNumber();
225 fRun = run;
226 for (Int_t i=0;i<n;++i) {
227 AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
228 AliESDtrack *track=fESD->GetTrack(i);
229 TObject *calibObject=0;
230 AliTPCseed *seed=0;
231 for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
232 if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
233 break;
234 if (track&&seed) FindMirror(track,seed);
235 //
236 }
263d466a 237
e9f38a4b 238 FitDriftV();
c18f4385 239 MakeDistHisto();
e9f38a4b 240 //
c18f4385 241 for (Int_t id=0; id<336; id++){
e9f38a4b 242 //
243 //
244 if (!fTracksEsdParam.At(id)) continue;
245 DumpLaser(id);
263d466a 246// RefitLaser(id);
247 RefitLaserJW(id);
248
c18f4385 249 }
263d466a 250// fEvent++;
c18f4385 251}
252
253void AliTPCcalibLaser::MakeDistHisto(){
254 //
255 //
256 //
257 for (Int_t id=0; id<336; id++){
258 //
259 //
263d466a 260 if (!fTracksEsdParam.At(id)) continue;
c18f4385 261 if (!AcceptLaser(id)) continue;
262 //
263 //
264 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
e5c59218 265 TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
266 TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
267
c18f4385 268 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
269 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
2b35e8e6 270 TH1F * hisSignal = (TH1F*)fSignals.At(id);
271
263d466a 272 if (!hisdz){
c03e3250 273 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
274 hisdz->SetDirectory(0);
c18f4385 275 fDeltaZ.AddAt(hisdz,id);
e5c59218 276
277 hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
278 hisP3->SetDirectory(0);
279 fDeltaP3.AddAt(hisP3,id);
280 //
281 hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
282 hisP4->SetDirectory(0);
283 fDeltaP4.AddAt(hisP4,id);
284
c18f4385 285 //
c03e3250 286 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
287 hisdphi->SetDirectory(0);
c18f4385 288 fDeltaPhi.AddAt(hisdphi,id);
289 //
c03e3250 290 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
291 hisdphiP->SetDirectory(0);
c18f4385 292 fDeltaPhiP.AddAt(hisdphiP,id);
e5c59218 293 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
c03e3250 294 hisSignal->SetDirectory(0);
2b35e8e6 295 fSignals.AddAt(hisSignal,id);
c18f4385 296 }
297
298 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
299 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
2b35e8e6 300 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
c03e3250 301 if (!param) return;
302 if (!ltrp) return;
303 if (!track) return;
c18f4385 304 Double_t xyz[3];
305 Double_t pxyz[3];
306 Double_t lxyz[3];
307 Double_t lpxyz[3];
308 param->GetXYZ(xyz);
309 param->GetPxPyPz(pxyz);
310 ltrp->GetXYZ(lxyz);
311 ltrp->GetPxPyPz(lpxyz);
312 //
313 Float_t dz = param->GetZ()-ltrp->GetZ();
314 Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
315 Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
c03e3250 316 if (hisdz) hisdz->Fill(dz);
e5c59218 317 if (hisP3) hisP3->Fill(param->GetParameter()[3]);
318 if (hisP4) hisP4->Fill(param->GetParameter()[4]);
c03e3250 319 if (hisdphi) hisdphi->Fill(dphi);
263d466a 320 if (hisdphiP) hisdphiP->Fill(dphiP);
e5c59218 321 if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
c18f4385 322 }
323}
324
325void AliTPCcalibLaser::FitDriftV(){
326 //
327 // Fit drift velocity - linear approximation in the z and global y
328 //
329 static TLinearFitter fdriftA(3,"hyp2");
330 static TLinearFitter fdriftC(3,"hyp2");
331 fdriftA.ClearPoints();
332 fdriftC.ClearPoints();
333 //
334 for (Int_t id=0; id<336; id++){
263d466a 335 if (!fTracksEsdParam.At(id)) continue;
c18f4385 336 if (!AcceptLaser(id)) continue;
337 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
338 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
339 Double_t xyz[3];
340 Double_t pxyz[3];
341 Double_t lxyz[3];
342 Double_t lpxyz[3];
343 param->GetXYZ(xyz);
344 param->GetPxPyPz(pxyz);
345 ltrp->GetXYZ(lxyz);
346 ltrp->GetPxPyPz(lpxyz);
347 Double_t xxx[2] = {lxyz[2],lxyz[1]};
348 if (ltrp->GetSide()==0){
349 fdriftA.AddPoint(xxx,xyz[2],1);
350 }else{
351 fdriftC.AddPoint(xxx,xyz[2],1);
352 }
353 }
354 Float_t chi2A = 0;
355 Float_t chi2C = 0;
356 Int_t npointsA=0;
357 Int_t npointsC=0;
358 //
359 if (fdriftA.GetNpoints()>10){
360 fdriftA.Eval();
361 fdriftA.EvalRobust(0.8);
362 fdriftA.GetParameters(*fFitAside);
363 npointsA= fdriftA.GetNpoints();
364 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
365 }
366 if (fdriftC.GetNpoints()>10){
367 fdriftC.Eval();
368 fdriftC.EvalRobust(0.8);
369 fdriftC.GetParameters(*fFitCside);
370 npointsC= fdriftC.GetNpoints();
371 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
e9f38a4b 372 }
263d466a 373
c18f4385 374 if (fStreamLevel>0){
375 TTreeSRedirector *cstream = GetDebugStreamer();
376 Int_t time = fESD->GetTimeStamp();
377 if (cstream){
378 (*cstream)<<"driftv"<<
379 "driftA.="<<fFitAside<<
380 "driftC.="<<fFitCside<<
381 "chi2A="<<chi2A<<
382 "chi2C="<<chi2C<<
383 "nA="<<npointsA<<
384 "nC="<<npointsC<<
385 "time="<<time<<
386 "\n";
387 }
388 }
389 //
390}
391
392
393Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
394 //
395 //
396 //
397 /*
e5c59218 398 TCut cutP0("cutP0","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03");
399 TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
400 TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
401 TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
402 TCut cutP4("cutPt","abs(Tr.fP[4])<0.1");
403
404 TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
c18f4385 405 */
406 AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
407 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
408 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
e5c59218 409 Double_t xyz[3];
410 Double_t lxyz[3];
411 param->GetXYZ(xyz);
412 ltrp->GetXYZ(lxyz);
413 if (TMath::Abs(TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))>0.03) return kFALSE; //cut y- P0
414 if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE; // cutZ -P1
415 if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE; // cut -P2
416 if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE; // cut Tl -P3
417 if (TMath::Abs(param->GetParameter()[4])>0.1) return kFALSE; // cut Pt -P4
c18f4385 418 //
c18f4385 419 //
e5c59218 420
263d466a 421 return kTRUE;
e9f38a4b 422}
423
424Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
425 //
426 // Find corresponding mirror
427 // add the corresponding tracks
c6914c83 428 //
1bd37cc0 429 Float_t kRadius0 = 252;
c18f4385 430 Float_t kRadius = 253.4;
e9f38a4b 431 if (!track->GetOuterParam()) return -1;
c6914c83 432 AliExternalTrackParam param(*(track->GetOuterParam()));
1bd37cc0 433 AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
434 AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
7b18d067 435 AliTPCLaserTrack ltr;
436 AliTPCLaserTrack *ltrp=0x0;
5732d662 437 // AliTPCLaserTrack *ltrpjw=0x0;
1bd37cc0 438 //
263d466a 439 Int_t id = AliTPCLaserTrack::IdentifyTrack(&param);
440 // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(&param);
441 //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
442
443 if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
1bd37cc0 444 ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
263d466a 445 else
1bd37cc0 446 ltrp=&ltr;
263d466a 447
448
e9f38a4b 449 if (id>=0){
c18f4385 450 //
451 //
452 Float_t radius=TMath::Abs(ltrp->GetX());
453 AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
454 //
e9f38a4b 455 if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
456 fTracksEsdParam.AddAt(param.Clone(),id);
457 fTracksEsd.AddAt(track,id);
458 fTracksTPC.AddAt(seed,id);
c18f4385 459 //
e9f38a4b 460 }
c18f4385 461
e9f38a4b 462 return id;
463}
464
465
466
467void AliTPCcalibLaser::DumpLaser(Int_t id) {
468 //
469 // Dump Laser info to the tree
470 //
471 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
472 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
473 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1bd37cc0 474 //
475 // Fast laser ID
476 //
477 Double_t xyz[3];
478 Double_t pxyz[3];
c18f4385 479 Double_t lxyz[3];
480 Double_t lpxyz[3];
e9f38a4b 481 param->GetXYZ(xyz);
482 param->GetPxPyPz(pxyz);
c18f4385 483 ltrp->GetXYZ(lxyz);
484 ltrp->GetPxPyPz(lpxyz);
485
c6914c83 486 if (fStreamLevel>0){
487 TTreeSRedirector *cstream = GetDebugStreamer();
c18f4385 488 Int_t time = fESD->GetTimeStamp();
489 Bool_t accept = AcceptLaser(id);
c6914c83 490 if (cstream){
491 (*cstream)<<"Track"<<
e9f38a4b 492 "run="<<fRun<<
7b18d067 493 "id="<<id<<
c18f4385 494 "accept="<<accept<<
495 "driftA.="<<fFitAside<<
496 "driftC.="<<fFitCside<<
497 "time="<<time<<
7b18d067 498 //
499 "LTr.="<<ltrp<<
500 "Esd.="<<track<<
e9f38a4b 501 "Tr.="<<param<<
c6914c83 502 "x0="<<xyz[0]<<
503 "x1="<<xyz[1]<<
504 "x2="<<xyz[2]<<
505 "px0="<<pxyz[0]<<
506 "px1="<<pxyz[1]<<
507 "px2="<<pxyz[2]<<
c18f4385 508 //
509 "lx0="<<lxyz[0]<<
510 "lx1="<<lxyz[1]<<
511 "lx2="<<lxyz[2]<<
512 "lpx0="<<lpxyz[0]<<
513 "lpx1="<<lpxyz[1]<<
514 "lpx2="<<lpxyz[2]<<
c6914c83 515 "\n";
516 }
517 }
518}
519
263d466a 520void AliTPCcalibLaser::RefitLaserJW(Int_t id){
521 //
522 // Refit the track with different tracklet models:
523 // 1. Per ROC using the kalman filter, different edge cuts
524 // 2. Per ROC linear in y and z
525 // 3. Per ROC quadratic in y and z
526 // 4. Per track offset for each sector, linear for each sector, common quadratic
527 // store x, y, z information for all models and the cluster to calculate the residuals
528 //
529 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
530 AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
531 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
532
533 AliTPCclusterMI dummyCl;
534
535 //two tracklets
536 Int_t kMaxTracklets=2;
0728a4f6 537 //=============================================//
538 // Linear Fitters for the Different Approaches //
539 //=============================================//
540 //linear fit model in y and z; inner - outer sector
541 static TLinearFitter fy1I(2,"hyp1");
542 static TLinearFitter fy1O(2,"hyp1");
543 static TLinearFitter fz1I(2,"hyp1");
544 static TLinearFitter fz1O(2,"hyp1");
545 //quadratic fit model in y and z; inner - sector
546 static TLinearFitter fy2I(3,"hyp2");
547 static TLinearFitter fy2O(3,"hyp2");
548 static TLinearFitter fz2I(3,"hyp2");
549 static TLinearFitter fz2O(3,"hyp2");
263d466a 550 //common quadratic fit for IROC and OROC in y and z
551 static TLinearFitter fy4(5,"hyp4");
552 static TLinearFitter fz4(5,"hyp4");
553
554
555 //set standard cuts
556 if ( fNcuts==0 ){
557 fNcuts=1;
558 fEdgeXcuts[0]=4;
559 fEdgeYcuts[0]=3;
560 fNClCuts[0]=20;
561 }
0728a4f6 562 //=============================//
563 // Loop over all Tracklet Cuts //
564 //=============================//
263d466a 565 for (Int_t icut=0; icut<fNcuts; icut++){
566 AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
567 //cut parameters
568 Double_t edgeCutX = fEdgeXcuts[icut];
569 Double_t edgeCutY = fEdgeYcuts[icut];
e5c59218 570 Int_t nclCut = (Int_t)fNClCuts[icut];
0728a4f6 571 //=========================//
572 // Parameters to calculate //
573 //=========================//
574 //fit parameter inner, outer tracklet and combined fit
575 TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
576 TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
263d466a 577 //
0728a4f6 578 TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
579 TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
263d466a 580 TVectorD vecy4res(5),vecz4res(5);
581 // cluster and track positions for each row - used for residuals
582 TVectorD vecX(159); // x is the same for all (row center)
583 TVectorD vecYkalman(159); // y from kalman fit
584 TVectorD vecZkalman(159); // z from kalman fit
585 TVectorD vecY1(159); // y from pol1 fit per ROC
586 TVectorD vecZ1(159); // z from pol1 fit per ROC
587 TVectorD vecY2(159); // y from pol2 fit per ROC
588 TVectorD vecZ2(159); // z from pol2 fit per ROC
589 TVectorD vecY4(159); // y from sector fit
590 TVectorD vecZ4(159); // z from sector fit
591 TVectorD vecClY(159); // y cluster position
592 TVectorD vecClZ(159); // z cluster position
593 TVectorD vecSec(159); // sector for each row
0728a4f6 594 //chi2 of fits
263d466a 595 Double_t chi2I1z=-1; // chi2 of pol1 fit in z (inner)
596 Double_t chi2I1y=-1; // chi2 of pol1 fit in y (inner)
597 Double_t chi2O1z=-1; // chi2 of pol1 fit in z (outer)
598 Double_t chi2O1y=-1; // chi2 of pol1 fit in y (outer)
599 Double_t chi2I2z=-1; // chi2 of pol2 fit in z (inner)
600 Double_t chi2I2y=-1; // chi2 of pol2 fit in y (inner)
601 Double_t chi2O2z=-1; // chi2 of pol2 fit in z (outer)
602 Double_t chi2O2y=-1; // chi2 of pol2 fit in y (outer)
603 Double_t chi2IOz=-1; // chi2 of hyp4 fit in z (inner+outer)
604 Double_t chi2IOy=-1; // chi2 of hyp4 fit in y (inner+outer)
0728a4f6 605 //more
606 Int_t innerSector = -1; // number of inner sector
607 Int_t outerSector = -1; // number of outer sector
263d466a 608 Int_t nclI=0; // number of clusters (inner)
609 Int_t nclO=0; // number of clusters (outer)
0728a4f6 610 //=================================================//
611 // Perform the Kalman Fit using the Tracklet Class //
612 //=================================================//
263d466a 613 AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
614 TObjArray tracklets=
615 AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
616 kFALSE,nclCut,kMaxTracklets);
263d466a 617 // tracklet pointers
263d466a 618 AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
619 AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
0728a4f6 620 AliTPCTracklet *tr=0x0;
621 AliTPCTracklet dummy;
622 //continue if we didn't find a tracklet
263d466a 623 if ( !trInner && !trOuter ) continue;
0728a4f6 624 //================================================================================//
625 // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
626 //================================================================================//
263d466a 627 if ( trInner && trOuter ){
628 if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
629 if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
630 tr = trInner;
631 trInner=trOuter;
632 trOuter=tr;
633 AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
634 }
635 } else {
636 if ( trInner ){
637 if ( !trInner->GetInner() ) continue;
638 trOuter=&dummy;
639 if ( trInner->GetSector()>35 ){
640 trOuter=trInner;
641 trInner=&dummy;
642 }
643 } else { //trOuter
644 if ( !trOuter->GetInner() ) continue;
645 trInner=&dummy;
646 if ( trOuter->GetSector()<36 ){
647 trInner=trOuter;
648 trOuter=&dummy;
649 }
650 }
651 }
652 innerSector = trInner->GetSector();
653 if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
654 outerSector = trOuter->GetSector();
655 if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
656
263d466a 657 // array of clusters
658 TClonesArray arrCl("AliTPCclusterMI",159);
659 arrCl.ExpandCreateFast(159);
0728a4f6 660 //=======================================//
661 // fill fitters with cluster information //
662 //=======================================//
663 AliDebug(3,"Filing Cluster Information");
263d466a 664 for (Int_t irow=158;irow>-1;--irow) {
665 AliTPCclusterMI *c=track->GetClusterPointer(irow);
0728a4f6 666 AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
667 cl=dummyCl;
668 //
263d466a 669 vecSec[irow]=-1;
670 if (!c) continue;
e5c59218 671 Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
672 Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
673
0728a4f6 674 //
675 Int_t roc = static_cast<Int_t>(c->GetDetector());
676 if ( roc!=innerSector && roc!=outerSector ) continue;
677 vecSec[irow]=roc;
263d466a 678 //store clusters in clones array
679 cl=*c;
680 //cluster position
681 vecX[irow] = c->GetX();
682 vecClY[irow] = c->GetY();
683 vecClZ[irow] = c->GetZ();
263d466a 684 //
0728a4f6 685 Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
686 Double_t y = vecClY[irow];
687 Double_t z = vecClZ[irow];
688 //
689 Double_t x2[2]={x,x*x}; //linear and parabolic parameters
690 Double_t x4[4]={0,0,0,0}; //hyp4 parameters
263d466a 691 if ( roc == innerSector ) {
0728a4f6 692 x4[0]=1; //offset inner - outer sector
693 x4[1]=x; //slope parameter inner sector
263d466a 694 } else {
0728a4f6 695 x4[2]=x; //slope parameter outer sector
263d466a 696 }
0728a4f6 697 x4[3]=x*x; //common parabolic shape
e5c59218 698 if (pedgeX < fEdgeXcuts[icut]) continue;
699 if (pedgeY < fEdgeYcuts[icut]) continue;
263d466a 700 //
0728a4f6 701 if ( roc==innerSector ){
702 fy1I.AddPoint(x2,y);
703 fz1I.AddPoint(x2,z);
704 fy2I.AddPoint(x2,y);
705 fz2I.AddPoint(x2,z);
706 ++nclI;
707 }
708 if ( roc==outerSector ){
709 fy1O.AddPoint(x2,y);
710 fz1O.AddPoint(x2,z);
711 fy2O.AddPoint(x2,y);
712 fz2O.AddPoint(x2,z);
713 ++nclO;
714 }
263d466a 715 fy4.AddPoint(x4,y);
716 fz4.AddPoint(x4,z);
717 }
0728a4f6 718 //======================================//
719 // Evaluate and retrieve fit parameters //
720 //======================================//
721 AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
722 //inner sector
723 if ( (innerSector>-1) && (fy1I.GetNpoints()>0) ){
724 fy1I.Eval();
725 fz1I.Eval();
726 fy2I.Eval();
727 fz2I.Eval();
728 fy1I.GetParameters(vecy1resInner);
729 fz1I.GetParameters(vecz1resInner);
730 fy2I.GetParameters(vecy2resInner);
731 fz2I.GetParameters(vecz2resInner);
732 chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
733 chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
734 chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
735 chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
736 }
737 //outer sector
738 if ( (outerSector>-1) && (fy1O.GetNpoints()>0) ){
739 fy1O.Eval();
740 fz1O.Eval();
741 fy2O.Eval();
742 fz2O.Eval();
743 fy1O.GetParameters(vecy1resOuter);
744 fz1O.GetParameters(vecz1resOuter);
745 fy2O.GetParameters(vecy2resOuter);
746 fz2O.GetParameters(vecz2resOuter);
747 chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
748 chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
749 chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
750 chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
751 }
752 //combined hyp4 fit
263d466a 753 if ( innerSector>0 && outerSector>0 ){
0728a4f6 754 if (fy4.GetNpoints()>0) {
755 fy4.Eval();
756 fy4.GetParameters(vecy4res);
757 chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
758 }
759 if (fz4.GetNpoints()>0) {
760 fz4.Eval();
761 fz4.GetParameters(vecz4res);
762 chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
763 }
263d466a 764 }
0728a4f6 765 //clear points
766 fy4.ClearPoints(); fz4.ClearPoints();
767 fy1I.ClearPoints(); fy1O.ClearPoints();
768 fz1I.ClearPoints(); fz1O.ClearPoints();
769 fy2I.ClearPoints(); fy2O.ClearPoints();
770 fz2I.ClearPoints(); fz2O.ClearPoints();
771 //==============================//
772 // calculate tracklet positions //
773 //==============================//
263d466a 774 AliDebug(4,"Calculate tracklet positions");
775 for (Int_t irow=158;irow>-1;--irow) {
0728a4f6 776 if ( vecSec[irow]==-1 ) continue; //no cluster info
263d466a 777 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
778 tr=&dummy;
779 Double_t x=vecX[irow];
0728a4f6 780 Double_t xref=x-133.4;
781 //
263d466a 782 Double_t yoffInner=0;
783 Double_t zoffInner=0;
784 Double_t yslopeInner=0;
785 Double_t yslopeOuter=0;
786 Double_t zslopeInner=0;
0728a4f6 787 Double_t zslopeOuter=0;
788 //positions of hyperplane fits
263d466a 789 if ( vecSec[irow] == outerSector ) {
790 tr=trOuter;
791 vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
792 vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
793 vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
794 vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
795 yslopeOuter=vecy4res[3];
796 zslopeOuter=vecz4res[3];
263d466a 797 } else {
798 tr=trInner;
799 vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
800 vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
801 vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
802 vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
803 yoffInner=vecy4res[1];
804 zoffInner=vecz4res[1];
805 yslopeInner=vecy4res[2];
806 zslopeInner=vecz4res[2];
263d466a 807 }
808 vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
809 vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
0728a4f6 810 //positions of kalman fits
263d466a 811 Double_t gxyz[3],xyz[3];
812 AliExternalTrackParam *param = 0x0;
0728a4f6 813 //
263d466a 814 param=tr->GetInner();
815 if (param){
816 param->GetXYZ(gxyz);
817 Float_t bz = AliTracker::GetBz(gxyz);
818 param->GetYAt(x, bz, xyz[1]);
819 param->GetZAt(x, bz, xyz[2]);
820 vecYkalman[irow]=xyz[1];
821 vecZkalman[irow]=xyz[2];
822 }
263d466a 823 }
0728a4f6 824 //=====================================================================//
825 // write results from the different tracklet fits with debug streamers //
826 //=====================================================================//
263d466a 827 if (fStreamLevel>4){
828 TTreeSRedirector *cstream = GetDebugStreamer();
829 if (cstream){
830 Float_t dedx = track->GetdEdx();
831 (*cstream)<<"FitModels"<<
832 "cutNr=" << icut <<
833 "edgeCutX=" << edgeCutX <<
834 "edgeCutY=" << edgeCutY <<
835 "nclCut=" << nclCut <<
836 "innerSector="<< innerSector <<
837 "outerSector="<< outerSector <<
838 "dEdx=" << dedx <<
839 "LTr.=" << ltrp <<
840 "Tr.=" << extparam <<
841 "yPol1In.=" << &vecy1resInner <<
842 "zPol1In.=" << &vecz1resInner <<
843 "yPol2In.=" << &vecy2resInner <<
844 "zPol2In.=" << &vecz2resInner <<
845 "yPol1Out.=" << &vecy1resOuter <<
846 "zPol1Out.=" << &vecz1resOuter <<
847 "yPol2Out.=" << &vecy2resOuter <<
848 "zPol2Out.=" << &vecz2resOuter <<
849 "yInOut.=" << &vecy4res <<
850 "zInOut.=" << &vecz4res <<
851 "chi2y1In=" << chi2I1y <<
852 "chi2z1In=" << chi2I1z <<
853 "chi2y1Out=" << chi2O1y <<
854 "chi2z1Out=" << chi2O1z <<
855 "chi2y2In=" << chi2I2y <<
856 "chi2z2In=" << chi2I2z <<
857 "chi2y2Out=" << chi2O2y <<
858 "chi2z2Out=" << chi2O2z <<
859 "chi2yInOut=" << chi2IOy <<
860 "chi2zInOut=" << chi2IOz <<
861 "trletIn.=" << trInner <<
862 "trletOut.=" << trOuter <<
863 "nclI=" << nclI <<
864 "nclO=" << nclO <<
865 "\n";
866 }
867 }
868
869 // wirte residuals information
870 if (fStreamLevel>5){
871 TTreeSRedirector *cstream = GetDebugStreamer();
872 if (cstream){
873 Float_t dedx = track->GetdEdx();
874 (*cstream)<<"Residuals"<<
875 "cutNr=" << icut <<
876 "edgeCutX=" << edgeCutX <<
877 "edgeCutY=" << edgeCutY <<
878 "nclCut=" << nclCut <<
879 "LTr.=" << ltrp <<
880 "Tr.=" << extparam<<
881 "dEdx=" << dedx <<
882 "Cl.=" << &arrCl <<
883 "TrX.=" << &vecX <<
884 "TrYpol1.=" << &vecY1 <<
885 "TrZpol1.=" << &vecZ1 <<
886 "TrYpol2.=" << &vecY2 <<
887 "TrZpol2.=" << &vecZ2 <<
888 "TrYInOut.=" << &vecY4 <<
889 "TrZInOut.=" << &vecZ4 <<
890 "ClY.=" << &vecClY <<
891 "ClZ.=" << &vecClZ <<
0728a4f6 892 "sec.=" << &vecSec <<
263d466a 893 "nclI=" << nclI <<
894 "nclO=" << nclO <<
895 "yInOut.=" << &vecy4res <<
896 "zInOut.=" << &vecz4res <<
897 "chi2y1In=" << chi2I1y <<
898 "chi2z1In=" << chi2I1z <<
899 "chi2y1Out=" << chi2O1y <<
900 "chi2z1Out=" << chi2O1z <<
901 "chi2y2In=" << chi2I2y <<
902 "chi2z2In=" << chi2I2z <<
903 "chi2y2Out=" << chi2O2y <<
904 "chi2z2Out=" << chi2O2z <<
905 "chi2yInOut=" << chi2IOy <<
906 "chi2zInOut=" << chi2IOz <<
907 "\n";
908
909 }
910 }
95a0e09b 911 //==========================//
912 // Fill Residual Histograms //
913 //==========================//
914 TProfile *profy = (TProfile*)fDeltaYres.UncheckedAt(id);
915 TProfile *profz = (TProfile*)fDeltaZres.UncheckedAt(id);
916 if (!profy){
917 profy=new TProfile(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d",id),115,80,250);
5732d662 918 profy->SetDirectory(0);
95a0e09b 919 fDeltaYres.AddAt(profy,id);
920 }
921 if (!profz){
922 profz=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),115,80,250);
5732d662 923 profz->SetDirectory(0);
95a0e09b 924 fDeltaZres.AddAt(profz,id);
925 }
926 for (Int_t irow=158;irow>-1;--irow) {
927 if (vecSec[irow]==-1)continue; //no cluster info
928 Double_t x = vecX[irow];
929 Double_t ycl = vecClY[irow];
930 Double_t yfit = vecY1[irow];
931 Double_t zcl = vecClZ[irow];
932 Double_t zfit = vecZ1[irow];
5732d662 933 if (profy)
934 if (profy->GetEntries()<1000000)
935 profy->Fill(x,yfit-ycl);
936 if (profz)
937 if (profz->GetEntries()<1000000)
938 profz->Fill(x,zfit-zcl);
95a0e09b 939 }
ecb7e705 940 //===============================//
941 // Fill Fit Parameter Histograms //
942 //===============================//
943 TH1F *pol2InnerY =(TH1F*)fPol2Par2InY.UncheckedAt(id);
944 TH1F *diff1InnerY=(TH1F*)fDiffPar1InY.UncheckedAt(id);
945 TH1F *pol2OuterY =(TH1F*)fPol2Par2OutY.UncheckedAt(id);
946 TH1F *diff1OuterY=(TH1F*)fDiffPar1OutY.UncheckedAt(id);
947 TH1F *pol2InnerZ =(TH1F*)fPol2Par2InZ.UncheckedAt(id);
948 TH1F *diff1InnerZ=(TH1F*)fDiffPar1InZ.UncheckedAt(id);
949 TH1F *pol2OuterZ =(TH1F*)fPol2Par2OutZ.UncheckedAt(id);
950 TH1F *diff1OuterZ=(TH1F*)fDiffPar1OutZ.UncheckedAt(id);
951 //create histograms if the do not already exist
952 if (!pol2InnerY){
953 pol2InnerY =new TH1F(Form("pol2par2inY%03d",id),
954 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (inner sector)",id),
674ee45c 955 500,-.005,.005);
5732d662 956 pol2InnerY->SetDirectory(0);
ecb7e705 957 pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
958 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
674ee45c 959 500,0.01,.01);
5732d662 960 pol2OuterY->SetDirectory(0);
961
ecb7e705 962 diff1InnerY=new TH1F(Form("diff1inY%03d",id),
963 Form("diff of 1st derivative from pol1 and pol2 in Y fit for Laser beam %03d (inner sector)",id),
674ee45c 964 500,-.5,.5);
5732d662 965 diff1InnerY->SetDirectory(0);
966
ecb7e705 967 diff1OuterY=new TH1F(Form("diff1outY%03d",id),
968 Form("diff of 1st derivative from pol1 and pol2 in Yfit for Laser beam %03d (outer sector)",id),
674ee45c 969 500,-1,1);
5732d662 970 diff1OuterY->SetDirectory(0);
971
ecb7e705 972 pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
973 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
674ee45c 974 500,-.002,.002);
5732d662 975 pol2InnerZ->SetDirectory(0);
976
ecb7e705 977 pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
978 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
674ee45c 979 500,-.005,.005);
5732d662 980 pol2OuterZ->SetDirectory(0);
ecb7e705 981 diff1InnerZ=new TH1F(Form("diff1inZ%03d",id),
982 Form("diff of 1st derivative from pol1 and pol2 in Z fit for Laser beam %03d (inner sector)",id),
674ee45c 983 500,-.02,.02);
5732d662 984 diff1InnerZ->SetDirectory(0);
ecb7e705 985 diff1OuterZ=new TH1F(Form("diff1outZ%03d",id),
986 Form("diff of 1st derivative from pol1 and pol2 in Zfit for Laser beam %03d (outer sector)",id),
674ee45c 987 500,-.03,.03);
5732d662 988 diff1OuterZ->SetDirectory(0);
ecb7e705 989 //add
990 fPol2Par2InY.AddAt(pol2InnerY,id);
991 fDiffPar1InY.AddAt(diff1InnerY,id);
992 fPol2Par2OutY.AddAt(pol2OuterY,id);
993 fDiffPar1OutY.AddAt(diff1OuterY,id);
994 fPol2Par2InZ.AddAt(pol2InnerZ,id);
995 fDiffPar1InZ.AddAt(diff1InnerZ,id);
996 fPol2Par2OutZ.AddAt(pol2OuterZ,id);
997 fDiffPar1OutZ.AddAt(diff1OuterZ,id);
998 }
999 //fill histograms
1000 pol2InnerY ->Fill(vecy2resInner[2]);
1001 pol2OuterY ->Fill(vecy2resOuter[2]);
1002 diff1InnerY->Fill(vecy2resInner[1]-vecy1resInner[1]);
1003 diff1OuterY->Fill(vecy2resOuter[1]-vecy1resOuter[1]);
1004 pol2InnerZ ->Fill(vecz2resInner[2]);
1005 pol2OuterZ ->Fill(vecz2resOuter[2]);
1006 diff1InnerZ->Fill(vecz2resInner[1]-vecz1resInner[1]);
1007 diff1OuterZ->Fill(vecz2resOuter[1]-vecz1resOuter[1]);
95a0e09b 1008
263d466a 1009
ecb7e705 1010
1011 }
263d466a 1012}
1013
e9f38a4b 1014
1015void AliTPCcalibLaser::RefitLaser(Int_t id){
1016 //
1017 // Refit the track store residuals
1018 //
1019
1020 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
1021 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1022 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
263d466a 1023
1024 //linear fit model in y and z per sector
e9f38a4b 1025 static TLinearFitter fy1(2,"hyp1");
1026 static TLinearFitter fz1(2,"hyp1");
263d466a 1027 //quadratic fit model in y and z per sector
1028 static TLinearFitter fy2(3,"hyp2");
1029 static TLinearFitter fz2(3,"hyp2");
e9f38a4b 1030 static TVectorD vecy2,vecz2,vecy1,vecz1;
1031
1032 const Int_t kMinClusters=20;
263d466a 1033 Int_t nclusters[72];
e9f38a4b 1034 //
1035 for (Int_t i=0;i<72;++i) nclusters[i]=0;
1036
263d466a 1037 for (Int_t i=0;i<160;++i) {
e9f38a4b 1038 AliTPCclusterMI *c=track->GetClusterPointer(i);
1039 if (c) nclusters[c->GetDetector()]++;
1040 }
263d466a 1041
e9f38a4b 1042 for (Int_t isec=0; isec<72;isec++){
1043 if (nclusters[isec]<kMinClusters) continue;
1044 fy2.ClearPoints();
1045 fz2.ClearPoints();
1046 fy1.ClearPoints();
1047 fz1.ClearPoints();
1048 //
263d466a 1049 for (Int_t irow=0;irow<160;++irow) {
e9f38a4b 1050 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1051 //if (c && RejectCluster(c)) continue;
263d466a 1052 Double_t xd = c->GetX()-133.4; // reference x is beteen iroc and oroc
e9f38a4b 1053 if (c&&c->GetDetector()==isec) {
e9f38a4b 1054 Double_t x[2]={xd,xd*xd};
1055 fy2.AddPoint(x,c->GetY());
1056 fz2.AddPoint(x,c->GetZ());
1057 //
1058 fy1.AddPoint(x,c->GetY());
1059 fz1.AddPoint(x,c->GetZ());
1060 }
1061 }
1062 fy2.Eval();
1063 fz2.Eval();
1064 fy1.Eval();
1065 fz1.Eval();
1066 fy1.GetParameters(vecy1);
1067 fy2.GetParameters(vecy2);
1068 fz1.GetParameters(vecz1);
1069 fz2.GetParameters(vecz2);
263d466a 1070
e9f38a4b 1071 if (fStreamLevel>0){
1072 TTreeSRedirector *cstream = GetDebugStreamer();
1073 if (cstream){
1074 Float_t dedx = track->GetdEdx();
1075 (*cstream)<<"Tracklet"<<
1076 "LTr.="<<ltrp<<
1077 "Tr.="<<param<<
1078 "sec="<<isec<<
1079 "ncl="<<nclusters[isec]<<
1080 "dedx="<<dedx<<
1081 "dedx="<<dedx<<
1082 "vy1.="<<&vecy1<<
1083 "vy2.="<<&vecy2<<
1084 "vz1.="<<&vecz1<<
1085 "vz2.="<<&vecz2<<
1086 "\n";
1087 }
1088 }
1089 }
1090 //
1091 //
1092 //
263d466a 1093 // for (Int_t irow=0;irow<160;++irow) {
e9f38a4b 1094 // AliTPCclusterMI *c=track->GetClusterPointer(irow);
1095 // if (c && RejectCluster(c)) continue;
1096 // if (c&&c->GetDetector()==isec) {
1097 // Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
1098 // fy2.AddPoint(&x,c->GetY());
1099 // fz2.AddPoint(&x,c->GetZ());
1100 // //
1101 // fy1.AddPoint(&x,c->GetY());
1102 // fz1.AddPoint(&x,c->GetZ());
1103 // }
263d466a 1104 // }
1105
e9f38a4b 1106}
1107
1108
e5c59218 1109void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run, Int_t minEntries){
2b35e8e6 1110 //
1fd56785 1111 // Dump information about laser beams
1112 // isOK variable indicates usability of the beam
1113 // Beam is not usable if:
1114 // a. No entries in range (krmsCut0)
1115 // b. Big sperad (krmscut1)
1116 // c. RMSto fit sigma bigger then (kmultiCut)
1117 // d. Too big angular spread
1118 //
1119
1120 const Float_t krmsCut0=0.001;
1121 const Float_t krmsCut1=0.16;
1122 const Float_t kmultiCut=2;
1123 const Float_t kcutP0=0.002;
2b35e8e6 1124 //
1125 AliTPCcalibLaser *laser = this;
1126 TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
1127 TF1 fg("fg","gaus");
1128 //
1129 //
1130 for (Int_t id=0; id<336; id++){
1fd56785 1131 Bool_t isOK=kTRUE;
2b35e8e6 1132 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
1133 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
1134 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
e5c59218 1135 TH1F * hisP3 = (TH1F*)laser->fDeltaP3.At(id);
1136 TH1F * hisP4 = (TH1F*)laser->fDeltaP4.At(id);
1fd56785 1137 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
2b35e8e6 1138 if (!hisphi) continue;;
1139 Double_t entries = hisphi->GetEntries();
1fd56785 1140 if (entries<minEntries) continue;
2b35e8e6 1141 //
1142 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1143 if (!ltrp) {
1144 AliTPCLaserTrack::LoadTracks();
1145 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1146 }
e5c59218 1147 pcstream->GetFile()->cd();
1148 if (hisphi) hisphi->Write();
1149 if (hisphiP) hisphiP->Write();
1150 if (hisZ) hisZ->Write();
1151 if (hisP3) hisP3->Write();
1152 if (hisP4) hisP4->Write();
1153
2b35e8e6 1154 Float_t meanphi = hisphi->GetMean();
1155 Float_t rmsphi = hisphi->GetRMS();
1fd56785 1156 //
2b35e8e6 1157 Float_t meanphiP = hisphiP->GetMean();
1158 Float_t rmsphiP = hisphiP->GetRMS();
1159 Float_t meanZ = hisZ->GetMean();
1160 Float_t rmsZ = hisZ->GetRMS();
1161 hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
263d466a 1162 Double_t gphi1 = fg.GetParameter(1);
1163 Double_t gphi2 = fg.GetParameter(2);
2b35e8e6 1164 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
263d466a 1165 Double_t gphiP1 = fg.GetParameter(1);
1166 Double_t gphiP2 = fg.GetParameter(2);
e5c59218 1167 //
2b35e8e6 1168 hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
263d466a 1169 Double_t gz1 = fg.GetParameter(1);
1170 Double_t gz2 = fg.GetParameter(2);
1fd56785 1171 //
e5c59218 1172 hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
1173 Double_t gp31 = fg.GetParameter(1);
1174 Double_t gp32 = fg.GetParameter(2);
1175 //
1176 hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
1177 Double_t gp41 = fg.GetParameter(1);
1178 Double_t gp42 = fg.GetParameter(2);
1179 //
1fd56785 1180 Float_t meanS=hisS->GetMean();
2b35e8e6 1181 //
1182 Double_t lxyz[3];
1183 Double_t lpxyz[3];
1184 ltrp->GetXYZ(lxyz);
1185 ltrp->GetPxPyPz(lpxyz);
1fd56785 1186
1187 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
1188 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
1189 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
1190 if (gphiP2>kcutP0) isOK=kFALSE;
2b35e8e6 1191 //
e5c59218 1192 if (run<=0) run=fRun;
2b35e8e6 1193 (*pcstream)<<"Mean"<<
e5c59218 1194 "run="<<run<< //
1fd56785 1195 "isOK="<<isOK<<
2b35e8e6 1196 "entries="<<entries<< // number of entries
1197 "bz="<<bfield<< // bfield
1198 "LTr.="<<ltrp<< // refernece track
1fd56785 1199 //
2b35e8e6 1200 "lx0="<<lxyz[0]<< // reference x
1201 "lx1="<<lxyz[1]<< // reference y
263d466a 1202 "lx2="<<lxyz[2]<< // refernece z
2b35e8e6 1203 "lpx0="<<lpxyz[0]<< // reference x
1204 "lpx1="<<lpxyz[1]<< // reference y
263d466a 1205 "lpx2="<<lpxyz[2]<< // refernece z
1fd56785 1206 //
1207 "msig="<<meanS<<
1208 //
2b35e8e6 1209 "mphi="<<meanphi<< //
1210 "rmsphi="<<rmsphi<< //
1211 "gphi1="<<gphi1<<
1212 "gphi2="<<gphi2<<
1fd56785 1213 //
2b35e8e6 1214 "mphiP="<<meanphiP<< //
1215 "rmsphiP="<<rmsphiP<< //
1216 "gphiP1="<<gphiP1<<
1217 "gphiP2="<<gphiP2<<
1fd56785 1218 //
2b35e8e6 1219 "meanZ="<<meanZ<<
1220 "rmsZ="<<rmsZ<<
1221 "gz1="<<gz1<<
1222 "gz2="<<gz2<<
e5c59218 1223 //
1224 "gp31="<<gp31<< //gaus mean - tgl
1225 "gp32="<<gp32<< //gaus rms -tgl
1226 "gp41="<<gp41<< //gaus mean - P4
1227 "gp42="<<gp42<< //gaus rms - P4
2b35e8e6 1228
1229 "\n";
1230 }
1231 delete pcstream;
1232}
1233
1fd56785 1234
1235
1236void AliTPCcalibLaser::DumpScanInfo(TTree * chain){
1237 //
1238 //
1239 //
1240 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
1241 TFile * f = pcstream->GetFile();
1242 f->mkdir("dirphi");
1243 f->mkdir("dirphiP");
1244 f->mkdir("dirZ");
1245 TF1 fp("p1","pol1");
1246 //
1247 //
1248 char cut[1000];
1249 char grname[1000];
1250 char grnamefull[1000];
1251
1252 Double_t mphi[100];
1253 Double_t mphiP[100];
1254 Double_t smphi[100];
1255 Double_t smphiP[100];
1256 Double_t mZ[100];
1257 Double_t smZ[100];
1258 Double_t bz[100];
1259 Double_t sbz[100];
1260 // fit parameters
1261 Double_t pphi[3];
1262 Double_t pphiP[3];
1263 Double_t pmZ[3];
1264 //
1265 for (Int_t id=0; id<336; id++){
1266 // id =205;
1267 sprintf(cut,"isOK&&fId==%d",id);
1268 Int_t entries = chain->Draw("bz",cut,"goff");
1269 if (entries<3) continue;
e5c59218 1270 AliTPCLaserTrack *ltrp = 0;
1fd56785 1271 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
1272 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1273 Double_t lxyz[3];
1274 Double_t lpxyz[3];
1275 ltrp->GetXYZ(lxyz);
1276 ltrp->GetPxPyPz(lpxyz);
1277
1278 chain->Draw("bz",cut,"goff");
1279 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
1280 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
1281 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
1282 //
1283 chain->Draw("gphi1",cut,"goff");
1284 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
1285 chain->Draw("0.05*abs(mphi)+gphi2",cut,"goff");
1286 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
1287 //
1288 chain->Draw("gphiP1",cut,"goff");
1289 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
1290 chain->Draw("0.05*abs(mphiP)+gphiP2",cut,"goff");
1291 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
1292 //
1293 chain->Draw("gz1",cut,"goff");
1294 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
1295 chain->Draw("0.01*abs(meanZ)+gz2",cut,"goff");
1296 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
1297 //
1298 //
1299 sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
1300 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
1301 // store data
1302 // phi
1303 f->cd("dirphi");
1304 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
1305 grphi->Draw("a*");
1306 grphi->Fit(&fp);
1307 pphi[0] = fp.GetParameter(0); // offset
1308 pphi[1] = fp.GetParameter(1); // slope
1309 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1310 sprintf(grname,"phi_id%d",id);
1311 grphi->SetName(grname); grphi->SetTitle(grnamefull);
1312 grphi->GetXaxis()->SetTitle("b_{z} (T)");
1313 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
ae69f16f 1314 grphi->SetMaximum(1.2);
1315 grphi->SetMinimum(-1.2);
1fd56785 1316 grphi->Draw("a*");
1317
1318 grphi->Write();
1319 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
1320 // phiP
1321 f->cd("dirphiP)");
1322 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
1323 grphiP->Draw("a*");
1324 grphiP->Fit(&fp);
1325 pphiP[0] = fp.GetParameter(0); // offset
1326 pphiP[1] = fp.GetParameter(1); // slope
1327 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1328 sprintf(grname,"phiP_id%d",id);
1329 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
1330 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
1331 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
ae69f16f 1332 grphiP->SetMaximum(pphiP[0]+0.005);
1333 grphiP->SetMinimum(pphiP[0]-0.005);
1fd56785 1334
1335 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
1336 grphiP->Write();
1337 //
1338 //Z
1339 f->cd("dirZ");
1340 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
1341 grmZ->Draw("a*");
1342 grmZ->Fit(&fp);
1343 pmZ[0] = fp.GetParameter(0); // offset
1344 pmZ[1] = fp.GetParameter(1); // slope
1345 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1346 sprintf(grname,"mZ_id%d",id);
1347 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
1348 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
1349 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
1350
1351 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
1352 grmZ->Write();
1353
1354
1355 for (Int_t ientry=0; ientry<entries; ientry++){
1356 (*pcstream)<<"Mean"<<
1357 "id="<<id<<
1358 "LTr.="<<ltrp<<
1359 "entries="<<entries<<
1360 "bz="<<bz[ientry]<<
1361 "lx0="<<lxyz[0]<< // reference x
1362 "lx1="<<lxyz[1]<< // reference y
1363 "lx2="<<lxyz[2]<< // refernece z
1364 "lpx0="<<lpxyz[0]<< // reference x
1365 "lpx1="<<lpxyz[1]<< // reference y
1366 "lpx2="<<lpxyz[2]<< // refernece z
1367 //values
1368 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
1369 "pphi0="<<pphi[0]<< // offset
1370 "pphi1="<<pphi[1]<< // mean
1371 "pphi2="<<pphi[2]<< // norm chi2
1372 //
1373 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
1374 "pphiP0="<<pphiP[0]<< // offset
1375 "pphiP1="<<pphiP[1]<< // mean
1376 "pphiP2="<<pphiP[2]<< // norm chi2
1377 //
1378 "gz1="<<mZ[ientry]<<
1379 "pmZ0="<<pmZ[0]<< // offset
1380 "pmZ1="<<pmZ[1]<< // mean
1381 "pmZ2="<<pmZ[2]<< // norm chi2
1382 "\n";
1383 }
1384 }
1385
1386 delete pcstream;
1387
1388}
1389
1390
c6914c83 1391void AliTPCcalibLaser::Analyze(){
1392 //
1393 //
1394 //
1395}
1396
1397
c03e3250 1398Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
c6914c83 1399
c03e3250 1400 TIterator* iter = li->MakeIterator();
1401 AliTPCcalibLaser* cal = 0;
c6914c83 1402
c03e3250 1403 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
1404 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
1405 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1406 return -1;
1407 }
1408 //
1409 // fHistNTracks->Add(cal->fHistNTracks);
1410// fClusters->Add(cal->fClusters);
1411// fModules->Add(cal->fModules);
1412// fHistPt->Add(cal->fHistPt);
1413// fPtResolution->Add(cal->fPtResolution);
1414// fDeDx->Add(cal->fDeDx);
263d466a 1415
c03e3250 1416
1417 TH1F *h=0x0;
1418 TH1F *hm=0x0;
95a0e09b 1419 TProfile *hp=0x0;
1420 TProfile *hpm=0x0;
c03e3250 1421
1422 for (Int_t id=0; id<336; id++){
1423 // merge fDeltaZ histograms
263d466a 1424 hm = (TH1F*)cal->fDeltaZ.At(id);
c03e3250 1425 h = (TH1F*)fDeltaZ.At(id);
1426 if (!h) {
1427 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
5732d662 1428 h->SetDirectory(0);
c03e3250 1429 fDeltaZ.AddAt(h,id);
1430 }
1431 if (hm) h->Add(hm);
e5c59218 1432 // merge fP3 histograms
1433 hm = (TH1F*)cal->fDeltaP3.At(id);
1434 h = (TH1F*)fDeltaP3.At(id);
1435 if (!h) {
1436 h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
1437 h->SetDirectory(0);
1438 fDeltaP3.AddAt(h,id);
1439 }
1440 if (hm) h->Add(hm);
1441 // merge fP4 histograms
1442 hm = (TH1F*)cal->fDeltaP4.At(id);
1443 h = (TH1F*)fDeltaP4.At(id);
1444 if (!h) {
1445 h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
1446 h->SetDirectory(0);
1447 fDeltaP4.AddAt(h,id);
1448 }
1449 if (hm) h->Add(hm);
1450
1451 //
c03e3250 1452 // merge fDeltaPhi histograms
263d466a 1453 hm = (TH1F*)cal->fDeltaPhi.At(id);
c03e3250 1454 h = (TH1F*)fDeltaPhi.At(id);
1455 if (!h) {
1456 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
5732d662 1457 h->SetDirectory(0);
c03e3250 1458 fDeltaPhi.AddAt(h,id);
1459 }
1460 if (hm) h->Add(hm);
1461 // merge fDeltaPhiP histograms
263d466a 1462 hm = (TH1F*)cal->fDeltaPhiP.At(id);
c03e3250 1463 h = (TH1F*)fDeltaPhiP.At(id);
1464 if (!h) {
1465 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
5732d662 1466 h->SetDirectory(0);
c03e3250 1467 fDeltaPhiP.AddAt(h,id);
1468 }
1469 if (hm) h->Add(hm);
1470 // merge fSignals histograms
263d466a 1471 hm = (TH1F*)cal->fSignals.At(id);
c03e3250 1472 h = (TH1F*)fSignals.At(id);
1473 if (!h) {
e5c59218 1474 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
5732d662 1475 h->SetDirectory(0);
c03e3250 1476 fSignals.AddAt(h,id);
1477 }
263d466a 1478 if (hm) h->Add(hm);
95a0e09b 1479 //
1480 //
1481 // merge ProfileY histograms
1482 hpm = (TProfile*)cal->fDeltaYres.At(id);
1483 hp = (TProfile*)fDeltaYres.At(id);
1484 if (!hp) {
5732d662 1485 hp=new TProfile(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d",id),115,80,250);
1486 hp->SetDirectory(0);
95a0e09b 1487 fDeltaYres.AddAt(hp,id);
1488 }
1489 if (hpm) hp->Add(hpm);
1490 //
1491 hpm = (TProfile*)cal->fDeltaZres.At(id);
1492 hp = (TProfile*)fDeltaZres.At(id);
1493 if (!hp) {
5732d662 1494 hp=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),115,80,250);
1495 hp->SetDirectory(0);
95a0e09b 1496 fDeltaZres.AddAt(hp,id);
1497 }
1498 if (hpm) hp->Add(hpm);
1499 //
1500 //
ecb7e705 1501 //merge fit param histograms
1502 //local hists
1503 TH1F *pol2InnerY =(TH1F*)fPol2Par2InY.UncheckedAt(id);
1504 TH1F *diff1InnerY=(TH1F*)fDiffPar1InY.UncheckedAt(id);
1505 TH1F *pol2OuterY =(TH1F*)fPol2Par2OutY.UncheckedAt(id);
1506 TH1F *diff1OuterY=(TH1F*)fDiffPar1OutY.UncheckedAt(id);
1507 TH1F *pol2InnerZ =(TH1F*)fPol2Par2InZ.UncheckedAt(id);
1508 TH1F *diff1InnerZ=(TH1F*)fDiffPar1InZ.UncheckedAt(id);
1509 TH1F *pol2OuterZ =(TH1F*)fPol2Par2OutZ.UncheckedAt(id);
1510 TH1F *diff1OuterZ=(TH1F*)fDiffPar1OutZ.UncheckedAt(id);
1511 //hists to merge
1512 TH1F *pol2InnerYm =(TH1F*)cal->fPol2Par2InY.UncheckedAt(id);
1513 TH1F *diff1InnerYm=(TH1F*)cal->fDiffPar1InY.UncheckedAt(id);
1514 TH1F *pol2OuterYm =(TH1F*)cal->fPol2Par2OutY.UncheckedAt(id);
1515 TH1F *diff1OuterYm=(TH1F*)cal->fDiffPar1OutY.UncheckedAt(id);
1516 TH1F *pol2InnerZm =(TH1F*)cal->fPol2Par2InZ.UncheckedAt(id);
1517 TH1F *diff1InnerZm=(TH1F*)cal->fDiffPar1InZ.UncheckedAt(id);
1518 TH1F *pol2OuterZm =(TH1F*)cal->fPol2Par2OutZ.UncheckedAt(id);
1519 TH1F *diff1OuterZm=(TH1F*)cal->fDiffPar1OutZ.UncheckedAt(id);
1520 //create histos if they do not exist
1521 if (!pol2InnerY){
1522 pol2InnerY =new TH1F(Form("pol2par2inY%03d",id),
1523 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (inner sector)",id),
674ee45c 1524 500,-.005,.005);
5732d662 1525 pol2InnerY->SetDirectory(0);
1526
ecb7e705 1527 pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
1528 Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
674ee45c 1529 500,0.01,.01);
5732d662 1530 pol2OuterY->SetDirectory(0);
ecb7e705 1531 diff1InnerY=new TH1F(Form("diff1inY%03d",id),
1532 Form("diff of 1st derivative from pol1 and pol2 in Y fit for Laser beam %03d (inner sector)",id),
674ee45c 1533 500,-.5,.5);
ecb7e705 1534 diff1OuterY=new TH1F(Form("diff1outY%03d",id),
1535 Form("diff of 1st derivative from pol1 and pol2 in Yfit for Laser beam %03d (outer sector)",id),
674ee45c 1536 500,-1,1);
5732d662 1537 diff1InnerY->SetDirectory(0);
1538
1539
ecb7e705 1540 pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
1541 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
674ee45c 1542 500,-.002,.002);
5732d662 1543 pol2InnerZ->SetDirectory(0);
1544
ecb7e705 1545 pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
1546 Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
674ee45c 1547 500,-.005,.005);
5732d662 1548 pol2OuterZ->SetDirectory(0);
ecb7e705 1549 diff1InnerZ=new TH1F(Form("diff1inZ%03d",id),
1550 Form("diff of 1st derivative from pol1 and pol2 in Z fit for Laser beam %03d (inner sector)",id),
674ee45c 1551 500,-.02,.02);
5732d662 1552 diff1InnerZ->SetDirectory(0);
ecb7e705 1553 diff1OuterZ=new TH1F(Form("diff1outZ%03d",id),
1554 Form("diff of 1st derivative from pol1 and pol2 in Zfit for Laser beam %03d (outer sector)",id),
674ee45c 1555 500,-.03,.03);
5732d662 1556 diff1OuterZ->SetDirectory(0);
ecb7e705 1557 }
1558 pol2InnerYm =(TH1F*)cal->fPol2Par2InY.UncheckedAt(id);
1559 diff1InnerYm=(TH1F*)cal->fDiffPar1InY.UncheckedAt(id);
1560 pol2OuterYm =(TH1F*)cal->fPol2Par2OutY.UncheckedAt(id);
1561 diff1OuterYm=(TH1F*)cal->fDiffPar1OutY.UncheckedAt(id);
1562 pol2InnerZm =(TH1F*)cal->fPol2Par2InZ.UncheckedAt(id);
1563 diff1InnerZm=(TH1F*)cal->fDiffPar1InZ.UncheckedAt(id);
1564 pol2OuterZm =(TH1F*)cal->fPol2Par2OutZ.UncheckedAt(id);
1565 diff1OuterZm=(TH1F*)cal->fDiffPar1OutZ.UncheckedAt(id);
c03e3250 1566 }
1567 }
1568 return 0;
1569}
1570
e5c59218 1571void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
1572 //
1573 // Dump fit information - collect information from the streamers
1574 //
1575 /*
1576 TChain * chainFit=0;
1577 TChain * chainTrack=0;
1578 TChain * chain=0;
1579 //
1580 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
1581 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
1582 AliXRDPROOFtoolkit tool;
1583 chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
1584 chainTrack->Lookup();
1585 chainTrack->SetProof(kTRUE);
1586 chain = tool.MakeChain("laser.txt","Residuals",0,10200);
1587 chain->Lookup();
1588 chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
1589 chainFit->Lookup();
1590 chainFit->SetProof(kTRUE);
1591 chain->SetProof(kTRUE);
1592 AliTPCLaserTrack::LoadTracks();
1593 //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
1594
1595 */
1596 //
1597 // Fit cuts
1598 //
1599 TCut cutP3("abs(Tr.fP[3])<0.1");
1600 TCut cutP4("abs(Tr.fP[4])<0.5");
1601 TCut cutPx = cutP3+cutP4;
1602 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
1603 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
1604 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
1605 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
1606 //
1607 TCut cutdEdx("sqrt(dEdx)>3");
1608 TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
1609 TCut cutN("nclO>20&&nclI>20");
1610 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx;
1611 //
1612 // Cluster cuts
1613 //
1614 TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
1615 TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
1616 TCut cutClX("abs(Cl[].fX)>10");
1617 TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
1618 TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
1619 TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
1620 TCut cutQ("sqrt(Cl[].fMax)>4");
1621 TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
1622
1623
1624 TH1F * phisAl = 0;
1625 TH1F * phisAccept = 0;
1626 TH1F * phisOut = 0;
1627 TProfile * pdEdx = 0;
1628
1629 TProfile * pP0 = 0;
1630 TProfile * pP1 = 0;
1631 TProfile * pP2 = 0;
1632 TProfile * pP3 = 0;
1633 TProfile * pP4 = 0;
1634 //
1635 TProfile * pNclI = 0;
1636 TProfile * pNclO = 0;
1637 //
1638 TProfile * pchi2YIn =0;
1639 TProfile * pchi2ZIn =0;
1640 TProfile * pchi2YOut =0;
1641 TProfile * pchi2ZOut =0;
1642 TProfile * pchi2YInOut =0;
1643 TProfile * pchi2ZInOut =0;;
1644 // laser counters
1645 chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
1646 phisAl = (TH1F*)gROOT->FindObject("hisAl");
1647 chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
1648 phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
1649 chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
1650 phisOut = (TH1F*)gROOT->FindObject("hisOut");
1651 //
1652 chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
1653 pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
1654 // track param
1655 //
1656 chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
1657 pP0 = (TProfile*)gROOT->FindObject("hP0");
1658 chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
1659 pP1 = (TProfile*)gROOT->FindObject("hP1");
1660 chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
1661 pP2 = (TProfile*)gROOT->FindObject("hP2");
1662 chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
1663 pP3 = (TProfile*)gROOT->FindObject("hP3");
1664 chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
1665 pP4 = (TProfile*)gROOT->FindObject("hP4");
1666 //
1667 chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
1668 pNclI = (TProfile*)gROOT->FindObject("hNclI");
1669 chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
1670 pNclO = (TProfile*)gROOT->FindObject("hNclO");
1671 //
1672 //
1673 chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
1674 pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
1675 chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
1676 pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
1677 chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
1678 pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
1679 chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
1680 pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
1681 chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
1682 pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
1683 chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
1684 pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
1685 //
1686 // second derivatives
1687 //
1688 TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
1689 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
1690 TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
1691 chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
1692 TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
1693 chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
1694 //
1695 phisPy2In->FitSlicesY();
1696 TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
1697 TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
1698 TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
1699 //
1700 phisPy2Out->FitSlicesY();
1701 TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
1702 TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
1703 TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
1704 //
1705 phisPy2InOut->FitSlicesY();
1706 TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
1707 TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
1708 TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
1709 //
1710 TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
1711 chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
1712 TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
1713 chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
1714 TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
1715 chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
1716 //
1717 phisPz2In->FitSlicesY();
1718 TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
1719 TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
1720 TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
1721 //
1722 phisPz2Out->FitSlicesY();
1723 TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
1724 TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
1725 TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
1726 //
1727 phisPz2InOut->FitSlicesY();
1728 TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
1729 TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
1730 TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
1731 //
1732 //
1733 //
1734
1735
1736 {
1737 TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
1738 for (Int_t ilaser=0; ilaser<336; ilaser++){
1739 Float_t all=phisAl->GetBinContent(ilaser+1);
1740 Float_t accept=phisAccept->GetBinContent(ilaser+1);
1741 Float_t out=phisOut->GetBinContent(ilaser+1);
1742 Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
1743 Float_t mP0 = pP0->GetBinContent(ilaser+1);
1744 Float_t mP1 = pP1->GetBinContent(ilaser+1);
1745 Float_t mP2 = pP2->GetBinContent(ilaser+1);
1746 Float_t mP3 = pP3->GetBinContent(ilaser+1);
1747 Float_t mP4 = pP4->GetBinContent(ilaser+1);
1748
1749
1750 Float_t nclI = pNclI->GetBinContent(ilaser+1);
1751 Float_t nclO = pNclO->GetBinContent(ilaser+1);
1752 //
1753 Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
1754 Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
1755 Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
1756 Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
1757 Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
1758 Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
1759 //
1760 Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
1761 Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
1762 Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
1763 //
1764 Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
1765 Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
1766 Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
1767 //
1768 Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
1769 Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
1770 Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
1771 //
1772 Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
1773 Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
1774 Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
1775 //
1776 Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
1777 Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
1778 Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
1779 //
1780 Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
1781 Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
1782 Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
1783
1784 AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
1785 (*pcstream)<<"Scan"<<
1786 "Run="<<id<<
1787 "LTr.="<<ltrp<<
1788 "all="<<all<<
1789 "accept="<<accept<<
1790 "out="<<out<<
1791 "sdedx="<<sdedx<<
1792 "mP0="<<mP0<<
1793 "mP1="<<mP1<<
1794 "mP2="<<mP2<<
1795 "mP3="<<mP3<<
1796 "mP4="<<mP4<<
1797 "nclI="<<nclI<<
1798 "nclO="<<nclO<<
1799 "chi2YIn="<<chi2YIn<<
1800 "chi2YOut="<<chi2YOut<<
1801 "chi2YInOut="<<chi2YInOut<<
1802 "chi2ZIn="<<chi2ZIn<<
1803 "chi2ZOut="<<chi2ZOut<<
1804 "chi2ZInOut="<<chi2ZInOut<<
1805 //
1806 "nPy2In="<<entriesPy2In<<
1807 "mPy2In="<<meanPy2In<<
1808 "sPy2In="<<sigmaPy2In<<
1809 //
1810 "nPy2Out="<<entriesPy2Out<<
1811 "mPy2Out="<<meanPy2Out<<
1812 "sPy2Out="<<sigmaPy2Out<<
1813 //
1814 "nPy2InOut="<<entriesPy2InOut<<
1815 "mPy2InOut="<<meanPy2InOut<<
1816 "sPy2InOut="<<sigmaPy2InOut<<
1817 //
1818 "nPz2In="<<entriesPz2In<<
1819 "mPz2In="<<meanPz2In<<
1820 "sPz2In="<<sigmaPz2In<<
1821 //
1822 "nPz2Out="<<entriesPz2Out<<
1823 "mPz2Out="<<meanPz2Out<<
1824 "sPz2Out="<<sigmaPz2Out<<
1825 //
1826 "nPz2InOut="<<entriesPz2InOut<<
1827 "mPz2InOut="<<meanPz2InOut<<
1828 "sPz2InOut="<<sigmaPz2InOut<<
1829 "\n";
1830 }
1831
1832 delete pcstream;
1833 }
1834 /*
1835 TFile f("vscan.root");
1836 */
1837
1838 /*
1839 pad binning effect
1840 chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
1841 //
1842 chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
1843 //
1844
1845chain->Draw("Cl.fY-TrYpol1.fElements-AliTPCClusterParam::SPosCorrection(0,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
1846
1847
1848chain->Draw("Cl[].fZ-TrZpol1.fElements-0*AliTPCClusterParam::SPosCorrection(1,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fTimeBin-int(Cl[].fTimeBin)",cutA+cutCl+"Cl[].fZ>0","prof",10000)
1849
1850 */
1851
1852
1853
1854
1855
1856 /*
1857 // check edge effects
1858 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
1859 //
1860 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
1861
1862 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
1863
1864
1865
1866 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
1867 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
1868
1869*/
1870
1871
1872
1873
1874
1875 /*
1876 Edge y effect
1877
1878 dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
1879
1880
1881 chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):pi/18-abs(Cl.fY/Cl.fX)>>hisYdphi(100,0,0.03)",""+cutA+cutCl,"prof",10000)
1882
1883 chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdy(100,0,5)",""+cutA+cutCl,"prof",10000)
1884
1885
1886
1887
1888
1889 chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyIROC(100,0,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
1890
1891 chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyOROC(100,0,5)","Cl.fDetector>36"+cutA+cutCl,"prof",100000)
1892
1893
1894
1895 chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>his(100,-5,5)",""+cutA+cutCl,"prof",100000)
1896
1897 chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>hisdyInner(100,-5,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
1898
1899
1900
1901*/
1902
1903
1904/*
1905
1906chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutDY,"prof")
1907
1908chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&&LTr.fP[1]<0"+cutA+cutDY,"prof")
1909
1910
1911
1912chainFit->Draw("LTr.fId","nclI>10",100000)
1913
1914chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
1915
1916chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
1917
1918TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
1919
1920*/
1921
1922}
1923
1924
1925
c03e3250 1926
1927
263d466a 1928
c03e3250 1929/*
1930 gSystem->Load("libSTAT.so")
1931 TStatToolkit toolkit;
1932 Double_t chi2;
1933 TVectorD fitParam;
1934 TMatrixD covMatrix;
1935 Int_t npoints;
e5c59218 1936
ae69f16f 1937 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
c03e3250 1938
1939
1940TString fstring="";
1fd56785 1941//
1942fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
1943fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
1944fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
1945fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
1946//
1947fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
1948fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
1949fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
1950fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
1951//
1952fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
1953fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
1954fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
1955fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
c5595838 1956
c5595838 1957
c5595838 1958
c5595838 1959
1fd56785 1960 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c03e3250 1961
1fd56785 1962 treeT->SetAlias("fit",strq0->Data());
1963
c03e3250 1964
1fd56785 1965 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c5595838 1966
1fd56785 1967 treeT->SetAlias("fitP",strqP->Data());
c5595838 1968
1969
ae69f16f 1970 TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
1971 treeT->SetAlias("fitD",strqDrift->Data());
c5595838 1972
c03e3250 1973
ae69f16f 1974treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
1975{
1976for (Int_t i=0; i<6;i++){
1977treeT->SetLineColor(i+2);
1978treeT->SetMarkerSize(1);
1979treeT->SetMarkerStyle(22+i);
1980treeT->SetMarkerColor(i+2);
c03e3250 1981
ae69f16f 1982treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
1983}
1984}
c03e3250 1985 */
68ff0583 1986
1987
1988
1989/*
1990 TTree * tree = (TTree*)f.Get("FitModels");
1991
1992 TEventList listLFit0("listLFit0","listLFit0");
1993 TEventList listLFit1("listLFit1","listLFit1");
68ff0583 1994 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
1995 tree->SetEventList(&listLFit0);
1996
1997
95a0e09b 1998
1999
2000 gSystem->Load("libSTAT.so")
2001 TStatToolkit toolkit;
2002 Double_t chi2;
2003 TVectorD fitParam;
2004 TMatrixD covMatrix;
2005 Int_t npoints;
2006
2007 chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
2008 chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
2009
2010
2011 TString fstring="";
2012 fstring+="cos(dp)++";
2013 fstring+="sin(dp)++";
2014 fstring+="cos(dt)++";
2015 fstring+="sin(dt)++";
2016
2017 TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
2018
2019
e5c59218 2020
2021*/
2022
2023
2024
2025/*
2026 Edge effects
2027 //
2028 //
2029 //
2030 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2031 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2032 AliXRDPROOFtoolkit tool;
2033 TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
2034 chainTrack->Lookup();
2035 chainTrack->SetProof(kTRUE);
2036
2037 TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
2038 chain->Lookup();
2039 TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
2040 chainFit->Lookup();
2041 chainFit->SetProof(kTRUE);
2042 chain->SetProof(kTRUE);
2043 //
2044 // Fit cuts
2045 //
2046 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
2047 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
2048 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
2049 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
2050 //
2051 TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
2052 TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
2053 TCut cutN("nclO>20&&nclI>20");
2054 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
2055 //
2056 // Cluster cuts
2057 //
2058 TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
2059 TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
2060 TCut cutClX("abs(Cl.fX)>10");
2061 TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
2062 TCut cutCl=cutClY+cutClZ+cutClX;
2063
2064
2065 // check edge effects
2066 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
2067 //
2068 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
2069
2070 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
2071
2072
2073
2074 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
2075 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
2076
68ff0583 2077*/
e5c59218 2078