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   //
  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)
  25   //
  26   // 3. The drift velocity and jitter is calculated event by event
  27   //    (see function drift velocity)
  28   //
  29   //
  30   //
  31   gSystem->Load("libANALYSIS");
  32   gSystem->Load("libTPCcalib");
  33   TFile fcalib("CalibObjects.root");
  34   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  35   AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)array->FindObject("laserTPC");
  36   laser->DumpMeanInfo(-0.4)
  37   TFile fmean("laserMean.root")
  38
  39   //  laser track clasification;
  40
  41   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
  42   TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
  43   TCut cutN("cutN","fTPCncls>70");
  44   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
  45   TCut cutA = cutT+cutPt+cutP;
  46
  47   TFile f("laserTPCDebug.root");
  48   TTree * treeT = (TTree*)f.Get("Track");
  49
  50
  51   // LASER scan
  52   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  53   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  54   AliXRDPROOFtoolkit tool;
  55   TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
  56   chain->Lookup();
  57
  58
  59   treeT->Draw("(atan2(x1,x0)-atan2(lx1,lx0))*250.:fBundle","fSide==1&&fRod==0"+cutA,"prof")
  60
  61   gSystem->Load("libSTAT.so")
  62   TStatToolkit toolkit;
  63   Double_t chi2;
  64   TVectorD fitParam;
  65   TMatrixD covMatrix;
  66   Int_t npoints;
  67   TString *strq0 = toolkit.FitPlane(treeT,"Tr.fP[1]-LTr.fP[1]","lx1++lx2", "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
  68
  69 */
  70
  71
  72
  73 #include "TLinearFitter.h"
  74 #include "AliTPCcalibLaser.h"
  75 #include "AliExternalTrackParam.h"
  76 #include "AliESDEvent.h"
  77 #include "AliESDfriend.h"
  78 #include "AliESDtrack.h"
  79 #include "AliTPCTracklet.h"
  80 #include "TH1D.h"
  81 #include "TVectorD.h"
  82 #include "TTreeStream.h"
  83 #include "TFile.h"
  84 #include "TF1.h"
  85 #include "TGraphErrors.h"
  86 #include "AliTPCclusterMI.h"
  87 #include "AliTPCseed.h"
  88 #include "AliTracker.h"
  89 #include "TClonesArray.h"
  90
  91
  92 #include "TTreeStream.h"
  93 #include <iostream>
  94 #include <sstream>
  95 #include "AliTPCLaserTrack.h"
  96
  97 using namespace std;
  98
  99 ClassImp(AliTPCcalibLaser)
 100
 101 AliTPCcalibLaser::AliTPCcalibLaser():
 102   AliTPCcalibBase(),
 103   fESD(0),
 104   fESDfriend(),
 105   fTracksMirror(336),
 106   fTracksEsd(336),
 107   fTracksEsdParam(336),
 108   fTracksTPC(336),
 109   fDeltaZ(336),          // array of histograms of delta z for each track
 110   fDeltaPhi(336),          // array of histograms of delta z for each track
 111   fDeltaPhiP(336),          // array of histograms of delta z for each track
 112   fSignals(336),           // array of dedx signals
 113   fFitAside(new TVectorD(3)),        // drift fit - A side
 114   fFitCside(new TVectorD(3)),        // drift fit - C- side
 115   fRun(0)
 116 {
 117   //
 118   // Constructor
 119   //
 120   fTracksEsdParam.SetOwner(kTRUE);
 121 }
 122
 123 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
 124   AliTPCcalibBase(),
 125   fESD(0),
 126   fESDfriend(0),
 127   fTracksMirror(336),
 128   fTracksEsd(336),
 129   fTracksEsdParam(336),
 130   fTracksTPC(336),
 131   fDeltaZ(336),          // array of histograms of delta z for each track
 132   fDeltaPhi(336),          // array of histograms of delta z for each track
 133   fDeltaPhiP(336),          // array of histograms of delta z for each track
 134   fSignals(336),           // array of dedx signals
 135   fFitAside(new TVectorD(3)),        // drift fit - A side
 136   fFitCside(new TVectorD(3)),        // drift fit - C- side
 137   fRun(0)
 138 {
 139   SetName(name);
 140   SetTitle(title);
 141   //
 142   // Constructor
 143   //
 144   fTracksEsdParam.SetOwner(kTRUE);
 145 }
 146
 147 AliTPCcalibLaser::~AliTPCcalibLaser() {
 148   //
 149   // destructor
 150   //
 151 }
 152
 153
 154
 155 void AliTPCcalibLaser::Process(AliESDEvent * event) {
 156   //
 157   //
 158   // Loop over tracks and call  Process function
 159   //
 160   fESD = event;
 161   if (!fESD) {
 162     return;
 163   }
 164   fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 165   if (!fESDfriend) {
 166     return;
 167   }
 168   fTracksTPC.Clear();
 169   fTracksEsd.Clear();
 170   fTracksEsdParam.Delete();
 171   //
 172   Int_t n=fESD->GetNumberOfTracks();
 173   Int_t run = fESD->GetRunNumber();
 174   fRun = run;
 175   for (Int_t i=0;i<n;++i) {
 176     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
 177     AliESDtrack *track=fESD->GetTrack(i);
 178     TObject *calibObject=0;
 179     AliTPCseed *seed=0;
 180     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
 181       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
 182         break;
 183     if (track&&seed) FindMirror(track,seed);
 184     //
 185   }
 186
 187   FitDriftV();
 188   MakeDistHisto();
 189   //
 190   for (Int_t id=0; id<336; id++){
 191     //
 192     //
 193     if (!fTracksEsdParam.At(id)) continue;
 194     DumpLaser(id);
 195     RefitLaser(id);
 196
 197   }
 198 }
 199
 200 void AliTPCcalibLaser::MakeDistHisto(){
 201   //
 202   //
 203   //
 204   for (Int_t id=0; id<336; id++){
 205     //
 206     //
 207     if (!fTracksEsdParam.At(id)) continue;
 208     if (!AcceptLaser(id)) continue;
 209     //
 210     //
 211     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
 212     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
 213     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
 214     TH1F * hisSignal = (TH1F*)fSignals.At(id);
 215
 216     if (!hisdz){
 217       hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
 218       hisdz->SetDirectory(0);
 219       fDeltaZ.AddAt(hisdz,id);
 220       //
 221       hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
 222       hisdphi->SetDirectory(0);
 223       fDeltaPhi.AddAt(hisdphi,id);
 224       //
 225       hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
 226       hisdphiP->SetDirectory(0);
 227       fDeltaPhiP.AddAt(hisdphiP,id);
 228       hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
 229       hisSignal->SetDirectory(0);
 230       fSignals.AddAt(hisSignal,id);
 231     }
 232
 233     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 234     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 235     AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 236     if (!param) return;
 237     if (!ltrp) return;
 238     if (!track) return;
 239     Double_t xyz[3];
 240     Double_t pxyz[3];
 241     Double_t lxyz[3];
 242     Double_t lpxyz[3];
 243     param->GetXYZ(xyz);
 244     param->GetPxPyPz(pxyz);
 245     ltrp->GetXYZ(lxyz);
 246     ltrp->GetPxPyPz(lpxyz);
 247     //
 248     Float_t dz   = param->GetZ()-ltrp->GetZ();
 249     Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
 250     Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
 251     if (hisdz) hisdz->Fill(dz);
 252     if (hisdphi) hisdphi->Fill(dphi);
 253     if (hisdphiP) hisdphiP->Fill(dphiP);
 254     if (hisSignal) hisSignal->Fill(track->GetTPCsignal());
 255   }
 256 }
 257
 258 void AliTPCcalibLaser::FitDriftV(){
 259   //
 260   // Fit drift velocity - linear approximation in the z and global y
 261   //
 262   static TLinearFitter fdriftA(3,"hyp2");
 263   static TLinearFitter fdriftC(3,"hyp2");
 264   fdriftA.ClearPoints();
 265   fdriftC.ClearPoints();
 266   //
 267   for (Int_t id=0; id<336; id++){
 268     if (!fTracksEsdParam.At(id)) continue;
 269     if (!AcceptLaser(id)) continue;
 270     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 271     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 272     Double_t xyz[3];
 273     Double_t pxyz[3];
 274     Double_t lxyz[3];
 275     Double_t lpxyz[3];
 276     param->GetXYZ(xyz);
 277     param->GetPxPyPz(pxyz);
 278     ltrp->GetXYZ(lxyz);
 279     ltrp->GetPxPyPz(lpxyz);
 280     Double_t xxx[2] = {lxyz[2],lxyz[1]};
 281     if (ltrp->GetSide()==0){
 282       fdriftA.AddPoint(xxx,xyz[2],1);
 283     }else{
 284       fdriftC.AddPoint(xxx,xyz[2],1);
 285     }
 286   }
 287   Float_t chi2A = 0;
 288   Float_t chi2C = 0;
 289   Int_t npointsA=0;
 290   Int_t npointsC=0;
 291   //
 292   if (fdriftA.GetNpoints()>10){
 293     fdriftA.Eval();
 294     fdriftA.EvalRobust(0.8);
 295     fdriftA.GetParameters(*fFitAside);
 296     npointsA= fdriftA.GetNpoints();
 297     chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
 298   }
 299   if (fdriftC.GetNpoints()>10){
 300     fdriftC.Eval();
 301     fdriftC.EvalRobust(0.8);
 302     fdriftC.GetParameters(*fFitCside);
 303     npointsC= fdriftC.GetNpoints();
 304     chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
 305   }
 306
 307   if (fStreamLevel>0){
 308     TTreeSRedirector *cstream = GetDebugStreamer();
 309     Int_t time = fESD->GetTimeStamp();
 310     if (cstream){
 311       (*cstream)<<"driftv"<<
 312         "driftA.="<<fFitAside<<
 313         "driftC.="<<fFitCside<<
 314         "chi2A="<<chi2A<<
 315         "chi2C="<<chi2C<<
 316         "nA="<<npointsA<<
 317         "nC="<<npointsC<<
 318         "time="<<time<<
 319         "\n";
 320     }
 321   }
 322   //
 323 }
 324
 325
 326 Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
 327   //
 328   //
 329   //
 330   /*
 331   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
 332   TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
 333   TCut cutN("cutN","fTPCncls>70");
 334   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
 335   TCut cutA = cutT+cutPt+cutP;
 336   */
 337   AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
 338   AliTPCLaserTrack *ltrp  = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 339   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 340
 341   if (TMath::Abs(param->GetParameter()[4])>0.03) return kFALSE;
 342   if (TMath::Abs(param->GetParameter()[3])>0.06) return kFALSE;
 343   if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.06) return kFALSE;
 344   if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>10) return kFALSE;
 345   //
 346   // dedx cut
 347   //
 348   if (TMath::Abs(track->GetTPCsignal())<20) return kFALSE;
 349   if (TMath::Abs(track->GetTPCsignal())>800) return kFALSE;
 350   //
 351   return kTRUE;
 352 }
 353
 354 Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
 355   //
 356   // Find corresponding mirror
 357   // add the corresponding tracks
 358   //
 359   Float_t kRadius0 = 252;
 360   Float_t kRadius  = 253.4;
 361   if (!track->GetOuterParam()) return -1;
 362   AliExternalTrackParam param(*(track->GetOuterParam()));
 363   AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
 364   AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
 365   AliTPCLaserTrack ltr;
 366   AliTPCLaserTrack *ltrp=0x0;
 367   //
 368   Int_t id = AliTPCLaserTrack::IdentifyTrack(&param);
 369   if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
 370     ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
 371   else
 372     ltrp=&ltr;
 373
 374   if (id>=0){
 375     //
 376     //
 377     Float_t radius=TMath::Abs(ltrp->GetX());
 378     AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
 379     //
 380     if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
 381     fTracksEsdParam.AddAt(param.Clone(),id);
 382     fTracksEsd.AddAt(track,id);
 383     fTracksTPC.AddAt(seed,id);
 384     //
 385   }
 386
 387   return id;
 388 }
 389
 390
 391
 392 void AliTPCcalibLaser::DumpLaser(Int_t id) {
 393   //
 394   //  Dump Laser info to the tree
 395   //
 396   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 397   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 398   AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 399   //
 400   // Fast laser ID
 401   //
 402   Double_t xyz[3];
 403   Double_t pxyz[3];
 404   Double_t lxyz[3];
 405   Double_t lpxyz[3];
 406   param->GetXYZ(xyz);
 407   param->GetPxPyPz(pxyz);
 408   ltrp->GetXYZ(lxyz);
 409   ltrp->GetPxPyPz(lpxyz);
 410
 411   if (fStreamLevel>0){
 412     TTreeSRedirector *cstream = GetDebugStreamer();
 413     Int_t time = fESD->GetTimeStamp();
 414     Bool_t accept = AcceptLaser(id);
 415     if (cstream){
 416       (*cstream)<<"Track"<<
 417         "run="<<fRun<<
 418         "id="<<id<<
 419         "accept="<<accept<<
 420         "driftA.="<<fFitAside<<
 421         "driftC.="<<fFitCside<<
 422         "time="<<time<<
 423         //
 424         "LTr.="<<ltrp<<
 425         "Esd.="<<track<<
 426         "Tr.="<<param<<
 427         "x0="<<xyz[0]<<
 428         "x1="<<xyz[1]<<
 429         "x2="<<xyz[2]<<
 430         "px0="<<pxyz[0]<<
 431         "px1="<<pxyz[1]<<
 432         "px2="<<pxyz[2]<<
 433         //
 434         "lx0="<<lxyz[0]<<
 435         "lx1="<<lxyz[1]<<
 436         "lx2="<<lxyz[2]<<
 437         "lpx0="<<lpxyz[0]<<
 438         "lpx1="<<lpxyz[1]<<
 439         "lpx2="<<lpxyz[2]<<
 440         "\n";
 441     }
 442   }
 443 }
 444
 445
 446 void AliTPCcalibLaser::RefitLaser(Int_t id){
 447   //
 448   // Refit the track store residuals
 449   //
 450
 451   AliTPCseed *track    = (AliTPCseed*)fTracksTPC.At(id);
 452   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 453   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
 454
 455   //
 456   static TLinearFitter fy2(3,"hyp2");
 457   static TLinearFitter fz2(3,"hyp2");
 458   static TLinearFitter fy1(2,"hyp1");
 459   static TLinearFitter fz1(2,"hyp1");
 460   static TVectorD vecy2,vecz2,vecy1,vecz1;
 461
 462   const Int_t kMinClusters=20;
 463   Int_t nclusters[72];
 464   //
 465   for (Int_t i=0;i<72;++i) nclusters[i]=0;
 466
 467   for (Int_t i=0;i<160;++i) {
 468     AliTPCclusterMI *c=track->GetClusterPointer(i);
 469     if (c) nclusters[c->GetDetector()]++;
 470   }
 471
 472   for (Int_t isec=0; isec<72;isec++){
 473     if (nclusters[isec]<kMinClusters) continue;
 474     fy2.ClearPoints();
 475     fz2.ClearPoints();
 476     fy1.ClearPoints();
 477     fz1.ClearPoints();
 478     //
 479     for (Int_t irow=0;irow<160;++irow) {
 480       AliTPCclusterMI *c=track->GetClusterPointer(irow);
 481       //if (c && RejectCluster(c)) continue;
 482       if (c&&c->GetDetector()==isec) {
 483         Double_t xd = c->GetX()-120;;
 484         Double_t x[2]={xd,xd*xd};
 485         fy2.AddPoint(x,c->GetY());
 486         fz2.AddPoint(x,c->GetZ());
 487         //
 488         fy1.AddPoint(x,c->GetY());
 489         fz1.AddPoint(x,c->GetZ());
 490       }
 491     }
 492     fy2.Eval();
 493     fz2.Eval();
 494     fy1.Eval();
 495     fz1.Eval();
 496     fy1.GetParameters(vecy1);
 497     fy2.GetParameters(vecy2);
 498     fz1.GetParameters(vecz1);
 499     fz2.GetParameters(vecz2);
 500
 501     if (fStreamLevel>0){
 502       TTreeSRedirector *cstream = GetDebugStreamer();
 503       if (cstream){
 504         Float_t dedx = track->GetdEdx();
 505         (*cstream)<<"Tracklet"<<
 506           "LTr.="<<ltrp<<
 507           "Tr.="<<param<<
 508           "sec="<<isec<<
 509           "ncl="<<nclusters[isec]<<
 510           "dedx="<<dedx<<
 511           "dedx="<<dedx<<
 512           "vy1.="<<&vecy1<<
 513           "vy2.="<<&vecy2<<
 514           "vz1.="<<&vecz1<<
 515           "vz2.="<<&vecz2<<
 516           "\n";
 517       }
 518     }
 519   }
 520   //
 521   //
 522   //
 523   //   for (Int_t irow=0;irow<160;++irow) {
 524   //       AliTPCclusterMI *c=track->GetClusterPointer(irow);
 525   //       if (c && RejectCluster(c)) continue;
 526   //       if (c&&c->GetDetector()==isec) {
 527   //    Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
 528   //    fy2.AddPoint(&x,c->GetY());
 529   //    fz2.AddPoint(&x,c->GetZ());
 530   //    //
 531   //    fy1.AddPoint(&x,c->GetY());
 532   //    fz1.AddPoint(&x,c->GetZ());
 533   //       }
 534   //     }
 535
 536 }
 537
 538
 539 void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield){
 540   //
 541   //
 542   //
 543   AliTPCcalibLaser *laser = this;
 544   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
 545   TF1 fg("fg","gaus");
 546   //
 547   //
 548   for (Int_t id=0; id<336; id++){
 549     TH1F * hisphi  = (TH1F*)laser->fDeltaPhi.At(id);
 550     TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
 551     TH1F * hisZ    = (TH1F*)laser->fDeltaZ.At(id);
 552     //    TH1F * hisS    = laser->fDeltaZ->At(id);
 553     if (!hisphi) continue;;
 554     Double_t entries = hisphi->GetEntries();
 555     if (entries<30) continue;
 556     //
 557     AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
 558     if (!ltrp) {
 559      AliTPCLaserTrack::LoadTracks();
 560       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
 561     }
 562     Float_t meanphi = hisphi->GetMean();
 563     Float_t rmsphi = hisphi->GetRMS();
 564     Float_t meanphiP = hisphiP->GetMean();
 565     Float_t rmsphiP = hisphiP->GetRMS();
 566     Float_t meanZ = hisZ->GetMean();
 567     Float_t rmsZ = hisZ->GetRMS();
 568     hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
 569     Double_t gphi1 = fg.GetParameter(1);
 570     Double_t gphi2 = fg.GetParameter(2);
 571     hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
 572     Double_t gphiP1 = fg.GetParameter(1);
 573     Double_t gphiP2 = fg.GetParameter(2);
 574     hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
 575     Double_t gz1 = fg.GetParameter(1);
 576     Double_t gz2 = fg.GetParameter(2);
 577
 578     //
 579     Double_t lxyz[3];
 580     Double_t lpxyz[3];
 581     ltrp->GetXYZ(lxyz);
 582     ltrp->GetPxPyPz(lpxyz);
 583     //
 584     (*pcstream)<<"Mean"<<
 585       "entries="<<entries<<      // number of entries
 586       "bz="<<bfield<<            // bfield
 587       "LTr.="<<ltrp<<             // refernece track
 588       "lx0="<<lxyz[0]<<          // reference x
 589       "lx1="<<lxyz[1]<<          // reference y
 590       "lx2="<<lxyz[2]<<          // refernece z
 591       "lpx0="<<lpxyz[0]<<          // reference x
 592       "lpx1="<<lpxyz[1]<<          // reference y
 593       "lpx2="<<lpxyz[2]<<          // refernece z
 594       "mphi="<<meanphi<<         //
 595       "rmsphi="<<rmsphi<<        //
 596       "gphi1="<<gphi1<<
 597       "gphi2="<<gphi2<<
 598       "mphiP="<<meanphiP<<       //
 599       "rmsphiP="<<rmsphiP<<      //
 600       "gphiP1="<<gphiP1<<
 601       "gphiP2="<<gphiP2<<
 602       "meanZ="<<meanZ<<
 603       "rmsZ="<<rmsZ<<
 604       "gz1="<<gz1<<
 605       "gz2="<<gz2<<
 606
 607       "\n";
 608   }
 609   delete pcstream;
 610 }
 611
 612 void AliTPCcalibLaser::Analyze(){
 613   //
 614   //
 615   //
 616 }
 617
 618
 619 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
 620
 621   TIterator* iter = li->MakeIterator();
 622   AliTPCcalibLaser* cal = 0;
 623
 624   while ((cal = (AliTPCcalibLaser*)iter->Next())) {
 625     if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
 626       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 627       return -1;
 628     }
 629     //
 630    //  fHistNTracks->Add(cal->fHistNTracks);
 631 //     fClusters->Add(cal->fClusters);
 632 //     fModules->Add(cal->fModules);
 633 //     fHistPt->Add(cal->fHistPt);
 634 //     fPtResolution->Add(cal->fPtResolution);
 635 //     fDeDx->Add(cal->fDeDx);
 636
 637
 638     TH1F *h=0x0;
 639     TH1F *hm=0x0;
 640
 641     for (Int_t id=0; id<336; id++){
 642       // merge fDeltaZ histograms
 643       hm = (TH1F*)cal->fDeltaZ.At(id);
 644       h  = (TH1F*)fDeltaZ.At(id);
 645       if (!h) {
 646         h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
 647         fDeltaZ.AddAt(h,id);
 648       }
 649       if (hm) h->Add(hm);
 650       // merge fDeltaPhi histograms
 651       hm = (TH1F*)cal->fDeltaPhi.At(id);
 652       h  = (TH1F*)fDeltaPhi.At(id);
 653       if (!h) {
 654         h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
 655         fDeltaPhi.AddAt(h,id);
 656       }
 657       if (hm) h->Add(hm);
 658       // merge fDeltaPhiP histograms
 659       hm = (TH1F*)cal->fDeltaPhiP.At(id);
 660       h  = (TH1F*)fDeltaPhiP.At(id);
 661       if (!h) {
 662         h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
 663         fDeltaPhiP.AddAt(h,id);
 664       }
 665       if (hm) h->Add(hm);
 666       // merge fSignals histograms
 667       hm = (TH1F*)cal->fSignals.At(id);
 668       h  = (TH1F*)fSignals.At(id);
 669       if (!h) {
 670         h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
 671         fSignals.AddAt(h,id);
 672       }
 673       if (hm) h->Add(hm);
 674     }
 675   }
 676   return 0;
 677 }
 678
 679
 680
 681 /*
 682  gSystem->Load("libSTAT.so")
 683  TStatToolkit toolkit;
 684  Double_t chi2;
 685  TVectorD fitParam;
 686  TMatrixD covMatrix;
 687  Int_t npoints;
 688  TCut cutA("gphi2<0.2&&abs(mphi-gphi1)<0.1&&entries>60")
 689
 690 //  TString *strq0 = toolkit.FitPlane(treeT,"Tr.fP[1]-LTr.fP[1]","lx1++lx2", "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
 691
 692
 693 TString fstring="";
 694 fstring+="LTr.fP[2]++"                        // 1
 695 fstring+="LTr.fP[2]*cos(atan2(lx1,lx0))++"    // 2
 696 fstring+="LTr.fP[2]*sin(atan2(lx1,lx0))++"    // 3
 697 fstring+="cos(atan2(lx1,lx0))++"              // 4
 698 fstring+="sin(atan2(lx1,lx0))++"              // 5
 699 fstring+="gphiP1++"
 700
 701 fstring+="bz++";
 702 fstring+="bz*LTr.fP[2]++";                      // 7
 703 fstring+="bz*cos(atan2(lx1,lx0))++"           // 8
 704 fstring+="bz*sin(atan2(lx1,lx0))++"           // 9
 705 fstring+="bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" // 10
 706 fstring+="bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"   // 11
 707
 708 Int_t entries = chain->Draw("fId","fSide==1&&fBundle==3&&bz==0"+cutA)
 709 {
 710 for (Int_t i=2;i<entries;i++){
 711 fstring+=Form("(abs(fId-%f)<0.1)++",chain->GetV1()[i]);
 712 }
 713 }
 714
 715
 716 TString *strq0 = toolkit.FitPlane(chain,"gphi1",fstring->Data(), "fSide==1&&fBundle==3"+cutA, chi2,npoints,fitParam,covMatrix);
 717
 718 chain->SetAlias("fit",strq0->Data());
 719
 720 gSystem->Load("libSTAT.so")
 721 TStatToolkit toolkit;
 722 Double_t chi2;
 723  TVectorD fitParam;
 724  TMatrixD covMatrix;
 725  Int_t npoints;
 726  TCut cutA("gphi2<0.2&&abs(mphi-gphi1)<0.1&&entries>60")
 727
 728 //  TString *strq0 = toolkit.FitPlane(treeT,"Tr.fP[1]-LTr.fP[1]","lx1++lx2", "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
 729
 730
 731 TString fstring="";
 732
 733 fstring+="LTr.fP[2]++"                        // 1
 734 fstring+="LTr.fP[2]*cos(atan2(lx1,lx0))++"    // 2
 735 fstring+="LTr.fP[2]*sin(atan2(lx1,lx0))++"    // 3
 736 fstring+="cos(atan2(lx1,lx0))++"              // 4
 737 fstring+="sin(atan2(lx1,lx0))++"              // 5
 738 fstring+="gphiP1++"
 739 //
 740 //
 741 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++";
 742 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++";                     // 7
 743 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++"            // 8
 744 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++"            // 9
 745 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  // 10
 746 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  // 11
 747 //
 748 fstring+="(abs(LTr.fP[1]/250)^2-1)*bz++";
 749 fstring+="(abs(LTr.fP[1]/250)^2-1)*bz*LTr.fP[2]++";                     // 7
 750 fstring+="(abs(LTr.fP[1]/250)^2-1)*bz*cos(atan2(lx1,lx0))++"            // 8
 751 fstring+="(abs(LTr.fP[1]/250)^2-1)*bz*sin(atan2(lx1,lx0))++"            // 9
 752 fstring+="(abs(LTr.fP[1]/250)^2-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  // 10
 753 fstring+="(abs(LTr.fP[1]/250)^2-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  // 11
 754 //
 755 fstring+="(abs(LTr.fP[1]/250)-1)*bz++";
 756 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++";                     // 7
 757 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++"            // 8
 758 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++"            // 9
 759 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  // 10
 760 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  // 11
 761
 762
 763 Int_t entries = chain->Draw("fId","fSide==1&&bz==0"+cutA)
 764 {
 765 for (Int_t i=2;i<entries;i++){
 766 fstring+=Form("(abs(fId-%f)<0.1)++",chain->GetV1()[i]);
 767 }
 768 }
 769
 770
 771  TString *strq0 = toolkit.FitPlane(chain,"gphi1",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
 772
 773  chain->SetAlias("fit",strq0->Data());
 774
 775
 776  */