1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
20 // 1. The laser track is associated with the mirror
21 // see function FindMirror
23 // 2. The laser track is accepted for the analysis under certain condition
24 // (see function Accpet laser)
26 // 3. The drift velocity and jitter is calculated event by event
27 // (see function drift velocity)
31 // To make laser scan the user interaction neccessary
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")
42 // laser track clasification;
44 TCut cutT("cutT","abs(Tr.fP[3])<0.06");
45 TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
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");
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);
60 AliTPCcalibLaser::DumpScanInfo(chain)
61 TFile fscan("laserScan.root")
62 TTree * treeT = (TTree*)fscan.Get("Mean")
68 #include "TLinearFitter.h"
69 #include "AliTPCcalibLaser.h"
70 #include "AliExternalTrackParam.h"
71 #include "AliESDEvent.h"
72 #include "AliESDfriend.h"
73 #include "AliESDtrack.h"
74 #include "AliTPCTracklet.h"
77 #include "TTreeStream.h"
80 #include "TGraphErrors.h"
81 #include "AliTPCclusterMI.h"
82 #include "AliTPCseed.h"
83 #include "AliTracker.h"
85 #include "TClonesArray.h"
89 #include "TTreeStream.h"
92 #include "AliTPCLaserTrack.h"
96 ClassImp(AliTPCcalibLaser)
98 AliTPCcalibLaser::AliTPCcalibLaser():
104 fTracksEsdParam(336),
110 fFitAside(new TVectorD(3)),
111 fFitCside(new TVectorD(3)),
122 fTracksEsdParam.SetOwner(kTRUE);
125 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
131 fTracksEsdParam(336),
133 fDeltaZ(336), // array of histograms of delta z for each track
134 fDeltaPhi(336), // array of histograms of delta z for each track
135 fDeltaPhiP(336), // array of histograms of delta z for each track
136 fSignals(336), // array of dedx signals
137 fFitAside(new TVectorD(3)), // drift fit - A side
138 fFitCside(new TVectorD(3)), // drift fit - C- side
139 fEdgeXcuts(5), // cuts in local x direction; used in the refit of the laser tracks
140 fEdgeYcuts(5), // cuts in local y direction; used in the refit of the laser tracks
141 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
142 fNcuts(0), // number of cuts
151 fTracksEsdParam.SetOwner(kTRUE);
154 AliTPCcalibLaser::~AliTPCcalibLaser() {
162 void AliTPCcalibLaser::Process(AliESDEvent * event) {
165 // Loop over tracks and call Process function
171 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
175 AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
178 fTracksEsdParam.Delete();
180 Int_t n=fESD->GetNumberOfTracks();
181 Int_t run = fESD->GetRunNumber();
183 for (Int_t i=0;i<n;++i) {
184 AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
185 AliESDtrack *track=fESD->GetTrack(i);
186 TObject *calibObject=0;
188 for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
189 if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
191 if (track&&seed) FindMirror(track,seed);
198 for (Int_t id=0; id<336; id++){
201 if (!fTracksEsdParam.At(id)) continue;
210 void AliTPCcalibLaser::MakeDistHisto(){
214 for (Int_t id=0; id<336; id++){
217 if (!fTracksEsdParam.At(id)) continue;
218 if (!AcceptLaser(id)) continue;
221 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
222 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
223 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
224 TH1F * hisSignal = (TH1F*)fSignals.At(id);
227 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
228 hisdz->SetDirectory(0);
229 fDeltaZ.AddAt(hisdz,id);
231 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
232 hisdphi->SetDirectory(0);
233 fDeltaPhi.AddAt(hisdphi,id);
235 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
236 hisdphiP->SetDirectory(0);
237 fDeltaPhiP.AddAt(hisdphiP,id);
238 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
239 hisSignal->SetDirectory(0);
240 fSignals.AddAt(hisSignal,id);
243 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
244 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
245 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
254 param->GetPxPyPz(pxyz);
256 ltrp->GetPxPyPz(lpxyz);
258 Float_t dz = param->GetZ()-ltrp->GetZ();
259 Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
260 Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
261 if (hisdz) hisdz->Fill(dz);
262 if (hisdphi) hisdphi->Fill(dphi);
263 if (hisdphiP) hisdphiP->Fill(dphiP);
264 if (hisSignal) hisSignal->Fill(track->GetTPCsignal());
268 void AliTPCcalibLaser::FitDriftV(){
270 // Fit drift velocity - linear approximation in the z and global y
272 static TLinearFitter fdriftA(3,"hyp2");
273 static TLinearFitter fdriftC(3,"hyp2");
274 fdriftA.ClearPoints();
275 fdriftC.ClearPoints();
277 for (Int_t id=0; id<336; id++){
278 if (!fTracksEsdParam.At(id)) continue;
279 if (!AcceptLaser(id)) continue;
280 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
281 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
287 param->GetPxPyPz(pxyz);
289 ltrp->GetPxPyPz(lpxyz);
290 Double_t xxx[2] = {lxyz[2],lxyz[1]};
291 if (ltrp->GetSide()==0){
292 fdriftA.AddPoint(xxx,xyz[2],1);
294 fdriftC.AddPoint(xxx,xyz[2],1);
302 if (fdriftA.GetNpoints()>10){
304 fdriftA.EvalRobust(0.8);
305 fdriftA.GetParameters(*fFitAside);
306 npointsA= fdriftA.GetNpoints();
307 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
309 if (fdriftC.GetNpoints()>10){
311 fdriftC.EvalRobust(0.8);
312 fdriftC.GetParameters(*fFitCside);
313 npointsC= fdriftC.GetNpoints();
314 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
318 TTreeSRedirector *cstream = GetDebugStreamer();
319 Int_t time = fESD->GetTimeStamp();
321 (*cstream)<<"driftv"<<
322 "driftA.="<<fFitAside<<
323 "driftC.="<<fFitCside<<
336 Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
341 TCut cutT("cutT","abs(Tr.fP[3])<0.06");
342 TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
343 TCut cutN("cutN","fTPCncls>70");
344 TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
345 TCut cutA = cutT+cutPt+cutP;
347 AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
348 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
349 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
351 if (TMath::Abs(param->GetParameter()[4])>0.03) return kFALSE;
352 if (TMath::Abs(param->GetParameter()[3])>0.06) return kFALSE;
353 if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.06) return kFALSE;
354 if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>10) return kFALSE;
358 if (TMath::Abs(track->GetTPCsignal())<20) return kFALSE;
359 if (TMath::Abs(track->GetTPCsignal())>800) return kFALSE;
364 Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
366 // Find corresponding mirror
367 // add the corresponding tracks
369 Float_t kRadius0 = 252;
370 Float_t kRadius = 253.4;
371 if (!track->GetOuterParam()) return -1;
372 AliExternalTrackParam param(*(track->GetOuterParam()));
373 AliTracker::PropagateTrackTo(¶m,kRadius0,0.10566,3,kTRUE);
374 AliTracker::PropagateTrackTo(¶m,kRadius,0.10566,0.1,kTRUE);
375 AliTPCLaserTrack ltr;
376 AliTPCLaserTrack *ltrp=0x0;
377 AliTPCLaserTrack *ltrpjw=0x0;
379 Int_t id = AliTPCLaserTrack::IdentifyTrack(¶m);
380 // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(¶m);
381 //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
383 if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
384 ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
388 if (idjw!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(idjw)))
389 ltrpjw=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(idjw);
395 TTreeSRedirector *cstream = GetDebugStreamer();
397 (*cstream)<<"idcmp"<<
414 Float_t radius=TMath::Abs(ltrp->GetX());
415 AliTracker::PropagateTrackTo(¶m,radius,0.10566,0.01,kTRUE);
417 if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
418 fTracksEsdParam.AddAt(param.Clone(),id);
419 fTracksEsd.AddAt(track,id);
420 fTracksTPC.AddAt(seed,id);
429 void AliTPCcalibLaser::DumpLaser(Int_t id) {
431 // Dump Laser info to the tree
433 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
434 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
435 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
444 param->GetPxPyPz(pxyz);
446 ltrp->GetPxPyPz(lpxyz);
449 TTreeSRedirector *cstream = GetDebugStreamer();
450 Int_t time = fESD->GetTimeStamp();
451 Bool_t accept = AcceptLaser(id);
453 (*cstream)<<"Track"<<
457 "driftA.="<<fFitAside<<
458 "driftC.="<<fFitCside<<
482 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
484 // Refit the track with different tracklet models:
485 // 1. Per ROC using the kalman filter, different edge cuts
486 // 2. Per ROC linear in y and z
487 // 3. Per ROC quadratic in y and z
488 // 4. Per track offset for each sector, linear for each sector, common quadratic
489 // store x, y, z information for all models and the cluster to calculate the residuals
491 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
492 AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
493 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
495 AliTPCclusterMI dummyCl;
498 Int_t kMaxTracklets=2;
499 //=============================================//
500 // Linear Fitters for the Different Approaches //
501 //=============================================//
502 //linear fit model in y and z; inner - outer sector
503 static TLinearFitter fy1I(2,"hyp1");
504 static TLinearFitter fy1O(2,"hyp1");
505 static TLinearFitter fz1I(2,"hyp1");
506 static TLinearFitter fz1O(2,"hyp1");
507 //quadratic fit model in y and z; inner - sector
508 static TLinearFitter fy2I(3,"hyp2");
509 static TLinearFitter fy2O(3,"hyp2");
510 static TLinearFitter fz2I(3,"hyp2");
511 static TLinearFitter fz2O(3,"hyp2");
512 //common quadratic fit for IROC and OROC in y and z
513 static TLinearFitter fy4(5,"hyp4");
514 static TLinearFitter fz4(5,"hyp4");
524 //=============================//
525 // Loop over all Tracklet Cuts //
526 //=============================//
527 for (Int_t icut=0; icut<fNcuts; icut++){
528 AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
530 Double_t edgeCutX = fEdgeXcuts[icut];
531 Double_t edgeCutY = fEdgeYcuts[icut];
532 Int_t nclCut = fNClCuts[icut];
533 //=========================//
534 // Parameters to calculate //
535 //=========================//
536 //fit parameter inner, outer tracklet and combined fit
537 TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
538 TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
540 TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
541 TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
542 TVectorD vecy4res(5),vecz4res(5);
543 // cluster and track positions for each row - used for residuals
544 TVectorD vecX(159); // x is the same for all (row center)
545 TVectorD vecYkalman(159); // y from kalman fit
546 TVectorD vecZkalman(159); // z from kalman fit
547 TVectorD vecY1(159); // y from pol1 fit per ROC
548 TVectorD vecZ1(159); // z from pol1 fit per ROC
549 TVectorD vecY2(159); // y from pol2 fit per ROC
550 TVectorD vecZ2(159); // z from pol2 fit per ROC
551 TVectorD vecY4(159); // y from sector fit
552 TVectorD vecZ4(159); // z from sector fit
553 TVectorD vecClY(159); // y cluster position
554 TVectorD vecClZ(159); // z cluster position
555 TVectorD vecSec(159); // sector for each row
557 Double_t chi2I1z=-1; // chi2 of pol1 fit in z (inner)
558 Double_t chi2I1y=-1; // chi2 of pol1 fit in y (inner)
559 Double_t chi2O1z=-1; // chi2 of pol1 fit in z (outer)
560 Double_t chi2O1y=-1; // chi2 of pol1 fit in y (outer)
561 Double_t chi2I2z=-1; // chi2 of pol2 fit in z (inner)
562 Double_t chi2I2y=-1; // chi2 of pol2 fit in y (inner)
563 Double_t chi2O2z=-1; // chi2 of pol2 fit in z (outer)
564 Double_t chi2O2y=-1; // chi2 of pol2 fit in y (outer)
565 Double_t chi2IOz=-1; // chi2 of hyp4 fit in z (inner+outer)
566 Double_t chi2IOy=-1; // chi2 of hyp4 fit in y (inner+outer)
568 Int_t innerSector = -1; // number of inner sector
569 Int_t outerSector = -1; // number of outer sector
570 Int_t nclI=0; // number of clusters (inner)
571 Int_t nclO=0; // number of clusters (outer)
572 //=================================================//
573 // Perform the Kalman Fit using the Tracklet Class //
574 //=================================================//
575 AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
577 AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
578 kFALSE,nclCut,kMaxTracklets);
580 AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
581 AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
582 AliTPCTracklet *tr=0x0;
583 AliTPCTracklet dummy;
584 //continue if we didn't find a tracklet
585 if ( !trInner && !trOuter ) continue;
586 //================================================================================//
587 // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
588 //================================================================================//
589 if ( trInner && trOuter ){
590 if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
591 if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
595 AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
599 if ( !trInner->GetInner() ) continue;
601 if ( trInner->GetSector()>35 ){
606 if ( !trOuter->GetInner() ) continue;
608 if ( trOuter->GetSector()<36 ){
614 innerSector = trInner->GetSector();
615 if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
616 outerSector = trOuter->GetSector();
617 if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
620 TClonesArray arrCl("AliTPCclusterMI",159);
621 arrCl.ExpandCreateFast(159);
622 //=======================================//
623 // fill fitters with cluster information //
624 //=======================================//
625 AliDebug(3,"Filing Cluster Information");
626 for (Int_t irow=158;irow>-1;--irow) {
627 AliTPCclusterMI *c=track->GetClusterPointer(irow);
628 AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
634 Int_t roc = static_cast<Int_t>(c->GetDetector());
635 if ( roc!=innerSector && roc!=outerSector ) continue;
637 //store clusters in clones array
640 vecX[irow] = c->GetX();
641 vecClY[irow] = c->GetY();
642 vecClZ[irow] = c->GetZ();
644 Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
645 Double_t y = vecClY[irow];
646 Double_t z = vecClZ[irow];
648 Double_t x2[2]={x,x*x}; //linear and parabolic parameters
649 Double_t x4[4]={0,0,0,0}; //hyp4 parameters
650 if ( roc == innerSector ) {
651 x4[0]=1; //offset inner - outer sector
652 x4[1]=x; //slope parameter inner sector
654 x4[2]=x; //slope parameter outer sector
656 x4[3]=x*x; //common parabolic shape
658 if ( roc==innerSector ){
665 if ( roc==outerSector ){
675 //======================================//
676 // Evaluate and retrieve fit parameters //
677 //======================================//
678 AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
680 if ( (innerSector>-1) && (fy1I.GetNpoints()>0) ){
685 fy1I.GetParameters(vecy1resInner);
686 fz1I.GetParameters(vecz1resInner);
687 fy2I.GetParameters(vecy2resInner);
688 fz2I.GetParameters(vecz2resInner);
689 chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
690 chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
691 chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
692 chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
695 if ( (outerSector>-1) && (fy1O.GetNpoints()>0) ){
700 fy1O.GetParameters(vecy1resOuter);
701 fz1O.GetParameters(vecz1resOuter);
702 fy2O.GetParameters(vecy2resOuter);
703 fz2O.GetParameters(vecz2resOuter);
704 chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
705 chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
706 chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
707 chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
710 if ( innerSector>0 && outerSector>0 ){
711 if (fy4.GetNpoints()>0) {
713 fy4.GetParameters(vecy4res);
714 chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
716 if (fz4.GetNpoints()>0) {
718 fz4.GetParameters(vecz4res);
719 chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
723 fy4.ClearPoints(); fz4.ClearPoints();
724 fy1I.ClearPoints(); fy1O.ClearPoints();
725 fz1I.ClearPoints(); fz1O.ClearPoints();
726 fy2I.ClearPoints(); fy2O.ClearPoints();
727 fz2I.ClearPoints(); fz2O.ClearPoints();
728 //==============================//
729 // calculate tracklet positions //
730 //==============================//
731 AliDebug(4,"Calculate tracklet positions");
732 for (Int_t irow=158;irow>-1;--irow) {
733 if ( vecSec[irow]==-1 ) continue; //no cluster info
734 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
736 Double_t x=vecX[irow];
737 Double_t xref=x-133.4;
739 Double_t yoffInner=0;
740 Double_t zoffInner=0;
741 Double_t yslopeInner=0;
742 Double_t yslopeOuter=0;
743 Double_t zslopeInner=0;
744 Double_t zslopeOuter=0;
745 //positions of hyperplane fits
746 if ( vecSec[irow] == outerSector ) {
748 vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
749 vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
750 vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
751 vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
752 yslopeOuter=vecy4res[3];
753 zslopeOuter=vecz4res[3];
756 vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
757 vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
758 vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
759 vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
760 yoffInner=vecy4res[1];
761 zoffInner=vecz4res[1];
762 yslopeInner=vecy4res[2];
763 zslopeInner=vecz4res[2];
765 vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
766 vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
767 //positions of kalman fits
768 Double_t gxyz[3],xyz[3];
769 AliExternalTrackParam *param = 0x0;
771 param=tr->GetInner();
774 Float_t bz = AliTracker::GetBz(gxyz);
775 param->GetYAt(x, bz, xyz[1]);
776 param->GetZAt(x, bz, xyz[2]);
777 vecYkalman[irow]=xyz[1];
778 vecZkalman[irow]=xyz[2];
781 //=====================================================================//
782 // write results from the different tracklet fits with debug streamers //
783 //=====================================================================//
785 TTreeSRedirector *cstream = GetDebugStreamer();
787 Float_t dedx = track->GetdEdx();
788 (*cstream)<<"FitModels"<<
790 "edgeCutX=" << edgeCutX <<
791 "edgeCutY=" << edgeCutY <<
792 "nclCut=" << nclCut <<
793 "innerSector="<< innerSector <<
794 "outerSector="<< outerSector <<
797 "Tr.=" << extparam <<
798 "yPol1In.=" << &vecy1resInner <<
799 "zPol1In.=" << &vecz1resInner <<
800 "yPol2In.=" << &vecy2resInner <<
801 "zPol2In.=" << &vecz2resInner <<
802 "yPol1Out.=" << &vecy1resOuter <<
803 "zPol1Out.=" << &vecz1resOuter <<
804 "yPol2Out.=" << &vecy2resOuter <<
805 "zPol2Out.=" << &vecz2resOuter <<
806 "yInOut.=" << &vecy4res <<
807 "zInOut.=" << &vecz4res <<
808 "chi2y1In=" << chi2I1y <<
809 "chi2z1In=" << chi2I1z <<
810 "chi2y1Out=" << chi2O1y <<
811 "chi2z1Out=" << chi2O1z <<
812 "chi2y2In=" << chi2I2y <<
813 "chi2z2In=" << chi2I2z <<
814 "chi2y2Out=" << chi2O2y <<
815 "chi2z2Out=" << chi2O2z <<
816 "chi2yInOut=" << chi2IOy <<
817 "chi2zInOut=" << chi2IOz <<
818 "trletIn.=" << trInner <<
819 "trletOut.=" << trOuter <<
826 // wirte residuals information
828 TTreeSRedirector *cstream = GetDebugStreamer();
830 Float_t dedx = track->GetdEdx();
831 (*cstream)<<"Residuals"<<
833 "edgeCutX=" << edgeCutX <<
834 "edgeCutY=" << edgeCutY <<
835 "nclCut=" << nclCut <<
841 "TrYpol1.=" << &vecY1 <<
842 "TrZpol1.=" << &vecZ1 <<
843 "TrYpol2.=" << &vecY2 <<
844 "TrZpol2.=" << &vecZ2 <<
845 "TrYInOut.=" << &vecY4 <<
846 "TrZInOut.=" << &vecZ4 <<
847 "ClY.=" << &vecClY <<
848 "ClZ.=" << &vecClZ <<
849 "sec.=" << &vecSec <<
852 "yInOut.=" << &vecy4res <<
853 "zInOut.=" << &vecz4res <<
854 "chi2y1In=" << chi2I1y <<
855 "chi2z1In=" << chi2I1z <<
856 "chi2y1Out=" << chi2O1y <<
857 "chi2z1Out=" << chi2O1z <<
858 "chi2y2In=" << chi2I2y <<
859 "chi2z2In=" << chi2I2z <<
860 "chi2y2Out=" << chi2O2y <<
861 "chi2z2Out=" << chi2O2z <<
862 "chi2yInOut=" << chi2IOy <<
863 "chi2zInOut=" << chi2IOz <<
871 Int_t indexMaxCut[kMaxTracklets];
873 Float_t xMinCut[kMaxTracklets];
874 Float_t xMedCut[kMaxTracklets];
875 Float_t xMaxCut[kMaxTracklets];
877 for (Int_t i=0; i<kMaxTracklets; i++){
878 trMinCut = (AliTPCTracklet*)trackletsMinCuts->At(i);
879 trMedCut = (AliTPCTracklet*)trackletsMedCuts->At(i);
880 trMaxCut = (AliTPCTracklet*)trackletsMaxCuts->At(i);
881 if (!trMinCut ) trMinCut=&dummy;
882 if (!trMedCut ) trMedCut=&dummy;
883 if (!trMaxCut ) trMaxCut=&dummy;
884 xMinCut[i]=trMinCut->GetInner()->GetX();
885 xMedCut[i]=trMedCut->GetInner()->GetX();
886 xMaxCut[i]=trMaxCut->GetInner()->GetX();
888 TMath::Sort(kMaxTracklets, xMinCut, indexMinCut);
889 TMath::Sort(kMaxTracklets, xMedCut, indexMedCut);
890 TMath::Sort(kMaxTracklets, xMaxCut, indexMaxCut);
896 void AliTPCcalibLaser::RefitLaser(Int_t id){
898 // Refit the track store residuals
901 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
902 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
903 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
905 //linear fit model in y and z per sector
906 static TLinearFitter fy1(2,"hyp1");
907 static TLinearFitter fz1(2,"hyp1");
908 //quadratic fit model in y and z per sector
909 static TLinearFitter fy2(3,"hyp2");
910 static TLinearFitter fz2(3,"hyp2");
911 static TVectorD vecy2,vecz2,vecy1,vecz1;
913 const Int_t kMinClusters=20;
916 for (Int_t i=0;i<72;++i) nclusters[i]=0;
918 for (Int_t i=0;i<160;++i) {
919 AliTPCclusterMI *c=track->GetClusterPointer(i);
920 if (c) nclusters[c->GetDetector()]++;
923 for (Int_t isec=0; isec<72;isec++){
924 if (nclusters[isec]<kMinClusters) continue;
930 for (Int_t irow=0;irow<160;++irow) {
931 AliTPCclusterMI *c=track->GetClusterPointer(irow);
932 //if (c && RejectCluster(c)) continue;
933 Double_t xd = c->GetX()-133.4; // reference x is beteen iroc and oroc
934 if (c&&c->GetDetector()==isec) {
935 Double_t x[2]={xd,xd*xd};
936 fy2.AddPoint(x,c->GetY());
937 fz2.AddPoint(x,c->GetZ());
939 fy1.AddPoint(x,c->GetY());
940 fz1.AddPoint(x,c->GetZ());
947 fy1.GetParameters(vecy1);
948 fy2.GetParameters(vecy2);
949 fz1.GetParameters(vecz1);
950 fz2.GetParameters(vecz2);
953 TTreeSRedirector *cstream = GetDebugStreamer();
955 Float_t dedx = track->GetdEdx();
956 (*cstream)<<"Tracklet"<<
960 "ncl="<<nclusters[isec]<<
974 // for (Int_t irow=0;irow<160;++irow) {
975 // AliTPCclusterMI *c=track->GetClusterPointer(irow);
976 // if (c && RejectCluster(c)) continue;
977 // if (c&&c->GetDetector()==isec) {
978 // Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
979 // fy2.AddPoint(&x,c->GetY());
980 // fz2.AddPoint(&x,c->GetZ());
982 // fy1.AddPoint(&x,c->GetY());
983 // fz1.AddPoint(&x,c->GetZ());
990 void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield,Int_t minEntries){
992 // Dump information about laser beams
993 // isOK variable indicates usability of the beam
994 // Beam is not usable if:
995 // a. No entries in range (krmsCut0)
996 // b. Big sperad (krmscut1)
997 // c. RMSto fit sigma bigger then (kmultiCut)
998 // d. Too big angular spread
1001 const Float_t krmsCut0=0.001;
1002 const Float_t krmsCut1=0.16;
1003 const Float_t kmultiCut=2;
1004 const Float_t kcutP0=0.002;
1006 AliTPCcalibLaser *laser = this;
1007 TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
1008 TF1 fg("fg","gaus");
1011 for (Int_t id=0; id<336; id++){
1013 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
1014 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
1015 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
1016 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
1017 if (!hisphi) continue;;
1018 Double_t entries = hisphi->GetEntries();
1019 if (entries<minEntries) continue;
1021 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1023 AliTPCLaserTrack::LoadTracks();
1024 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1026 Float_t meanphi = hisphi->GetMean();
1027 Float_t rmsphi = hisphi->GetRMS();
1029 Float_t meanphiP = hisphiP->GetMean();
1030 Float_t rmsphiP = hisphiP->GetRMS();
1031 Float_t meanZ = hisZ->GetMean();
1032 Float_t rmsZ = hisZ->GetRMS();
1033 hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
1034 Double_t gphi1 = fg.GetParameter(1);
1035 Double_t gphi2 = fg.GetParameter(2);
1036 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
1037 Double_t gphiP1 = fg.GetParameter(1);
1038 Double_t gphiP2 = fg.GetParameter(2);
1039 hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
1040 Double_t gz1 = fg.GetParameter(1);
1041 Double_t gz2 = fg.GetParameter(2);
1043 Float_t meanS=hisS->GetMean();
1048 ltrp->GetPxPyPz(lpxyz);
1050 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
1051 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
1052 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
1053 if (gphiP2>kcutP0) isOK=kFALSE;
1055 (*pcstream)<<"Mean"<<
1057 "entries="<<entries<< // number of entries
1058 "bz="<<bfield<< // bfield
1059 "LTr.="<<ltrp<< // refernece track
1061 "lx0="<<lxyz[0]<< // reference x
1062 "lx1="<<lxyz[1]<< // reference y
1063 "lx2="<<lxyz[2]<< // refernece z
1064 "lpx0="<<lpxyz[0]<< // reference x
1065 "lpx1="<<lpxyz[1]<< // reference y
1066 "lpx2="<<lpxyz[2]<< // refernece z
1070 "mphi="<<meanphi<< //
1071 "rmsphi="<<rmsphi<< //
1075 "mphiP="<<meanphiP<< //
1076 "rmsphiP="<<rmsphiP<< //
1092 void AliTPCcalibLaser::DumpScanInfo(TTree * chain){
1096 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
1097 TFile * f = pcstream->GetFile();
1099 f->mkdir("dirphiP");
1101 TF1 fp("p1","pol1");
1106 char grnamefull[1000];
1109 Double_t mphiP[100];
1110 Double_t smphi[100];
1111 Double_t smphiP[100];
1121 for (Int_t id=0; id<336; id++){
1123 sprintf(cut,"isOK&&fId==%d",id);
1124 Int_t entries = chain->Draw("bz",cut,"goff");
1125 if (entries<3) continue;
1126 AliTPCLaserTrack *ltrp = 0;;
1127 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
1128 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1132 ltrp->GetPxPyPz(lpxyz);
1134 chain->Draw("bz",cut,"goff");
1135 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
1136 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
1137 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
1139 chain->Draw("gphi1",cut,"goff");
1140 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
1141 chain->Draw("0.05*abs(mphi)+gphi2",cut,"goff");
1142 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
1144 chain->Draw("gphiP1",cut,"goff");
1145 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
1146 chain->Draw("0.05*abs(mphiP)+gphiP2",cut,"goff");
1147 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
1149 chain->Draw("gz1",cut,"goff");
1150 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
1151 chain->Draw("0.01*abs(meanZ)+gz2",cut,"goff");
1152 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
1155 sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
1156 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
1160 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
1163 pphi[0] = fp.GetParameter(0); // offset
1164 pphi[1] = fp.GetParameter(1); // slope
1165 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1166 sprintf(grname,"phi_id%d",id);
1167 grphi->SetName(grname); grphi->SetTitle(grnamefull);
1168 grphi->GetXaxis()->SetTitle("b_{z} (T)");
1169 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
1170 grphi->SetMaximum(1.2);
1171 grphi->SetMinimum(-1.2);
1175 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
1178 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
1181 pphiP[0] = fp.GetParameter(0); // offset
1182 pphiP[1] = fp.GetParameter(1); // slope
1183 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1184 sprintf(grname,"phiP_id%d",id);
1185 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
1186 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
1187 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
1188 grphiP->SetMaximum(pphiP[0]+0.005);
1189 grphiP->SetMinimum(pphiP[0]-0.005);
1191 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
1196 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
1199 pmZ[0] = fp.GetParameter(0); // offset
1200 pmZ[1] = fp.GetParameter(1); // slope
1201 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
1202 sprintf(grname,"mZ_id%d",id);
1203 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
1204 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
1205 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
1207 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
1211 for (Int_t ientry=0; ientry<entries; ientry++){
1212 (*pcstream)<<"Mean"<<
1215 "entries="<<entries<<
1217 "lx0="<<lxyz[0]<< // reference x
1218 "lx1="<<lxyz[1]<< // reference y
1219 "lx2="<<lxyz[2]<< // refernece z
1220 "lpx0="<<lpxyz[0]<< // reference x
1221 "lpx1="<<lpxyz[1]<< // reference y
1222 "lpx2="<<lpxyz[2]<< // refernece z
1224 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
1225 "pphi0="<<pphi[0]<< // offset
1226 "pphi1="<<pphi[1]<< // mean
1227 "pphi2="<<pphi[2]<< // norm chi2
1229 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
1230 "pphiP0="<<pphiP[0]<< // offset
1231 "pphiP1="<<pphiP[1]<< // mean
1232 "pphiP2="<<pphiP[2]<< // norm chi2
1234 "gz1="<<mZ[ientry]<<
1235 "pmZ0="<<pmZ[0]<< // offset
1236 "pmZ1="<<pmZ[1]<< // mean
1237 "pmZ2="<<pmZ[2]<< // norm chi2
1247 void AliTPCcalibLaser::Analyze(){
1254 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
1256 TIterator* iter = li->MakeIterator();
1257 AliTPCcalibLaser* cal = 0;
1259 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
1260 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
1261 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1265 // fHistNTracks->Add(cal->fHistNTracks);
1266 // fClusters->Add(cal->fClusters);
1267 // fModules->Add(cal->fModules);
1268 // fHistPt->Add(cal->fHistPt);
1269 // fPtResolution->Add(cal->fPtResolution);
1270 // fDeDx->Add(cal->fDeDx);
1276 for (Int_t id=0; id<336; id++){
1277 // merge fDeltaZ histograms
1278 hm = (TH1F*)cal->fDeltaZ.At(id);
1279 h = (TH1F*)fDeltaZ.At(id);
1281 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
1282 fDeltaZ.AddAt(h,id);
1285 // merge fDeltaPhi histograms
1286 hm = (TH1F*)cal->fDeltaPhi.At(id);
1287 h = (TH1F*)fDeltaPhi.At(id);
1289 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
1290 fDeltaPhi.AddAt(h,id);
1293 // merge fDeltaPhiP histograms
1294 hm = (TH1F*)cal->fDeltaPhiP.At(id);
1295 h = (TH1F*)fDeltaPhiP.At(id);
1297 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
1298 fDeltaPhiP.AddAt(h,id);
1301 // merge fSignals histograms
1302 hm = (TH1F*)cal->fSignals.At(id);
1303 h = (TH1F*)fSignals.At(id);
1305 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
1306 fSignals.AddAt(h,id);
1318 gSystem->Load("libSTAT.so")
1319 TStatToolkit toolkit;
1324 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
1329 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
1330 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
1331 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
1332 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
1334 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
1335 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
1336 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
1337 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
1339 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
1340 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
1341 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
1342 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
1347 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
1349 treeT->SetAlias("fit",strq0->Data());
1352 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
1354 treeT->SetAlias("fitP",strqP->Data());
1357 TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
1358 treeT->SetAlias("fitD",strqDrift->Data());
1361 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
1363 for (Int_t i=0; i<6;i++){
1364 treeT->SetLineColor(i+2);
1365 treeT->SetMarkerSize(1);
1366 treeT->SetMarkerStyle(22+i);
1367 treeT->SetMarkerColor(i+2);
1369 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
1377 TTree * tree = (TTree*)f.Get("FitModels");
1379 TEventList listLFit0("listLFit0","listLFit0");
1380 TEventList listLFit1("listLFit1","listLFit1");
1382 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
1383 tree->SetEventList(&listLFit0);