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     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  39     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  40     AliXRDPROOFtoolkit tool;
  41     TChain * chainCosmic = tool.MakeChainRandom("cosmic.txt","Track0",0,10000);
  42
  43
  44 */
  45
  46
  47
  48 #include "Riostream.h"
  49 #include "TChain.h"
  50 #include "TTree.h"
  51 #include "TH1F.h"
  52 #include "TH2F.h"
  53 #include "TList.h"
  54 #include "TMath.h"
  55 #include "TCanvas.h"
  56 #include "TFile.h"
  57 #include "TF1.h"
  58 #include "THnSparse.h"
  59
  60 #include "AliTPCclusterMI.h"
  61 #include "AliTPCseed.h"
  62 #include "AliESDVertex.h"
  63 #include "AliESDEvent.h"
  64 #include "AliESDfriend.h"
  65 #include "AliESDInputHandler.h"
  66 #include "AliAnalysisManager.h"
  67
  68 #include "AliTracker.h"
  69 #include "AliMagF.h"
  70 #include "AliTPCCalROC.h"
  71
  72 #include "AliLog.h"
  73
  74 #include "AliTPCcalibCosmic.h"
  75 #include "TTreeStream.h"
  76 #include "AliTPCTracklet.h"
  77 #include "AliESDcosmic.h"
  78
  79
  80 ClassImp(AliTPCcalibCosmic)
  81
  82
  83 AliTPCcalibCosmic::AliTPCcalibCosmic()
  84   :AliTPCcalibBase(),
  85    fHistNTracks(0),
  86    fClusters(0),
  87    fModules(0),
  88    fHistPt(0),
  89    fDeDx(0),
  90    fDeDxMIP(0),
  91    fMIPvalue(1),
  92    fCutMaxD(5),        // maximal distance in rfi ditection
  93    fCutMaxDz(40),      // maximal distance in z ditection
  94    fCutTheta(0.03),    // maximal distan theta
  95    fCutMinDir(-0.99)   // direction vector products
  96 {
  97   AliInfo("Default Constructor");
  98   for (Int_t ihis=0; ihis<6;ihis++){
  99     fHistoDelta[ihis]=0;
 100     fHistoPull[ihis]=0;
 101   }
 102   for (Int_t ihis=0; ihis<4;ihis++){
 103     fHistodEdxMax[ihis]    =0;
 104     fHistodEdxTot[ihis]    =0;
 105   }
 106 }
 107
 108
 109 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
 110   :AliTPCcalibBase(),
 111    fHistNTracks(0),
 112    fClusters(0),
 113    fModules(0),
 114    fHistPt(0),
 115    fDeDx(0),
 116    fDeDxMIP(0),
 117    fMIPvalue(1),
 118    fCutMaxD(5),        // maximal distance in rfi ditection
 119    fCutMaxDz(40),      // maximal distance in z ditection
 120    fCutTheta(0.03),    // maximal distan theta
 121    fCutMinDir(-0.99)   // direction vector products
 122 {
 123   SetName(name);
 124   SetTitle(title);
 125
 126   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
 127   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
 128   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
 129   fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
 130   fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
 131   BinLogX(fDeDx);
 132   fDeDxMIP =  new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
 133   Init();
 134   AliInfo("Non Default Constructor");
 135   //
 136 }
 137
 138 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
 139   //
 140   //
 141   //
 142   for (Int_t ihis=0; ihis<6;ihis++){
 143     delete fHistoDelta[ihis];
 144     delete fHistoPull[ihis];
 145   }
 146   for (Int_t ihis=0; ihis<4;ihis++){
 147     delete fHistodEdxTot[ihis];
 148     delete fHistodEdxMax[ihis];
 149   }
 150
 151   delete fHistNTracks;            //  histogram showing number of ESD tracks per event
 152   delete fClusters;               //  histogram showing the number of clusters per track
 153   delete fModules;                //  2d histogram of tracks which are propagated to the ACORDE scintillator array
 154   delete fHistPt;                 //  Pt histogram of reconstructed tracks
 155   delete fDeDx;                   //  dEdx spectrum showing the different particle types
 156   delete fDeDxMIP;                //  TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV
 157 }
 158
 159
 160 void AliTPCcalibCosmic::Init(){
 161   //
 162   // init component
 163   // Make performance histograms
 164   //
 165
 166   // tracking performance bins
 167   // 0 - delta of interest
 168   // 1 - min (track0, track1) number of clusters
 169   // 2 - R  - vertex radius
 170   // 3 - P1 - mean z
 171   // 4 - P2 - snp(phi)    at inner wall of TPC
 172   // 5 - P3 - tan(theta)  at inner wall of TPC
 173   // 6 - P4 - 1/pt mean
 174   // 7 - pt - pt mean
 175   // 8 - alpha
 176
 177   Double_t xminTrack[9], xmaxTrack[9];
 178   Int_t binsTrack[9];
 179   TString axisName[9];
 180   //
 181   binsTrack[0] =100;
 182   axisName[0]  ="#Delta";
 183   //
 184   binsTrack[1] =8;
 185   xminTrack[1] =80; xmaxTrack[1]=160;
 186   axisName[1]  ="N_{cl}";
 187   //
 188   binsTrack[2] =10;
 189   xminTrack[2] =0; xmaxTrack[2]=90;  //
 190   axisName[2]  ="dca_{r} (cm)";
 191   //
 192   binsTrack[3] =25;
 193   xminTrack[3] =-250; xmaxTrack[3]=250;  //
 194   axisName[3]  ="z (cm)";
 195   //
 196   binsTrack[4] =10;
 197   xminTrack[4] =-0.8; xmaxTrack[4]=0.8;  //
 198   axisName[4]  ="sin(#phi)";
 199   //
 200   binsTrack[5] =10;
 201   xminTrack[5] =-1; xmaxTrack[5]=1;  //
 202   axisName[5]  ="tan(#theta)";
 203   //
 204   binsTrack[6] =10;
 205   xminTrack[6] =0; xmaxTrack[6]=2;  //
 206   axisName[6]  ="1/pt (1/GeV)";
 207   //
 208   binsTrack[7] =40;
 209   xminTrack[7] =0.2; xmaxTrack[7]=50;  //
 210   axisName[7]  ="pt (GeV)";
 211   //
 212   binsTrack[8] =32;
 213   xminTrack[8] =0; xmaxTrack[8]=TMath::Pi();  //
 214   axisName[8]  ="alpha";
 215   //
 216   // delta y
 217   xminTrack[0] =-1; xmaxTrack[0]=1;  //
 218   fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
 219   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 220   fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 221   //
 222   // delta z
 223   xminTrack[0] =-1; xmaxTrack[0]=1;  //
 224   fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
 225   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 226   fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 227   //
 228   // delta P2
 229   xminTrack[0] =-10; xmaxTrack[0]=10;  //
 230   fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
 231   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 232   fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 233   //
 234   // delta P3
 235   xminTrack[0] =-10; xmaxTrack[0]=10;  //
 236   fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
 237   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 238   fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 239   //
 240   // delta P4
 241   xminTrack[0] =-0.2; xmaxTrack[0]=0.2;  //
 242   fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", 9, binsTrack,xminTrack, xmaxTrack);
 243   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 244   fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 245
 246   //
 247   // delta Pt
 248   xminTrack[0] =-0.5; xmaxTrack[0]=0.5;  //
 249   fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", 9, binsTrack,xminTrack, xmaxTrack);
 250   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 251   fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 252   //
 253
 254   for (Int_t idedx=0;idedx<4;idedx++){
 255     xminTrack[0] =0.5; xmaxTrack[0]=1.5;  //
 256     binsTrack[1] =40;
 257     xminTrack[1] =10; xmaxTrack[1]=160;
 258
 259     fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
 260                                           Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
 261                                           9, binsTrack,xminTrack, xmaxTrack);
 262     fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
 263                                           Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
 264                                           9, binsTrack,xminTrack, xmaxTrack);
 265   }
 266
 267
 268
 269   for (Int_t ivar=0;ivar<6;ivar++){
 270     for (Int_t ivar2=0;ivar2<9;ivar2++){
 271       fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 272       fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 273       fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 274       fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 275       BinLogX(fHistoDelta[ivar],7);
 276       BinLogX(fHistoPull[ivar],7);
 277       if (ivar<4){
 278         fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 279         fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 280         fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 281         fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 282         BinLogX(fHistodEdxMax[ivar],7);
 283         BinLogX(fHistodEdxTot[ivar],7);
 284       }
 285     }
 286   }
 287 }
 288
 289
 290 void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){
 291   //
 292   //
 293   //
 294   for (Int_t ivar=0; ivar<6;ivar++){
 295     if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){
 296       fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]);
 297     }
 298     if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){
 299       fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]);
 300     }
 301   }
 302   for (Int_t ivar=0; ivar<4;ivar++){
 303     if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){
 304       fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]);
 305     }
 306     if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){
 307       fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
 308     }
 309   }
 310 }
 311
 312
 313
 314
 315 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
 316   //
 317   //
 318   //
 319   if (!event) {
 320     Printf("ERROR: ESD not available");
 321     return;
 322   }
 323   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 324   if (!ESDfriend) {
 325    Printf("ERROR: ESDfriend not available");
 326    return;
 327   }
 328
 329
 330   FindPairs(event); // nearly everything takes place in find pairs...
 331
 332   if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
 333   Int_t ntracks=event->GetNumberOfTracks();
 334   fHistNTracks->Fill(ntracks);
 335   if (ntracks==0) return;
 336  //  AliESDcosmic cosmicESD;
 337 //   TTreeSRedirector * cstream =  GetDebugStreamer();
 338 //   cosmicESD.SetDebugStreamer(cstream);
 339 //   cosmicESD.ProcessEvent(event);
 340 //   if (cstream) cosmicESD.DumpToTree();
 341
 342
 343 }
 344
 345
 346 void AliTPCcalibCosmic::FillHistoPerformance(AliExternalTrackParam *par0, AliExternalTrackParam *par1, AliExternalTrackParam *inner0, AliExternalTrackParam *inner1, AliTPCseed *seed0,  AliTPCseed *seed1){
 347   //
 348   //
 349   //
 350   Int_t kMinCldEdx =20;
 351   Int_t ncl0 = seed0->GetNumberOfClusters();
 352   Int_t ncl1 = seed1->GetNumberOfClusters();
 353
 354   const Double_t kpullCut    = 10;
 355   Double_t x[9];
 356   Double_t xyz0[3];
 357   Double_t xyz1[3];
 358   par0->GetXYZ(xyz0);
 359   par1->GetXYZ(xyz1);
 360   Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
 361   Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
 362   inner0->GetXYZ(xyz0);
 363   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 364   // bin parameters
 365   x[1] = TMath::Min(ncl0,ncl1);
 366   x[2] = (radius0+radius1)*0.5;
 367   x[3] = (inner0->GetZ()+inner1->GetZ())*0.5;
 368   x[4] = (inner0->GetSnp()-inner1->GetSnp())*0.5;
 369   x[5] = (inner0->GetTgl()-inner1->GetTgl())*0.5;
 370   x[6] = (1/par0->Pt()+1/par1->Pt())*0.5;
 371   x[7] = (par0->Pt()+par1->Pt())*0.5;
 372   x[8] = alpha;
 373   // deltas
 374   Double_t delta[6];
 375   Double_t sigma[6];
 376   delta[0] = (par0->GetY()+par1->GetY());
 377   delta[1] = (par0->GetZ()-par1->GetZ());
 378   delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi());
 379   delta[3] = (par0->GetTgl()+par1->GetTgl());
 380   delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]);
 381   delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5);
 382   //
 383   sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2());
 384   sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2());
 385   sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2());
 386   sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2());
 387   sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2());
 388   sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5);
 389   //
 390   Bool_t isOK = kTRUE;
 391   for (Int_t ivar=0;ivar<6;ivar++){
 392     if (sigma[ivar]==0) isOK=kFALSE;
 393     x[0]= delta[ivar]/sigma[ivar];
 394     if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE;
 395   }
 396   //
 397
 398   if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
 399     x[0]= delta[ivar]/TMath::Sqrt(2);
 400     if (ivar==2 || ivar ==3) x[0]*=1000;
 401     fHistoDelta[ivar]->Fill(x);
 402     if (sigma[ivar]>0){
 403       x[0]= delta[ivar]/sigma[ivar];
 404       fHistoPull[ivar]->Fill(x);
 405     }
 406   }
 407
 408   //
 409   // Fill dedx performance
 410   //
 411   for (Int_t ipad=0; ipad<4;ipad++){
 412     //
 413     //
 414     //
 415     Int_t row0=0;
 416     Int_t row1=160;
 417     if (ipad==0) row1=63;
 418     if (ipad==1) {row0=63; row1=63+64;}
 419     if (ipad==2) {row0=128;}
 420     Int_t   nclUp       = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
 421     Int_t   nclDown     = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
 422     Int_t   minCl       = TMath::Min(nclUp,nclDown);
 423     if (minCl<kMinCldEdx) continue;
 424     x[1] = minCl;
 425     //
 426     Float_t dEdxTotUp   = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
 427     Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
 428     Float_t dEdxMaxUp   = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
 429     Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
 430     //
 431     if (dEdxTotDown<=0) continue;
 432     if (dEdxMaxDown<=0) continue;
 433     x[0]=dEdxTotUp/dEdxTotDown;
 434     fHistodEdxTot[ipad]->Fill(x);
 435     x[0]=dEdxMaxUp/dEdxMaxDown;
 436     fHistodEdxMax[ipad]->Fill(x);
 437   }
 438
 439
 440
 441 }
 442
 443
 444
 445 void AliTPCcalibCosmic::Analyze() {
 446
 447   fMIPvalue = CalculateMIPvalue(fDeDxMIP);
 448
 449   return;
 450
 451 }
 452
 453
 454
 455 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
 456   //
 457   // Find cosmic pairs
 458   //
 459   // Track0 is choosen in upper TPC part
 460   // Track1 is choosen in lower TPC part
 461   //
 462   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
 463   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 464   Int_t ntracks=event->GetNumberOfTracks();
 465   TObjArray  tpcSeeds(ntracks);
 466   if (ntracks==0) return;
 467   Double_t vtxx[3]={0,0,0};
 468   Double_t svtxx[3]={0.000001,0.000001,100.};
 469   AliESDVertex vtx(vtxx,svtxx);
 470   //
 471   //track loop
 472   //
 473   for (Int_t i=0;i<ntracks;++i) {
 474    AliESDtrack *track = event->GetTrack(i);
 475    fClusters->Fill(track->GetTPCNcls());
 476
 477    const AliExternalTrackParam * trackIn = track->GetInnerParam();
 478    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 479    if (!trackIn) continue;
 480    if (!trackOut) continue;
 481    if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue;  // filter laser
 482
 483
 484    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 485    TObject *calibObject;
 486    AliTPCseed *seed = 0;
 487    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 488      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 489    }
 490    if (seed) tpcSeeds.AddAt(seed,i);
 491
 492    Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 493    if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
 494      fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
 495      //
 496      if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
 497      //
 498      if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
 499        TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
 500        if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
 501        if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
 502      }
 503
 504    }
 505
 506   }
 507
 508   if (ntracks<2) return;
 509   //
 510   // Find pairs
 511   //
 512   for (Int_t i=0;i<ntracks;++i) {
 513     AliESDtrack *track0 = event->GetTrack(i);
 514     // track0 - choosen upper part
 515     if (!track0) continue;
 516     if (!track0->GetOuterParam()) continue;
 517     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 518     Double_t dir0[3];
 519     track0->GetDirection(dir0);
 520     for (Int_t j=0;j<ntracks;++j) {
 521       if (i==j) continue;
 522       AliESDtrack *track1 = event->GetTrack(j);
 523       //track 1 lower part
 524       if (!track1) continue;
 525       if (!track1->GetOuterParam()) continue;
 526       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 527       //
 528       Double_t dir1[3];
 529       track1->GetDirection(dir1);
 530
 531       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 532       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 533       if (! seed0) continue;
 534       if (! seed1) continue;
 535       Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
 536       Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
 537       //
 538       Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
 539       Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
 540       //
 541       Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
 542       Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
 543       //
 544       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
 545       Float_t d0  = track0->GetLinearD(0,0);
 546       Float_t d1  = track1->GetLinearD(0,0);
 547       //
 548       // conservative cuts - convergence to be guarantied
 549       // applying before track propagation
 550       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 551       if (dir>fCutMinDir) continue;               // direction vector product
 552       Float_t bz = AliTracker::GetBz();
 553       Float_t dvertex0[2];   //distance to 0,0
 554       Float_t dvertex1[2];   //distance to 0,0
 555       track0->GetDZ(0,0,0,bz,dvertex0);
 556       track1->GetDZ(0,0,0,bz,dvertex1);
 557       if (TMath::Abs(dvertex0[1])>250) continue;
 558       if (TMath::Abs(dvertex1[1])>250) continue;
 559       //
 560       //
 561       //
 562       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 563       AliExternalTrackParam param0(*track0);
 564       AliExternalTrackParam param1(*track1);
 565       //
 566       // Propagate using Magnetic field and correct fo material budget
 567       //
 568       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 569       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 570       //
 571       // Propagate rest to the 0,0 DCA - z should be ignored
 572       //
 573       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
 574       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
 575       //
 576       param0.GetDZ(0,0,0,bz,dvertex0);
 577       param1.GetDZ(0,0,0,bz,dvertex1);
 578       if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
 579       //
 580       Double_t xyz0[3];//,pxyz0[3];
 581       Double_t xyz1[3];//,pxyz1[3];
 582       param0.GetXYZ(xyz0);
 583       param1.GetXYZ(xyz1);
 584       Bool_t isPair = IsPair(&param0,&param1);
 585       //
 586       if (isPair) FillAcordeHist(track0);
 587       //
 588       // combined track params
 589       //
 590       AliExternalTrackParam *par0U=MakeCombinedTrack(&param0,&param1);
 591       AliExternalTrackParam *par1U=MakeCombinedTrack(&param1,&param0);
 592
 593
 594       //
 595       if (fStreamLevel>0){
 596         TTreeSRedirector * cstream =  GetDebugStreamer();
 597         //printf("My stream=%p\n",(void*)cstream);
 598         AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
 599         AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
 600         AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
 601         AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
 602         Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
 603         Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
 604         Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
 605         Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
 606         //
 607         //
 608         //
 609         FillHistoPerformance(&param0, &param1, ip0, ip1, seed0, seed1);
 610
 611         if (cstream) {
 612           (*cstream) << "Track0" <<
 613             "run="<<fRun<<              //  run number
 614             "event="<<fEvent<<          //  event number
 615             "time="<<fTime<<            //  time stamp of event
 616             "trigger="<<fTrigger<<      //  trigger
 617             "triggerClass="<<&fTriggerClass<<      //  trigger
 618             "mag="<<fMagF<<             //  magnetic field
 619             "dir="<<dir<<               //  direction
 620             "OK="<<isPair<<             //  will be accepted
 621             "b0="<<b0<<                 //  propagate status
 622             "b1="<<b1<<                 //  propagate status
 623             "crossI="<<isCrossI<<       //  cross inner
 624             "crossO="<<isCrossO<<       //  cross outer
 625             //
 626             "Orig0.=" << track0 <<      //  original track  0
 627             "Orig1.=" << track1 <<      //  original track  1
 628             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
 629             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0
 630             "Ip0.="<<ip0<<              //  inner param - upper
 631             "Ip1.="<<ip1<<              //  inner param - lower
 632             "Op0.="<<op0<<              //  outer param - upper
 633             "Op1.="<<op1<<              //  outer param - lower
 634             "Up0.="<<par0U<<           //  combined track 0
 635             "Up1.="<<par1U<<           //  combined track 1
 636             //
 637             "v00="<<dvertex0[0]<<       //  distance using kalman
 638             "v01="<<dvertex0[1]<<       //
 639             "v10="<<dvertex1[0]<<       //
 640             "v11="<<dvertex1[1]<<       //
 641             "d0="<<d0<<                 //  linear distance to 0,0
 642             "d1="<<d1<<                 //  linear distance to 0,0
 643             //
 644             //
 645             //
 646             "x00="<<xyz0[0]<<           // global position close to vertex
 647             "x01="<<xyz0[1]<<
 648             "x02="<<xyz0[2]<<
 649             //
 650             "x10="<<xyz1[0]<<           // global position close to vertex
 651             "x11="<<xyz1[1]<<
 652             "x12="<<xyz1[2]<<
 653             //
 654             "alpha0="<<alpha0<<
 655             "alpha1="<<alpha1<<
 656             "dir00="<<dir0[0]<<           // direction upper
 657             "dir01="<<dir0[1]<<
 658             "dir02="<<dir0[2]<<
 659             //
 660             "dir10="<<dir1[0]<<           // direction lower
 661             "dir11="<<dir1[1]<<
 662             "dir12="<<dir1[2]<<
 663             //
 664             //
 665             "Seed0.=" << seed0 <<       //  original seed 0
 666             "Seed1.=" << seed1 <<       //  original seed 1
 667             //
 668             "dedx0="<<dedx0<<           //  dedx0 - all
 669             "dedx1="<<dedx1<<           //  dedx1 - all
 670             //
 671             "dedx0I="<<dedx0I<<         //  dedx0 - inner ROC
 672             "dedx1I="<<dedx1I<<         //  dedx1 - inner ROC
 673             //
 674             "dedx0O="<<dedx0O<<         //  dedx0 - outer ROC
 675             "dedx1O="<<dedx1O<<         //  dedx1 - outer ROC
 676             "\n";
 677         }
 678       }
 679       delete par0U;
 680       delete par1U;
 681     }
 682   }
 683 }
 684
 685
 686
 687
 688 void  AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
 689
 690   // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
 691   if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
 692
 693   const Double_t AcordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
 694   const Double_t roof = 210.5;       // distance from x =0 to end of magnet roof
 695
 696   Double_t r[3];
 697   upperTrack->GetXYZ(r);
 698   Double_t d[3];
 699   upperTrack->GetDirection(d);
 700   Double_t x,z;
 701   z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
 702   x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
 703
 704   if (x > roof) {
 705     x = r[0] + (d[0]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
 706     z = r[2] + (d[2]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
 707   }
 708   if (x < -roof) {
 709     x = r[0] + (d[0]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
 710     z = r[2] + (d[2]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
 711   }
 712
 713   fModules->Fill(z, x);
 714
 715 }
 716
 717
 718
 719 Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
 720
 721   TIterator* iter = li->MakeIterator();
 722   AliTPCcalibCosmic* cal = 0;
 723
 724   while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
 725     if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
 726       //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 727       return -1;
 728     }
 729
 730     fHistNTracks->Add(cal->GetHistNTracks());
 731     fClusters->Add(cal-> GetHistClusters());
 732     fModules->Add(cal->GetHistAcorde());
 733     fHistPt->Add(cal->GetHistPt());
 734     fDeDx->Add(cal->GetHistDeDx());
 735     fDeDxMIP->Add(cal->GetHistMIP());
 736     Add(cal);
 737   }
 738   return 0;
 739
 740 }
 741
 742
 743 Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 744   //
 745   //
 746   /*
 747   // 0. Same direction - OPOSITE  - cutDir +cutT
 748   TCut cutDir("cutDir","dir<-0.99")
 749   // 1.
 750   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 751   //
 752   // 2. The same rphi
 753   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 754   //
 755   //
 756   //
 757   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 758   // 1/Pt diff cut
 759   */
 760   const Double_t *p0 = tr0->GetParameter();
 761   const Double_t *p1 = tr1->GetParameter();
 762   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 763   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
 764   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 765
 766   Double_t d0[3], d1[3];
 767   tr0->GetDirection(d0);
 768   tr1->GetDirection(d1);
 769   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 770   //
 771   return kTRUE;
 772 }
 773
 774
 775
 776 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
 777
 778   TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
 779   funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
 780                                 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
 781   hist->Fit(funcDoubleGaus);
 782   Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
 783
 784   delete funcDoubleGaus;
 785
 786   return MIPvalue;
 787
 788 }
 789
 790
 791
 792
 793 void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
 794   //
 795   // Not implemented yet
 796   //
 797   return;
 798
 799 }
 800
 801
 802 void AliTPCcalibCosmic::BinLogX(THnSparse *h, Int_t axisDim) {
 803
 804   // Method for the correct logarithmic binning of histograms
 805
 806   TAxis *axis = h->GetAxis(axisDim);
 807   int bins = axis->GetNbins();
 808
 809   Double_t from = axis->GetXmin();
 810   Double_t to = axis->GetXmax();
 811   Double_t *new_bins = new Double_t[bins + 1];
 812
 813   new_bins[0] = from;
 814   Double_t factor = pow(to/from, 1./bins);
 815
 816   for (int i = 1; i <= bins; i++) {
 817    new_bins[i] = factor * new_bins[i-1];
 818   }
 819   axis->Set(bins, new_bins);
 820   delete new_bins;
 821
 822 }
 823
 824
 825 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
 826
 827   // Method for the correct logarithmic binning of histograms
 828
 829   TAxis *axis = h->GetXaxis();
 830   int bins = axis->GetNbins();
 831
 832   Double_t from = axis->GetXmin();
 833   Double_t to = axis->GetXmax();
 834   Double_t *new_bins = new Double_t[bins + 1];
 835
 836   new_bins[0] = from;
 837   Double_t factor = pow(to/from, 1./bins);
 838
 839   for (int i = 1; i <= bins; i++) {
 840    new_bins[i] = factor * new_bins[i-1];
 841   }
 842   axis->Set(bins, new_bins);
 843   delete new_bins;
 844
 845 }
 846
 847
 848 AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
 849   //
 850   //
 851   //
 852   AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
 853   par1R->Rotate(track0->GetAlpha());
 854   par1R->PropagateTo(track0->GetX(),AliTracker::GetBz());
 855   //
 856   //
 857   Double_t * param = (Double_t*)par1R->GetParameter();
 858   Double_t * covar = (Double_t*)par1R->GetCovariance();
 859
 860   param[0]*=1;  //OK
 861   param[1]*=1;  //OK
 862   param[2]*=1;  //?
 863   param[3]*=-1; //OK
 864   param[4]*=-1; //OK
 865   //
 866   covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
 867   //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
 868   covar[13]*=-1.;
 869   return par1R;
 870 }
 871
 872 AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
 873   //
 874   // Make combined track
 875   //
 876   //
 877   AliExternalTrackParam * par1T = MakeTrack(track0,track1);
 878   AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
 879   //
 880   UpdateTrack(*par0U,*par1T);
 881   delete par1T;
 882   return par0U;
 883 }
 884
 885
 886 void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
 887   //
 888   // Update track 1 with track 2
 889   //
 890   //
 891   //
 892   TMatrixD vecXk(5,1);    // X vector
 893   TMatrixD covXk(5,5);    // X covariance
 894   TMatrixD matHk(5,5);    // vector to mesurement
 895   TMatrixD measR(5,5);    // measurement error
 896   TMatrixD vecZk(5,1);    // measurement
 897   //
 898   TMatrixD vecYk(5,1);    // Innovation or measurement residual
 899   TMatrixD matHkT(5,5);
 900   TMatrixD matSk(5,5);    // Innovation (or residual) covariance
 901   TMatrixD matKk(5,5);    // Optimal Kalman gain
 902   TMatrixD mat1(5,5);     // update covariance matrix
 903   TMatrixD covXk2(5,5);   //
 904   TMatrixD covOut(5,5);
 905   //
 906   Double_t *param1=(Double_t*) track1.GetParameter();
 907   Double_t *covar1=(Double_t*) track1.GetCovariance();
 908   Double_t *param2=(Double_t*) track2.GetParameter();
 909   Double_t *covar2=(Double_t*) track2.GetCovariance();
 910   //
 911   // copy data to the matrix
 912   for (Int_t ipar=0; ipar<5; ipar++){
 913     for (Int_t jpar=0; jpar<5; jpar++){
 914       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
 915       measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
 916       matHk(ipar,jpar)=0;
 917       mat1(ipar,jpar)=0;
 918     }
 919     vecXk(ipar,0) = param1[ipar];
 920     vecZk(ipar,0) = param2[ipar];
 921     matHk(ipar,ipar)=1;
 922     mat1(ipar,ipar)=0;
 923   }
 924   //
 925   //
 926   //
 927   //
 928   //
 929   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
 930   matHkT=matHk.T(); matHk.T();
 931   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
 932   matSk.Invert();
 933   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
 934   vecXk += matKk*vecYk;                      //  updated vector
 935   covXk2 = (mat1-(matKk*matHk));
 936   covOut =  covXk2*covXk;
 937   //
 938   //
 939   //
 940   // copy from matrix to parameters
 941   if (0) {
 942     vecXk.Print();
 943     vecZk.Print();
 944     //
 945     measR.Print();
 946     covXk.Print();
 947     covOut.Print();
 948     //
 949     track1.Print();
 950     track2.Print();
 951   }
 952
 953   for (Int_t ipar=0; ipar<5; ipar++){
 954     param1[ipar]= vecXk(ipar,0) ;
 955     for (Int_t jpar=0; jpar<5; jpar++){
 956       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
 957     }
 958   }
 959 }
 960
 961
 962
 963