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    fCutTheta(0.03),    // maximal distan theta
  86    fCutMinDir(-0.99)   // direction vector products
  87 {
  88   AliInfo("Default Constructor");
  89 }
  90
  91
  92 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
  93   :AliTPCcalibBase(),
  94    fGainMap(0),
  95    fHistNTracks(0),
  96    fClusters(0),
  97    fModules(0),
  98    fHistPt(0),
  99    fDeDx(0),
 100    fDeDxMIP(0),
 101    fMIPvalue(1),
 102    fCutMaxD(5),        // maximal distance in rfi ditection
 103    fCutTheta(0.03),    // maximal distan theta
 104    fCutMinDir(-0.99)   // direction vector products
 105 {
 106   SetName(name);
 107   SetTitle(title);
 108   AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
 109   AliTracker::SetFieldMap(field, kTRUE);
 110
 111   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
 112   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
 113   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
 114   fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
 115   fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
 116   BinLogX(fDeDx);
 117   fDeDxMIP =  new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
 118
 119   AliInfo("Non Default Constructor");
 120   //
 121 }
 122
 123 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
 124   //
 125   //
 126   //
 127 }
 128
 129
 130
 131
 132
 133 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
 134   //
 135   //
 136   //
 137   if (!event) {
 138     Printf("ERROR: ESD not available");
 139     return;
 140   }
 141   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 142   if (!ESDfriend) {
 143    Printf("ERROR: ESDfriend not available");
 144    return;
 145   }
 146
 147   FindPairs(event); // nearly everything takes place in find pairs...
 148
 149   if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
 150   Int_t ntracks=event->GetNumberOfTracks();
 151   fHistNTracks->Fill(ntracks);
 152   if (ntracks==0) return;
 153
 154 }
 155
 156
 157
 158 void AliTPCcalibCosmic::Analyze() {
 159
 160   fMIPvalue = CalculateMIPvalue(fDeDxMIP);
 161
 162   return;
 163
 164 }
 165
 166
 167
 168 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
 169   //
 170   // Find cosmic pairs
 171   //
 172   // Track0 is choosen in upper TPC part
 173   // Track1 is choosen in lower TPC part
 174   //
 175   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
 176   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 177   Int_t ntracks=event->GetNumberOfTracks();
 178   TObjArray  tpcSeeds(ntracks);
 179   if (ntracks==0) return;
 180   Double_t vtxx[3]={0,0,0};
 181   Double_t svtxx[3]={0.000001,0.000001,100.};
 182   AliESDVertex vtx(vtxx,svtxx);
 183   //
 184   //track loop
 185   //
 186   for (Int_t i=0;i<ntracks;++i) {
 187    AliESDtrack *track = event->GetTrack(i);
 188    fClusters->Fill(track->GetTPCNcls());
 189
 190    const AliExternalTrackParam * trackIn = track->GetInnerParam();
 191    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 192    if (!trackIn) continue;
 193    if (!trackOut) continue;
 194
 195    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 196    TObject *calibObject;
 197    AliTPCseed *seed = 0;
 198    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 199      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 200    }
 201    if (seed) tpcSeeds.AddAt(seed,i);
 202
 203    Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 204    if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
 205      fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
 206      //
 207      if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
 208      //
 209      if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap)>300) {
 210        TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
 211        if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
 212        if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
 213      }
 214
 215    }
 216
 217   }
 218
 219   if (ntracks<2) return;
 220   //
 221   // Find pairs
 222   //
 223   for (Int_t i=0;i<ntracks;++i) {
 224     AliESDtrack *track0 = event->GetTrack(i);
 225     // track0 - choosen upper part
 226     if (!track0) continue;
 227     if (!track0->GetOuterParam()) continue;
 228     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 229     Double_t dir0[3];
 230     track0->GetDirection(dir0);
 231     for (Int_t j=0;j<ntracks;++j) {
 232       if (i==j) continue;
 233       AliESDtrack *track1 = event->GetTrack(j);
 234       //track 1 lower part
 235       if (!track1) continue;
 236       if (!track1->GetOuterParam()) continue;
 237       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 238       //
 239       Double_t dir1[3];
 240       track1->GetDirection(dir1);
 241
 242       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 243       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 244       if (! seed0) continue;
 245       if (! seed1) continue;
 246       Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
 247       Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
 248       //
 249       Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
 250       Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
 251       //
 252       Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
 253       Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
 254       //
 255       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
 256       Float_t d0  = track0->GetLinearD(0,0);
 257       Float_t d1  = track1->GetLinearD(0,0);
 258       //
 259       // conservative cuts - convergence to be guarantied
 260       // applying before track propagation
 261       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 262       if (dir>fCutMinDir) continue;               // direction vector product
 263       Float_t bz = AliTracker::GetBz();
 264       Float_t dvertex0[2];   //distance to 0,0
 265       Float_t dvertex1[2];   //distance to 0,0
 266       track0->GetDZ(0,0,0,bz,dvertex0);
 267       track1->GetDZ(0,0,0,bz,dvertex1);
 268       if (TMath::Abs(dvertex0[1])>250) continue;
 269       if (TMath::Abs(dvertex1[1])>250) continue;
 270       //
 271       //
 272       //
 273       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 274       AliExternalTrackParam param0(*track0);
 275       AliExternalTrackParam param1(*track1);
 276       //
 277       // Propagate using Magnetic field and correct fo material budget
 278       //
 279       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 280       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 281       //
 282       // Propagate rest to the 0,0 DCA - z should be ignored
 283       //
 284       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
 285       Bool_t b1 = 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       if (isPair) FillAcordeHist(track0);
 297       //
 298       if (fStreamLevel>0){
 299         TTreeSRedirector * cstream =  GetDebugStreamer();
 300         printf("My stream=%p\n",(void*)cstream);
 301         AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
 302         AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
 303         AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
 304         AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
 305         Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
 306         Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
 307         Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
 308         Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
 309         Int_t runNr = event->GetRunNumber();
 310         UInt_t time = event->GetTimeStamp();
 311         if (cstream) {
 312           (*cstream) << "Track0" <<
 313             "runNr="<<runNr<<           //  run number
 314             "time="<<time<<             //  time stamp of event
 315             "dir="<<dir<<               //  direction
 316             "OK="<<isPair<<             //  will be accepted
 317             "b0="<<b0<<                 //  propagate status
 318             "b1="<<b1<<                 //  propagate status
 319             "crossI="<<isCrossI<<       //  cross inner
 320             "crossO="<<isCrossO<<       //  cross outer
 321             //
 322             "Orig0.=" << track0 <<      //  original track  0
 323             "Orig1.=" << track1 <<      //  original track  1
 324             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
 325             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0
 326             "Ip0.="<<ip0<<              //  inner param - upper
 327             "Ip1.="<<ip1<<              //  inner param - lower
 328             "Op0.="<<op0<<              //  outer param - upper
 329             "Op1.="<<op1<<              //  outer param - lower
 330             //
 331             "v00="<<dvertex0[0]<<       //  distance using kalman
 332             "v01="<<dvertex0[1]<<       //
 333             "v10="<<dvertex1[0]<<       //
 334             "v11="<<dvertex1[1]<<       //
 335             "d0="<<d0<<                 //  linear distance to 0,0
 336             "d1="<<d1<<                 //  linear distance to 0,0
 337             //
 338             //
 339             //
 340             "x00="<<xyz0[0]<<           // global position close to vertex
 341             "x01="<<xyz0[1]<<
 342             "x02="<<xyz0[2]<<
 343             //
 344             "x10="<<xyz1[0]<<           // global position close to vertex
 345             "x11="<<xyz1[1]<<
 346             "x12="<<xyz1[2]<<
 347             //
 348             "alpha0="<<alpha0<<
 349             "alpha1="<<alpha1<<
 350             "dir00="<<dir0[0]<<           // direction upper
 351             "dir01="<<dir0[1]<<
 352             "dir02="<<dir0[2]<<
 353             //
 354             "dir10="<<dir1[0]<<           // direction lower
 355             "dir11="<<dir1[1]<<
 356             "dir12="<<dir1[2]<<
 357             //
 358             //
 359             "Seed0.=" << seed0 <<       //  original seed 0
 360             "Seed1.=" << seed1 <<       //  original seed 1
 361             //
 362             "dedx0="<<dedx0<<           //  dedx0 - all
 363             "dedx1="<<dedx1<<           //  dedx1 - all
 364             //
 365             "dedx0I="<<dedx0I<<         //  dedx0 - inner ROC
 366             "dedx1I="<<dedx1I<<         //  dedx1 - inner ROC
 367             //
 368             "dedx0O="<<dedx0O<<         //  dedx0 - outer ROC
 369             "dedx1O="<<dedx1O<<         //  dedx1 - outer ROC
 370             "\n";
 371         }
 372       }
 373     }
 374   }
 375 }
 376
 377
 378
 379
 380 void  AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
 381
 382   // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
 383   if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
 384
 385   const Double_t AcordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
 386   const Double_t roof = 210.5;       // distance from x =0 to end of magnet roof
 387
 388   Double_t r[3];
 389   upperTrack->GetXYZ(r);
 390   Double_t d[3];
 391   upperTrack->GetDirection(d);
 392   Double_t x,z;
 393   z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
 394   x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
 395
 396   if (x > roof) {
 397     x = r[0] + (d[0]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
 398     z = r[2] + (d[2]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
 399   }
 400   if (x < -roof) {
 401     x = r[0] + (d[0]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
 402     z = r[2] + (d[2]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
 403   }
 404
 405   fModules->Fill(z, x);
 406
 407 }
 408
 409
 410
 411 Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
 412
 413   TIterator* iter = li->MakeIterator();
 414   AliTPCcalibCosmic* cal = 0;
 415
 416   while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
 417     if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
 418       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 419       return -1;
 420     }
 421
 422     fHistNTracks->Add(cal->GetHistNTracks());
 423     fClusters->Add(cal-> GetHistClusters());
 424     fModules->Add(cal->GetHistAcorde());
 425     fHistPt->Add(cal->GetHistPt());
 426     fDeDx->Add(cal->GetHistDeDx());
 427     fDeDxMIP->Add(cal->GetHistMIP());
 428
 429   }
 430
 431   return 0;
 432
 433 }
 434
 435
 436 Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 437   //
 438   //
 439   /*
 440   // 0. Same direction - OPOSITE  - cutDir +cutT
 441   TCut cutDir("cutDir","dir<-0.99")
 442   // 1.
 443   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 444   //
 445   // 2. The same rphi
 446   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 447   //
 448   //
 449   //
 450   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 451   // 1/Pt diff cut
 452   */
 453   const Double_t *p0 = tr0->GetParameter();
 454   const Double_t *p1 = tr1->GetParameter();
 455   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 456   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 457   Double_t d0[3], d1[3];
 458   tr0->GetDirection(d0);
 459   tr1->GetDirection(d1);
 460   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 461   //
 462   return kTRUE;
 463 }
 464
 465
 466
 467 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
 468
 469   TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
 470   funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
 471                                 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
 472   hist->Fit(funcDoubleGaus);
 473   Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
 474
 475   delete funcDoubleGaus;
 476
 477   return MIPvalue;
 478
 479 }
 480
 481
 482
 483
 484 void AliTPCcalibCosmic::CalculateBetheParams(TH2F *hist, Double_t * initialParam) {
 485
 486   return;
 487
 488 }
 489
 490
 491 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
 492
 493   // Method for the correct logarithmic binning of histograms
 494
 495   TAxis *axis = h->GetXaxis();
 496   int bins = axis->GetNbins();
 497
 498   Double_t from = axis->GetXmin();
 499   Double_t to = axis->GetXmax();
 500   Double_t *new_bins = new Double_t[bins + 1];
 501
 502   new_bins[0] = from;
 503   Double_t factor = pow(to/from, 1./bins);
 504
 505   for (int i = 1; i <= bins; i++) {
 506    new_bins[i] = factor * new_bins[i-1];
 507   }
 508   axis->Set(bins, new_bins);
 509   delete new_bins;
 510
 511 }
 512
 513 AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
 514 {
 515   //
 516   // Invert paramerameter  - not covariance yet
 517   //
 518   AliExternalTrackParam *output = new AliExternalTrackParam(*input);
 519   Double_t * param = (Double_t*)output->GetParameter();
 520   param[0]*=-1;
 521   param[3]*=-1;
 522   param[4]*=-1;
 523   //
 524   return output;
 525 }
 526
 527 /*
 528
 529 void Init(){
 530
 531 .x ~/UliStyle.C
 532 .x ~/rootlogon.C
 533 gSystem->Load("libSTAT.so");
 534 gSystem->Load("libANALYSIS");
 535 gSystem->Load("libTPCcalib");
 536 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 537
 538 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
 539 AliXRDPROOFtoolkit tool;
 540 TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
 541 chain->Lookup();
 542
 543 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01");  // OK
 544 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2");     // OK
 545 TCut cutP1("cutP1","abs(Tr0.fP[1]-Tr1.fP[1])<4");   // OK
 546 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 547 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
 548 TCut cutA=cutT+cutD+cutPt+cutN+cutP1;
 549
 550 TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
 551
 552 //
 553 chain->Draw(">>listEL",cutA,"entryList");
 554 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
 555 chain->SetEntryList(elist);
 556 //
 557 chain->Draw(">>listV40Z100","abs(d0)<40&&abs(v01)<100","entryList");
 558 TEntryList *elistV40Z100 = (TEntryList*)gDirectory->Get("listV40Z100");
 559 chain->SetEntryList(elistV40Z100);
 560
 561 //
 562 // Aliases
 563 //
 564 chain->SetAlias("side","(-1+(Tr0.fP[1]>0)*2)");
 565 chain->SetAlias("hpt","abs(Tr0.fP[4])<0.2");
 566 chain->SetAlias("dtheta","(Tr0.fP[3]+Tr1.fP[3])");
 567 chain->SetAlias("mtheta","(Tr0.fP[3]-Tr1.fP[3])*0.5")
 568 chain->SetAlias("sa","(sin(Tr0.fAlpha+0.))");
 569 chain->SetAlias("ca","(cos(Tr0.fAlpha+0.))");
 570
 571 }
 572
 573
 574 */
 575
 576
 577 /*
 578 void MatchTheta(){
 579
 580 TStatToolkit toolkit;
 581 Double_t chi2=0;
 582 Int_t    npoints=0;
 583 TVectorD fitParamA0;
 584 TVectorD fitParamA1;
 585 TVectorD fitParamC0;
 586 TVectorD fitParamC1;
 587 TMatrixD covMatrix;
 588
 589
 590 TString fstring="";
 591 //
 592 fstring+="ca++";
 593 fstring+="sa++";
 594 fstring+="ca*mtheta++";
 595 fstring+="sa*mtheta++";
 596 fstring+="side++";
 597 fstring+="side*ca++";
 598 fstring+="side*sa++";
 599 fstring+="side*ca*mtheta++";
 600 fstring+="side*sa*mtheta++";
 601
 602
 603 TString *strTheta0 = toolkit.FitPlane(chain,"dtheta",fstring->Data(), "hpt"+cutS, chi2,npoints,fitParamA0,covMatrix,0.8);
 604 chain->SetAlias("dtheta0",strTheta0.Data());
 605 strTheta0->Tokenize("+")->Print();
 606
 607
 608 //fstring+="mtheta++";
 609 //fstring+="mtheta^2++";
 610 //fstring+="ca*mtheta^2++";
 611 //fstring+="sa*mtheta^2++";
 612
 613
 614
 615 }
 616
 617 */
 618
 619
 620
 621
 622 /*
 623  void PosCorrection()
 624
 625
 626
 627
 628  TStatToolkit toolkit;
 629  Double_t chi2=0;
 630  Int_t    npoints=0;
 631  TVectorD fitParam;
 632  TMatrixD covMatrix;
 633
 634  //Theta
 635 chain->SetAlias("dthe","(Tr0.fP[3]+Tr1.fP[3])");
 636 chain->SetAlias("sign","(-1+(Tr0.fP[1]>0)*2)");
 637 chain->SetAlias("di","(1-abs(Tr0.fP[1])/250)");
 638 //
 639 //
 640 TString strFit="";
 641 //
 642 strFit+="sign++";                              //1
 643 strFit+="Tr0.fP[3]++";                         //2
 644 //
 645 strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]>0)++";     //3
 646 strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]<0)++";     //4
 647 //
 648 strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]>0)++";     //5
 649 strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]<0)++";     //6
 650 //
 651 strFit+="sin(Tr0.fAlpha)*Tr0.fP[3]++";         //7
 652 strFit+="cos(Tr0.fAlpha)*Tr0.fP[3]++";         //8
 653
 654
 655  //
 656  TString * thetaParam = toolkit.FitPlane(chain,"dthe", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix,0.8,0,10000)
 657
 658  chain->SetAlias("corTheta",thetaParam->Data());
 659  chain->Draw("dthe:Tr0.fP[1]","","",50000);
 660
 661
 662
 663 */
 664
 665
 666
 667 /*
 668
 669 void AliTPCcalibCosmic::dEdxCorrection(){
 670   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01");  // OK
 671   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2");     // OK
 672   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 673   TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
 674   TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
 675   TCut cutA=cutT+cutD+cutPt+cutN+cutS;
 676
 677
 678  .x ~/UliStyle.C
 679   gSystem->Load("libANALYSIS");
 680   gSystem->Load("libTPCcalib");
 681   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 682   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
 683  AliXRDPROOFtoolkit tool;
 684   TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
 685   chain->Lookup();
 686
 687   chain->Draw(">>listEL",cutA,"entryList");
 688   TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
 689   chain->SetEntryList(elist);
 690
 691   .x ~/rootlogon.C
 692    gSystem->Load("libSTAT.so");
 693    TStatToolkit toolkit;
 694   Double_t chi2=0;
 695   Int_t    npoints=0;
 696   TVectorD fitParam;
 697   TMatrixD covMatrix;
 698
 699   chain->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
 700
 701   TString strFit;
 702   strFit+="(Tr0.fP[1]/250)++";
 703   strFit+="(Tr0.fP[1]/250)^2++";
 704   strFit+="(Tr0.fP[3])++";
 705   strFit+="(Tr0.fP[3])^2++";
 706
 707   TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix)
 708
 709
 710
 711 */
 712
 713
 714 /*
 715 gSystem->Load("libANALYSIS");
 716 gSystem->Load("libSTAT");
 717 gSystem->Load("libTPCcalib");
 718
 719 TStatToolkit toolkit;
 720 Double_t chi2;
 721 TVectorD fitParam;
 722 TMatrixD covMatrix;
 723 Int_t npoints;
 724 //
 725 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03");  // OK
 726 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5");     // OK
 727 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 728 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
 729 TCut cutA="abs(norm-1)<0.3"+cutT+cutD+cutPt+cutN;
 730
 731
 732
 733 TTree * chain = Track0;
 734
 735
 736 chain->SetAlias("norm","signalTot0.fElements[3]/signalTot1.fElements[3]");
 737 //
 738 chain->SetAlias("dr1","(signalTot0.fElements[1]/signalTot0.fElements[3])");
 739 chain->SetAlias("dr2","(signalTot0.fElements[2]/signalTot0.fElements[3])");
 740 chain->SetAlias("dr4","(signalTot0.fElements[4]/signalTot0.fElements[3])");
 741 chain->SetAlias("dr5","(signalTot0.fElements[5]/signalTot0.fElements[3])");
 742
 743 TString fstring="";
 744 fstring+="dr1++";
 745 fstring+="dr2++";
 746 fstring+="dr4++";
 747 fstring+="dr5++";
 748 //
 749 fstring+="dr1*dr2++";
 750 fstring+="dr1*dr4++";
 751 fstring+="dr1*dr5++";
 752 fstring+="dr2*dr4++";
 753 fstring+="dr2*dr5++";
 754 fstring+="dr4*dr5++";
 755
 756
 757
 758 TString *strqdedx = toolkit.FitPlane(chain,"norm",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
 759
 760 chain->SetAlias("corQT",strqdedx->Data());
 761
 762 */
 763
 764