TPC/AliTPCcalibTime.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     Comments to be written here:
  18     1. What do we calibrate.
  19
  20 Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
  21
  22   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
  23
  24
  25     2. How to interpret results
  26     3. Simple example
  27
  28 a.) determine the required time range:
  29
  30 AliXRDPROOFtoolkit tool;
  31 TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
  32 chain->Draw("GetTimeStamp()")
  33
  34 b) analyse calibration object on Proof in calibration train
  35
  36 AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift, IntegrationTimeDeDx);
  37 s
  38 c) plot results
  39  .x ~/NimStyle.C
  40   gSystem->Load("libANALYSIS");
  41   gSystem->Load("libTPCcalib");
  42
  43 TFile f("CalibObjects.root");
  44 AliTPCcalibTime *cal = (AliTPCcalibTime *)f->Get("TPCCalib")->FindObject("calibTime");
  45 cal->GetHistVdrift()->Projection(1,0)->Draw()
  46
  47     4. Analysis using debug streamers.
  48
  49     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  50     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  51     AliXRDPROOFtoolkit tool;
  52     TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
  53     chainTime->Lookup();
  54
  55
  56 */
  57
  58
  59 #include "Riostream.h"
  60 #include "TChain.h"
  61 #include "TTree.h"
  62 #include "TH1F.h"
  63 #include "TH2F.h"
  64 #include "TH3F.h"
  65 #include "THnSparse.h"
  66 #include "TList.h"
  67 #include "TMath.h"
  68 #include "TCanvas.h"
  69 #include "TFile.h"
  70 #include "TF1.h"
  71 #include "TVectorD.h"
  72 #include "TProfile.h"
  73 #include "TGraphErrors.h"
  74 #include "TCanvas.h"
  75
  76 #include "AliTPCclusterMI.h"
  77 #include "AliTPCseed.h"
  78 #include "AliESDVertex.h"
  79 #include "AliESDEvent.h"
  80 #include "AliESDfriend.h"
  81 #include "AliESDInputHandler.h"
  82 #include "AliAnalysisManager.h"
  83
  84 #include "AliTracker.h"
  85 #include "AliMagF.h"
  86 #include "AliTPCCalROC.h"
  87
  88 #include "AliLog.h"
  89
  90 #include "AliTPCcalibTime.h"
  91
  92 #include "TTreeStream.h"
  93 #include "AliTPCTracklet.h"
  94 #include "TTimeStamp.h"
  95 #include "AliTPCcalibDB.h"
  96 #include "AliTPCcalibLaser.h"
  97 #include "AliDCSSensorArray.h"
  98 #include "AliDCSSensor.h"
  99
 100 ClassImp(AliTPCcalibTime)
 101
 102
 103 AliTPCcalibTime::AliTPCcalibTime()
 104   :AliTPCcalibBase(),
 105    fHistDeDxTgl(0),
 106    fHistDeDx(0),
 107    fHistVdrift(0),
 108    fIntegrationTimeDeDx(0),
 109    fIntegrationTimeVdrift(0),
 110    fLaser(0),       // pointer to laser calibration
 111    fDz(0),          // current delta z
 112    fdEdx(0),        // current dEdx
 113    fdEdxRatio(0),   // current dEdx ratio
 114    fTl(0),          // current tan(lambda)
 115    fCutMaxD(5),        // maximal distance in rfi ditection
 116    fCutMaxDz(20),      // maximal distance in rfi ditection
 117    fCutTheta(0.03),    // maximal distan theta
 118    fCutMinDir(-0.99),  // direction vector products
 119    fMapDz(0),          //NEW! Tmap of V drifts for different triggers
 120    fTimeBins(0),
 121    fTimeStart(0),
 122    fTimeEnd(0),
 123    fPtBins(0),
 124    fPtStart(0),
 125    fPtEnd(0),
 126    fVdriftBins(0),
 127    fVdriftStart(0),
 128    fVdriftEnd(0),
 129    fRunBins(0),
 130    fRunStart(0),
 131    fRunEnd(0)
 132 //   fBinsVdrift(fTimeBins,fPtBins,fVdriftBins),
 133 //   fXminVdrift(fTimeStart,fPtStart,fVdriftStart),
 134 //   fXmaxVdrift(fTimeEnd,fPtEnd,fVdriftEnd)
 135
 136 {
 137   AliInfo("Default Constructor");
 138   for (Int_t i=0;i<3;i++) {
 139     fHistVdriftLaserA[i]=0;
 140     fHistVdriftLaserC[i]=0;
 141   }
 142 }
 143
 144
 145 AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeDeDx, Int_t deltaIntegrationTimeVdrift)
 146   :AliTPCcalibBase(),
 147    fHistDeDxTgl(0),
 148    fHistDeDx(0),
 149    fHistVdrift(0),
 150    fIntegrationTimeDeDx(0),
 151    fIntegrationTimeVdrift(0),
 152    fLaser(0),       // pointer to laser calibration
 153    fDz(0),          // current delta z
 154    fdEdx(0),        // current dEdx
 155    fdEdxRatio(0),   // current dEdx ratio
 156    fTl(0),          // current tan(lambda)
 157    fCutMaxD(5),        // maximal distance in rfi ditection
 158    fCutMaxDz(20),      // maximal distance in rfi ditection
 159    fCutTheta(0.03),    // maximal distan theta
 160    fCutMinDir(-0.99),  // direction vector products
 161    fMapDz(0),          //NEW! Tmap of V drifts for different triggers
 162    fTimeBins(0),
 163    fTimeStart(0),
 164    fTimeEnd(0),
 165    fPtBins(0),
 166    fPtStart(0),
 167    fPtEnd(0),
 168    fVdriftBins(0),
 169    fVdriftStart(0),
 170    fVdriftEnd(0),
 171    fRunBins(0),
 172    fRunStart(0),
 173    fRunEnd(0)
 174 {
 175
 176   SetName(name);
 177   SetTitle(title);
 178   for (Int_t i=0;i<3;i++) {
 179     fHistVdriftLaserA[i]=0;
 180     fHistVdriftLaserC[i]=0;
 181   }
 182
 183   AliInfo("Non Default Constructor");
 184
 185
 186   fIntegrationTimeDeDx = deltaIntegrationTimeDeDx;
 187   fIntegrationTimeVdrift = deltaIntegrationTimeVdrift;
 188
 189   Double_t deltaTime = EndTime - StartTime;
 190
 191   Int_t binsVdrift[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeVdrift), 100};
 192   Double_t xminVdrift[2] = {StartTime, -20};
 193   Double_t xmaxVdrift[2] = {EndTime, 20};
 194   fHistVdrift = new THnSparseF("HistVdrift","vDrift; time;#Delta z",2,binsVdrift,xminVdrift,xmaxVdrift);
 195
 196   Int_t binsDeDxTgl[4] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),30,100,10000};
 197   Double_t xminDeDxTgl[4] = {StartTime,-1,0.7,-0.5};
 198   Double_t xmaxDeDxTgl[4] = {EndTime,1,1.3,999.5};
 199   fHistDeDxTgl = new THnSparseF("HistDeDxTgl","dEdx vs tgl;time;tgl;dEdxUp/dEdxLow",3,binsDeDxTgl,xminDeDxTgl,xmaxDeDxTgl);
 200
 201   Int_t binsDeDx[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),100};
 202   Double_t xminDeDx[2] = {StartTime,1};
 203   Double_t xmaxDeDx[2] = {EndTime,100};
 204   fHistDeDx = new THnSparseF("HistDeDx","dEdx l;time;dEdx",2,binsDeDx,xminDeDx,xmaxDeDx);
 205
 206 //  fLaser =  new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
 207
 208   fTimeBins   =TMath::Nint(deltaTime/deltaIntegrationTimeVdrift);
 209   fTimeStart  =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi();
 210   fTimeEnd    =EndTime;   //(((TObjString*)(mapGRP->GetValue("fAliceStopTime")))->GetString()).Atoi();
 211   fPtBins     = 200;
 212   fPtStart    = -0.04;
 213   fPtEnd      =  0.04;
 214   fVdriftBins = 200;
 215   fVdriftStart= -20.0/500.0;
 216   fVdriftEnd  =  20.0/500.0;
 217   fRunBins    = 10000;
 218   fRunStart   = -0.5;
 219   fRunEnd     = 0.5;
 220
 221   Int_t    binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins };
 222   Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart  , fRunStart};
 223   Double_t xmaxVdriftLaser[4] = {fTimeEnd  , fPtEnd  , fVdriftEnd    , fRunEnd  };
 224
 225   for (Int_t i=0;i<3;i++) {
 226     fHistVdriftLaserA[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
 227     fHistVdriftLaserC[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
 228   }
 229   fBinsVdrift[0] = fTimeBins;
 230   fBinsVdrift[1] = fPtBins;
 231   fBinsVdrift[2] = fVdriftBins;
 232   fBinsVdrift[3] = fRunBins;
 233   fXminVdrift[0] = fTimeStart;
 234   fXminVdrift[1] = fPtStart;
 235   fXminVdrift[2] = fVdriftStart;
 236   fXminVdrift[3] = fRunStart;
 237   fXmaxVdrift[0] = fTimeEnd;
 238   fXmaxVdrift[1] = fPtEnd;
 239   fXmaxVdrift[2] = fVdriftEnd;
 240   fXmaxVdrift[3] = fRunEnd;
 241
 242   fMapDz=new TMap();
 243 }
 244
 245
 246
 247 AliTPCcalibTime::~AliTPCcalibTime(){
 248   //
 249   // Destructor
 250   //
 251   delete fHistDeDxTgl;
 252   delete fHistDeDx;
 253   delete fHistVdrift;
 254   for (Int_t i=0;i<3;i++){
 255     delete fHistVdriftLaserA[i];
 256     delete fHistVdriftLaserC[i];
 257   }
 258   fMapDz->SetOwner();
 259   fMapDz->Delete();
 260   delete fMapDz;
 261 }
 262
 263
 264 void AliTPCcalibTime::ResetCurrent(){
 265   //
 266   // reset current values
 267   //
 268   fDz=0;          // current delta z
 269   fdEdx=0;        // current dEdx
 270   fdEdxRatio=0;   // current dEdx ratio
 271   fTl=0;          // current tan(lambda)
 272
 273 }
 274
 275 Bool_t AliTPCcalibTime::IsLaser(AliESDEvent *event){
 276   return ((event->GetFiredTriggerClasses()).Contains("0LSR")==1);
 277 }
 278
 279 void AliTPCcalibTime::Process(AliESDEvent *event){
 280   if(!event) return;
 281   if (event->GetNumberOfTracks()<2) return;
 282   ResetCurrent();
 283
 284 //  if(IsLaser(event))
 285   ProcessLaser (event);
 286 //  else
 287   ProcessCosmic(event);
 288 }
 289
 290 void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){
 291   //
 292   // Fit drift velocity using laser
 293   //
 294   // 0. cuts
 295   const Int_t    kMinTracks     = 20;    // minimal number of laser tracks
 296   const Int_t    kMinTracksSide = 10;    // minimal number of tracks per side
 297   const Float_t  kMaxRMS        = 0.1;   // maximal RMS of tracks
 298   const Float_t  kMaxDeltaZ     = 3.;    // maximal deltaZ A-C side
 299   const Float_t  kMaxDeltaV     = 0.01;  // maximal deltaV A-C side
 300   const Float_t  kMaxDeltaY     = 2.;    // maximal deltaY A-C side
 301   /*
 302     TCut cutRMS("sqrt(laserA.fElements[4])<0.1&&sqrt(laserC.fElements[4])<0.1");
 303     TCut cutZ("abs(laserA.fElements[0]-laserC.fElements[0])<3");
 304     TCut cutV("abs(laserA.fElements[1]-laserC.fElements[1])<0.01");
 305     TCut cutY("abs(laserA.fElements[2]-laserC.fElements[2])<2");
 306     TCut cutAll = cutRMS+cutZ+cutV+cutY;
 307   */
 308   if (event->GetNumberOfTracks()<kMinTracks) return;
 309   //
 310   if(!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
 311   fLaser->Process(event);
 312   if (fLaser->GetNtracks()<kMinTracks) return;   // small amount of tracks cut
 313   if (fLaser->fFitAside->GetNrows()==0) return;  // no fit A side
 314   if (fLaser->fFitCside->GetNrows()==0) return;  // no fit C side
 315   //
 316   // debug streamer  - activate stream level
 317   // Use it for tuning of the cuts
 318   //
 319   if (fStreamLevel>0){
 320     printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
 321
 322     TTreeSRedirector *cstream = GetDebugStreamer();
 323     if (cstream){
 324       TTimeStamp tstamp(fTime);
 325       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 326       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 327       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 328       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 329       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 330       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 331       TVectorD vecGoofie(20);
 332       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 333       if (goofieArray){
 334         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 335           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 336           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 337         }
 338       }
 339       (*cstream)<<"laserInfo"<<
 340         "run="<<fRun<<              //  run number
 341         "event="<<fEvent<<          //  event number
 342         "time="<<fTime<<            //  time stamp of event
 343         "trigger="<<fTrigger<<      //  trigger
 344         "mag="<<fMagF<<             //  magnetic field
 345         // Environment values
 346         "press0="<<valuePressure0<<
 347         "press1="<<valuePressure1<<
 348         "pt0="<<ptrelative0<<
 349         "pt1="<<ptrelative1<<
 350         "temp0="<<temp0<<
 351         "temp1="<<temp1<<
 352         "vecGoofie.=<<"<<&vecGoofie<<
 353         //laser
 354         "laserA.="<<fLaser->fFitAside<<
 355         "laserC.="<<fLaser->fFitCside<<
 356         "laserAC.="<<fLaser->fFitACside<<
 357         "trigger="<<event->GetFiredTriggerClasses()<<
 358         "\n";
 359     }
 360   }
 361   //
 362   // Apply custs
 363   //
 364   if ((*fLaser->fFitAside)[3] <kMinTracksSide) return; // enough tracks A side
 365   if ((*fLaser->fFitCside)[3]<kMinTracksSide) return;  // enough tracks C side
 366   //
 367   if (TMath::Abs((*fLaser->fFitAside)[0]-(*fLaser->fFitCside)[0])>kMaxDeltaZ) return;
 368   if (TMath::Abs((*fLaser->fFitAside)[2]-(*fLaser->fFitCside)[2])>kMaxDeltaY) return;
 369   if (TMath::Abs((*fLaser->fFitAside)[1]-(*fLaser->fFitCside)[1])>kMaxDeltaV) return;
 370   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitAside)[4]))>kMaxRMS) return;
 371   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitCside)[4]))>kMaxRMS) return;
 372   //
 373   // fill histos
 374   //
 375   TVectorD vdriftA(5), vdriftC(5),vdriftAC(5);
 376   vdriftA=*(fLaser->fFitAside);
 377   vdriftC=*(fLaser->fFitCside);
 378   vdriftAC=*(fLaser->fFitACside);
 379   Int_t npointsA=0, npointsC=0;
 380   Float_t chi2A=0, chi2C=0;
 381   npointsA= TMath::Nint(vdriftA[3]);
 382   chi2A= vdriftA[4];
 383   npointsC= TMath::Nint(vdriftC[3]);
 384   chi2C= vdriftC[4];
 385   //
 386   //
 387
 388   //
 389   TTimeStamp tstamp(fTime);
 390   Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 391   Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 392   //
 393   Double_t vecVdriftLaserA[4]={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitAside))[1])-1,event->GetRunNumber()};
 394   Double_t vecVdriftLaserC[4]={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitCside))[1])-1,event->GetRunNumber()};
 395
 396   for (Int_t i=0;i<3;i++){
 397     if (i==0){ //z0 shift
 398       vecVdriftLaserA[3]=(*(fLaser->fFitAside))[0]/250.;
 399       vecVdriftLaserA[3]=(*(fLaser->fFitCside))[0]/250.;
 400     }
 401     if (i==1){ //vdrel shift
 402       vecVdriftLaserA[3]=1./(*(fLaser->fFitAside))[1]-1.;
 403       vecVdriftLaserA[3]=1./(*(fLaser->fFitCside))[1]-1.;
 404     }
 405     if (i==2){ //gy shift - full gy - full drift
 406       vecVdriftLaserA[3]=(*(fLaser->fFitAside))[2]/250.;
 407       vecVdriftLaserA[3]=(*(fLaser->fFitCside))[2]/250.;
 408     }
 409     fHistVdriftLaserA[i]->Fill(vecVdriftLaserA);
 410     fHistVdriftLaserC[i]->Fill(vecVdriftLaserC);
 411   }
 412   //
 413   //
 414   //
 415
 416
 417 }
 418
 419 void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){
 420
 421   if (!event) {
 422     Printf("ERROR: ESD not available");
 423     return;
 424   }
 425   if (event->GetTimeStamp() == 0 ) {
 426     Printf("no time stamp!");
 427     return;
 428   }
 429
 430
 431   UInt_t time = event->GetTimeStamp();
 432
 433   //fd
 434   // Find cosmic pairs
 435   //
 436   // Track0 is choosen in upper TPC part
 437   // Track1 is choosen in lower TPC part
 438   //
 439   Int_t ntracks=event->GetNumberOfTracks();
 440   if (ntracks==0) return;
 441   if (ntracks > 10) return; // temporary debug to remove LASER events
 442
 443
 444   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 445   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 446
 447   TObjArray  tpcSeeds(ntracks);
 448   Double_t vtxx[3]={0,0,0};
 449   Double_t svtxx[3]={0.000001,0.000001,100.};
 450   AliESDVertex vtx(vtxx,svtxx);
 451   //
 452   // track loop
 453   //
 454   for (Int_t i=0;i<ntracks;++i) {
 455    AliESDtrack *track = event->GetTrack(i);
 456
 457    const AliExternalTrackParam * trackIn = track->GetInnerParam();
 458    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 459    if (!trackIn) continue;
 460    if (!trackOut) continue;
 461
 462    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 463    TObject *calibObject;
 464    AliTPCseed *seed = 0;
 465    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 466      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 467    }
 468    if (seed) {
 469      tpcSeeds.AddAt(seed,i);
 470      if (track->GetTPCNcls() > 50) {
 471        Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 472        Double_t TPCsignal = seed->CookdEdxNorm(0.0,0.6,1,0,159,0x0,kTRUE,kTRUE);
 473        Double_t vecDeDx[2] = {time, TPCsignal};
 474        if (meanP > 12) {
 475          fdEdx = TPCsignal;
 476          fHistDeDx->Fill(vecDeDx);
 477        }
 478      }
 479    }
 480   }
 481
 482   if (ntracks<2) return;
 483   //
 484   // Find pairs
 485   //
 486   for (Int_t i=0;i<ntracks;++i) {
 487     AliESDtrack *track0 = event->GetTrack(i);
 488     // track0 - choosen upper part
 489     if (!track0) continue;
 490     if (!track0->GetOuterParam()) continue;
 491     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 492     Double_t d1[3];
 493     track0->GetDirection(d1);
 494     for (Int_t j=0;j<ntracks;++j) {
 495       if (i==j) continue;
 496        AliESDtrack *track1 = event->GetTrack(j);
 497        //track 1 lower part
 498        if (!track1) continue;
 499        if (!track1->GetOuterParam()) continue;
 500        if (track1->GetOuterParam()->GetAlpha()>0) continue;
 501        //
 502        Double_t d2[3];
 503        track1->GetDirection(d2);
 504
 505        AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 506        AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 507        if (! seed0) continue;
 508        if (! seed1) continue;
 509        Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
 510        Float_t dist0  = track0->GetLinearD(0,0);
 511        Float_t dist1  = track1->GetLinearD(0,0);
 512        //
 513        // conservative cuts - convergence to be guarantied
 514        // applying before track propagation
 515        if (TMath::Abs(dist0+dist1)>fCutMaxD) continue;   // distance to the 0,0
 516        if (dir>fCutMinDir) continue;               // direction vector product
 517        Float_t bz = AliTracker::GetBz();
 518        Float_t dvertex0[2];   //distance to 0,0
 519        Float_t dvertex1[2];   //distance to 0,0
 520        track0->GetDZ(0,0,0,bz,dvertex0);
 521        track1->GetDZ(0,0,0,bz,dvertex1);
 522        if (TMath::Abs(dvertex0[1])>250) continue;
 523        if (TMath::Abs(dvertex1[1])>250) continue;
 524        //
 525        //
 526        //
 527        Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
 528        AliExternalTrackParam param0(*track0);
 529        AliExternalTrackParam param1(*track1);
 530        //
 531        // Propagate using Magnetic field and correct fo material budget
 532        //
 533        AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 534        AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 535        //
 536        // Propagate rest to the 0,0 DCA - z should be ignored
 537        //
 538        //Bool_t b0 = ;
 539        param0.PropagateToDCA(&vtx,bz,1000);
 540        //Bool_t b1 =
 541        param1.PropagateToDCA(&vtx,bz,1000);
 542        //
 543        param0.GetDZ(0,0,0,bz,dvertex0);
 544        param1.GetDZ(0,0,0,bz,dvertex1);
 545        //
 546        Double_t xyz0[3];//,pxyz0[3];
 547        Double_t xyz1[3];//,pxyz1[3];
 548        param0.GetXYZ(xyz0);
 549        param1.GetXYZ(xyz1);
 550        Bool_t isPair = IsPair(&param0,&param1);
 551
 552        Double_t z0 = track0->GetOuterParam()->GetZ();
 553        Double_t z1 = track1->GetOuterParam()->GetZ();
 554
 555        Double_t z0inner = track0->GetInnerParam()->GetZ();
 556        Double_t z1inner = track1->GetInnerParam()->GetZ();
 557
 558        if (isPair && z0>0 && z0inner>0 && z1<0 && z1inner<0) {
 559          Double_t vecVdrift[2] = {time, param0.GetZ() - param1.GetZ()};
 560
 561          if (track0->GetTPCNcls() > 80) {
 562            fDz = param0.GetZ() - param1.GetZ();
 563            TTimeStamp tstamp(fTime);
 564            Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 565            Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 566            THnSparse* curHist   =(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
 567            if(!curHist){
 568              curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 569              fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
 570            }
 571            Double_t vecVdrift2[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
 572            curHist->Fill(vecVdrift2);
 573
 574            fHistVdrift->Fill(vecVdrift);
 575          }
 576        }
 577        if (isPair && z0 > 0 && z1 > 0) {
 578          if (track1->GetTPCNcls()> 110 && track0->GetTPCNcls()> 110 && seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE) > 0) {
 579            Float_t dedxratio = seed0->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)/seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE);
 580            Double_t vecDeDxTgl[3] = {time, track0->GetTgl(), seed0->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)/seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)};
 581            fHistDeDxTgl->Fill(vecDeDxTgl);
 582            fdEdxRatio = dedxratio;
 583            fTl = track0->GetTgl() ;
 584          }
 585        }
 586
 587     } // end 2nd order loop
 588   } // end 1st order loop
 589
 590   if (fStreamLevel>0){
 591     TTreeSRedirector *cstream = GetDebugStreamer();
 592     if (cstream){
 593       TTimeStamp tstamp(fTime);
 594       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 595       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 596       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 597       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 598       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 599       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 600       TVectorD vecGoofie(20);
 601       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 602       if (goofieArray){
 603         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 604           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 605           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 606         }
 607       }
 608       (*cstream)<<"timeInfo"<<
 609         "run="<<fRun<<              //  run number
 610         "event="<<fEvent<<          //  event number
 611         "time="<<fTime<<            //  time stamp of event
 612         "trigger="<<fTrigger<<      //  trigger
 613         "mag="<<fMagF<<             //  magnetic field
 614         // Environment values
 615         "press0="<<valuePressure0<<
 616         "press1="<<valuePressure1<<
 617         "pt0="<<ptrelative0<<
 618         "pt1="<<ptrelative1<<
 619         "temp0="<<temp0<<
 620         "temp1="<<temp1<<
 621         "vecGoofie.=<<"<<&vecGoofie<<
 622         //
 623         // accumulated values
 624         //
 625         "fDz="<<fDz<<          //! current delta z
 626         "fdEdx="<<fdEdx<<        //! current dEdx
 627         "fdEdxRatio="<<fdEdxRatio<<   //! current dEdx ratio
 628         "fTl="<<fTl<<          //! current tan(lambda)
 629         "trigger="<<event->GetFiredTriggerClasses()<<
 630         "\n";
 631     }
 632   }
 633   printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
 634 }
 635
 636 void AliTPCcalibTime::Analyze() {
 637   //
 638   //
 639   //
 640   TH2D * hVdrift = GetHistVdrift()->Projection(1,0);
 641   hVdrift->Draw();
 642 }
 643
 644 Long64_t AliTPCcalibTime::Merge(TCollection *li) {
 645
 646   TIterator* iter = li->MakeIterator();
 647   AliTPCcalibTime* cal = 0;
 648
 649   while ((cal = (AliTPCcalibTime*)iter->Next())) {
 650     if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
 651       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 652       return -1;
 653     }
 654
 655     // add histograms here...
 656     fHistDeDxTgl->Add(cal->GetHistDeDxVsTgl());
 657     fHistVdrift->Add(cal->GetHistVdrift());
 658     fHistDeDx->Add(cal->GetHistDeDx());
 659
 660     for (Int_t imeas=0; imeas<3; imeas++){
 661       if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
 662         fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
 663         fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
 664       }
 665     }
 666     TMap * addMap=cal->GetMapDz();
 667     if(!addMap) return 0;
 668     TIterator* iterator = addMap->MakeIterator();
 669     iterator->Reset();
 670     TPair* addPair=0;
 671     while((addPair=(TPair *)(addMap->FindObject(iterator->Next())))){
 672       THnSparse* addHist=dynamic_cast<THnSparseF*>(addPair->Value());
 673       if (!addHist) continue;
 674       addHist->Print();
 675       THnSparse* localHist=dynamic_cast<THnSparseF*>(fMapDz->GetValue(addHist->GetName()));
 676       if(!localHist){
 677         localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 678         fMapDz->Add(new TObjString(addHist->GetName()),localHist);
 679       }
 680       localHist->Add(addHist);
 681     }
 682   }
 683   return 0;
 684 }
 685
 686 Bool_t  AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 687   //
 688   //
 689   /*
 690   // 0. Same direction - OPOSITE  - cutDir +cutT
 691   TCut cutDir("cutDir","dir<-0.99")
 692   // 1.
 693   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 694   //
 695   // 2. The same rphi
 696   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 697   //
 698   //
 699   //
 700   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 701   // 1/Pt diff cut
 702   */
 703   const Double_t *p0 = tr0->GetParameter();
 704   const Double_t *p1 = tr1->GetParameter();
 705   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 706   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 707   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz)  return kFALSE;
 708   Double_t d0[3], d1[3];
 709   tr0->GetDirection(d0);
 710   tr1->GetDirection(d1);
 711   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 712   //
 713   return kTRUE;
 714 }
 715
 716