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
  40 TFile f("CalibObjects.root");
  41 AliTPCcalibTime *cal = (AliTPCcalibTime *)f->Get("TPCCalib")->FindObject("cosmicTime");
  42 cal->GetHistVdrift()->Projection(1,0)->Draw()
  43
  44     4. Analysis using debug streamers.
  45
  46     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  47     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  48     AliXRDPROOFtoolkit tool;
  49     TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
  50     chainTime->Lookup();
  51
  52
  53 */
  54
  55
  56 #include "Riostream.h"
  57 #include "TChain.h"
  58 #include "TTree.h"
  59 #include "TH1F.h"
  60 #include "TH2F.h"
  61 #include "TH3F.h"
  62 #include "THnSparse.h"
  63 #include "TList.h"
  64 #include "TMath.h"
  65 #include "TCanvas.h"
  66 #include "TFile.h"
  67 #include "TF1.h"
  68 #include "TVectorD.h"
  69 #include "TProfile.h"
  70 #include "TGraphErrors.h"
  71 #include "TCanvas.h"
  72
  73 #include "AliTPCclusterMI.h"
  74 #include "AliTPCseed.h"
  75 #include "AliESDVertex.h"
  76 #include "AliESDEvent.h"
  77 #include "AliESDfriend.h"
  78 #include "AliESDInputHandler.h"
  79 #include "AliAnalysisManager.h"
  80
  81 #include "AliTracker.h"
  82 #include "AliMagFMaps.h"
  83 #include "AliTPCCalROC.h"
  84
  85 #include "AliLog.h"
  86
  87 #include "AliTPCcalibTime.h"
  88
  89 #include "TTreeStream.h"
  90 #include "AliTPCTracklet.h"
  91 #include "TTimeStamp.h"
  92 #include "AliTPCcalibDB.h"
  93 #include "AliTPCcalibLaser.h"
  94
  95 ClassImp(AliTPCcalibTime)
  96
  97
  98 AliTPCcalibTime::AliTPCcalibTime()
  99   :AliTPCcalibBase(),
 100    fTriggerMask(0),
 101    fHistDeDxTgl(0),
 102    fHistDeDx(0),
 103    fHistVdrift(0),
 104    fIntegrationTimeDeDx(0),
 105    fIntegrationTimeVdrift(0),
 106    fLaser(0),       // pointer to laser calibration
 107    fDz(0),          // current delta z
 108    fdEdx(0),        // current dEdx
 109    fdEdxRatio(0),   // current dEdx ratio
 110    fTl(0),          // current tan(lambda)
 111    fCutMaxD(5),        // maximal distance in rfi ditection
 112    fCutMaxDz(20),        // maximal distance in rfi ditection
 113    fCutTheta(0.03),    // maximal distan theta
 114    fCutMinDir(-0.99)   // direction vector products
 115
 116 {
 117   AliInfo("Default Constructor");
 118 }
 119
 120
 121 AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, ULong64_t TriggerMask, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeDeDx, Int_t deltaIntegrationTimeVdrift)
 122   :AliTPCcalibBase(),
 123    fTriggerMask(0),
 124    fHistDeDxTgl(0),
 125    fHistDeDx(0),
 126    fHistVdrift(0),
 127    fIntegrationTimeDeDx(0),
 128    fIntegrationTimeVdrift(0),
 129    fLaser(0),       // pointer to laser calibration
 130    fDz(0),          // current delta z
 131    fdEdx(0),        // current dEdx
 132    fdEdxRatio(0),   // current dEdx ratio
 133    fTl(0),          // current tan(lambda)
 134    fCutMaxD(5),        // maximal distance in rfi ditection
 135    fCutMaxDz(20),        // maximal distance in rfi ditection
 136    fCutTheta(0.03),    // maximal distan theta
 137    fCutMinDir(-0.99)   // direction vector products
 138 {
 139
 140   SetName(name);
 141   SetTitle(title);
 142
 143   AliInfo("Non Default Constructor");
 144
 145   fTriggerMask = TriggerMask;
 146
 147   fIntegrationTimeDeDx = deltaIntegrationTimeDeDx;
 148   fIntegrationTimeVdrift = deltaIntegrationTimeVdrift;
 149
 150   Double_t deltaTime = EndTime - StartTime;
 151
 152   Int_t binsVdrift[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeVdrift), 100};
 153   Double_t xminVdrift[2] = {StartTime, -20};
 154   Double_t xmaxVdrift[2] = {EndTime, 20};
 155   fHistVdrift = new THnSparseF("HistVdrift","vDrift; time;#Delta z",2,binsVdrift,xminVdrift,xmaxVdrift);
 156
 157   Int_t binsDeDxTgl[3] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),30,100};
 158   Double_t xminDeDxTgl[3] = {StartTime,-1,0.7};
 159   Double_t xmaxDeDxTgl[3] = {EndTime,1,1.3};
 160   fHistDeDxTgl = new THnSparseF("HistDeDxTgl","dEdx vs tgl;time;tgl;dEdxUp/dEdxLow",3,binsDeDxTgl,xminDeDxTgl,xmaxDeDxTgl);
 161
 162   Int_t binsDeDx[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),100};
 163   Double_t xminDeDx[2] = {StartTime,1};
 164   Double_t xmaxDeDx[2] = {EndTime,100};
 165   fHistDeDx = new THnSparseF("HistDeDx","dEdx l;time;dEdx",2,binsDeDx,xminDeDx,xmaxDeDx);
 166
 167 }
 168
 169
 170
 171 AliTPCcalibTime::~AliTPCcalibTime(){
 172   //
 173   //
 174   //
 175 }
 176 void AliTPCcalibTime::ResetCurrent(){
 177   //
 178   // reset current values
 179   //
 180   fDz=0;          // current delta z
 181   fdEdx=0;        // current dEdx
 182   fdEdxRatio=0;   // current dEdx ratio
 183   fTl=0;          // current tan(lambda)
 184
 185 }
 186
 187
 188 void AliTPCcalibTime::Process(AliESDEvent *event) {
 189   //
 190   //
 191   //
 192   Int_t ntracks=event->GetNumberOfTracks();
 193   if (ntracks<2) return;
 194   ResetCurrent();
 195   //
 196   ProcessCosmic(event);
 197   if (fTrigger==16){
 198     if (!fLaser) fLaser =  new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
 199     fLaser->Process(event);
 200   }
 201   //
 202   // fill debug streamer
 203   if (fStreamLevel>0 && fDz!=0){
 204     TTreeSRedirector *cstream = GetDebugStreamer();
 205     if (cstream){
 206       TTimeStamp tstamp(fTime);
 207       Float_t valuePressure0  = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 208       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 209       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 210       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 211       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 212       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 213       TVectorD vdriftA, vdriftC,vdriftAC;
 214       if (fLaser && fTrigger==16) {
 215         if (fLaser->fFitAside)  vdriftA=*(fLaser->fFitAside);
 216         if (fLaser->fFitCside)  vdriftC=*(fLaser->fFitCside);
 217         if (fLaser->fFitACside) vdriftAC=*(fLaser->fFitACside);
 218       }
 219       (*cstream)<<"timeInfo"<<
 220         "run="<<fRun<<              //  run number
 221         "event="<<fEvent<<          //  event number
 222         "time="<<fTime<<            //  time stamp of event
 223         "trigger="<<fTrigger<<      //  trigger
 224         "mag="<<fMagF<<             //  magnetic field
 225         // Environment values
 226         "press0="<<valuePressure0<<
 227         "press1="<<valuePressure1<<
 228         "pt0="<<ptrelative0<<
 229         "pt1="<<ptrelative1<<
 230         "temp0="<<temp0<<
 231         "temp1="<<temp1<<
 232         //
 233         // accumulated values
 234         //
 235         "fDz="<<fDz<<          //! current delta z
 236         "fdEdx="<<fdEdx<<        //! current dEdx
 237         "fdEdxRatio="<<fdEdxRatio<<   //! current dEdx ratio
 238         "fTl="<<fTl<<          //! current tan(lambda)
 239         //
 240         //laser
 241         //
 242         "laserA.="<<&vdriftA<<
 243         "laserC.="<<&vdriftC<<
 244         "laserAC.="<<&vdriftAC<<
 245         "\n";
 246     }
 247   }
 248 }
 249
 250
 251
 252 void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event) {
 253
 254   if (!event) {
 255     Printf("ERROR: ESD not available");
 256     return;
 257   }
 258   if (event->GetTimeStamp() == 0 ) {
 259     Printf("no time stamp!");
 260     return;
 261   }
 262
 263   if (fTriggerMask != 0 && event->GetTriggerMask() != fTriggerMask) return;
 264
 265   UInt_t time = event->GetTimeStamp();
 266
 267   //
 268   // Find cosmic pairs
 269   //
 270   // Track0 is choosen in upper TPC part
 271   // Track1 is choosen in lower TPC part
 272   //
 273   Int_t ntracks=event->GetNumberOfTracks();
 274   if (ntracks==0) return;
 275   if (ntracks > 10) return; // temporary debug to remove LASER events
 276
 277
 278   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 279   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 280
 281   TObjArray  tpcSeeds(ntracks);
 282   Double_t vtxx[3]={0,0,0};
 283   Double_t svtxx[3]={0.000001,0.000001,100.};
 284   AliESDVertex vtx(vtxx,svtxx);
 285   //
 286   // track loop
 287   //
 288   for (Int_t i=0;i<ntracks;++i) {
 289    AliESDtrack *track = event->GetTrack(i);
 290
 291    const AliExternalTrackParam * trackIn = track->GetInnerParam();
 292    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 293    if (!trackIn) continue;
 294    if (!trackOut) continue;
 295
 296    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 297    TObject *calibObject;
 298    AliTPCseed *seed = 0;
 299    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 300      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 301    }
 302    if (seed) {
 303      tpcSeeds.AddAt(seed,i);
 304      if (track->GetTPCNcls() > 50) {
 305        Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 306        Double_t TPCsignal = seed->CookdEdxNorm(0.0,0.6,1,0,159,0x0,kTRUE,kTRUE);
 307        Double_t vecDeDx[2] = {time, TPCsignal};
 308        if (meanP > 12) {
 309          fdEdx = TPCsignal;
 310          fHistDeDx->Fill(vecDeDx);
 311        }
 312      }
 313    }
 314   }
 315
 316   if (ntracks<2) return;
 317   //
 318   // Find pairs
 319   //
 320   for (Int_t i=0;i<ntracks;++i) {
 321     AliESDtrack *track0 = event->GetTrack(i);
 322     // track0 - choosen upper part
 323     if (!track0) continue;
 324     if (!track0->GetOuterParam()) continue;
 325     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 326     Double_t d1[3];
 327     track0->GetDirection(d1);
 328     for (Int_t j=0;j<ntracks;++j) {
 329       if (i==j) continue;
 330        AliESDtrack *track1 = event->GetTrack(j);
 331        //track 1 lower part
 332        if (!track1) continue;
 333        if (!track1->GetOuterParam()) continue;
 334        if (track1->GetOuterParam()->GetAlpha()>0) continue;
 335        //
 336        Double_t d2[3];
 337        track1->GetDirection(d2);
 338
 339        AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 340        AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 341        if (! seed0) continue;
 342        if (! seed1) continue;
 343        Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
 344        Float_t dist0  = track0->GetLinearD(0,0);
 345        Float_t dist1  = track1->GetLinearD(0,0);
 346        //
 347        // conservative cuts - convergence to be guarantied
 348        // applying before track propagation
 349        if (TMath::Abs(dist0+dist1)>fCutMaxD) continue;   // distance to the 0,0
 350        if (dir>fCutMinDir) continue;               // direction vector product
 351        Float_t bz = AliTracker::GetBz();
 352        Float_t dvertex0[2];   //distance to 0,0
 353        Float_t dvertex1[2];   //distance to 0,0
 354        track0->GetDZ(0,0,0,bz,dvertex0);
 355        track1->GetDZ(0,0,0,bz,dvertex1);
 356        if (TMath::Abs(dvertex0[1])>250) continue;
 357        if (TMath::Abs(dvertex1[1])>250) continue;
 358        //
 359        //
 360        //
 361        Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
 362        AliExternalTrackParam param0(*track0);
 363        AliExternalTrackParam param1(*track1);
 364        //
 365        // Propagate using Magnetic field and correct fo material budget
 366        //
 367        AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 368        AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 369        //
 370        // Propagate rest to the 0,0 DCA - z should be ignored
 371        //
 372        //Bool_t b0 = ;
 373        param0.PropagateToDCA(&vtx,bz,1000);
 374        //Bool_t b1 =
 375        param1.PropagateToDCA(&vtx,bz,1000);
 376        //
 377        param0.GetDZ(0,0,0,bz,dvertex0);
 378        param1.GetDZ(0,0,0,bz,dvertex1);
 379        //
 380        Double_t xyz0[3];//,pxyz0[3];
 381        Double_t xyz1[3];//,pxyz1[3];
 382        param0.GetXYZ(xyz0);
 383        param1.GetXYZ(xyz1);
 384        Bool_t isPair = IsPair(&param0,&param1);
 385
 386        Double_t z0 = track0->GetOuterParam()->GetZ();
 387        Double_t z1 = track1->GetOuterParam()->GetZ();
 388
 389        Double_t z0inner = track0->GetInnerParam()->GetZ();
 390        Double_t z1inner = track1->GetInnerParam()->GetZ();
 391
 392        if (isPair && z0>0 && z0inner>0 && z1<0 && z1inner<0) {
 393          Double_t vecVdrift[2] = {time, param0.GetZ() - param1.GetZ()};
 394
 395          if (track0->GetTPCNcls() > 80) {
 396            fHistVdrift->Fill(vecVdrift);
 397            fDz = param0.GetZ() - param1.GetZ();
 398          }
 399        }
 400        if (isPair && z0 > 0 && z1 > 0) {
 401          if (track1->GetTPCNcls()> 110 && track0->GetTPCNcls()> 110 && seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE) > 0) {
 402            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);
 403            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)};
 404            fHistDeDxTgl->Fill(vecDeDxTgl);
 405            fdEdxRatio = dedxratio;
 406            fTl = track0->GetTgl() ;
 407          }
 408        }
 409
 410     } // end 2nd order loop
 411   } // end 1st order loop
 412
 413 }
 414
 415
 416 void AliTPCcalibTime::Analyze() {
 417   //
 418   //
 419   //
 420   TH2D * hVdrift = GetHistVdrift()->Projection(1,0);
 421   hVdrift->Draw();
 422 }
 423
 424
 425 Long64_t AliTPCcalibTime::Merge(TCollection *li) {
 426
 427   TIterator* iter = li->MakeIterator();
 428   AliTPCcalibTime* cal = 0;
 429
 430   while ((cal = (AliTPCcalibTime*)iter->Next())) {
 431     if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
 432       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 433       return -1;
 434     }
 435
 436     // add histograms here...
 437     fHistDeDxTgl->Add(cal->GetHistDeDxVsTgl());
 438     fHistVdrift->Add(cal->GetHistVdrift());
 439     fHistDeDx->Add(cal->GetHistDeDx());
 440
 441   }
 442
 443   return 0;
 444
 445 }
 446
 447
 448
 449 Bool_t  AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 450   //
 451   //
 452   /*
 453   // 0. Same direction - OPOSITE  - cutDir +cutT
 454   TCut cutDir("cutDir","dir<-0.99")
 455   // 1.
 456   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 457   //
 458   // 2. The same rphi
 459   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 460   //
 461   //
 462   //
 463   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 464   // 1/Pt diff cut
 465   */
 466   const Double_t *p0 = tr0->GetParameter();
 467   const Double_t *p1 = tr1->GetParameter();
 468   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 469   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 470   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz)  return kFALSE;
 471   Double_t d0[3], d1[3];
 472   tr0->GetDirection(d0);
 473   tr1->GetDirection(d1);
 474   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 475   //
 476   return kTRUE;
 477 }