TPC/AliTPCcalibLaser.cxx

   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 /*
  17   laser track clasification;
  18   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
  19   TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
  20   TCut cutN("cutN","fTPCncls>70");
  21   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
  22   TCut cutA = cutT+cutPt+cutP;
  23   TFile f("laserTPCDebug.root");
  24   TTree * treeT = (TTree*)f.Get("Track");
  25
  26
  27   treeT->Draw("(atan2(x1,x0)-atan2(lx1,lx0))*250.:fBundle","fSide==1&&fRod==0"+cutA,"prof")
  28
  29   gSystem->Load("libSTAT.so")
  30   TStatToolkit toolkit;
  31   Double_t chi2;
  32   TVectorD fitParam;
  33   TMatrixD covMatrix;
  34   Int_t npoints;
  35   TString *strq0 = toolkit.FitPlane(treeT,"Tr.fP[1]-LTr.fP[1]","lx1++lx2", "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
  36
  37 */
  38
  39
  40
  41 #include "TLinearFitter.h"
  42 #include "AliTPCcalibLaser.h"
  43 #include "AliExternalTrackParam.h"
  44 #include "AliESDEvent.h"
  45 #include "AliESDfriend.h"
  46 #include "AliESDtrack.h"
  47 #include "AliTPCTracklet.h"
  48 #include "TH1D.h"
  49 #include "TVectorD.h"
  50 #include "TTreeStream.h"
  51 #include "TFile.h"
  52 #include "TF1.h"
  53 #include "TGraphErrors.h"
  54 #include "AliTPCclusterMI.h"
  55 #include "AliTPCseed.h"
  56 #include "AliTracker.h"
  57 #include "TClonesArray.h"
  58
  59
  60 #include "TTreeStream.h"
  61 #include <iostream>
  62 #include <sstream>
  63 #include "AliTPCLaserTrack.h"
  64
  65 using namespace std;
  66
  67 ClassImp(AliTPCcalibLaser)
  68
  69 AliTPCcalibLaser::AliTPCcalibLaser():
  70   AliTPCcalibBase(),
  71   fESD(0),
  72   fESDfriend(),
  73   fTracksMirror(336),
  74   fTracksEsd(336),
  75   fTracksEsdParam(336),
  76   fTracksTPC(336),
  77   fDeltaZ(336),          // array of histograms of delta z for each track
  78   fDeltaPhi(336),          // array of histograms of delta z for each track
  79   fDeltaPhiP(336),          // array of histograms of delta z for each track
  80   fFitAside(new TVectorD(3)),        // drift fit - A side
  81   fFitCside(new TVectorD(3)),        // drift fit - C- side
  82   fRun(0)
  83 {
  84   //
  85   // Constructor
  86   //
  87   fTracksEsdParam.SetOwner(kTRUE);
  88 }
  89
  90 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
  91   AliTPCcalibBase(),
  92   fESD(0),
  93   fESDfriend(0),
  94   fTracksMirror(336),
  95   fTracksEsd(336),
  96   fTracksEsdParam(336),
  97   fTracksTPC(336),
  98   fDeltaZ(336),          // array of histograms of delta z for each track
  99   fDeltaPhi(336),          // array of histograms of delta z for each track
 100   fDeltaPhiP(336),          // array of histograms of delta z for each track
 101   fFitAside(new TVectorD(3)),        // drift fit - A side
 102   fFitCside(new TVectorD(3)),        // drift fit - C- side
 103   fRun(0)
 104 {
 105   SetName(name);
 106   SetTitle(title);
 107   //
 108   // Constructor
 109   //
 110   fTracksEsdParam.SetOwner(kTRUE);
 111 }
 112
 113 AliTPCcalibLaser::~AliTPCcalibLaser() {
 114   //
 115   // destructor
 116   //
 117 }
 118
 119
 120
 121 void AliTPCcalibLaser::Process(AliESDEvent * event) {
 122   //
 123   //
 124   // Loop over tracks and call  Process function
 125   //
 126   fESD = event;
 127   if (!fESD) {
 128     return;
 129   }
 130   fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 131   if (!fESDfriend) {
 132     return;
 133   }
 134   fTracksTPC.Clear();
 135   fTracksEsd.Clear();
 136   fTracksEsdParam.Delete();
 137   //
 138   Int_t n=fESD->GetNumberOfTracks();
 139   Int_t run = fESD->GetRunNumber();
 140   fRun = run;
 141   for (Int_t i=0;i<n;++i) {
 142     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
 143     AliESDtrack *track=fESD->GetTrack(i);
 144     TObject *calibObject=0;
 145     AliTPCseed *seed=0;
 146     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
 147       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
 148         break;
 149     if (track&&seed) FindMirror(track,seed);
 150     //
 151   }
 152
 153   FitDriftV();
 154   MakeDistHisto();
 155   //
 156   for (Int_t id=0; id<336; id++){
 157     //
 158     //
 159     if (!fTracksEsdParam.At(id)) continue;
 160     DumpLaser(id);
 161     RefitLaser(id);
 162
 163   }
 164 }
 165
 166 void AliTPCcalibLaser::MakeDistHisto(){
 167   //
 168   //
 169   //
 170   for (Int_t id=0; id<336; id++){
 171     //
 172     //
 173     if (!fTracksEsdParam.At(id)) continue;
 174     if (!AcceptLaser(id)) continue;
 175     //
 176     //
 177     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
 178     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
 179     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
 180     if (!hisdz){
 181       hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),500,-10,10);
 182       fDeltaZ.AddAt(hisdz,id);
 183       //
 184       hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),500,-2,2);
 185       fDeltaPhi.AddAt(hisdphi,id);
 186       //
 187       hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),500,-0.02,0.02);
 188       fDeltaPhiP.AddAt(hisdphiP,id);
 189     }
 190
 191     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 192     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 193     Double_t xyz[3];
 194     Double_t pxyz[3];
 195     Double_t lxyz[3];
 196     Double_t lpxyz[3];
 197     param->GetXYZ(xyz);
 198     param->GetPxPyPz(pxyz);
 199     ltrp->GetXYZ(lxyz);
 200     ltrp->GetPxPyPz(lpxyz);
 201     //
 202     Float_t dz   = param->GetZ()-ltrp->GetZ();
 203     Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
 204     Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
 205     hisdz->Fill(dz);
 206     hisdphi->Fill(dphi);
 207     hisdphiP->Fill(dphiP);
 208   }
 209 }
 210
 211 void AliTPCcalibLaser::FitDriftV(){
 212   //
 213   // Fit drift velocity - linear approximation in the z and global y
 214   //
 215   static TLinearFitter fdriftA(3,"hyp2");
 216   static TLinearFitter fdriftC(3,"hyp2");
 217   fdriftA.ClearPoints();
 218   fdriftC.ClearPoints();
 219   //
 220   for (Int_t id=0; id<336; id++){
 221     if (!fTracksEsdParam.At(id)) continue;
 222     if (!AcceptLaser(id)) continue;
 223     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 224     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 225     Double_t xyz[3];
 226     Double_t pxyz[3];
 227     Double_t lxyz[3];
 228     Double_t lpxyz[3];
 229     param->GetXYZ(xyz);
 230     param->GetPxPyPz(pxyz);
 231     ltrp->GetXYZ(lxyz);
 232     ltrp->GetPxPyPz(lpxyz);
 233     Double_t xxx[2] = {lxyz[2],lxyz[1]};
 234     if (ltrp->GetSide()==0){
 235       fdriftA.AddPoint(xxx,xyz[2],1);
 236     }else{
 237       fdriftC.AddPoint(xxx,xyz[2],1);
 238     }
 239   }
 240   Float_t chi2A = 0;
 241   Float_t chi2C = 0;
 242   Int_t npointsA=0;
 243   Int_t npointsC=0;
 244   //
 245   if (fdriftA.GetNpoints()>10){
 246     fdriftA.Eval();
 247     fdriftA.EvalRobust(0.8);
 248     fdriftA.GetParameters(*fFitAside);
 249     npointsA= fdriftA.GetNpoints();
 250     chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
 251   }
 252   if (fdriftC.GetNpoints()>10){
 253     fdriftC.Eval();
 254     fdriftC.EvalRobust(0.8);
 255     fdriftC.GetParameters(*fFitCside);
 256     npointsC= fdriftC.GetNpoints();
 257     chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
 258   }
 259
 260   if (fStreamLevel>0){
 261     TTreeSRedirector *cstream = GetDebugStreamer();
 262     Int_t time = fESD->GetTimeStamp();
 263     if (cstream){
 264       (*cstream)<<"driftv"<<
 265         "driftA.="<<fFitAside<<
 266         "driftC.="<<fFitCside<<
 267         "chi2A="<<chi2A<<
 268         "chi2C="<<chi2C<<
 269         "nA="<<npointsA<<
 270         "nC="<<npointsC<<
 271         "time="<<time<<
 272         "\n";
 273     }
 274   }
 275   //
 276 }
 277
 278
 279 Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
 280   //
 281   //
 282   //
 283   /*
 284   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
 285   TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
 286   TCut cutN("cutN","fTPCncls>70");
 287   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
 288   TCut cutA = cutT+cutPt+cutP;
 289   */
 290   AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
 291   AliTPCLaserTrack *ltrp  = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 292   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 293
 294   if (TMath::Abs(param->GetParameter()[4])>0.03) return kFALSE;
 295   if (TMath::Abs(param->GetParameter()[3])>0.06) return kFALSE;
 296   if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.06) return kFALSE;
 297   if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>10) return kFALSE;
 298   //
 299   // dedx cut
 300   //
 301   if (TMath::Abs(track->GetTPCsignal())<20) return kFALSE;
 302   if (TMath::Abs(track->GetTPCsignal())>800) return kFALSE;
 303   //
 304   return kTRUE;
 305 }
 306
 307 Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
 308   //
 309   // Find corresponding mirror
 310   // add the corresponding tracks
 311   //
 312   Float_t kRadius0 = 252;
 313   Float_t kRadius  = 253.4;
 314   if (!track->GetOuterParam()) return -1;
 315   AliExternalTrackParam param(*(track->GetOuterParam()));
 316   AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
 317   AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
 318   AliTPCLaserTrack ltr;
 319   AliTPCLaserTrack *ltrp=0x0;
 320   //
 321   Int_t id = AliTPCLaserTrack::IdentifyTrack(&param);
 322   if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
 323     ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
 324   else
 325     ltrp=&ltr;
 326
 327   if (id>=0){
 328     //
 329     //
 330     Float_t radius=TMath::Abs(ltrp->GetX());
 331     AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
 332     //
 333     if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
 334     fTracksEsdParam.AddAt(param.Clone(),id);
 335     fTracksEsd.AddAt(track,id);
 336     fTracksTPC.AddAt(seed,id);
 337     //
 338   }
 339
 340   return id;
 341 }
 342
 343
 344
 345 void AliTPCcalibLaser::DumpLaser(Int_t id) {
 346   //
 347   //  Dump Laser info to the tree
 348   //
 349   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 350   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 351   AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 352   //
 353   // Fast laser ID
 354   //
 355   Double_t xyz[3];
 356   Double_t pxyz[3];
 357   Double_t lxyz[3];
 358   Double_t lpxyz[3];
 359   param->GetXYZ(xyz);
 360   param->GetPxPyPz(pxyz);
 361   ltrp->GetXYZ(lxyz);
 362   ltrp->GetPxPyPz(lpxyz);
 363
 364   if (fStreamLevel>0){
 365     TTreeSRedirector *cstream = GetDebugStreamer();
 366     Int_t time = fESD->GetTimeStamp();
 367     Bool_t accept = AcceptLaser(id);
 368     if (cstream){
 369       (*cstream)<<"Track"<<
 370         "run="<<fRun<<
 371         "id="<<id<<
 372         "accept="<<accept<<
 373         "driftA.="<<fFitAside<<
 374         "driftC.="<<fFitCside<<
 375         "time="<<time<<
 376         //
 377         "LTr.="<<ltrp<<
 378         "Esd.="<<track<<
 379         "Tr.="<<param<<
 380         "x0="<<xyz[0]<<
 381         "x1="<<xyz[1]<<
 382         "x2="<<xyz[2]<<
 383         "px0="<<pxyz[0]<<
 384         "px1="<<pxyz[1]<<
 385         "px2="<<pxyz[2]<<
 386         //
 387         "lx0="<<lxyz[0]<<
 388         "lx1="<<lxyz[1]<<
 389         "lx2="<<lxyz[2]<<
 390         "lpx0="<<lpxyz[0]<<
 391         "lpx1="<<lpxyz[1]<<
 392         "lpx2="<<lpxyz[2]<<
 393         "\n";
 394     }
 395   }
 396 }
 397
 398
 399 void AliTPCcalibLaser::RefitLaser(Int_t id){
 400   //
 401   // Refit the track store residuals
 402   //
 403
 404   AliTPCseed *track    = (AliTPCseed*)fTracksTPC.At(id);
 405   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 406   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
 407
 408   //
 409   static TLinearFitter fy2(3,"hyp2");
 410   static TLinearFitter fz2(3,"hyp2");
 411   static TLinearFitter fy1(2,"hyp1");
 412   static TLinearFitter fz1(2,"hyp1");
 413   static TVectorD vecy2,vecz2,vecy1,vecz1;
 414
 415   const Int_t kMinClusters=20;
 416   Int_t nclusters[72];
 417   //
 418   for (Int_t i=0;i<72;++i) nclusters[i]=0;
 419
 420   for (Int_t i=0;i<160;++i) {
 421     AliTPCclusterMI *c=track->GetClusterPointer(i);
 422     if (c) nclusters[c->GetDetector()]++;
 423   }
 424
 425   for (Int_t isec=0; isec<72;isec++){
 426     if (nclusters[isec]<kMinClusters) continue;
 427     fy2.ClearPoints();
 428     fz2.ClearPoints();
 429     fy1.ClearPoints();
 430     fz1.ClearPoints();
 431     //
 432     for (Int_t irow=0;irow<160;++irow) {
 433       AliTPCclusterMI *c=track->GetClusterPointer(irow);
 434       //if (c && RejectCluster(c)) continue;
 435       if (c&&c->GetDetector()==isec) {
 436         Double_t xd = c->GetX()-120;;
 437         Double_t x[2]={xd,xd*xd};
 438         fy2.AddPoint(x,c->GetY());
 439         fz2.AddPoint(x,c->GetZ());
 440         //
 441         fy1.AddPoint(x,c->GetY());
 442         fz1.AddPoint(x,c->GetZ());
 443       }
 444     }
 445     fy2.Eval();
 446     fz2.Eval();
 447     fy1.Eval();
 448     fz1.Eval();
 449     fy1.GetParameters(vecy1);
 450     fy2.GetParameters(vecy2);
 451     fz1.GetParameters(vecz1);
 452     fz2.GetParameters(vecz2);
 453
 454     if (fStreamLevel>0){
 455       TTreeSRedirector *cstream = GetDebugStreamer();
 456       if (cstream){
 457         Float_t dedx = track->GetdEdx();
 458         (*cstream)<<"Tracklet"<<
 459           "LTr.="<<ltrp<<
 460           "Tr.="<<param<<
 461           "sec="<<isec<<
 462           "ncl="<<nclusters[isec]<<
 463           "dedx="<<dedx<<
 464           "dedx="<<dedx<<
 465           "vy1.="<<&vecy1<<
 466           "vy2.="<<&vecy2<<
 467           "vz1.="<<&vecz1<<
 468           "vz2.="<<&vecz2<<
 469           "\n";
 470       }
 471     }
 472   }
 473   //
 474   //
 475   //
 476   //   for (Int_t irow=0;irow<160;++irow) {
 477   //       AliTPCclusterMI *c=track->GetClusterPointer(irow);
 478   //       if (c && RejectCluster(c)) continue;
 479   //       if (c&&c->GetDetector()==isec) {
 480   //    Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
 481   //    fy2.AddPoint(&x,c->GetY());
 482   //    fz2.AddPoint(&x,c->GetZ());
 483   //    //
 484   //    fy1.AddPoint(&x,c->GetY());
 485   //    fz1.AddPoint(&x,c->GetZ());
 486   //       }
 487   //     }
 488
 489 }
 490
 491
 492 void AliTPCcalibLaser::Analyze(){
 493   //
 494   //
 495   //
 496 }
 497
 498
 499
 500