TPC/AliTPCcalibCosmic.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 #include "Riostream.h"
  17 #include "TChain.h"
  18 #include "TTree.h"
  19 #include "TH1F.h"
  20 #include "TH2F.h"
  21 #include "TList.h"
  22 #include "TMath.h"
  23 #include "TCanvas.h"
  24 #include "TFile.h"
  25
  26 #include "AliTPCseed.h"
  27 #include "AliESDVertex.h"
  28 #include "AliESDEvent.h"
  29 #include "AliESDfriend.h"
  30 #include "AliESDInputHandler.h"
  31
  32 #include "AliTracker.h"
  33 #include "AliMagFMaps.h"
  34
  35 #include "AliLog.h"
  36
  37 #include "AliTPCcalibCosmic.h"
  38
  39 #include "TTreeStream.h"
  40 #include "AliTPCTracklet.h"
  41
  42 ClassImp(AliTPCcalibCosmic)
  43
  44
  45 AliTPCcalibCosmic::AliTPCcalibCosmic()
  46   :AliTPCcalibBase(),
  47    fGainMap(0),
  48    fHistNTracks(0),
  49    fClusters(0),
  50    fModules(0),
  51    fHistPt(0),
  52    fPtResolution(0),
  53    fDeDx(0),
  54    fCutMaxD(5),        // maximal distance in rfi ditection
  55    fCutTheta(0.03),    // maximal distan theta
  56    fCutMinDir(-0.99)   // direction vector products
  57 {
  58   AliInfo("Defualt Constructor");
  59   TFile f("/u/miranov/calibKr.root");
  60   AliTPCCalPad *gainMap =  (AliTPCCalPad *)f.Get("spectrMean");
  61   gainMap->Multiply(1/gainMap->GetMedian());
  62   fGainMap =gainMap;
  63 }
  64
  65
  66 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
  67   :AliTPCcalibBase(),
  68    fGainMap(0),
  69    fHistNTracks(0),
  70    fClusters(0),
  71    fModules(0),
  72    fHistPt(0),
  73    fPtResolution(0),
  74    fDeDx(0),
  75    fCutMaxD(5),        // maximal distance in rfi ditection
  76    fCutTheta(0.03),    // maximal distan theta
  77    fCutMinDir(-0.99)   // direction vector products
  78 {
  79   SetName(name);
  80   SetTitle(title);
  81   AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
  82   AliTracker::SetFieldMap(field, kTRUE);
  83   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event",501,-0.5,500.5);
  84   fClusters = new TH1F("signal","Number of Clusters per track",160,0,160);
  85   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-1200,1200,600,-1000,1000);
  86   fHistPt = new TH1F("Pt","Pt distribution",2000,0,50);
  87   fPtResolution = new TH1F("PtResolution","Pt resolution",100,-50,50);
  88   fDeDx = new TH2F("DeDx","dEdx",500,0.01,20.,500,0.,500);
  89   BinLogX(fDeDx);
  90   AliInfo("Non Default Constructor");
  91   //
  92   TFile f("/u/miranov/calibKr.root");
  93   AliTPCCalPad *gainMap = (AliTPCCalPad *)f.Get("spectrMean");
  94   gainMap->Multiply(1/gainMap->GetMedian());
  95   fGainMap =gainMap;
  96 }
  97
  98 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
  99   //
 100   //
 101   //
 102 }
 103
 104
 105
 106
 107
 108 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
 109   //
 110   //
 111   //
 112   if (!event) {
 113     Printf("ERROR: ESD not available");
 114     return;
 115   }
 116   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 117   if (!ESDfriend) {
 118    Printf("ERROR: ESDfriend not available");
 119    return;
 120   }
 121   FindPairs(event);
 122
 123   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 124   Int_t ntracks=event->GetNumberOfTracks();
 125   fHistNTracks->Fill(ntracks);
 126   TObjArray  tpcSeeds(ntracks);
 127   if (ntracks==0) return;
 128   //
 129   //track loop
 130   //
 131   for (Int_t i=0;i<ntracks;++i) {
 132    AliESDtrack *track = event->GetTrack(i);
 133    fClusters->Fill(track->GetTPCNcls());
 134    AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
 135
 136    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 137    TObject *calibObject;
 138    AliTPCseed *seed = 0;
 139    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 140      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 141    }
 142    if (seed) tpcSeeds.AddAt(seed,i);
 143    if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0,0,159,fGainMap));
 144   }
 145   if (ntracks<2) return;
 146
 147
 148   // dE/dx,pt and ACORDE study --> studies which need the pair selection
 149   for (Int_t i=0;i<ntracks;++i) {
 150     AliESDtrack *track1 = event->GetTrack(i);
 151
 152     Double_t d1[3];
 153     track1->GetDirection(d1);
 154
 155     for (Int_t j=i+1;j<ntracks;++j) {
 156      AliESDtrack *track2 = event->GetTrack(j);
 157      Double_t d2[3];
 158      track2->GetDirection(d2);
 159
 160      if (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2] < -0.999) {
 161
 162       /*___________________________________ Pt resolution ________________________________________*/
 163       if (track1->Pt() != 0 && track1->GetTPCNcls() > 80 && track2->GetTPCNcls() > 80) {
 164        Double_t res = (track1->Pt() - track2->Pt());
 165        res = res/(2*(track1->Pt() + track2->Pt()));
 166        fPtResolution->Fill(100*res);
 167       }
 168
 169       /*_______________________________ Propagation to ACORDE ___________________________________*/
 170       const Double_t AcordePlane = 850.; //distance of the central Acorde detectors to the beam line at y =0
 171       const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
 172
 173       if (d1[1] > 0 && d2[1] < 0 && track1->GetTPCNcls() > 50) {
 174        Double_t r[3];
 175        track1->GetXYZ(r);
 176        Double_t x,z;
 177        z = r[2] + (d1[2]/d1[1])*(AcordePlane - r[1]);
 178        x = r[0] + (d1[0]/d1[1])*(AcordePlane - r[1]);
 179
 180        if (x > roof) {
 181         x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
 182         z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
 183        }
 184        if (x < -roof) {
 185         x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
 186         z = z -  TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
 187        }
 188
 189        fModules->Fill(z, x);
 190       }
 191
 192       if (d2[1] > 0 && d1[1] < 0 && track2->GetTPCNcls() > 50) {
 193        Double_t r[3];
 194        track2->GetXYZ(r);
 195        Double_t x,z;
 196        z = r[2] + (d2[2]/d2[1])*(AcordePlane - r[1]);
 197        x = r[0] + (d2[0]/d2[1])*(AcordePlane - r[1]);
 198
 199        if (x > roof) {
 200         x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
 201         z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
 202        }
 203        if (x < -roof) {
 204         x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
 205         z = z -  TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
 206        }
 207
 208        fModules->Fill(z, x);
 209       }
 210
 211   //     AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
 212 //       AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
 213 //       delete trackOut;
 214
 215
 216
 217
 218
 219       break;
 220      }
 221     }
 222   }
 223
 224
 225
 226
 227 }
 228
 229 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
 230   //
 231   // Find cosmic pairs
 232   //
 233   // Track0 is choosen in upper TPC part
 234   // Track1 is choosen in lower TPC part
 235   //
 236   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 237   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 238   Int_t ntracks=event->GetNumberOfTracks();
 239   TObjArray  tpcSeeds(ntracks);
 240   if (ntracks==0) return;
 241   Double_t vtxx[3]={0,0,0};
 242   Double_t svtxx[3]={0.000001,0.000001,100.};
 243   AliESDVertex vtx(vtxx,svtxx);
 244   //
 245   //track loop
 246   //
 247   for (Int_t i=0;i<ntracks;++i) {
 248    AliESDtrack *track = event->GetTrack(i);
 249    fClusters->Fill(track->GetTPCNcls());
 250    if (!track->GetInnerParam()) continue;
 251    AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
 252
 253    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 254    TObject *calibObject;
 255    AliTPCseed *seed = 0;
 256    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 257      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 258    }
 259    if (seed) tpcSeeds.AddAt(seed,i);
 260    if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0,0,159,fGainMap));
 261   }
 262   if (ntracks<2) return;
 263   //
 264   // Find pairs
 265   //
 266   for (Int_t i=0;i<ntracks;++i) {
 267     AliESDtrack *track0 = event->GetTrack(i);
 268     // track0 - choosen upper part
 269     if (!track0) continue;
 270     if (!track0->GetOuterParam()) continue;
 271     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 272     Double_t d1[3];
 273     track0->GetDirection(d1);
 274     for (Int_t j=0;j<ntracks;++j) {
 275       if (i==j) continue;
 276       AliESDtrack *track1 = event->GetTrack(j);
 277       //track 1 lower part
 278       if (!track1) continue;
 279       if (!track1->GetOuterParam()) continue;
 280       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 281       //
 282       Double_t d2[3];
 283       track1->GetDirection(d2);
 284       printf("My stream level=%d\n",fStreamLevel);
 285       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 286       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 287       if (! seed0) continue;
 288       if (! seed1) continue;
 289       Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
 290       Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
 291       //
 292       Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
 293       Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
 294       //
 295       Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
 296       Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
 297       //
 298       Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
 299       Float_t d0  = track0->GetLinearD(0,0);
 300       Float_t d1  = track1->GetLinearD(0,0);
 301       //
 302       // conservative cuts - convergence to be guarantied
 303       // applying before track propagation
 304       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 305       if (dir>fCutMinDir) continue;               // direction vector product
 306       Float_t bz = AliTracker::GetBz();
 307       Float_t dvertex0[2];   //distance to 0,0
 308       Float_t dvertex1[2];   //distance to 0,0
 309       track0->GetDZ(0,0,0,bz,dvertex0);
 310       track1->GetDZ(0,0,0,bz,dvertex1);
 311       if (TMath::Abs(dvertex0[1])>250) continue;
 312       if (TMath::Abs(dvertex1[1])>250) continue;
 313       //
 314       //
 315       //
 316       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 317       AliExternalTrackParam param0(*track0);
 318       AliExternalTrackParam param1(*track1);
 319       //
 320       // Propagate using Magnetic field and correct fo material budget
 321       //
 322       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 323       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 324       //
 325       // Propagate rest to the 0,0 DCA - z should be ignored
 326       //
 327       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
 328       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
 329       //
 330       param0.GetDZ(0,0,0,bz,dvertex0);
 331       param1.GetDZ(0,0,0,bz,dvertex1);
 332       //
 333       Double_t xyz0[3];//,pxyz0[3];
 334       Double_t xyz1[3];//,pxyz1[3];
 335       param0.GetXYZ(xyz0);
 336       param1.GetXYZ(xyz1);
 337       Bool_t isPair = IsPair(&param0,&param1);
 338       //
 339       if (fStreamLevel>0){
 340         TTreeSRedirector * cstream =  GetDebugStreamer();
 341         printf("My stream=%p\n",(void*)cstream);
 342         if (cstream) {
 343           (*cstream) << "Track0" <<
 344             "dir="<<dir<<               //  direction
 345             "OK="<<isPair<<             // will be accepted
 346             "b0="<<b0<<                 //  propagate status
 347             "b1="<<b1<<                 //  propagate status
 348             "Orig0.=" << track0 <<      //  original track  0
 349             "Orig1.=" << track1 <<      //  original track  1
 350             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
 351             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0
 352             "v00="<<dvertex0[0]<<       //  distance using kalman
 353             "v01="<<dvertex0[1]<<       //
 354             "v10="<<dvertex1[0]<<       //
 355             "v11="<<dvertex1[1]<<       //
 356             "d0="<<d0<<                 //  linear distance to 0,0
 357             "d1="<<d1<<                 //  linear distance to 0,0
 358             //
 359             "x00="<<xyz0[0]<<
 360             "x01="<<xyz0[1]<<
 361             "x02="<<xyz0[2]<<
 362             //
 363             "x10="<<xyz1[0]<<
 364             "x11="<<xyz1[1]<<
 365             "x12="<<xyz1[2]<<
 366             //
 367             "Seed0.=" << track0 <<      //  original seed 0
 368             "Seed1.=" << track1 <<      //  original seed 1
 369             //
 370             "dedx0="<<dedx0<<           //  dedx0 - all
 371             "dedx1="<<dedx1<<           //  dedx1 - all
 372             //
 373             "dedx0I="<<dedx0I<<           //  dedx0 - inner
 374             "dedx1I="<<dedx1I<<           //  dedx1 - inner
 375             //
 376             "dedx0O="<<dedx0O<<           //  dedx0 - outer
 377             "dedx1O="<<dedx1O<<           //  dedx1 -outer
 378             "\n";
 379         }
 380       }
 381     }
 382   }
 383 }
 384
 385
 386
 387
 388 Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
 389
 390   TIterator* iter = li->MakeIterator();
 391   AliTPCcalibCosmic* cal = 0;
 392
 393   while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
 394     if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
 395       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 396       return -1;
 397     }
 398
 399     fHistNTracks->Add(cal->fHistNTracks);
 400     fClusters->Add(cal->fClusters);
 401     fModules->Add(cal->fModules);
 402     fHistPt->Add(cal->fHistPt);
 403     fPtResolution->Add(cal->fPtResolution);
 404     fDeDx->Add(cal->fDeDx);
 405
 406   }
 407
 408   return 0;
 409
 410 }
 411
 412 Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 413   //
 414   //
 415   /*
 416   // 0. Same direction - OPOSITE  - cutDir +cutT
 417   TCut cutDir("cutDir","dir<-0.99")
 418   // 1.
 419   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 420   //
 421   // 2. The same rphi
 422   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 423   //
 424   //
 425   //
 426   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 427   // 1/Pt diff cut
 428   */
 429   const Double_t *p0 = tr0->GetParameter();
 430   const Double_t *p1 = tr1->GetParameter();
 431   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 432   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 433   Double_t d0[3], d1[3];
 434   tr0->GetDirection(d0);
 435   tr1->GetDirection(d1);
 436   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 437   //
 438   return kTRUE;
 439 }
 440
 441
 442
 443 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
 444
 445   // Method for the correct logarithmic binning of histograms
 446
 447   TAxis *axis = h->GetXaxis();
 448   int bins = axis->GetNbins();
 449
 450   Double_t from = axis->GetXmin();
 451   Double_t to = axis->GetXmax();
 452   Double_t *new_bins = new Double_t[bins + 1];
 453
 454   new_bins[0] = from;
 455   Double_t factor = pow(to/from, 1./bins);
 456
 457   for (int i = 1; i <= bins; i++) {
 458    new_bins[i] = factor * new_bins[i-1];
 459   }
 460   axis->Set(bins, new_bins);
 461   delete new_bins;
 462
 463 }
 464
 465
 466
 467 /*
 468
 469
 470
 471
 472 void AliTPCcalibCosmic::dEdxCorrection(){
 473   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03");
 474   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5");
 475   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 476   TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
 477   TCut cutA=cutT+cutD+cutPt+cutN;
 478
 479   TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
 480   chain->Lookup();
 481
 482   .x ~/rootlogon.C
 483    gSystem->Load("libSTAT.so");
 484
 485   Double_t chi2=0;
 486   Int_t    npoints=0;
 487   TVectorD fitParam;
 488   TMatrixD covMatrix;
 489
 490   chain->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
 491
 492   TString strFit;
 493   strFit+="(Tr0.fP[1]/250)++";
 494   strFit+="(Tr0.fP[1]/250)^2++";
 495   strFit+="(Tr0.fP[3])++";
 496   strFit+="(Tr0.fP[3])^2++";
 497
 498   TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),cutA, chi2,npoints,fitParam,covMatrix)
 499
 500 strFit+="(Tr0.fP[1]/250)++";
 501 strFit+="(Tr0.fP[1]/250)^2++";
 502 strFit+="(Tr0.fP[3])++";
 503 strFit+="(Tr0.fP[3])^2++";
 504 strFit+="(Tr0.fP[1]/250)^2*Tr0.fP[3]++";
 505 strFit+="(Tr0.fP[1]/250)^2*Tr0.fP[3]^2++";
 506 //
 507
 508 strFit+="sign(Tr0.fP[1])++"
 509 strFit+="sign(Tr0.fP[1])*(1-abs(Tr0.fP[1]/250))"
 510
 511 TString * thetaParam = TStatToolkit::FitPlane(chain,"Tr0.fP[3]+Tr1.fP[3]", strFit.Data(),cutA, chi2,npoints,fitParam,covMatrix)
 512
 513
 514
 515 (-0.009263+(Tr0.fP[1]/250)*(-0.009693)+(Tr0.fP[1]/250)^2*(0.009062)+(Tr0.fP[3])*(0.002256)+(Tr0.fP[3])^2*(-0.052775)+(Tr0.fP[1]/250)^2*Tr0.fP[3]*(-0.020627)+(Tr0.fP[1]/250)^2*Tr0.fP[3]^2*(0.049030))
 516
 517 */
 518
 519