TPC/AliTPCcalibTime.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /*
  18 Comments to be written here:
  19
  20 1. What do we calibrate.
  21
  22   Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
  23
  24   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
  25
  26 2. How to interpret results
  27
  28 3. Simple example
  29
  30   a) determine the required time range:
  31
  32   AliXRDPROOFtoolkit tool;
  33   TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
  34   chain->Draw("GetTimeStamp()")
  35
  36   b) analyse calibration object on Proof in calibration train
  37
  38   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift);
  39
  40   c) plot results
  41   .x ~/NimStyle.C
  42   gSystem->Load("libANALYSIS");
  43   gSystem->Load("libTPCcalib");
  44
  45   TFile f("CalibObjects.root");
  46   AliTPCcalibTime *cal = (AliTPCcalibTime *)f->Get("TPCCalib")->FindObject("calibTime");
  47   cal->GetHistoDrift("all")->Projection(2,0)->Draw()
  48   cal->GetFitDrift("all")->Draw("lp")
  49
  50 4. Analysis using debug streamers.
  51
  52   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  53   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  54   AliXRDPROOFtoolkit tool;
  55   TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
  56   chainTime->Lookup();
  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    fLaser(0),       // pointer to laser calibration
 106    fDz(0),          // current delta z
 107    fCutMaxD(3),        // maximal distance in rfi ditection
 108    fCutMaxDz(25),      // maximal distance in rfi ditection
 109    fCutTheta(0.03),    // maximal distan theta
 110    fCutMinDir(-0.99),  // direction vector products
 111    fCutTracks(10),
 112    fMapDz(0),          //NEW! Tmap of V drifts for different triggers
 113    fTimeBins(0),
 114    fTimeStart(0),
 115    fTimeEnd(0),
 116    fPtBins(0),
 117    fPtStart(0),
 118    fPtEnd(0),
 119    fVdriftBins(0),
 120    fVdriftStart(0),
 121    fVdriftEnd(0),
 122    fRunBins(0),
 123    fRunStart(0),
 124    fRunEnd(0)
 125 //   fBinsVdrift(fTimeBins,fPtBins,fVdriftBins),
 126 //   fXminVdrift(fTimeStart,fPtStart,fVdriftStart),
 127 //   fXmaxVdrift(fTimeEnd,fPtEnd,fVdriftEnd)
 128
 129 {
 130   AliInfo("Default Constructor");
 131   for (Int_t i=0;i<3;i++) {
 132     fHistVdriftLaserA[i]=0;
 133     fHistVdriftLaserC[i]=0;
 134   }
 135   for (Int_t i=0;i<10;i++) {
 136     fCosmiMatchingHisto[i]=0;
 137   }
 138 }
 139
 140
 141 AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeVdrift)
 142   :AliTPCcalibBase(),
 143    fLaser(0),            // pointer to laser calibration
 144    fDz(0),               // current delta z
 145    fCutMaxD(5*0.5356),   // maximal distance in rfi ditection
 146    fCutMaxDz(5*4.541),   // maximal distance in rfi ditection
 147    fCutTheta(5*0.004644),// maximal distan theta
 148    fCutMinDir(-0.99),    // direction vector products
 149    fCutTracks(10),
 150    fMapDz(0),            //Tmap of V drifts for different triggers
 151    fTimeBins(0),
 152    fTimeStart(0),
 153    fTimeEnd(0),
 154    fPtBins(0),
 155    fPtStart(0),
 156    fPtEnd(0),
 157    fVdriftBins(0),
 158    fVdriftStart(0),
 159    fVdriftEnd(0),
 160    fRunBins(0),
 161    fRunStart(0),
 162    fRunEnd(0)
 163 {
 164
 165   SetName(name);
 166   SetTitle(title);
 167   for (Int_t i=0;i<3;i++) {
 168     fHistVdriftLaserA[i]=0;
 169     fHistVdriftLaserC[i]=0;
 170   }
 171
 172   AliInfo("Non Default Constructor");
 173
 174   fTimeBins   =(EndTime-StartTime)/deltaIntegrationTimeVdrift;
 175   fTimeStart  =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi();
 176   fTimeEnd    =EndTime;   //(((TObjString*)(mapGRP->GetValue("fAliceStopTime")))->GetString()).Atoi();
 177   fPtBins     = 200;
 178   fPtStart    = -0.04;
 179   fPtEnd      =  0.04;
 180   fVdriftBins = 200;
 181   fVdriftStart= -20.0/500.0;
 182   fVdriftEnd  =  20.0/500.0;
 183   fRunBins    = 100000;
 184   fRunStart   = -0.5;
 185   fRunEnd     = 99999.5;
 186
 187   Int_t    binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins };
 188   Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart  , fRunStart};
 189   Double_t xmaxVdriftLaser[4] = {fTimeEnd  , fPtEnd  , fVdriftEnd    , fRunEnd  };
 190
 191   for (Int_t i=0;i<3;i++) {
 192     fHistVdriftLaserA[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
 193     fHistVdriftLaserC[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
 194   }
 195   fBinsVdrift[0] = fTimeBins;
 196   fBinsVdrift[1] = fPtBins;
 197   fBinsVdrift[2] = fVdriftBins;
 198   fBinsVdrift[3] = fRunBins;
 199   fXminVdrift[0] = fTimeStart;
 200   fXminVdrift[1] = fPtStart;
 201   fXminVdrift[2] = fVdriftStart;
 202   fXminVdrift[3] = fRunStart;
 203   fXmaxVdrift[0] = fTimeEnd;
 204   fXmaxVdrift[1] = fPtEnd;
 205   fXmaxVdrift[2] = fVdriftEnd;
 206   fXmaxVdrift[3] = fRunEnd;
 207
 208   fMapDz=new TMap();
 209
 210   fCosmiMatchingHisto[0]=new TH1F("Cosmics matching","p0-all"   ,100,-10*0.5356  ,10*0.5356  );
 211   fCosmiMatchingHisto[1]=new TH1F("Cosmics matching","p1-all"   ,100,-10*4.541   ,10*4.541   );
 212   fCosmiMatchingHisto[2]=new TH1F("Cosmics matching","p2-all"   ,100,-10*0.01134 ,10*0.01134 );
 213   fCosmiMatchingHisto[3]=new TH1F("Cosmics matching","p3-all"   ,100,-10*0.004644,10*0.004644);
 214   fCosmiMatchingHisto[4]=new TH1F("Cosmics matching","p4-all"   ,100,-10*0.03773 ,10*0.03773 );
 215   fCosmiMatchingHisto[5]=new TH1F("Cosmics matching","p0-isPair",100,-10*0.5356  ,10*0.5356  );
 216   fCosmiMatchingHisto[6]=new TH1F("Cosmics matching","p1-isPair",100,-10*4.541   ,10*4.541   );
 217   fCosmiMatchingHisto[7]=new TH1F("Cosmics matching","p2-isPair",100,-10*0.01134 ,10*0.01134 );
 218   fCosmiMatchingHisto[8]=new TH1F("Cosmics matching","p3-isPair",100,-10*0.004644,10*0.004644);
 219   fCosmiMatchingHisto[9]=new TH1F("Cosmics matching","p4-isPair",100,-10*0.03773 ,10*0.03773 );
 220 //  Char_t nameHisto[3]={'p','0','\n'};
 221 //  for (Int_t i=0;i<10;i++){
 222 //    fCosmiMatchingHisto[i]=new TH1F("Cosmics matching",nameHisto,8192,0,0);
 223 //    nameHisto[1]++;
 224 //    if(i==4) nameHisto[1]='0';
 225 //  }
 226 }
 227
 228 AliTPCcalibTime::~AliTPCcalibTime(){
 229   //
 230   // Destructor
 231   //
 232   for(Int_t i=0;i<3;i++){
 233     if(fHistVdriftLaserA[i]){
 234       delete fHistVdriftLaserA[i];
 235       fHistVdriftLaserA[i]=NULL;
 236     }
 237     if(fHistVdriftLaserC[i]){
 238       delete fHistVdriftLaserC[i];
 239       fHistVdriftLaserC[i]=NULL;
 240     }
 241   }
 242   if(fMapDz){
 243     fMapDz->SetOwner();
 244     fMapDz->Delete();
 245     delete fMapDz;
 246     fMapDz=NULL;
 247   }
 248   for(Int_t i=0;i<5;i++){
 249     if(fCosmiMatchingHisto[i]){
 250       delete fCosmiMatchingHisto[i];
 251       fCosmiMatchingHisto[i]=NULL;
 252     }
 253   }
 254 }
 255
 256 void AliTPCcalibTime::ResetCurrent(){
 257   //
 258   // reset current values
 259   //
 260   fDz=0;          // current delta z
 261 }
 262
 263 Bool_t AliTPCcalibTime::IsLaser(AliESDEvent *event){
 264   return ((event->GetFiredTriggerClasses()).Contains("0LSR")==1);
 265 }
 266
 267 void AliTPCcalibTime::Process(AliESDEvent *event){
 268   if(!event) return;
 269   if (event->GetNumberOfTracks()<2) return;
 270   ResetCurrent();
 271
 272 //  if(IsLaser(event))
 273   ProcessLaser (event);
 274 //  else
 275   ProcessCosmic(event);
 276 }
 277
 278 void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){
 279   //
 280   // Fit drift velocity using laser
 281   //
 282   // 0. cuts
 283   const Int_t    kMinTracks     = 20;    // minimal number of laser tracks
 284   const Int_t    kMinTracksSide = 10;    // minimal number of tracks per side
 285   const Float_t  kMaxRMS        = 0.1;   // maximal RMS of tracks
 286   const Float_t  kMaxDeltaZ     = 3.;    // maximal deltaZ A-C side
 287   const Float_t  kMaxDeltaV     = 0.01;  // maximal deltaV A-C side
 288   const Float_t  kMaxDeltaY     = 2.;    // maximal deltaY A-C side
 289   /*
 290     TCut cutRMS("sqrt(laserA.fElements[4])<0.1&&sqrt(laserC.fElements[4])<0.1");
 291     TCut cutZ("abs(laserA.fElements[0]-laserC.fElements[0])<3");
 292     TCut cutV("abs(laserA.fElements[1]-laserC.fElements[1])<0.01");
 293     TCut cutY("abs(laserA.fElements[2]-laserC.fElements[2])<2");
 294     TCut cutAll = cutRMS+cutZ+cutV+cutY;
 295   */
 296   if (event->GetNumberOfTracks()<kMinTracks) return;
 297   //
 298   if(!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
 299   fLaser->Process(event);
 300   if (fLaser->GetNtracks()<kMinTracks) return;   // small amount of tracks cut
 301   if (fLaser->fFitAside->GetNrows()==0) return;  // no fit A side
 302   if (fLaser->fFitCside->GetNrows()==0) return;  // no fit C side
 303   //
 304   // debug streamer  - activate stream level
 305   // Use it for tuning of the cuts
 306   //
 307   if (fStreamLevel>0){
 308     printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
 309
 310     TTreeSRedirector *cstream = GetDebugStreamer();
 311     if (cstream){
 312       TTimeStamp tstamp(fTime);
 313       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 314       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 315       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 316       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 317       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 318       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 319       TVectorD vecGoofie(20);
 320       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 321       if (goofieArray){
 322         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 323           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 324           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 325         }
 326       }
 327       (*cstream)<<"laserInfo"<<
 328         "run="<<fRun<<              //  run number
 329         "event="<<fEvent<<          //  event number
 330         "time="<<fTime<<            //  time stamp of event
 331         "trigger="<<fTrigger<<      //  trigger
 332         "mag="<<fMagF<<             //  magnetic field
 333         // Environment values
 334         "press0="<<valuePressure0<<
 335         "press1="<<valuePressure1<<
 336         "pt0="<<ptrelative0<<
 337         "pt1="<<ptrelative1<<
 338         "temp0="<<temp0<<
 339         "temp1="<<temp1<<
 340         "vecGoofie.=<<"<<&vecGoofie<<
 341         //laser
 342         "laserA.="<<fLaser->fFitAside<<
 343         "laserC.="<<fLaser->fFitCside<<
 344         "laserAC.="<<fLaser->fFitACside<<
 345         "trigger="<<event->GetFiredTriggerClasses()<<
 346         "\n";
 347     }
 348   }
 349   //
 350   // Apply custs
 351   //
 352   if ((*fLaser->fFitAside)[3] <kMinTracksSide) return; // enough tracks A side
 353   if ((*fLaser->fFitCside)[3]<kMinTracksSide) return;  // enough tracks C side
 354   //
 355   if (TMath::Abs((*fLaser->fFitAside)[0]-(*fLaser->fFitCside)[0])>kMaxDeltaZ) return;
 356   if (TMath::Abs((*fLaser->fFitAside)[2]-(*fLaser->fFitCside)[2])>kMaxDeltaY) return;
 357   if (TMath::Abs((*fLaser->fFitAside)[1]-(*fLaser->fFitCside)[1])>kMaxDeltaV) return;
 358   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitAside)[4]))>kMaxRMS) return;
 359   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitCside)[4]))>kMaxRMS) return;
 360   //
 361   // fill histos
 362   //
 363   TVectorD vdriftA(5), vdriftC(5),vdriftAC(5);
 364   vdriftA=*(fLaser->fFitAside);
 365   vdriftC=*(fLaser->fFitCside);
 366   vdriftAC=*(fLaser->fFitACside);
 367   Int_t npointsA=0, npointsC=0;
 368   Float_t chi2A=0, chi2C=0;
 369   npointsA= TMath::Nint(vdriftA[3]);
 370   chi2A= vdriftA[4];
 371   npointsC= TMath::Nint(vdriftC[3]);
 372   chi2C= vdriftC[4];
 373
 374   TTimeStamp tstamp(fTime);
 375   Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 376   Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 377
 378   Double_t vecVdriftLaserA[4]={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitAside))[1])-1,event->GetRunNumber()};
 379   Double_t vecVdriftLaserC[4]={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitCside))[1])-1,event->GetRunNumber()};
 380
 381   for (Int_t i=0;i<3;i++){
 382     if (i==0){ //z0 shift
 383       vecVdriftLaserA[3]=(*(fLaser->fFitAside))[0]/250.;
 384       vecVdriftLaserA[3]=(*(fLaser->fFitCside))[0]/250.;
 385     }
 386     if (i==1){ //vdrel shift
 387       vecVdriftLaserA[3]=1./(*(fLaser->fFitAside))[1]-1.;
 388       vecVdriftLaserA[3]=1./(*(fLaser->fFitCside))[1]-1.;
 389     }
 390     if (i==2){ //gy shift - full gy - full drift
 391       vecVdriftLaserA[3]=(*(fLaser->fFitAside))[2]/250.;
 392       vecVdriftLaserA[3]=(*(fLaser->fFitCside))[2]/250.;
 393     }
 394     fHistVdriftLaserA[i]->Fill(vecVdriftLaserA);
 395     fHistVdriftLaserC[i]->Fill(vecVdriftLaserC);
 396   }
 397 }
 398
 399 void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){
 400   if (!event) {
 401     Printf("ERROR: ESD not available");
 402     return;
 403   }
 404   if (event->GetTimeStamp() == 0 ) {
 405     Printf("no time stamp!");
 406     return;
 407   }
 408
 409   //fd
 410   // Find cosmic pairs
 411   //
 412   // Track0 is choosen in upper TPC part
 413   // Track1 is choosen in lower TPC part
 414   //
 415   Int_t ntracks=event->GetNumberOfTracks();
 416   if (ntracks==0) return;
 417   if (ntracks > fCutTracks) return;
 418
 419   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 420   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 421
 422   TObjArray  tpcSeeds(ntracks);
 423   Double_t vtxx[3]={0,0,0};
 424   Double_t svtxx[3]={0.000001,0.000001,100.};
 425   AliESDVertex vtx(vtxx,svtxx);
 426   //
 427   // track loop
 428   //
 429   for (Int_t i=0;i<ntracks;++i) {
 430     AliESDtrack *track = event->GetTrack(i);
 431
 432     const AliExternalTrackParam * trackIn = track->GetInnerParam();
 433     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 434     if (!trackIn) continue;
 435     if (!trackOut) continue;
 436
 437     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 438     TObject *calibObject;
 439     AliTPCseed *seed = 0;
 440     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 441     if (seed) tpcSeeds.AddAt(seed,i);
 442   }
 443   if (ntracks<2) return;
 444   //
 445   // Find pairs
 446   //
 447
 448   for (Int_t i=0;i<ntracks;++i) {
 449     AliESDtrack *track0 = event->GetTrack(i);
 450     // track0 - choosen upper part
 451     if (!track0) continue;
 452     if (!track0->GetOuterParam()) continue;
 453     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 454     Double_t d1[3];
 455     track0->GetDirection(d1);
 456     for (Int_t j=0;j<ntracks;++j) {
 457       if (i==j) continue;
 458       AliESDtrack *track1 = event->GetTrack(j);
 459       //track 1 lower part
 460       if (!track1) continue;
 461       if (!track1->GetOuterParam()) continue;
 462       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 463       //
 464       Double_t d2[3];
 465       track1->GetDirection(d2);
 466
 467       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 468       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 469       if (! seed0) continue;
 470       if (! seed1) continue;
 471       Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
 472       Float_t dist0  = track0->GetLinearD(0,0);
 473       Float_t dist1  = track1->GetLinearD(0,0);
 474       //
 475       // conservative cuts - convergence to be guarantied
 476       // applying before track propagation
 477       if (TMath::Abs(dist0+dist1)>fCutMaxD) continue;   // distance to the 0,0
 478       if (dir>fCutMinDir) continue;               // direction vector product
 479       Float_t bz = AliTracker::GetBz();
 480       Float_t dvertex0[2];   //distance to 0,0
 481       Float_t dvertex1[2];   //distance to 0,0
 482       track0->GetDZ(0,0,0,bz,dvertex0);
 483       track1->GetDZ(0,0,0,bz,dvertex1);
 484       if (TMath::Abs(dvertex0[1])>250) continue;
 485       if (TMath::Abs(dvertex1[1])>250) continue;
 486       //
 487       //
 488       //
 489       Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
 490       AliExternalTrackParam param0(*track0);
 491       AliExternalTrackParam param1(*track1);
 492       //
 493       // Propagate using Magnetic field and correct fo material budget
 494       //
 495       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 496       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 497       //
 498       // Propagate rest to the 0,0 DCA - z should be ignored
 499       //
 500       //Bool_t b0 = ;
 501       param0.PropagateToDCA(&vtx,bz,1000);
 502       //Bool_t b1 =
 503       param1.PropagateToDCA(&vtx,bz,1000);
 504       //
 505       param0.GetDZ(0,0,0,bz,dvertex0);
 506       param1.GetDZ(0,0,0,bz,dvertex1);
 507       //
 508       Double_t xyz0[3];//,pxyz0[3];
 509       Double_t xyz1[3];//,pxyz1[3];
 510       param0.GetXYZ(xyz0);
 511       param1.GetXYZ(xyz1);
 512       Bool_t isPair = IsPair(&param0,&param1);
 513       Bool_t isCross = IsCross(track0, track1);
 514
 515 //      Double_t z0 = track0->GetOuterParam()->GetZ();
 516 //      Double_t z1 = track1->GetOuterParam()->GetZ();
 517
 518 //      Double_t z0inner = track0->GetInnerParam()->GetZ();
 519 //      Double_t z1inner = track1->GetInnerParam()->GetZ();
 520
 521       if (isPair && isCross){
 522         if (track0->GetTPCNcls() > 80) {
 523           fDz = param0.GetZ() - param1.GetZ();
 524           if(track0->GetOuterParam()->GetZ()<0) fDz=-fDz;
 525           TTimeStamp tstamp(fTime);
 526           Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 527           Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 528           Double_t vecVdrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
 529           THnSparse* curHist=0;
 530
 531           curHist=(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
 532           if(!curHist){
 533             curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 534             fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
 535           }
 536           curHist->Fill(vecVdrift);
 537
 538           curHist=(THnSparseF*)(fMapDz->GetValue("all"));
 539           if(!curHist){
 540             curHist=new THnSparseF("all","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 541             fMapDz->Add(new TObjString("all"),curHist);
 542           }
 543           curHist->Fill(vecVdrift);
 544         }
 545       }
 546       TTreeSRedirector *cstream = GetDebugStreamer();
 547       if (fStreamLevel>0){
 548         if (cstream){
 549         (*cstream)<<"trackInfo"<<
 550           "tr0.="<<track0<<
 551           "tr1.="<<track1<<
 552           "p0.="<<&param0<<
 553           "p1.="<<&param1<<
 554           "isPair="<<isPair<<
 555           "isCross="<<isCross<<
 556           "fDz="<<fDz<<
 557           "fRun="<<fRun<<
 558           "fTime="<<fTime<<
 559           "\n";
 560         }
 561       }
 562     } // end 2nd order loop
 563   } // end 1st order loop
 564
 565   if (fStreamLevel>0){
 566     TTreeSRedirector *cstream = GetDebugStreamer();
 567     if (cstream){
 568       TTimeStamp tstamp(fTime);
 569       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 570       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 571       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 572       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 573       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 574       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 575       TVectorD vecGoofie(20);
 576       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 577       if (goofieArray){
 578         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 579           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 580           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 581         }
 582       }
 583       (*cstream)<<"timeInfo"<<
 584         "run="<<fRun<<              //  run number
 585         "event="<<fEvent<<          //  event number
 586         "time="<<fTime<<            //  time stamp of event
 587         "trigger="<<fTrigger<<      //  trigger
 588         "mag="<<fMagF<<             //  magnetic field
 589         // Environment values
 590         "press0="<<valuePressure0<<
 591         "press1="<<valuePressure1<<
 592         "pt0="<<ptrelative0<<
 593         "pt1="<<ptrelative1<<
 594         "temp0="<<temp0<<
 595         "temp1="<<temp1<<
 596         "vecGoofie.=<<"<<&vecGoofie<<
 597         //
 598         // accumulated values
 599         //
 600         "fDz="<<fDz<<          //! current delta z
 601         "trigger="<<event->GetFiredTriggerClasses()<<
 602         "\n";
 603     }
 604   }
 605   printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
 606 }
 607
 608 void AliTPCcalibTime::Analyze(){}
 609
 610 THnSparse* AliTPCcalibTime::GetHistoDrift(TObjString* name){
 611   return (THnSparseF*)(fMapDz->GetValue(name));
 612 }
 613 THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name){
 614   TObjString* objName=new TObjString(name);
 615   THnSparse* histoDrift=0;
 616   if(objName){
 617     histoDrift=GetHistoDrift(objName);
 618     delete objName;
 619     objName=0;
 620   }
 621   return histoDrift;
 622 }
 623
 624 TMap* AliTPCcalibTime::GetHistoDrift(){
 625   return fMapDz;
 626 }
 627
 628 TGraphErrors* AliTPCcalibTime::GetGraphDrift(TObjString* name){
 629   THnSparse* histoDrift=GetHistoDrift(name);
 630   TGraphErrors*    graphDrift=0;
 631   if(histoDrift) graphDrift=FitSlices(histoDrift,2,0,400,100,0.05,0.95, kTRUE);
 632   return graphDrift;
 633 }
 634
 635 TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){
 636   TObjString* objName=new TObjString(name);
 637   TGraphErrors* graphDrift=0;
 638   if(objName){
 639     graphDrift=GetGraphDrift(objName);
 640     delete objName;
 641     objName=0;
 642   }
 643   return graphDrift;
 644 }
 645
 646 TMap* AliTPCcalibTime::GetGraphDrift(){
 647   TMap* mapGraphDrift=new TMap();
 648   TIterator* iterator=fMapDz->MakeIterator();
 649   iterator->Reset();
 650   TPair* addPair=0;
 651   while((addPair=(TPair*)(fMapDz->FindObject(iterator->Next())))) mapGraphDrift->Add((TObjString*)addPair->Key(), GetGraphDrift((TObjString*)addPair->Key()));
 652   return mapGraphDrift;
 653 }
 654
 655 TGraph* AliTPCcalibTime::GetFitDrift(TObjString* name){
 656   TGraph* graphDrift=GetGraphDrift(name);
 657   TGraph* fitDrift=0;
 658   if(graphDrift && graphDrift->GetN()){
 659     AliSplineFit fit;
 660     fit.SetGraph(graphDrift);
 661     fit.SetMinPoints(graphDrift->GetN()+1);
 662     fit.InitKnots(graphDrift,2,0,0.001);
 663     fit.SplineFit(0);
 664     fitDrift=fit.MakeGraph(graphDrift->GetX()[0],graphDrift->GetX()[graphDrift->GetN()-1],50000,0);
 665     delete graphDrift;
 666     graphDrift=0;
 667   }
 668   return fitDrift;
 669 }
 670
 671 TGraph* AliTPCcalibTime::GetFitDrift(const char* name){
 672   TObjString* objName=new TObjString(name);
 673   TGraph* fitDrift=0;
 674   if(objName){
 675     fitDrift=GetFitDrift(objName);
 676     delete objName;
 677     objName=0;
 678   }
 679   return fitDrift;
 680 }
 681
 682 TMap* AliTPCcalibTime::GetFitDrift(){
 683   TMap* mapFitDrift=new TMap();
 684   TIterator* iterator = fMapDz->MakeIterator();
 685   iterator->Reset();
 686   TPair* addPair=0;
 687   while((addPair=(TPair*)(fMapDz->FindObject(iterator->Next())))) mapFitDrift->Add((TObjString*)addPair->Key(), GetFitDrift((TObjString*)addPair->Key()));
 688   return mapFitDrift;
 689 }
 690 Long64_t AliTPCcalibTime::Merge(TCollection *li) {
 691
 692   TIterator* iter = li->MakeIterator();
 693   AliTPCcalibTime* cal = 0;
 694
 695   while ((cal = (AliTPCcalibTime*)iter->Next())) {
 696     if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
 697       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 698       return -1;
 699     }
 700     for (Int_t imeas=0; imeas<3; imeas++){
 701       if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
 702         fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
 703         fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
 704       }
 705     }
 706     TMap * addMap=cal->GetHistoDrift();
 707     if(!addMap) return 0;
 708     TIterator* iterator = addMap->MakeIterator();
 709     iterator->Reset();
 710     TPair* addPair=0;
 711     while((addPair=(TPair *)(addMap->FindObject(iterator->Next())))){
 712       THnSparse* addHist=dynamic_cast<THnSparseF*>(addPair->Value());
 713       if (!addHist) continue;
 714       addHist->Print();
 715       THnSparse* localHist=dynamic_cast<THnSparseF*>(fMapDz->GetValue(addHist->GetName()));
 716       if(!localHist){
 717         localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 718         fMapDz->Add(new TObjString(addHist->GetName()),localHist);
 719       }
 720       localHist->Add(addHist);
 721     }
 722     for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
 723   }
 724   return 0;
 725 }
 726
 727 Bool_t  AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 728   /*
 729   // 0. Same direction - OPOSITE  - cutDir +cutT
 730   TCut cutDir("cutDir","dir<-0.99")
 731   // 1.
 732   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 733   //
 734   // 2. The same rphi
 735   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 736   //
 737   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 738   // 1/Pt diff cut
 739   */
 740   const Double_t *p0 = tr0->GetParameter();
 741   const Double_t *p1 = tr1->GetParameter();
 742   fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
 743   fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
 744   fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
 745   fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
 746   fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
 747
 748   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 749   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 750   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz)  return kFALSE;
 751   Double_t d0[3], d1[3];
 752   tr0->GetDirection(d0);
 753   tr1->GetDirection(d1);
 754   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 755
 756   fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
 757   fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
 758   fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
 759   fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
 760   fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
 761
 762   return kTRUE;
 763 }
 764 Bool_t AliTPCcalibTime::IsCross(AliESDtrack *tr0, AliESDtrack *tr1){
 765   return  tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0 && tr0->GetInnerParam()->GetZ()*tr1->GetInnerParam()->GetZ()<0 && tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0 && tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0;
 766 }
 767
 768 /*
 769 chainDrift->Draw("p0.fP[0]+p1.fP[0]","isPair");
 770   mean ~-0.02  ~-0.03913
 771   RMS  ~ 0.5   ~ 0.5356    --> 3    (fCutMaxD)
 772
 773 chainDrift->Draw("p0.fP[1]-p1.fP[1]","isPair");
 774   mean         ~ 0.1855
 775   RMS          ~ 4.541     -->25    (fCutMaxDz)
 776
 777 chainDrift->Draw("p1.fAlpha-p0.fAlpha+pi","isPair");
 778 //chainDrift->Draw("p1.fAlpha+p0.fAlpha","isPair");
 779 //chainDrift->Draw("p1.fP[2]-p0.fP[2]+pi","isPair");
 780 //chainDrift->Draw("p1.fP[2]+p0.fP[2]","isPair");
 781   mean ~ 0     ~ 0.001898
 782   RMS  ~ 0.009 ~ 0.01134   --> 0.06
 783
 784 chainDrift->Draw("p0.fP[3]+p1.fP[3]","isPair");
 785   mean ~ 0.0013 ~ 0.001539
 786   RMS  ~ 0.003  ~ 0.004644 --> 0.03  (fCutTheta)
 787
 788 chainDrift->Draw("p1.fP[4]+p0.fP[4]>>his(100,-0.2,0.2)","isPair")
 789   mean ~ 0.012  ~-0.0009729
 790   RMS  ~ 0.036  ~ 0.03773  --> 0.2
 791 */
 792