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    fHistNTracks(0),
  48    fClusters(0),
  49    fModules(0),
  50    fHistPt(0),
  51    fPtResolution(0),
  52    fDeDx(0),
  53    fCutMaxD(5),        // maximal distance in rfi ditection
  54    fCutTheta(0.03),    // maximal distan theta
  55    fCutMinDir(-0.99)   // direction vector products
  56 {
  57   AliInfo("Defualt Constructor");
  58 }
  59
  60
  61 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
  62   :AliTPCcalibBase(),
  63    fHistNTracks(0),
  64    fClusters(0),
  65    fModules(0),
  66    fHistPt(0),
  67    fPtResolution(0),
  68    fDeDx(0),
  69    fCutMaxD(5),        // maximal distance in rfi ditection
  70    fCutTheta(0.03),    // maximal distan theta
  71    fCutMinDir(-0.99)   // direction vector products
  72 {
  73   SetName(name);
  74   SetTitle(title);
  75   AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
  76   AliTracker::SetFieldMap(field, kTRUE);
  77   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event",501,-0.5,500.5);
  78   fClusters = new TH1F("signal","Number of Clusters per track",160,0,160);
  79   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-1200,1200,600,-1000,1000);
  80   fHistPt = new TH1F("Pt","Pt distribution",2000,0,50);
  81   fPtResolution = new TH1F("PtResolution","Pt resolution",100,-50,50);
  82   fDeDx = new TH2F("DeDx","dEdx",500,0.01,20.,500,0.,500);
  83   BinLogX(fDeDx);
  84   AliInfo("Non Default Constructor");
  85 }
  86
  87 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
  88   //
  89   //
  90   //
  91 }
  92
  93
  94
  95
  96
  97 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
  98   //
  99   //
 100   //
 101   if (!event) {
 102     Printf("ERROR: ESD not available");
 103     return;
 104   }
 105   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 106   if (!ESDfriend) {
 107    Printf("ERROR: ESDfriend not available");
 108    return;
 109   }
 110   FindPairs(event);
 111
 112   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 113   Int_t ntracks=event->GetNumberOfTracks();
 114   fHistNTracks->Fill(ntracks);
 115   TObjArray  tpcSeeds(ntracks);
 116   if (ntracks==0) return;
 117   //
 118   //track loop
 119   //
 120   for (Int_t i=0;i<ntracks;++i) {
 121    AliESDtrack *track = event->GetTrack(i);
 122    fClusters->Fill(track->GetTPCNcls());
 123    AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
 124
 125    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 126    TObject *calibObject;
 127    AliTPCseed *seed = 0;
 128    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 129      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 130    }
 131    if (seed) tpcSeeds.AddAt(seed,i);
 132    if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
 133   }
 134   if (ntracks<2) return;
 135
 136
 137   // dE/dx,pt and ACORDE study --> studies which need the pair selection
 138   for (Int_t i=0;i<ntracks;++i) {
 139     AliESDtrack *track1 = event->GetTrack(i);
 140
 141     Double_t d1[3];
 142     track1->GetDirection(d1);
 143
 144     for (Int_t j=i+1;j<ntracks;++j) {
 145      AliESDtrack *track2 = event->GetTrack(j);
 146      Double_t d2[3];
 147      track2->GetDirection(d2);
 148
 149      if (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2] < -0.999) {
 150
 151       /*___________________________________ Pt resolution ________________________________________*/
 152       if (track1->Pt() != 0 && track1->GetTPCNcls() > 80 && track2->GetTPCNcls() > 80) {
 153        Double_t res = (track1->Pt() - track2->Pt());
 154        res = res/(2*(track1->Pt() + track2->Pt()));
 155        fPtResolution->Fill(100*res);
 156       }
 157
 158       /*_______________________________ Propagation to ACORDE ___________________________________*/
 159       const Double_t AcordePlane = 850.; //distance of the central Acorde detectors to the beam line at y =0
 160       const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
 161
 162       if (d1[1] > 0 && d2[1] < 0 && track1->GetTPCNcls() > 50) {
 163        Double_t r[3];
 164        track1->GetXYZ(r);
 165        Double_t x,z;
 166        z = r[2] + (d1[2]/d1[1])*(AcordePlane - r[1]);
 167        x = r[0] + (d1[0]/d1[1])*(AcordePlane - r[1]);
 168
 169        if (x > roof) {
 170         x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
 171         z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
 172        }
 173        if (x < -roof) {
 174         x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
 175         z = z -  TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
 176        }
 177
 178        fModules->Fill(z, x);
 179       }
 180
 181       if (d2[1] > 0 && d1[1] < 0 && track2->GetTPCNcls() > 50) {
 182        Double_t r[3];
 183        track2->GetXYZ(r);
 184        Double_t x,z;
 185        z = r[2] + (d2[2]/d2[1])*(AcordePlane - r[1]);
 186        x = r[0] + (d2[0]/d2[1])*(AcordePlane - r[1]);
 187
 188        if (x > roof) {
 189         x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
 190         z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
 191        }
 192        if (x < -roof) {
 193         x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
 194         z = z -  TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
 195        }
 196
 197        fModules->Fill(z, x);
 198       }
 199
 200   //     AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
 201 //       AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
 202 //       delete trackOut;
 203
 204
 205
 206
 207
 208       break;
 209      }
 210     }
 211   }
 212
 213
 214
 215
 216 }
 217
 218 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
 219   //
 220   // Find cosmic pairs
 221   //
 222   // Track0 is choosen in upper TPC part
 223   // Track1 is choosen in lower TPC part
 224   //
 225   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 226   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 227   Int_t ntracks=event->GetNumberOfTracks();
 228   TObjArray  tpcSeeds(ntracks);
 229   if (ntracks==0) return;
 230   Double_t vtxx[3]={0,0,0};
 231   Double_t svtxx[3]={0.000001,0.000001,100.};
 232   AliESDVertex vtx(vtxx,svtxx);
 233   //
 234   //track loop
 235   //
 236   for (Int_t i=0;i<ntracks;++i) {
 237    AliESDtrack *track = event->GetTrack(i);
 238    fClusters->Fill(track->GetTPCNcls());
 239    AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
 240
 241    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 242    TObject *calibObject;
 243    AliTPCseed *seed = 0;
 244    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 245      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 246    }
 247    if (seed) tpcSeeds.AddAt(seed,i);
 248    if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
 249   }
 250   if (ntracks<2) return;
 251   //
 252   // Find pairs
 253   //
 254   for (Int_t i=0;i<ntracks;++i) {
 255     AliESDtrack *track0 = event->GetTrack(i);
 256     // track0 - choosen upper part
 257     if (!track0) continue;
 258     if (!track0->GetOuterParam()) continue;
 259     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 260     Double_t d1[3];
 261     track0->GetDirection(d1);
 262     for (Int_t j=0;j<ntracks;++j) {
 263       if (i==j) continue;
 264       AliESDtrack *track1 = event->GetTrack(j);
 265       //track 1 lower part
 266       if (!track1) continue;
 267       if (!track1->GetOuterParam()) continue;
 268       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 269       //
 270       Double_t d2[3];
 271       track1->GetDirection(d2);
 272       printf("My stream level=%d\n",fStreamLevel);
 273       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 274       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 275       if (! seed0) continue;
 276       if (! seed1) continue;
 277       Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0);
 278       Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0);
 279       Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
 280       Float_t d0  = track0->GetLinearD(0,0);
 281       Float_t d1  = track1->GetLinearD(0,0);
 282       //
 283       // conservative cuts - convergence to be guarantied
 284       // applying before track propagation
 285       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 286       if (dir>fCutMinDir) continue;               // direction vector product
 287       Float_t bz = AliTracker::GetBz();
 288       Float_t dvertex0[2];   //distance to 0,0
 289       Float_t dvertex1[2];   //distance to 0,0
 290       track0->GetDZ(0,0,0,bz,dvertex0);
 291       track1->GetDZ(0,0,0,bz,dvertex1);
 292       if (TMath::Abs(dvertex0[1])>250) continue;
 293       if (TMath::Abs(dvertex1[1])>250) continue;
 294       //
 295       //
 296       //
 297       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 298       AliExternalTrackParam param0(*track0);
 299       AliExternalTrackParam param1(*track1);
 300       //
 301       // Propagate using Magnetic field and correct fo material budget
 302       //
 303       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 304       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 305       //
 306       // Propagate rest to the 0,0 DCA - z should be ignored
 307       //
 308       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
 309       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
 310       //
 311       param0.GetDZ(0,0,0,bz,dvertex0);
 312       param1.GetDZ(0,0,0,bz,dvertex1);
 313       //
 314       Double_t xyz0[3];//,pxyz0[3];
 315       Double_t xyz1[3];//,pxyz1[3];
 316       param0.GetXYZ(xyz0);
 317       param1.GetXYZ(xyz1);
 318       Bool_t isPair = IsPair(&param0,&param1);
 319       //
 320       if (fStreamLevel>0){
 321         TTreeSRedirector * cstream =  GetDebugStreamer();
 322         printf("My stream=%p\n",(void*)cstream);
 323         if (cstream) {
 324           (*cstream) << "Track0" <<
 325             "dir="<<dir<<               //  direction
 326             "OK="<<isPair<<             // will be accepted
 327             "b0="<<b0<<                 //  propagate status
 328             "b1="<<b1<<                 //  propagate status
 329             "Orig0.=" << track0 <<      //  original track  0
 330             "Orig1.=" << track1 <<      //  original track  1
 331             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
 332             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0
 333             "v00="<<dvertex0[0]<<       //  distance using kalman
 334             "v01="<<dvertex0[1]<<       //
 335             "v10="<<dvertex1[0]<<       //
 336             "v11="<<dvertex1[1]<<       //
 337             "d0="<<d0<<                 //  linear distance to 0,0
 338             "d1="<<d1<<                 //  linear distance to 0,0
 339             //
 340             "x00="<<xyz0[0]<<
 341             "x01="<<xyz0[1]<<
 342             "x02="<<xyz0[2]<<
 343             //
 344             "x10="<<xyz1[0]<<
 345             "x11="<<xyz1[1]<<
 346             "x12="<<xyz1[2]<<
 347             //
 348             "Seed0.=" << track0 <<      //  original seed 0
 349             "Seed1.=" << track1 <<      //  original seed 1
 350             "dedx0="<<dedx0<<           //  dedx0
 351             "dedx1="<<dedx1<<           //  dedx1
 352             "\n";
 353         }
 354       }
 355     }
 356   }
 357 }
 358
 359 /*
 360
 361
 362 void AliTPCcalibCosmic::dEdxCorrection(){
 363   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 364   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 365   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 366   TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>70");
 367 }
 368
 369 */
 370
 371
 372
 373 Long64_t AliTPCcalibCosmic::Merge(TCollection */*li*/) {
 374
 375 }
 376
 377 Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 378   //
 379   //
 380   /*
 381   // 0. Same direction - OPOSITE  - cutDir +cutT
 382   TCut cutDir("cutDir","dir<-0.99")
 383   // 1.
 384   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 385   //
 386   // 2. The same rphi
 387   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 388   //
 389   //
 390   //
 391   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 392   // 1/Pt diff cut
 393   */
 394   const Double_t *p0 = tr0->GetParameter();
 395   const Double_t *p1 = tr1->GetParameter();
 396   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 397   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 398   Double_t d0[3], d1[3];
 399   tr0->GetDirection(d0);
 400   tr1->GetDirection(d1);
 401   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 402   //
 403   return kTRUE;
 404 }
 405
 406
 407
 408 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
 409
 410   // Method for the correct logarithmic binning of histograms
 411
 412   TAxis *axis = h->GetXaxis();
 413   int bins = axis->GetNbins();
 414
 415   Double_t from = axis->GetXmin();
 416   Double_t to = axis->GetXmax();
 417   Double_t *new_bins = new Double_t[bins + 1];
 418
 419   new_bins[0] = from;
 420   Double_t factor = pow(to/from, 1./bins);
 421
 422   for (int i = 1; i <= bins; i++) {
 423    new_bins[i] = factor * new_bins[i-1];
 424   }
 425   axis->Set(bins, new_bins);
 426   delete new_bins;
 427
 428 }
 429