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
  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 */
  47
  48
  49 #include "Riostream.h"
  50 #include "TChain.h"
  51 #include "TTree.h"
  52 #include "TH1F.h"
  53 #include "TH2F.h"
  54 #include "TH3F.h"
  55 #include "THnSparse.h"
  56 #include "TList.h"
  57 #include "TMath.h"
  58 #include "TCanvas.h"
  59 #include "TFile.h"
  60 #include "TF1.h"
  61 #include "TVectorD.h"
  62 #include "TProfile.h"
  63 #include "TGraphErrors.h"
  64 #include "TCanvas.h"
  65
  66 #include "AliTPCclusterMI.h"
  67 #include "AliTPCseed.h"
  68 #include "AliESDVertex.h"
  69 #include "AliESDEvent.h"
  70 #include "AliESDfriend.h"
  71 #include "AliESDInputHandler.h"
  72 #include "AliAnalysisManager.h"
  73
  74 #include "AliTracker.h"
  75 #include "AliMagFMaps.h"
  76 #include "AliTPCCalROC.h"
  77
  78 #include "AliLog.h"
  79
  80 #include "AliTPCcalibTime.h"
  81
  82 #include "TTreeStream.h"
  83 #include "AliTPCTracklet.h"
  84
  85 ClassImp(AliTPCcalibTime)
  86
  87
  88 AliTPCcalibTime::AliTPCcalibTime()
  89   :AliTPCcalibBase(),
  90    fHistDeDxTgl(0),
  91    fHistDeDx(0),
  92    fHistVdrift(0),
  93    fIntegrationTimeDeDx(0),
  94    fIntegrationTimeVdrift(0),
  95    fCutMaxD(5),        // maximal distance in rfi ditection
  96    fCutTheta(0.03),    // maximal distan theta
  97    fCutMinDir(-0.99)   // direction vector products
  98
  99 {
 100   AliInfo("Default Constructor");
 101 }
 102
 103
 104 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)
 105   :AliTPCcalibBase(),
 106    fTriggerMask(0),
 107    fHistDeDxTgl(0),
 108    fHistDeDx(0),
 109    fHistVdrift(0),
 110    fIntegrationTimeDeDx(0),
 111    fIntegrationTimeVdrift(0),
 112    fCutMaxD(5),        // maximal distance in rfi ditection
 113    fCutTheta(0.03),    // maximal distan theta
 114    fCutMinDir(-0.99)   // direction vector products
 115 {
 116
 117   SetName(name);
 118   SetTitle(title);
 119
 120   AliInfo("Non Default Constructor");
 121
 122   fTriggerMask = TriggerMask;
 123
 124   fIntegrationTimeDeDx = deltaIntegrationTimeDeDx;
 125   fIntegrationTimeVdrift = deltaIntegrationTimeVdrift;
 126
 127   Double_t deltaTime = EndTime - StartTime;
 128
 129   Int_t binsVdrift[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeVdrift), 100};
 130   Double_t xminVdrift[2] = {StartTime, -20};
 131   Double_t xmaxVdrift[2] = {EndTime, 20};
 132   fHistVdrift = new THnSparseF("HistVdrift","vDrift; time;#Delta z",2,binsVdrift,xminVdrift,xmaxVdrift);
 133
 134   Int_t binsDeDxTgl[3] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),30,100};
 135   Double_t xminDeDxTgl[3] = {StartTime,-1,0.7};
 136   Double_t xmaxDeDxTgl[3] = {EndTime,1,1.3};
 137   fHistDeDxTgl = new THnSparseF("HistDeDxTgl","dEdx vs tgl;time;tgl;dEdxUp/dEdxLow",3,binsDeDxTgl,xminDeDxTgl,xmaxDeDxTgl);
 138
 139   Int_t binsDeDx[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),100};
 140   Double_t xminDeDx[2] = {StartTime,1};
 141   Double_t xmaxDeDx[2] = {EndTime,100};
 142   fHistDeDx = new THnSparseF("HistDeDx","dEdx l;time;dEdx",2,binsDeDx,xminDeDx,xmaxDeDx);
 143
 144 }
 145
 146
 147
 148 AliTPCcalibTime::~AliTPCcalibTime(){
 149   //
 150   //
 151   //
 152 }
 153
 154 void AliTPCcalibTime::Process(AliESDEvent *event) {
 155   //
 156   //
 157   //
 158
 159   ProcessCosmic(event);
 160
 161 }
 162
 163
 164
 165 void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event) {
 166
 167   if (!event) {
 168     Printf("ERROR: ESD not available");
 169     return;
 170   }
 171   if (event->GetTimeStamp() == 0 ) {
 172     Printf("no time stamp!");
 173     return;
 174   }
 175
 176   if (fTriggerMask != 0 && event->GetTriggerMask() != fTriggerMask) return;
 177
 178   UInt_t time = event->GetTimeStamp();
 179
 180   //
 181   // Find cosmic pairs
 182   //
 183   // Track0 is choosen in upper TPC part
 184   // Track1 is choosen in lower TPC part
 185   //
 186   Int_t ntracks=event->GetNumberOfTracks();
 187   if (ntracks==0) return;
 188   if (ntracks > 10) return; // temporary debug to remove LASER events
 189
 190
 191   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 192   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 193
 194   TObjArray  tpcSeeds(ntracks);
 195   Double_t vtxx[3]={0,0,0};
 196   Double_t svtxx[3]={0.000001,0.000001,100.};
 197   AliESDVertex vtx(vtxx,svtxx);
 198   //
 199   // track loop
 200   //
 201   for (Int_t i=0;i<ntracks;++i) {
 202    AliESDtrack *track = event->GetTrack(i);
 203
 204    const AliExternalTrackParam * trackIn = track->GetInnerParam();
 205    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 206    if (!trackIn) continue;
 207    if (!trackOut) continue;
 208
 209    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 210    TObject *calibObject;
 211    AliTPCseed *seed = 0;
 212    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 213      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 214    }
 215    if (seed) {
 216      tpcSeeds.AddAt(seed,i);
 217      if (track->GetTPCNcls() > 50) {
 218        Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 219        Double_t TPCsignal = seed->CookdEdxNorm(0.0,0.6,1,0,159,0x0,kTRUE,kTRUE);
 220        Double_t vecDeDx[2] = {time, TPCsignal};
 221        if (meanP > 12) fHistDeDx->Fill(vecDeDx);
 222      }
 223    }
 224   }
 225
 226   if (ntracks<2) return;
 227   //
 228   // Find pairs
 229   //
 230   for (Int_t i=0;i<ntracks;++i) {
 231     AliESDtrack *track0 = event->GetTrack(i);
 232     // track0 - choosen upper part
 233     if (!track0) continue;
 234     if (!track0->GetOuterParam()) continue;
 235     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 236     Double_t d1[3];
 237     track0->GetDirection(d1);
 238     for (Int_t j=0;j<ntracks;++j) {
 239       if (i==j) continue;
 240        AliESDtrack *track1 = event->GetTrack(j);
 241        //track 1 lower part
 242        if (!track1) continue;
 243        if (!track1->GetOuterParam()) continue;
 244        if (track1->GetOuterParam()->GetAlpha()>0) continue;
 245        //
 246        Double_t d2[3];
 247        track1->GetDirection(d2);
 248
 249        AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 250        AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 251        if (! seed0) continue;
 252        if (! seed1) continue;
 253        Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
 254        Float_t d0  = track0->GetLinearD(0,0);
 255        Float_t d1  = track1->GetLinearD(0,0);
 256        //
 257        // conservative cuts - convergence to be guarantied
 258        // applying before track propagation
 259        if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 260        if (dir>fCutMinDir) continue;               // direction vector product
 261        Float_t bz = AliTracker::GetBz();
 262        Float_t dvertex0[2];   //distance to 0,0
 263        Float_t dvertex1[2];   //distance to 0,0
 264        track0->GetDZ(0,0,0,bz,dvertex0);
 265        track1->GetDZ(0,0,0,bz,dvertex1);
 266        if (TMath::Abs(dvertex0[1])>250) continue;
 267        if (TMath::Abs(dvertex1[1])>250) continue;
 268        //
 269        //
 270        //
 271        Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 272        AliExternalTrackParam param0(*track0);
 273        AliExternalTrackParam param1(*track1);
 274        //
 275        // Propagate using Magnetic field and correct fo material budget
 276        //
 277        AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 278        AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 279        //
 280        // Propagate rest to the 0,0 DCA - z should be ignored
 281        //
 282        //Bool_t b0 = ;
 283        param0.PropagateToDCA(&vtx,bz,1000);
 284        //Bool_t b1 =
 285        param1.PropagateToDCA(&vtx,bz,1000);
 286        //
 287        param0.GetDZ(0,0,0,bz,dvertex0);
 288        param1.GetDZ(0,0,0,bz,dvertex1);
 289        //
 290        Double_t xyz0[3];//,pxyz0[3];
 291        Double_t xyz1[3];//,pxyz1[3];
 292        param0.GetXYZ(xyz0);
 293        param1.GetXYZ(xyz1);
 294        Bool_t isPair = IsPair(&param0,&param1);
 295
 296        Double_t z0 = track0->GetOuterParam()->GetZ();
 297        Double_t z1 = track1->GetOuterParam()->GetZ();
 298
 299        Double_t z0inner = track0->GetInnerParam()->GetZ();
 300        Double_t z1inner = track1->GetInnerParam()->GetZ();
 301
 302        if (isPair && z0>0 && z0inner>0 && z1<0 && z1inner<0) {
 303          Double_t vecVdrift[2] = {time, param0.GetZ() - param1.GetZ()};
 304
 305          if (track0->GetTPCNcls() > 80) fHistVdrift->Fill(vecVdrift);
 306        }
 307        if (isPair && z0 > 0 && z1 > 0) {
 308          if (track1->GetTPCNcls()> 110 && track0->GetTPCNcls()> 110 && seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE) > 0) {
 309            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)};
 310            fHistDeDxTgl->Fill(vecDeDxTgl);
 311          }
 312        }
 313
 314     } // end 2nd order loop
 315   } // end 1st order loop
 316
 317 }
 318
 319
 320 void AliTPCcalibTime::Analyze() {
 321   //
 322   //
 323   //
 324   TH2D * hVdrift = GetHistVdrift()->Projection(1,0);
 325   hVdrift->Draw();
 326 }
 327
 328
 329 Long64_t AliTPCcalibTime::Merge(TCollection *li) {
 330
 331   TIterator* iter = li->MakeIterator();
 332   AliTPCcalibTime* cal = 0;
 333
 334   while ((cal = (AliTPCcalibTime*)iter->Next())) {
 335     if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
 336       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 337       return -1;
 338     }
 339
 340     // add histograms here...
 341     fHistDeDxTgl->Add(cal->GetHistDeDxVsTgl());
 342     fHistVdrift->Add(cal->GetHistVdrift());
 343     fHistDeDx->Add(cal->GetHistDeDx());
 344
 345   }
 346
 347   return 0;
 348
 349 }
 350
 351
 352
 353 Bool_t  AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 354   //
 355   //
 356   /*
 357   // 0. Same direction - OPOSITE  - cutDir +cutT
 358   TCut cutDir("cutDir","dir<-0.99")
 359   // 1.
 360   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 361   //
 362   // 2. The same rphi
 363   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 364   //
 365   //
 366   //
 367   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 368   // 1/Pt diff cut
 369   */
 370   const Double_t *p0 = tr0->GetParameter();
 371   const Double_t *p1 = tr1->GetParameter();
 372   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 373   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 374   Double_t d0[3], d1[3];
 375   tr0->GetDirection(d0);
 376   tr1->GetDirection(d1);
 377   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 378   //
 379   return kTRUE;
 380 }