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 /*
  17     Comments to be written here:
  18     1. What do we calibrate.
  19     2. How to interpret results
  20     3. Simple example
  21     4. Analysis using debug streamers.
  22
  23
  24
  25     3.Simple example
  26     // To make cosmic scan the user interaction neccessary
  27     //
  28     .x ~/UliStyle.C
  29     gSystem->Load("libANALYSIS");
  30     gSystem->Load("libTPCcalib");
  31     TFile fcalib("CalibObjects.root");
  32     TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  33     AliTPCcalibCosmic * cosmic = ( AliTPCcalibCosmic *)array->FindObject("cosmicTPC");
  34
  35
  36
  37 */
  38
  39
  40
  41 #include "Riostream.h"
  42 #include "TChain.h"
  43 #include "TTree.h"
  44 #include "TH1F.h"
  45 #include "TH2F.h"
  46 #include "TList.h"
  47 #include "TMath.h"
  48 #include "TCanvas.h"
  49 #include "TFile.h"
  50 #include "TF1.h"
  51
  52 #include "AliTPCclusterMI.h"
  53 #include "AliTPCseed.h"
  54 #include "AliESDVertex.h"
  55 #include "AliESDEvent.h"
  56 #include "AliESDfriend.h"
  57 #include "AliESDInputHandler.h"
  58 #include "AliAnalysisManager.h"
  59
  60 #include "AliTracker.h"
  61 #include "AliMagFMaps.h"
  62 #include "AliTPCCalROC.h"
  63
  64 #include "AliLog.h"
  65
  66 #include "AliTPCcalibCosmic.h"
  67
  68 #include "TTreeStream.h"
  69 #include "AliTPCTracklet.h"
  70
  71 ClassImp(AliTPCcalibCosmic)
  72
  73
  74 AliTPCcalibCosmic::AliTPCcalibCosmic()
  75   :AliTPCcalibBase(),
  76    fGainMap(0),
  77    fHistNTracks(0),
  78    fClusters(0),
  79    fModules(0),
  80    fHistPt(0),
  81    fDeDx(0),
  82    fDeDxMIP(0),
  83    fMIPvalue(1),
  84    fCutMaxD(5),        // maximal distance in rfi ditection
  85    fCutMaxDz(40),      // maximal distance in z ditection
  86    fCutTheta(0.03),    // maximal distan theta
  87    fCutMinDir(-0.99)   // direction vector products
  88 {
  89   AliInfo("Default Constructor");
  90 }
  91
  92
  93 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
  94   :AliTPCcalibBase(),
  95    fGainMap(0),
  96    fHistNTracks(0),
  97    fClusters(0),
  98    fModules(0),
  99    fHistPt(0),
 100    fDeDx(0),
 101    fDeDxMIP(0),
 102    fMIPvalue(1),
 103    fCutMaxD(5),        // maximal distance in rfi ditection
 104    fCutMaxDz(40),      // maximal distance in z ditection
 105    fCutTheta(0.03),    // maximal distan theta
 106    fCutMinDir(-0.99)   // direction vector products
 107 {
 108   SetName(name);
 109   SetTitle(title);
 110   AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
 111   AliTracker::SetFieldMap(field, kTRUE);
 112
 113   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
 114   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
 115   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
 116   fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
 117   fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
 118   BinLogX(fDeDx);
 119   fDeDxMIP =  new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
 120
 121   AliInfo("Non Default Constructor");
 122   //
 123 }
 124
 125 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
 126   //
 127   //
 128   //
 129 }
 130
 131
 132
 133
 134
 135 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
 136   //
 137   //
 138   //
 139   if (!event) {
 140     Printf("ERROR: ESD not available");
 141     return;
 142   }
 143   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 144   if (!ESDfriend) {
 145    Printf("ERROR: ESDfriend not available");
 146    return;
 147   }
 148
 149
 150   FindPairs(event); // nearly everything takes place in find pairs...
 151
 152   if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
 153   Int_t ntracks=event->GetNumberOfTracks();
 154   fHistNTracks->Fill(ntracks);
 155   if (ntracks==0) return;
 156
 157 }
 158
 159
 160
 161 void AliTPCcalibCosmic::Analyze() {
 162
 163   fMIPvalue = CalculateMIPvalue(fDeDxMIP);
 164
 165   return;
 166
 167 }
 168
 169
 170
 171 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
 172   //
 173   // Find cosmic pairs
 174   //
 175   // Track0 is choosen in upper TPC part
 176   // Track1 is choosen in lower TPC part
 177   //
 178   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
 179   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 180   Int_t ntracks=event->GetNumberOfTracks();
 181   TObjArray  tpcSeeds(ntracks);
 182   if (ntracks==0) return;
 183   Double_t vtxx[3]={0,0,0};
 184   Double_t svtxx[3]={0.000001,0.000001,100.};
 185   AliESDVertex vtx(vtxx,svtxx);
 186   //
 187   //track loop
 188   //
 189   for (Int_t i=0;i<ntracks;++i) {
 190    AliESDtrack *track = event->GetTrack(i);
 191    fClusters->Fill(track->GetTPCNcls());
 192
 193    const AliExternalTrackParam * trackIn = track->GetInnerParam();
 194    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 195    if (!trackIn) continue;
 196    if (!trackOut) continue;
 197
 198    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 199    TObject *calibObject;
 200    AliTPCseed *seed = 0;
 201    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 202      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 203    }
 204    if (seed) tpcSeeds.AddAt(seed,i);
 205
 206    Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 207    if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
 208      fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
 209      //
 210      if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
 211      //
 212      if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap)>300) {
 213        TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
 214        if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
 215        if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
 216      }
 217
 218    }
 219
 220   }
 221
 222   if (ntracks<2) return;
 223   //
 224   // Find pairs
 225   //
 226   for (Int_t i=0;i<ntracks;++i) {
 227     AliESDtrack *track0 = event->GetTrack(i);
 228     // track0 - choosen upper part
 229     if (!track0) continue;
 230     if (!track0->GetOuterParam()) continue;
 231     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 232     Double_t dir0[3];
 233     track0->GetDirection(dir0);
 234     for (Int_t j=0;j<ntracks;++j) {
 235       if (i==j) continue;
 236       AliESDtrack *track1 = event->GetTrack(j);
 237       //track 1 lower part
 238       if (!track1) continue;
 239       if (!track1->GetOuterParam()) continue;
 240       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 241       //
 242       Double_t dir1[3];
 243       track1->GetDirection(dir1);
 244
 245       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 246       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 247       if (! seed0) continue;
 248       if (! seed1) continue;
 249       Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
 250       Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
 251       //
 252       Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
 253       Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
 254       //
 255       Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
 256       Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
 257       //
 258       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
 259       Float_t d0  = track0->GetLinearD(0,0);
 260       Float_t d1  = track1->GetLinearD(0,0);
 261       //
 262       // conservative cuts - convergence to be guarantied
 263       // applying before track propagation
 264       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 265       if (dir>fCutMinDir) continue;               // direction vector product
 266       Float_t bz = AliTracker::GetBz();
 267       Float_t dvertex0[2];   //distance to 0,0
 268       Float_t dvertex1[2];   //distance to 0,0
 269       track0->GetDZ(0,0,0,bz,dvertex0);
 270       track1->GetDZ(0,0,0,bz,dvertex1);
 271       if (TMath::Abs(dvertex0[1])>250) continue;
 272       if (TMath::Abs(dvertex1[1])>250) continue;
 273       //
 274       //
 275       //
 276       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 277       AliExternalTrackParam param0(*track0);
 278       AliExternalTrackParam param1(*track1);
 279       //
 280       // Propagate using Magnetic field and correct fo material budget
 281       //
 282       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 283       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 284       //
 285       // Propagate rest to the 0,0 DCA - z should be ignored
 286       //
 287       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
 288       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
 289       //
 290       param0.GetDZ(0,0,0,bz,dvertex0);
 291       param1.GetDZ(0,0,0,bz,dvertex1);
 292       if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
 293       //
 294       Double_t xyz0[3];//,pxyz0[3];
 295       Double_t xyz1[3];//,pxyz1[3];
 296       param0.GetXYZ(xyz0);
 297       param1.GetXYZ(xyz1);
 298       Bool_t isPair = IsPair(&param0,&param1);
 299       //
 300       if (isPair) FillAcordeHist(track0);
 301       //
 302       if (fStreamLevel>0){
 303         TTreeSRedirector * cstream =  GetDebugStreamer();
 304         //printf("My stream=%p\n",(void*)cstream);
 305         AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
 306         AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
 307         AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
 308         AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
 309         Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
 310         Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
 311         Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
 312         Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
 313         if (cstream) {
 314           (*cstream) << "Track0" <<
 315             "run="<<fRun<<              //  run number
 316             "event="<<fEvent<<          //  event number
 317             "time="<<fTime<<            //  time stamp of event
 318             "trigger="<<fTrigger<<      //  trigger
 319             "mag="<<fMagF<<             //  magnetic field
 320             "dir="<<dir<<               //  direction
 321             "OK="<<isPair<<             //  will be accepted
 322             "b0="<<b0<<                 //  propagate status
 323             "b1="<<b1<<                 //  propagate status
 324             "crossI="<<isCrossI<<       //  cross inner
 325             "crossO="<<isCrossO<<       //  cross outer
 326             //
 327             "Orig0.=" << track0 <<      //  original track  0
 328             "Orig1.=" << track1 <<      //  original track  1
 329             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
 330             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0
 331             "Ip0.="<<ip0<<              //  inner param - upper
 332             "Ip1.="<<ip1<<              //  inner param - lower
 333             "Op0.="<<op0<<              //  outer param - upper
 334             "Op1.="<<op1<<              //  outer param - lower
 335             //
 336             "v00="<<dvertex0[0]<<       //  distance using kalman
 337             "v01="<<dvertex0[1]<<       //
 338             "v10="<<dvertex1[0]<<       //
 339             "v11="<<dvertex1[1]<<       //
 340             "d0="<<d0<<                 //  linear distance to 0,0
 341             "d1="<<d1<<                 //  linear distance to 0,0
 342             //
 343             //
 344             //
 345             "x00="<<xyz0[0]<<           // global position close to vertex
 346             "x01="<<xyz0[1]<<
 347             "x02="<<xyz0[2]<<
 348             //
 349             "x10="<<xyz1[0]<<           // global position close to vertex
 350             "x11="<<xyz1[1]<<
 351             "x12="<<xyz1[2]<<
 352             //
 353             "alpha0="<<alpha0<<
 354             "alpha1="<<alpha1<<
 355             "dir00="<<dir0[0]<<           // direction upper
 356             "dir01="<<dir0[1]<<
 357             "dir02="<<dir0[2]<<
 358             //
 359             "dir10="<<dir1[0]<<           // direction lower
 360             "dir11="<<dir1[1]<<
 361             "dir12="<<dir1[2]<<
 362             //
 363             //
 364             "Seed0.=" << seed0 <<       //  original seed 0
 365             "Seed1.=" << seed1 <<       //  original seed 1
 366             //
 367             "dedx0="<<dedx0<<           //  dedx0 - all
 368             "dedx1="<<dedx1<<           //  dedx1 - all
 369             //
 370             "dedx0I="<<dedx0I<<         //  dedx0 - inner ROC
 371             "dedx1I="<<dedx1I<<         //  dedx1 - inner ROC
 372             //
 373             "dedx0O="<<dedx0O<<         //  dedx0 - outer ROC
 374             "dedx1O="<<dedx1O<<         //  dedx1 - outer ROC
 375             "\n";
 376         }
 377       }
 378     }
 379   }
 380 }
 381
 382
 383
 384
 385 void  AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
 386
 387   // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
 388   if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
 389
 390   const Double_t AcordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
 391   const Double_t roof = 210.5;       // distance from x =0 to end of magnet roof
 392
 393   Double_t r[3];
 394   upperTrack->GetXYZ(r);
 395   Double_t d[3];
 396   upperTrack->GetDirection(d);
 397   Double_t x,z;
 398   z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
 399   x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
 400
 401   if (x > roof) {
 402     x = r[0] + (d[0]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
 403     z = r[2] + (d[2]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
 404   }
 405   if (x < -roof) {
 406     x = r[0] + (d[0]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
 407     z = r[2] + (d[2]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
 408   }
 409
 410   fModules->Fill(z, x);
 411
 412 }
 413
 414
 415
 416 Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
 417
 418   TIterator* iter = li->MakeIterator();
 419   AliTPCcalibCosmic* cal = 0;
 420
 421   while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
 422     if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
 423       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 424       return -1;
 425     }
 426
 427     fHistNTracks->Add(cal->GetHistNTracks());
 428     fClusters->Add(cal-> GetHistClusters());
 429     fModules->Add(cal->GetHistAcorde());
 430     fHistPt->Add(cal->GetHistPt());
 431     fDeDx->Add(cal->GetHistDeDx());
 432     fDeDxMIP->Add(cal->GetHistMIP());
 433
 434   }
 435
 436   return 0;
 437
 438 }
 439
 440
 441 Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 442   //
 443   //
 444   /*
 445   // 0. Same direction - OPOSITE  - cutDir +cutT
 446   TCut cutDir("cutDir","dir<-0.99")
 447   // 1.
 448   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 449   //
 450   // 2. The same rphi
 451   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 452   //
 453   //
 454   //
 455   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 456   // 1/Pt diff cut
 457   */
 458   const Double_t *p0 = tr0->GetParameter();
 459   const Double_t *p1 = tr1->GetParameter();
 460   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 461   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
 462   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 463
 464   Double_t d0[3], d1[3];
 465   tr0->GetDirection(d0);
 466   tr1->GetDirection(d1);
 467   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 468   //
 469   return kTRUE;
 470 }
 471
 472
 473
 474 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
 475
 476   TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
 477   funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
 478                                 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
 479   hist->Fit(funcDoubleGaus);
 480   Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
 481
 482   delete funcDoubleGaus;
 483
 484   return MIPvalue;
 485
 486 }
 487
 488
 489
 490
 491 void AliTPCcalibCosmic::CalculateBetheParams(TH2F *hist, Double_t * initialParam) {
 492
 493   return;
 494
 495 }
 496
 497
 498 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
 499
 500   // Method for the correct logarithmic binning of histograms
 501
 502   TAxis *axis = h->GetXaxis();
 503   int bins = axis->GetNbins();
 504
 505   Double_t from = axis->GetXmin();
 506   Double_t to = axis->GetXmax();
 507   Double_t *new_bins = new Double_t[bins + 1];
 508
 509   new_bins[0] = from;
 510   Double_t factor = pow(to/from, 1./bins);
 511
 512   for (int i = 1; i <= bins; i++) {
 513    new_bins[i] = factor * new_bins[i-1];
 514   }
 515   axis->Set(bins, new_bins);
 516   delete new_bins;
 517
 518 }
 519
 520 AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
 521 {
 522   //
 523   // Invert paramerameter  - not covariance yet
 524   //
 525   AliExternalTrackParam *output = new AliExternalTrackParam(*input);
 526   Double_t * param = (Double_t*)output->GetParameter();
 527   param[0]*=-1;
 528   param[3]*=-1;
 529   param[4]*=-1;
 530   //
 531   return output;
 532 }
 533
 534 /*
 535
 536 void Init(){
 537
 538 .x ~/UliStyle.C
 539 .x ~/rootlogon.C
 540 gSystem->Load("libSTAT.so");
 541 gSystem->Load("libANALYSIS");
 542 gSystem->Load("libTPCcalib");
 543 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 544
 545 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
 546 AliXRDPROOFtoolkit tool;
 547 TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
 548 chain->Lookup();
 549
 550 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01");  // OK
 551 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2");     // OK
 552 TCut cutP1("cutP1","abs(Tr0.fP[1]-Tr1.fP[1])<10");   // OK
 553 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 554 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
 555 TCut cutA=cutT+cutD+cutPt+cutN+cutP1;
 556
 557 TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
 558
 559 //
 560 chain->Draw(">>listEL",cutA,"entryList");
 561 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
 562 chain->SetEntryList(elist);
 563 //
 564 chain->Draw(">>listV40Z100","abs(d0)<40&&abs(v01)<100","entryList");
 565 TEntryList *elistV40Z100 = (TEntryList*)gDirectory->Get("listV40Z100");
 566 chain->SetEntryList(elistV40Z100);
 567
 568 //
 569 // Aliases
 570 //
 571 chain->SetAlias("side","(-1+(Tr0.fP[1]>0)*2)");
 572 chain->SetAlias("hpt","abs(Tr0.fP[4])<0.2");
 573 chain->SetAlias("signy","(-1+(Tr0.fP[0]>0)*2)");
 574
 575 chain->SetAlias("dy","Tr0.fP[0]+Tr1.fP[0]");
 576 chain->SetAlias("dz","Tr0.fP[1]-Tr1.fP[1]");
 577 chain->SetAlias("d1pt","Tr0.fP[4]+Tr1.fP[4]");
 578 chain->SetAlias("dtheta","(Tr0.fP[3]+Tr1.fP[3])");
 579 chain->SetAlias("dphi","(Tr0.fAlpha-Tr1.fAlpha-pi)");
 580
 581 chain->SetAlias("mtheta","(Tr0.fP[3]-Tr1.fP[3])*0.5")
 582 chain->SetAlias("sa","(sin(Tr0.fAlpha+0.))");
 583 chain->SetAlias("ca","(cos(Tr0.fAlpha+0.))");
 584
 585
 586
 587 chain->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
 588 hisdyA->FitSlicesY();
 589 hisdyA_2->SetXTitle("#sqrt{1/p_{t}}");
 590 hisdyA_2->SetYTitle("#sigma_{y}(cm)");
 591 hisdyA_2->SetTitle("Cosmic - Y matching");
 592 hisdyA_2->SetMaximum(0.5);
 593
 594
 595 chain->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
 596 hisdyC->FitSlicesY();
 597 hisdyC_2->SetXTitle("#sqrt{1/p_{t}}");
 598 hisdyC_2->SetYTitle("#sigma_{y}(cm)");
 599 hisdyC_2->SetTitle("Cosmic - Y matching");
 600 hisdyC_2->SetMaximum(1);
 601 hisdyC_2->SetMinimum(0);
 602 hisdyC_2->SetMarkerStyle(22);
 603 hisdyA_2->SetMarkerStyle(21);
 604 hisdyC_2->SetMarkerSize(1.5);
 605 hisdzA_2->SetMarkerSize(1.5);
 606 hisdyC_2->Draw();
 607 hisdyA_2->Draw("same");
 608 gPad->SaveAs("~/Calibration/Cosmic/pic/ymatching.gif")
 609
 610 chain->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
 611 hisdzA->FitSlicesY();
 612 hisdzA_2->SetXTitle("#sqrt{1/p_{t}}");
 613 hisdzA_2->SetYTitle("#sigma_{z}(cm)");
 614 hisdzA_2->SetTitle("Cosmic - Z matching - A side ");
 615 hisdzA_2->SetMaximum(0.5);
 616
 617 chain->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
 618 hisdzC->FitSlicesY();
 619 hisdzC_2->SetXTitle("#sqrt{1/p_{t}}");
 620 hisdzC_2->SetYTitle("#sigma_{z}(cm)");
 621 hisdzC_2->SetTitle("Cosmic - Z matching");
 622 hisdzC_2->SetMaximum(0.5);
 623 hisdzC_2->SetMarkerStyle(22);
 624 hisdzA_2->SetMarkerStyle(21);
 625 hisdzC_2->SetMarkerSize(1.5);
 626 hisdzA_2->SetMarkerSize(1.5);
 627
 628 hisdzC_2->Draw();
 629 hisdzA_2->Draw("same");
 630
 631
 632 //
 633 // PICTURE 1/pt
 634 //
 635 chain->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptA(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
 636 hisd1ptA->FitSlicesY();
 637 hisd1ptA_2->SetXTitle("#sqrt{1/p_{t}}");
 638 hisd1ptA_2->SetYTitle("#sigma_{z}(cm)");
 639 hisd1ptA_2->SetTitle("Cosmic - Z matching - A side ");
 640 hisd1ptA_2->SetMaximum(0.5);
 641
 642 chain->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptC(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
 643 hisd1ptC->FitSlicesY();
 644 hisd1ptC_2->SetXTitle("#sqrt{1/p_{t}}");
 645 hisd1ptC_2->SetYTitle("#sigma_{1/pt}(1/GeV)");
 646 hisd1ptC_2->SetTitle("Cosmic - 1/pt matching");
 647 hisd1ptC_2->SetMaximum(0.05);
 648 hisd1ptC_2->SetMarkerStyle(22);
 649 hisd1ptA_2->SetMarkerStyle(21);
 650 hisd1ptC_2->SetMarkerSize(1.5);
 651 hisd1ptA_2->SetMarkerSize(1.5);
 652
 653 hisd1ptC_2->Draw();
 654 hisd1ptA_2->Draw("same");
 655 gPad->SaveAs("~/Calibration/Cosmic/pic/1ptmatching.gif")
 656
 657 //
 658 // Theta
 659 //
 660 chain->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
 661 hisdthetaA->FitSlicesY();
 662 hisdthetaA_2->SetXTitle("#sqrt{1/p_{t}}");
 663 hisdthetaA_2->SetYTitle("#sigma_{#theta}(cm)");
 664 hisdthetaA_2->SetTitle("Cosmic - Z matching - A side ");
 665 hisdthetaA_2->SetMaximum(0.5);
 666
 667 chain->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
 668 hisdthetaC->FitSlicesY();
 669 hisdthetaC_2->SetXTitle("#sqrt{1/p_{t}}");
 670 hisdthetaC_2->SetYTitle("#sigma_{#theta}(rad)");
 671 hisdthetaC_2->SetTitle("Cosmic - Theta matching");
 672 hisdthetaC_2->SetMaximum(0.01);
 673 hisdthetaC_2->SetMinimum(0.0);
 674 hisdthetaC_2->SetMarkerStyle(22);
 675 hisdthetaA_2->SetMarkerStyle(21);
 676 hisdthetaC_2->SetMarkerSize(1.5);
 677 hisdthetaA_2->SetMarkerSize(1.5);
 678
 679 hisdthetaC_2->Draw();
 680 hisdthetaA_2->Draw("same");
 681 gPad->SaveAs("~/Calibration/Cosmic/pic/thetamatching.gif")
 682 //
 683 // Phi
 684 //
 685 chain->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
 686 hisdphiA->FitSlicesY();
 687 hisdphiA_2->SetXTitle("#sqrt{1/p_{t}}");
 688 hisdphiA_2->SetYTitle("#sigma_{#phi}(rad)");
 689 hisdphiA_2->SetTitle("Cosmic - Z matching - A side ");
 690 hisdphiA_2->SetMaximum(0.5);
 691
 692 chain->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
 693 hisdphiC->FitSlicesY();
 694 hisdphiC_2->SetXTitle("#sqrt{1/p_{t}}");
 695 hisdphiC_2->SetYTitle("#sigma_{#phi}(rad)");
 696 hisdphiC_2->SetTitle("Cosmic - Phi matching");
 697 hisdphiC_2->SetMaximum(0.01);
 698 hisdphiC_2->SetMinimum(0.0);
 699 hisdphiC_2->SetMarkerStyle(22);
 700 hisdphiA_2->SetMarkerStyle(21);
 701 hisdphiC_2->SetMarkerSize(1.5);
 702 hisdphiA_2->SetMarkerSize(1.5);
 703
 704 hisdphiC_2->Draw();
 705 hisdphiA_2->Draw("same");
 706 gPad->SaveAs("~/Calibration/Cosmic/pic/phimatching.gif")
 707
 708
 709
 710 }
 711
 712
 713 */
 714
 715
 716 /*
 717 void MatchTheta(){
 718
 719 TStatToolkit toolkit;
 720 Double_t chi2=0;
 721 Int_t    npoints=0;
 722 TVectorD fitParamA0;
 723 TVectorD fitParamA1;
 724 TVectorD fitParamC0;
 725 TVectorD fitParamC1;
 726 TMatrixD covMatrix;
 727
 728
 729 TString fstring="";
 730 //
 731 fstring+="mtheta++";
 732 fstring+="ca++";
 733 fstring+="sa++";
 734 fstring+="ca*mtheta++";
 735 fstring+="sa*mtheta++";
 736 //
 737 fstring+="side++";
 738 fstring+="side*mtheta++";
 739 fstring+="side*ca++";
 740 fstring+="side*sa++";
 741 fstring+="side*ca*mtheta++";
 742 fstring+="side*sa*mtheta++";
 743
 744
 745 TString *strTheta0 = toolkit.FitPlane(chain,"dtheta",fstring->Data(), "hpt&&!crossI&&!crossO", chi2,npoints,fitParamA0,covMatrix,0.8);
 746 chain->SetAlias("dtheta0",strTheta0.Data());
 747 strTheta0->Tokenize("+")->Print();
 748
 749
 750 //fstring+="mtheta++";
 751 //fstring+="mtheta^2++";
 752 //fstring+="ca*mtheta^2++";
 753 //fstring+="sa*mtheta^2++";
 754
 755
 756
 757 }
 758
 759 */
 760
 761
 762
 763
 764 /*
 765  void PosCorrection()
 766
 767
 768
 769
 770  TStatToolkit toolkit;
 771  Double_t chi2=0;
 772  Int_t    npoints=0;
 773  TVectorD fitParam;
 774  TMatrixD covMatrix;
 775
 776  //Theta
 777 chain->SetAlias("dthe","(Tr0.fP[3]+Tr1.fP[3])");
 778 chain->SetAlias("sign","(-1+(Tr0.fP[1]>0)*2)");
 779 chain->SetAlias("di","(1-abs(Tr0.fP[1])/250)");
 780 //
 781 //
 782 TString strFit="";
 783 //
 784 strFit+="sign++";                              //1
 785 strFit+="Tr0.fP[3]++";                         //2
 786 //
 787 strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]>0)++";     //3
 788 strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]<0)++";     //4
 789 //
 790 strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]>0)++";     //5
 791 strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]<0)++";     //6
 792 //
 793 strFit+="sin(Tr0.fAlpha)*Tr0.fP[3]++";         //7
 794 strFit+="cos(Tr0.fAlpha)*Tr0.fP[3]++";         //8
 795
 796
 797  //
 798  TString * thetaParam = toolkit.FitPlane(chain,"dthe", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix,0.8,0,10000)
 799
 800  chain->SetAlias("corTheta",thetaParam->Data());
 801  chain->Draw("dthe:Tr0.fP[1]","","",50000);
 802
 803
 804
 805 */
 806
 807
 808
 809 /*
 810
 811 void AliTPCcalibCosmic::dEdxCorrection(){
 812   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01");  // OK
 813   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2");     // OK
 814   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 815   TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
 816   TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
 817   TCut cutA=cutT+cutD+cutPt+cutN+cutS;
 818
 819
 820  .x ~/UliStyle.C
 821   gSystem->Load("libANALYSIS");
 822   gSystem->Load("libTPCcalib");
 823   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 824   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
 825  AliXRDPROOFtoolkit tool;
 826   TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
 827   chain->Lookup();
 828
 829   chain->Draw(">>listEL",cutA,"entryList");
 830   TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
 831   chain->SetEntryList(elist);
 832
 833   .x ~/rootlogon.C
 834    gSystem->Load("libSTAT.so");
 835    TStatToolkit toolkit;
 836   Double_t chi2=0;
 837   Int_t    npoints=0;
 838   TVectorD fitParam;
 839   TMatrixD covMatrix;
 840
 841   chain->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
 842
 843   TString strFit;
 844   strFit+="(Tr0.fP[1]/250)++";
 845   strFit+="(Tr0.fP[1]/250)^2++";
 846   strFit+="(Tr0.fP[3])++";
 847   strFit+="(Tr0.fP[3])^2++";
 848
 849   TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix)
 850
 851
 852
 853 */
 854
 855
 856 /*
 857 gSystem->Load("libANALYSIS");
 858 gSystem->Load("libSTAT");
 859 gSystem->Load("libTPCcalib");
 860
 861 TStatToolkit toolkit;
 862 Double_t chi2;
 863 TVectorD fitParam;
 864 TMatrixD covMatrix;
 865 Int_t npoints;
 866 //
 867 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03");  // OK
 868 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5");     // OK
 869 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 870 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
 871 TCut cutA=cutT+cutD+cutPt+cutN;
 872
 873
 874
 875 TTree * chain = Track0;
 876
 877
 878 chain->SetAlias("norm","signalTot0.fElements[3]/signalTot1.fElements[3]");
 879 //
 880 chain->SetAlias("dr1","(signalTot0.fElements[1]/signalTot0.fElements[3])");
 881 chain->SetAlias("dr2","(signalTot0.fElements[2]/signalTot0.fElements[3])");
 882 chain->SetAlias("dr4","(signalTot0.fElements[4]/signalTot0.fElements[3])");
 883 chain->SetAlias("dr5","(signalTot0.fElements[5]/signalTot0.fElements[3])");
 884
 885 TString fstring="";
 886 fstring+="dr1++";
 887 fstring+="dr2++";
 888 fstring+="dr4++";
 889 fstring+="dr5++";
 890 //
 891 fstring+="dr1*dr2++";
 892 fstring+="dr1*dr4++";
 893 fstring+="dr1*dr5++";
 894 fstring+="dr2*dr4++";
 895 fstring+="dr2*dr5++";
 896 fstring+="dr4*dr5++";
 897
 898
 899
 900 TString *strqdedx = toolkit.FitPlane(chain,"norm",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
 901
 902 chain->SetAlias("corQT",strqdedx->Data());
 903
 904 */
 905
 906
 907 /*
 908   chain->SetProof(kTRUE);
 909   chain->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE):Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)",""+cutA,"",1000);
 910
 911
 912 chain->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)/Seed1.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)>>his(100,0.5,1.5)","min(Orig0.fTPCncls,Orig1.fTPCncls)>130"+cutA,"",50000);
 913
 914 */
 915
 916