TPC/AliTPCcalibCosmic.cxx

   1
   2
   3 /**************************************************************************
   4  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   5  *                                                                        *
   6  * Author: The ALICE Off-line Project.                                    *
   7  * Contributors are mentioned in the code where appropriate.              *
   8  *                                                                        *
   9  * Permission to use, copy, modify and distribute this software and its   *
  10  * documentation strictly for non-commercial purposes is hereby granted   *
  11  * without fee, provided that the above copyright notice appears in all   *
  12  * copies and that both the copyright notice and this permission notice   *
  13  * appear in the supporting documentation. The authors make no claims     *
  14  * about the suitability of this software for any purpose. It is          *
  15  * provided "as is" without express or implied warranty.                  *
  16  **************************************************************************/
  17
  18 /*
  19     Comments to be written here:
  20     1. What do we calibrate.
  21     2. How to interpret results
  22     3. Simple example
  23     4. Analysis using debug streamers.
  24
  25
  26
  27     3.Simple example
  28     // To make cosmic scan the user interaction neccessary
  29     //
  30     .x ~/UliStyle.C
  31     gSystem->Load("libANALYSIS");
  32     gSystem->Load("libTPCcalib");
  33     TFile fcalib("CalibObjects.root");
  34     TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  35     AliTPCcalibCosmic * cosmic = ( AliTPCcalibCosmic *)array->FindObject("cosmicTPC");
  36
  37
  38     //
  39     //
  40     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  41     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  42     AliXRDPROOFtoolkit tool;
  43     TChain * chainCosmic = tool.MakeChainRandom("cosmicF.txt","Track0",0,10000);
  44     //
  45     TCut cutITSN="min(Orig0.fITSncls,Orig1.fITSncls)>2";
  46     TCut cutTPCN="min(Orig0.fTPCncls,Orig1.fTPCncls)>120";
  47
  48     chainCosmic->Draw(">>listITS",cutITSN+cutTPCN,"entryList");
  49     TEntryList *elistITS = (TEntryList*)gDirectory->Get("listITS");
  50     chainCosmic->SetEntryList(elistITS);
  51
  52   */
  53
  54
  55
  56 #include "Riostream.h"
  57 #include "TChain.h"
  58 #include "TTree.h"
  59 #include "TH1F.h"
  60 #include "TH2F.h"
  61 #include "TList.h"
  62 #include "TMath.h"
  63 #include "TCanvas.h"
  64 #include "TFile.h"
  65 #include "TF1.h"
  66 #include "THnSparse.h"
  67 #include "TDatabasePDG.h"
  68
  69 #include "AliTPCclusterMI.h"
  70 #include "AliTPCseed.h"
  71 #include "AliESDVertex.h"
  72 #include "AliESDEvent.h"
  73 #include "AliESDfriend.h"
  74 #include "AliESDInputHandler.h"
  75 #include "AliAnalysisManager.h"
  76
  77 #include "AliTracker.h"
  78 #include "AliMagF.h"
  79 #include "AliTPCCalROC.h"
  80 #include "AliTPCParam.h"
  81 #include "AliLog.h"
  82
  83 #include "AliTPCcalibCosmic.h"
  84 #include "TTreeStream.h"
  85 #include "AliTPCTracklet.h"
  86 //#include "AliESDcosmic.h"
  87 #include "AliRecoParam.h"
  88 #include "AliMultiplicity.h"
  89 #include "AliTPCTransform.h"
  90 #include "AliTPCcalibDB.h"
  91 #include "AliTPCseed.h"
  92 #include "AliGRPObject.h"
  93 #include "AliTPCCorrection.h"
  94 ClassImp(AliTPCcalibCosmic)
  95
  96
  97 AliTPCcalibCosmic::AliTPCcalibCosmic()
  98   :AliTPCcalibBase(),
  99    fHistNTracks(0),
 100    fClusters(0),
 101    fModules(0),
 102    fHistPt(0),
 103    fDeDx(0),
 104    fDeDxMIP(0),
 105    fMIPvalue(1),
 106    fCutMaxD(5),        // maximal distance in rfi ditection
 107    fCutMaxDz(40),      // maximal distance in z ditection
 108    fCutTheta(0.03),    // maximal distan theta
 109    fCutMinDir(-0.99),   // direction vector products
 110    fCosmicTree(0)      // tree with cosmic data
 111 {
 112   AliInfo("Default Constructor");
 113   for (Int_t ihis=0; ihis<6;ihis++){
 114     fHistoDelta[ihis]=0;
 115     fHistoPull[ihis]=0;
 116   }
 117   for (Int_t ihis=0; ihis<4;ihis++){
 118     fHistodEdxMax[ihis]    =0;
 119     fHistodEdxTot[ihis]    =0;
 120   }
 121 }
 122
 123
 124 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
 125   :AliTPCcalibBase(),
 126    fHistNTracks(0),
 127    fClusters(0),
 128    fModules(0),
 129    fHistPt(0),
 130    fDeDx(0),
 131    fDeDxMIP(0),
 132    fMIPvalue(1),
 133    fCutMaxD(5),        // maximal distance in rfi ditection
 134    fCutMaxDz(40),      // maximal distance in z ditection
 135    fCutTheta(0.03),    // maximal distan theta
 136    fCutMinDir(-0.99),  // direction vector products
 137    fCosmicTree(0)      // tree with cosmic data
 138 {
 139   SetName(name);
 140   SetTitle(title);
 141
 142   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
 143   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
 144   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
 145   fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
 146   fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
 147   BinLogX(fDeDx);
 148   fDeDxMIP =  new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
 149   Init();
 150   AliInfo("Non Default Constructor");
 151   //
 152 }
 153
 154 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
 155   //
 156   //
 157   //
 158   for (Int_t ihis=0; ihis<6;ihis++){
 159     delete fHistoDelta[ihis];
 160     delete fHistoPull[ihis];
 161   }
 162   for (Int_t ihis=0; ihis<4;ihis++){
 163     delete fHistodEdxTot[ihis];
 164     delete fHistodEdxMax[ihis];
 165   }
 166
 167   delete fHistNTracks;            //  histogram showing number of ESD tracks per event
 168   delete fClusters;               //  histogram showing the number of clusters per track
 169   delete fModules;                //  2d histogram of tracks which are propagated to the ACORDE scintillator array
 170   delete fHistPt;                 //  Pt histogram of reconstructed tracks
 171   delete fDeDx;                   //  dEdx spectrum showing the different particle types
 172   delete fDeDxMIP;                //  TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV
 173 }
 174
 175
 176 void AliTPCcalibCosmic::Init(){
 177   //
 178   // init component
 179   // Make performance histograms
 180   //
 181
 182   // tracking performance bins
 183   // 0 - delta of interest
 184   // 1 - min (track0, track1) number of clusters
 185   // 2 - R  - vertex radius
 186   // 3 - P1 - mean z
 187   // 4 - P2 - snp(phi)    at inner wall of TPC
 188   // 5 - P3 - tan(theta)  at inner wall of TPC
 189   // 6 - P4 - 1/pt mean
 190   // 7 - pt - pt mean
 191   // 8 - alpha
 192   // 9 - is corss indicator
 193   Int_t ndim=10;
 194   Double_t xminTrack[10], xmaxTrack[10];
 195   Int_t binsTrack[10];
 196   TString axisName[10];
 197   //
 198   binsTrack[0] =100;
 199   axisName[0]  ="#Delta";
 200   //
 201   binsTrack[1] =8;
 202   xminTrack[1] =80; xmaxTrack[1]=160;
 203   axisName[1]  ="N_{cl}";
 204   //
 205   binsTrack[2] =10;
 206   xminTrack[2] =0; xmaxTrack[2]=90;  //
 207   axisName[2]  ="dca_{r} (cm)";
 208   //
 209   binsTrack[3] =25;
 210   xminTrack[3] =-250; xmaxTrack[3]=250;  //
 211   axisName[3]  ="z (cm)";
 212   //
 213   binsTrack[4] =10;
 214   xminTrack[4] =-0.8; xmaxTrack[4]=0.8;  //
 215   axisName[4]  ="sin(#phi)";
 216   //
 217   binsTrack[5] =10;
 218   xminTrack[5] =-1; xmaxTrack[5]=1;  //
 219   axisName[5]  ="tan(#theta)";
 220   //
 221   binsTrack[6] =40;
 222   xminTrack[6] =-2; xmaxTrack[6]=2;  //
 223   axisName[6]  ="1/pt (1/GeV)";
 224   //
 225   binsTrack[7] =50;
 226   xminTrack[7] =1; xmaxTrack[7]=1000;  //
 227   axisName[7]  ="pt (GeV)";
 228   //
 229   binsTrack[8] =18;
 230   xminTrack[8] =0; xmaxTrack[8]=TMath::Pi();  //
 231   axisName[8]  ="alpha";
 232   //
 233   binsTrack[9] =3;
 234   xminTrack[9] =-0.1; xmaxTrack[9]=2.1;  //
 235   axisName[9]  ="cross";
 236   //
 237   // delta y
 238   xminTrack[0] =-1; xmaxTrack[0]=1;  //
 239   fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
 240   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 241   fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
 242   //
 243   // delta z
 244   xminTrack[0] =-1; xmaxTrack[0]=1;  //
 245   fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
 246   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 247   fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
 248   //
 249   // delta P2
 250   xminTrack[0] =-10; xmaxTrack[0]=10;  //
 251   fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
 252   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 253   fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
 254   //
 255   // delta P3
 256   xminTrack[0] =-10; xmaxTrack[0]=10;  //
 257   fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
 258   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 259   fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
 260   //
 261   // delta P4
 262   xminTrack[0] =-0.2; xmaxTrack[0]=0.2;  //
 263   fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", ndim, binsTrack,xminTrack, xmaxTrack);
 264   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 265   fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
 266
 267   //
 268   // delta Pt
 269   xminTrack[0] =-0.5; xmaxTrack[0]=0.5;  //
 270   fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", ndim, binsTrack,xminTrack, xmaxTrack);
 271   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 272   fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
 273   //
 274
 275   for (Int_t idedx=0;idedx<4;idedx++){
 276     xminTrack[0] =0.5; xmaxTrack[0]=1.5;  //
 277     binsTrack[1] =40;
 278     xminTrack[1] =10; xmaxTrack[1]=160;
 279
 280     fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
 281                                           Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
 282                                           ndim, binsTrack,xminTrack, xmaxTrack);
 283     fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
 284                                           Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
 285                                           ndim, binsTrack,xminTrack, xmaxTrack);
 286   }
 287
 288
 289
 290   for (Int_t ivar=0;ivar<6;ivar++){
 291     for (Int_t ivar2=0;ivar2<ndim;ivar2++){
 292       fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 293       fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 294       fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 295       fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 296       BinLogX(fHistoDelta[ivar],7);
 297       BinLogX(fHistoPull[ivar],7);
 298       if (ivar<4){
 299         fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 300         fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 301         fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 302         fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 303         BinLogX(fHistodEdxMax[ivar],7);
 304         BinLogX(fHistodEdxTot[ivar],7);
 305       }
 306     }
 307   }
 308 }
 309
 310
 311 void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){
 312   //
 313   //
 314   //
 315   for (Int_t ivar=0; ivar<6;ivar++){
 316     if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){
 317       fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]);
 318     }
 319     if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){
 320       fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]);
 321     }
 322   }
 323   for (Int_t ivar=0; ivar<4;ivar++){
 324     if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){
 325       fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]);
 326     }
 327     if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){
 328       fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
 329     }
 330   }
 331   if (cosmic->fCosmicTree){
 332     if (!fCosmicTree) {
 333       fCosmicTree = new TTree("pairs","pairs");
 334       fCosmicTree->SetDirectory(0);
 335     }
 336     AliTPCcalibCosmic::AddTree(fCosmicTree,cosmic->fCosmicTree);
 337   }
 338 }
 339
 340
 341
 342
 343 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
 344   //
 345   //
 346   //
 347   if (!event) {
 348     Printf("ERROR: ESD not available");
 349     return;
 350   }
 351   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 352   if (!esdFriend) {
 353    Printf("ERROR: esdFriend not available");
 354    return;
 355   }
 356
 357   //
 358   //Int_t isOK=kTRUE;
 359   // COSMIC not signed properly
 360   //  UInt_t specie = event->GetEventSpecie();  // select only cosmic events
 361   //if (specie==AliRecoParam::kCosmic || specie==AliRecoParam::kCalib) {
 362   //  isOK = kTRUE;
 363   //}
 364   //if (!isOK) return;
 365   // Work around
 366   FindCosmicPairs(event);
 367   return;
 368   const AliMultiplicity *multiplicity = event->GetMultiplicity();
 369   Int_t ntracklets = multiplicity->GetNumberOfTracklets();
 370   if (ntracklets>6) return; // filter out "normal" event with high multiplicity
 371   const TString &trigger = event->GetFiredTriggerClasses();
 372   if (trigger.Contains("C0OB0")==0) return;
 373
 374
 375   FindPairs(event); // nearly everything takes place in find pairs...
 376
 377   if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
 378   Int_t ntracks=event->GetNumberOfTracks();
 379   fHistNTracks->Fill(ntracks);
 380
 381 }
 382
 383
 384 void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0,  AliTPCseed *seed1, const AliExternalTrackParam *param0Combined , Int_t cross){
 385   //
 386   // par0,par1       - parameter of tracks at DCA 0
 387   // inner0,inner1   - parameter of tracks at the TPC entrance
 388   // seed0, seed1    - detailed track information
 389   // param0Combined  - Use combined track parameters for binning
 390   //
 391   Int_t kMinCldEdx =20;
 392   Int_t ncl0 = seed0->GetNumberOfClusters();
 393   Int_t ncl1 = seed1->GetNumberOfClusters();
 394   const Double_t kpullCut    = 10;
 395   Double_t x[10];
 396   Double_t xyz0[3];
 397   Double_t xyz1[3];
 398   par0->GetXYZ(xyz0);
 399   par1->GetXYZ(xyz1);
 400   Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
 401   Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
 402   inner0->GetXYZ(xyz0);
 403   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 404   // bin parameters
 405   x[1] = TMath::Min(ncl0,ncl1);
 406   x[2] = (radius0+radius1)*0.5;
 407   x[3] = param0Combined->GetZ();
 408   x[4] = inner0->GetSnp();
 409   x[5] = param0Combined->GetTgl();
 410   x[6] = param0Combined->GetSigned1Pt();
 411   x[7] = param0Combined->Pt();
 412   x[8] = alpha;
 413   x[9] = cross;
 414   // deltas
 415   Double_t delta[6];
 416   Double_t sigma[6];
 417   delta[0] = (par0->GetY()+par1->GetY());
 418   delta[1] = (par0->GetZ()-par1->GetZ());
 419   delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi());
 420   delta[3] = (par0->GetTgl()+par1->GetTgl());
 421   delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]);
 422   delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5);
 423   //
 424   sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2());
 425   sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2());
 426   sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2());
 427   sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2());
 428   sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2());
 429   sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5);
 430   //
 431   Bool_t isOK = kTRUE;
 432   for (Int_t ivar=0;ivar<6;ivar++){
 433     if (sigma[ivar]==0) isOK=kFALSE;
 434     x[0]= delta[ivar]/sigma[ivar];
 435     if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE;
 436   }
 437   //
 438
 439   if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
 440     x[0]= delta[ivar];    // Modifiation 10.10 use not normalized deltas
 441     if (ivar==2 || ivar ==3) x[0]*=1000;  // angles in radian
 442     fHistoDelta[ivar]->Fill(x);
 443     if (sigma[ivar]>0){
 444       x[0]= delta[ivar]/sigma[ivar];
 445       fHistoPull[ivar]->Fill(x);
 446     }
 447   }
 448
 449   //
 450   // Fill dedx performance
 451   //
 452   for (Int_t ipad=0; ipad<4;ipad++){
 453     //
 454     //
 455     //
 456     Int_t row0=0;
 457     Int_t row1=160;
 458     if (ipad==0) row1=63;
 459     if (ipad==1) {row0=63; row1=63+64;}
 460     if (ipad==2) {row0=128;}
 461     Int_t   nclUp       = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
 462     Int_t   nclDown     = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
 463     Int_t   minCl       = TMath::Min(nclUp,nclDown);
 464     if (minCl<kMinCldEdx) continue;
 465     x[1] = minCl;
 466     //
 467     Float_t dEdxTotUp   = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
 468     Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
 469     Float_t dEdxMaxUp   = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
 470     Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
 471     //
 472     if (dEdxTotDown<=0) continue;
 473     if (dEdxMaxDown<=0) continue;
 474     x[0]=dEdxTotUp/dEdxTotDown;
 475     fHistodEdxTot[ipad]->Fill(x);
 476     x[0]=dEdxMaxUp/dEdxMaxDown;
 477     fHistodEdxMax[ipad]->Fill(x);
 478   }
 479
 480
 481
 482 }
 483
 484 void AliTPCcalibCosmic::MaterialBudgetDump(AliExternalTrackParam *const par0, AliExternalTrackParam *const par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam *inner1, AliTPCseed *const seed0,  AliTPCseed *const seed1, AliExternalTrackParam *const param0Combined, AliExternalTrackParam *const param1Combined){
 485   //
 486   // matrial budget AOD dump
 487   //
 488   // par0,par1       - parameter of tracks at DCA 0
 489   // inner0,inner1   - parameter of tracks at the TPC entrance
 490   // seed0, seed1    - detailed track information
 491   // param0Combined  - Use combined track parameters for binning
 492   // param1Combined  -
 493   Double_t p0In = inner0->GetP();
 494   Double_t p1In = inner1->GetP();
 495   Double_t p0V  = par0->GetP();
 496   Double_t p1V  = par1->GetP();
 497   //
 498   Double_t pt0In = inner0->Pt();
 499   Double_t pt1In = inner1->Pt();
 500   Double_t pt0V  = par0->Pt();
 501   Double_t pt1V  = par1->Pt();
 502   Int_t ncl0 = seed0->GetNumberOfClusters();
 503   Int_t ncl1 = seed1->GetNumberOfClusters();
 504   Int_t nclmin=TMath::Min(ncl0,ncl1);
 505   Double_t sign = (param0Combined->GetSigned1Pt()>0) ? 1:-1.;
 506   //
 507   TTreeSRedirector * pcstream =  GetDebugStreamer();
 508   if (pcstream){
 509     (*pcstream)<<"material"<<
 510       "run="<<fRun<<              //  run number
 511       "event="<<fEvent<<          //  event number
 512       "time="<<fTime<<            //  time stamp of event
 513       "trigger="<<fTrigger<<      //  trigger
 514       "triggerClass="<<&fTriggerClass<<      //  trigger
 515       "mag="<<fMagF<<             //  magnetic field
 516       "sign="<<sign<<             // sign of the track
 517       //
 518       "ncl0="<<ncl0<<
 519       "ncl1="<<ncl1<<
 520       "nclmin="<<nclmin<<
 521       //
 522       "p0In="<<p0In<<
 523       "p1In="<<p1In<<
 524       "p0V="<<p0V<<
 525       "p1V="<<p1V<<
 526       "pt0In="<<pt0In<<
 527       "pt1In="<<pt1In<<
 528       "pt0V="<<pt0V<<
 529       "pt1V="<<pt1V<<
 530       "p0.="<<par0<<
 531       "p1.="<<par1<<
 532       "up0.="<<param0Combined<<
 533       "up1.="<<param1Combined<<
 534       "\n";
 535   }
 536
 537 }
 538
 539 void AliTPCcalibCosmic::Analyze() {
 540
 541   fMIPvalue = CalculateMIPvalue(fDeDxMIP);
 542
 543   return;
 544
 545 }
 546
 547
 548
 549 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
 550   //
 551   // Find cosmic pairs
 552   //
 553   // Track0 is choosen in upper TPC part
 554   // Track1 is choosen in lower TPC part
 555   //
 556   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
 557   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 558   Int_t ntracks=event->GetNumberOfTracks();
 559   TObjArray  tpcSeeds(ntracks);
 560   if (ntracks==0) return;
 561   Double_t vtxx[3]={0,0,0};
 562   Double_t svtxx[3]={0.000001,0.000001,100.};
 563   AliESDVertex vtx(vtxx,svtxx);
 564   //
 565   //track loop
 566   //
 567   for (Int_t i=0;i<ntracks;++i) {
 568    AliESDtrack *track = event->GetTrack(i);
 569    fClusters->Fill(track->GetTPCNcls());
 570
 571    const AliExternalTrackParam * trackIn = track->GetInnerParam();
 572    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 573    if (!trackIn) continue;
 574    if (!trackOut) continue;
 575    if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue;  // filter laser
 576
 577
 578    AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
 579    if (!friendTrack) continue;
 580    TObject *calibObject;
 581    AliTPCseed *seed = 0;
 582    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 583      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 584    }
 585    if (seed) tpcSeeds.AddAt(seed,i);
 586
 587    Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 588    if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
 589      fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
 590      //
 591      if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
 592      //
 593      // if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
 594 //        //TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
 595 //        //if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
 596 //        // if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
 597 //      }
 598
 599    }
 600
 601   }
 602
 603   if (ntracks<2) return;
 604   //
 605   // Find pairs
 606   //
 607   for (Int_t i=0;i<ntracks;++i) {
 608     AliESDtrack *track0 = event->GetTrack(i);
 609     // track0 - choosen upper part
 610     if (!track0) continue;
 611     if (!track0->GetOuterParam()) continue;
 612     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 613     Double_t dir0[3];
 614     track0->GetDirection(dir0);
 615     for (Int_t j=0;j<ntracks;++j) {
 616       if (i==j) continue;
 617       AliESDtrack *track1 = event->GetTrack(j);
 618       //track 1 lower part
 619       if (!track1) continue;
 620       if (!track1->GetOuterParam()) continue;
 621       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 622       //
 623       Double_t dir1[3];
 624       track1->GetDirection(dir1);
 625
 626       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 627       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 628       if (! seed0) continue;
 629       if (! seed1) continue;
 630       Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
 631       Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
 632       //
 633       Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
 634       Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
 635       //
 636       Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
 637       Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
 638       //
 639       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
 640       Float_t d0  = track0->GetLinearD(0,0);
 641       Float_t d1  = track1->GetLinearD(0,0);
 642       //
 643       // conservative cuts - convergence to be guarantied
 644       // applying before track propagation
 645       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 646       if (dir>fCutMinDir) continue;               // direction vector product
 647       Float_t bz = AliTracker::GetBz();
 648       Float_t dvertex0[2];   //distance to 0,0
 649       Float_t dvertex1[2];   //distance to 0,0
 650       track0->GetDZ(0,0,0,bz,dvertex0);
 651       track1->GetDZ(0,0,0,bz,dvertex1);
 652       if (TMath::Abs(dvertex0[1])>250) continue;
 653       if (TMath::Abs(dvertex1[1])>250) continue;
 654       //
 655       //
 656       //
 657       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 658       AliExternalTrackParam param0(*track0);
 659       AliExternalTrackParam param1(*track1);
 660       //
 661       // Propagate using Magnetic field and correct fo material budget
 662       //
 663       Double_t sign0=-1;
 664       Double_t sign1=1;
 665       Double_t maxsnp=0.90;
 666       AliTracker::PropagateTrackToBxByBz(&param0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign0);
 667       AliTracker::PropagateTrackToBxByBz(&param1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign1);
 668       //
 669       // Propagate rest to the 0,0 DCA - z should be ignored
 670       //
 671       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
 672       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
 673       //
 674       param0.GetDZ(0,0,0,bz,dvertex0);
 675       param1.GetDZ(0,0,0,bz,dvertex1);
 676       if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
 677       //
 678       Double_t xyz0[3];//,pxyz0[3];
 679       Double_t xyz1[3];//,pxyz1[3];
 680       param0.GetXYZ(xyz0);
 681       param1.GetXYZ(xyz1);
 682       Bool_t isPair = IsPair(&param0,&param1);
 683       //
 684       if (isPair) FillAcordeHist(track0);
 685       if (isPair &&param0.Pt()>1) {
 686         const TString &trigger = event->GetFiredTriggerClasses();
 687         UInt_t specie = event->GetEventSpecie();
 688         printf("COSMIC ?\t%s\t%d\t%f\t%f\n", trigger.Data(),specie, param0.GetZ(), param1.GetZ());
 689       }
 690       //
 691       // combined track params
 692       //
 693       AliExternalTrackParam *par0U=MakeCombinedTrack(&param0,&param1);
 694       AliExternalTrackParam *par1U=MakeCombinedTrack(&param1,&param0);
 695
 696       //
 697       if (fStreamLevel>0){
 698         TTreeSRedirector * cstream =  GetDebugStreamer();
 699         //printf("My stream=%p\n",(void*)cstream);
 700         AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
 701         AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
 702         AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
 703         AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
 704         Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
 705         Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
 706         Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
 707         Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
 708         //
 709         //
 710         //
 711         Int_t cross =0;  // 0 no cross, 2 cross on both sides
 712         if (isCrossI) cross+=1;
 713         if (isCrossO) cross+=1;
 714         FillHistoPerformance(&param0, &param1, ip0, ip1, seed0, seed1,par0U, cross);
 715         MaterialBudgetDump(&param0, &param1, ip0, ip1, seed0, seed1,par0U,par1U);
 716         if (cstream) {
 717           (*cstream) << "Track0" <<
 718             "run="<<fRun<<              //  run number
 719             "event="<<fEvent<<          //  event number
 720             "time="<<fTime<<            //  time stamp of event
 721             "trigger="<<fTrigger<<      //  trigger
 722             "triggerClass="<<&fTriggerClass<<      //  trigger
 723             "mag="<<fMagF<<             //  magnetic field
 724             "dir="<<dir<<               //  direction
 725             "OK="<<isPair<<             //  will be accepted
 726             "b0="<<b0<<                 //  propagate status
 727             "b1="<<b1<<                 //  propagate status
 728             "crossI="<<isCrossI<<       //  cross inner
 729             "crossO="<<isCrossO<<       //  cross outer
 730             //
 731             "Orig0.=" << track0 <<      //  original track  0
 732             "Orig1.=" << track1 <<      //  original track  1
 733             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
 734             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0
 735             "Ip0.="<<ip0<<              //  inner param - upper
 736             "Ip1.="<<ip1<<              //  inner param - lower
 737             "Op0.="<<op0<<              //  outer param - upper
 738             "Op1.="<<op1<<              //  outer param - lower
 739             "Up0.="<<par0U<<           //  combined track 0
 740             "Up1.="<<par1U<<           //  combined track 1
 741             //
 742             "v00="<<dvertex0[0]<<       //  distance using kalman
 743             "v01="<<dvertex0[1]<<       //
 744             "v10="<<dvertex1[0]<<       //
 745             "v11="<<dvertex1[1]<<       //
 746             "d0="<<d0<<                 //  linear distance to 0,0
 747             "d1="<<d1<<                 //  linear distance to 0,0
 748             //
 749             //
 750             //
 751             "x00="<<xyz0[0]<<           // global position close to vertex
 752             "x01="<<xyz0[1]<<
 753             "x02="<<xyz0[2]<<
 754             //
 755             "x10="<<xyz1[0]<<           // global position close to vertex
 756             "x11="<<xyz1[1]<<
 757             "x12="<<xyz1[2]<<
 758             //
 759             "alpha0="<<alpha0<<
 760             "alpha1="<<alpha1<<
 761             "dir00="<<dir0[0]<<           // direction upper
 762             "dir01="<<dir0[1]<<
 763             "dir02="<<dir0[2]<<
 764             //
 765             "dir10="<<dir1[0]<<           // direction lower
 766             "dir11="<<dir1[1]<<
 767             "dir12="<<dir1[2]<<
 768             //
 769             //
 770             "Seed0.=" << seed0 <<       //  original seed 0
 771             "Seed1.=" << seed1 <<       //  original seed 1
 772             //
 773             "dedx0="<<dedx0<<           //  dedx0 - all
 774             "dedx1="<<dedx1<<           //  dedx1 - all
 775             //
 776             "dedx0I="<<dedx0I<<         //  dedx0 - inner ROC
 777             "dedx1I="<<dedx1I<<         //  dedx1 - inner ROC
 778             //
 779             "dedx0O="<<dedx0O<<         //  dedx0 - outer ROC
 780             "dedx1O="<<dedx1O<<         //  dedx1 - outer ROC
 781             "\n";
 782         }
 783       }
 784       delete par0U;
 785       delete par1U;
 786     }
 787   }
 788 }
 789
 790
 791
 792
 793 void  AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
 794
 795   // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
 796   if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
 797
 798   const Double_t acordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
 799   const Double_t roof = 210.5;       // distance from x =0 to end of magnet roof
 800
 801   Double_t r[3];
 802   upperTrack->GetXYZ(r);
 803   Double_t d[3];
 804   upperTrack->GetDirection(d);
 805   Double_t x,z;
 806   z = r[2] + (d[2]/d[1])*(acordePlane - r[1]);
 807   x = r[0] + (d[0]/d[1])*(acordePlane - r[1]);
 808
 809   if (x > roof) {
 810     x = r[0] + (d[0]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
 811     z = r[2] + (d[2]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
 812   }
 813   if (x < -roof) {
 814     x = r[0] + (d[0]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
 815     z = r[2] + (d[2]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
 816   }
 817
 818   fModules->Fill(z, x);
 819
 820 }
 821
 822
 823
 824 Long64_t AliTPCcalibCosmic::Merge(TCollection *const li) {
 825
 826   TIterator* iter = li->MakeIterator();
 827   AliTPCcalibCosmic* cal = 0;
 828
 829   while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
 830     if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
 831       //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 832       return -1;
 833     }
 834
 835     fHistNTracks->Add(cal->GetHistNTracks());
 836     fClusters->Add(cal-> GetHistClusters());
 837     fModules->Add(cal->GetHistAcorde());
 838     fHistPt->Add(cal->GetHistPt());
 839     fDeDx->Add(cal->GetHistDeDx());
 840     fDeDxMIP->Add(cal->GetHistMIP());
 841     Add(cal);
 842   }
 843   return 0;
 844
 845 }
 846
 847
 848 Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 849   //
 850   //
 851   /*
 852   // 0. Same direction - OPOSITE  - cutDir +cutT
 853   TCut cutDir("cutDir","dir<-0.99")
 854   // 1.
 855   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 856   //
 857   // 2. The same rphi
 858   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 859   //
 860   //
 861   //
 862   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 863   // 1/Pt diff cut
 864   */
 865   const Double_t *p0 = tr0->GetParameter();
 866   const Double_t *p1 = tr1->GetParameter();
 867   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 868   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
 869   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 870
 871   Double_t d0[3], d1[3];
 872   tr0->GetDirection(d0);
 873   tr1->GetDirection(d1);
 874   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 875   //
 876   return kTRUE;
 877 }
 878
 879
 880
 881 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
 882
 883   TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
 884   funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
 885                                 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
 886   hist->Fit(funcDoubleGaus);
 887   Double_t aMIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
 888
 889   delete funcDoubleGaus;
 890
 891   return aMIPvalue;
 892
 893 }
 894
 895
 896
 897
 898 void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
 899   //
 900   // Not implemented yet
 901   //
 902   return;
 903
 904 }
 905
 906
 907 void AliTPCcalibCosmic::BinLogX(THnSparse *const h, Int_t axisDim) {
 908
 909   // Method for the correct logarithmic binning of histograms
 910
 911   TAxis *axis = h->GetAxis(axisDim);
 912   int bins = axis->GetNbins();
 913
 914   Double_t from = axis->GetXmin();
 915   Double_t to = axis->GetXmax();
 916   Double_t *newBins = new Double_t[bins + 1];
 917
 918   newBins[0] = from;
 919   Double_t factor = pow(to/from, 1./bins);
 920
 921   for (int i = 1; i <= bins; i++) {
 922    newBins[i] = factor * newBins[i-1];
 923   }
 924   axis->Set(bins, newBins);
 925   delete [] newBins;
 926
 927 }
 928
 929
 930 void AliTPCcalibCosmic::BinLogX(TH1 *const h) {
 931
 932   // Method for the correct logarithmic binning of histograms
 933
 934   TAxis *axis = h->GetXaxis();
 935   int bins = axis->GetNbins();
 936
 937   Double_t from = axis->GetXmin();
 938   Double_t to = axis->GetXmax();
 939   Double_t *newBins = new Double_t[bins + 1];
 940
 941   newBins[0] = from;
 942   Double_t factor = pow(to/from, 1./bins);
 943
 944   for (int i = 1; i <= bins; i++) {
 945    newBins[i] = factor * newBins[i-1];
 946   }
 947   axis->Set(bins, newBins);
 948   delete [] newBins;
 949
 950 }
 951
 952
 953 AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
 954   //
 955   //
 956   //
 957   AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
 958   par1R->Rotate(track0->GetAlpha());
 959   par1R->PropagateTo(track0->GetX(),AliTracker::GetBz());
 960   //
 961   //
 962   Double_t * param = (Double_t*)par1R->GetParameter();
 963   Double_t * covar = (Double_t*)par1R->GetCovariance();
 964
 965   param[0]*=1;  //OK
 966   param[1]*=1;  //OK
 967   param[2]*=1;  //?
 968   param[3]*=-1; //OK
 969   param[4]*=-1; //OK
 970   //
 971   covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
 972   //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
 973   covar[13]*=-1.;
 974   return par1R;
 975 }
 976
 977 AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
 978   //
 979   // Make combined track
 980   //
 981   //
 982   AliExternalTrackParam * par1T = MakeTrack(track0,track1);
 983   AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
 984   //
 985   UpdateTrack(*par0U,*par1T);
 986   delete par1T;
 987   return par0U;
 988 }
 989
 990
 991 void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
 992   //
 993   // Update track 1 with track 2
 994   //
 995   //
 996   //
 997   TMatrixD vecXk(5,1);    // X vector
 998   TMatrixD covXk(5,5);    // X covariance
 999   TMatrixD matHk(5,5);    // vector to mesurement
1000   TMatrixD measR(5,5);    // measurement error
1001   TMatrixD vecZk(5,1);    // measurement
1002   //
1003   TMatrixD vecYk(5,1);    // Innovation or measurement residual
1004   TMatrixD matHkT(5,5);
1005   TMatrixD matSk(5,5);    // Innovation (or residual) covariance
1006   TMatrixD matKk(5,5);    // Optimal Kalman gain
1007   TMatrixD mat1(5,5);     // update covariance matrix
1008   TMatrixD covXk2(5,5);   //
1009   TMatrixD covOut(5,5);
1010   //
1011   Double_t *param1=(Double_t*) track1.GetParameter();
1012   Double_t *covar1=(Double_t*) track1.GetCovariance();
1013   Double_t *param2=(Double_t*) track2.GetParameter();
1014   Double_t *covar2=(Double_t*) track2.GetCovariance();
1015   //
1016   // copy data to the matrix
1017   for (Int_t ipar=0; ipar<5; ipar++){
1018     for (Int_t jpar=0; jpar<5; jpar++){
1019       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
1020       measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
1021       matHk(ipar,jpar)=0;
1022       mat1(ipar,jpar)=0;
1023     }
1024     vecXk(ipar,0) = param1[ipar];
1025     vecZk(ipar,0) = param2[ipar];
1026     matHk(ipar,ipar)=1;
1027     mat1(ipar,ipar)=0;
1028   }
1029   //
1030   //
1031   //
1032   //
1033   //
1034   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
1035   matHkT=matHk.T(); matHk.T();
1036   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
1037   matSk.Invert();
1038   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
1039   vecXk += matKk*vecYk;                      //  updated vector
1040   covXk2 = (mat1-(matKk*matHk));
1041   covOut =  covXk2*covXk;
1042   //
1043   //
1044   //
1045   // copy from matrix to parameters
1046   if (0) {
1047     vecXk.Print();
1048     vecZk.Print();
1049     //
1050     measR.Print();
1051     covXk.Print();
1052     covOut.Print();
1053     //
1054     track1.Print();
1055     track2.Print();
1056   }
1057
1058   for (Int_t ipar=0; ipar<5; ipar++){
1059     param1[ipar]= vecXk(ipar,0) ;
1060     for (Int_t jpar=0; jpar<5; jpar++){
1061       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
1062     }
1063   }
1064 }
1065
1066
1067
1068 void AliTPCcalibCosmic::FindCosmicPairs(AliESDEvent * event){
1069   //
1070   // find cosmic pairs trigger by random trigger
1071   //
1072   //
1073   AliESDVertex *vertexSPD =  (AliESDVertex *)event->GetPrimaryVertexSPD();
1074   AliESDVertex *vertexTPC =  (AliESDVertex *)event->GetPrimaryVertexTPC();
1075   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1076   const Double_t kMinPt=1;
1077   const Double_t kMinPtMax=0.8;
1078   const Double_t kMinNcl=50;
1079   const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
1080   Int_t ntracks=event->GetNumberOfTracks();
1081   //  Float_t dcaTPC[2]={0,0};
1082   // Float_t covTPC[3]={0,0,0};
1083
1084   UInt_t specie = event->GetEventSpecie();  // skip laser events
1085   if (specie==AliRecoParam::kCalib) return;
1086
1087
1088
1089   for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
1090     AliESDtrack *track0 = event->GetTrack(itrack0);
1091     if (!track0) continue;
1092     if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
1093
1094     if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()<kMinPt) continue;
1095     if (track0->GetTPCncls()<kMinNcl) continue;
1096     if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue;
1097     if (track0->GetKinkIndex(0)>0) continue;
1098     const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
1099     //rm primaries
1100     //
1101     //track0->GetImpactParametersTPC(dcaTPC,covTPC);
1102     //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
1103     //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
1104     //    const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
1105     for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
1106       AliESDtrack *track1 = event->GetTrack(itrack1);
1107       if (!track1) continue;
1108       if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
1109       if (track1->GetKinkIndex(0)>0) continue;
1110       if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()<kMinPt) continue;
1111       if (track1->GetTPCncls()<kMinNcl) continue;
1112       if (TMath::Abs(AliTracker::GetBz())>1&&TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
1113       if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
1114       //track1->GetImpactParametersTPC(dcaTPC,covTPC);
1115       //      if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
1116       //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
1117       //
1118       const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
1119       //
1120       Bool_t isPair=kTRUE;
1121       for (Int_t ipar=0; ipar<5; ipar++){
1122         if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
1123         if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
1124       }
1125       if (!isPair) continue;
1126       if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
1127       //delta with correct sign
1128       /*
1129         TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
1130         TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
1131         TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
1132       */
1133       if  (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y   opposite sign
1134       if  (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
1135       if  (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
1136       if (!isPair) continue;
1137       //      const AliExternalTrackParam * trackIn1 = track1->GetInnerParam();
1138       //
1139       //
1140       TTreeSRedirector * pcstream =  GetDebugStreamer();
1141       Int_t ntracksSPD = vertexSPD->GetNContributors();
1142       Int_t ntracksTPC = vertexTPC->GetNContributors();
1143       //
1144       AliESDfriendTrack *friendTrack0 = esdFriend->GetTrack(itrack0);
1145       if (!friendTrack0) continue;
1146       AliESDfriendTrack *friendTrack1 = esdFriend->GetTrack(itrack1);
1147       if (!friendTrack1) continue;
1148       TObject *calibObject;
1149       AliTPCseed *seed0 = 0;
1150       AliTPCseed *seed1 = 0;
1151       //
1152       for (Int_t l=0;(calibObject=friendTrack0->GetCalibObject(l));++l) {
1153         if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1154       }
1155       for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) {
1156         if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1157       }
1158       //
1159       if (pcstream){
1160         (*pcstream)<<"pairs"<<
1161           "run="<<fRun<<              //  run number
1162           "event="<<fEvent<<          //  event number
1163           "time="<<fTime<<            //  time stamp of event
1164           "trigger="<<fTrigger<<      //  trigger
1165           "triggerClass="<<&fTriggerClass<<      //  trigger
1166           "mag="<<fMagF<<             //  magnetic field
1167           //
1168           "nSPD="<<ntracksSPD<<
1169           "nTPC="<<ntracksTPC<<
1170           "vSPD.="<<vertexSPD<<         //primary vertex -SPD
1171           "vTPC.="<<vertexTPC<<         //primary vertex -TPC
1172           "t0.="<<track0<<              //track0
1173           "t1.="<<track1<<              //track1
1174           "ft0.="<<friendTrack0<<       //track0
1175           "ft1.="<<friendTrack1<<       //track1
1176           "s0.="<<seed0<<               //track0
1177           "s1.="<<seed1<<               //track1
1178           "\n";
1179       }
1180       if (!fCosmicTree) {
1181         fCosmicTree = new TTree("pairs","pairs");
1182         fCosmicTree->SetDirectory(0);
1183       }
1184       if (fCosmicTree->GetEntries()==0){
1185         //
1186         fCosmicTree->SetDirectory(0);
1187         fCosmicTree->Branch("t0.",&track0);
1188         fCosmicTree->Branch("t1.",&track1);
1189         fCosmicTree->Branch("ft0.",&friendTrack0);
1190         fCosmicTree->Branch("ft1.",&friendTrack1);
1191       }else{
1192         fCosmicTree->SetBranchAddress("t0.",&track0);
1193         fCosmicTree->SetBranchAddress("t1.",&track1);
1194         fCosmicTree->SetBranchAddress("ft0.",&friendTrack0);
1195         fCosmicTree->SetBranchAddress("ft1.",&friendTrack1);
1196       }
1197       fCosmicTree->Fill();
1198     }
1199   }
1200 }
1201
1202
1203 void  AliTPCcalibCosmic::Terminate(){
1204   //
1205   // copy the cosmic tree to memory resident tree
1206   //
1207   static Int_t counter=0;
1208   printf("AliTPCcalibCosmic::Terminate\t%d\n",counter);
1209   counter++;
1210   AliTPCcalibBase::Terminate();
1211 }
1212
1213
1214 void AliTPCcalibCosmic::AddTree(TTree* treeOutput, TTree * treeInput){
1215   //
1216   // Add the content of tree:
1217   // Notice automatic copy of tree in ROOT does not work for such complicated tree
1218   //
1219   AliESDtrack *track0=new AliESDtrack;
1220   AliESDtrack *track1=new AliESDtrack;
1221   AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
1222   AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
1223   treeInput->SetBranchAddress("t0.",&track0);
1224   treeInput->SetBranchAddress("t1.",&track1);
1225   treeInput->SetBranchAddress("ft0.",&ftrack0);
1226   treeInput->SetBranchAddress("ft1.",&ftrack1);
1227   if (treeOutput->GetEntries()==0){
1228     //
1229     treeOutput->SetDirectory(0);
1230     treeOutput->Branch("t0.",&track0);
1231     treeOutput->Branch("t1.",&track1);
1232     treeOutput->Branch("ft0.",&ftrack0);
1233     treeOutput->Branch("ft1.",&ftrack1);
1234   }else{
1235      treeOutput->SetBranchAddress("t0.",&track0);
1236      treeOutput->SetBranchAddress("t1.",&track1);
1237      treeOutput->SetBranchAddress("ft0.",&ftrack0);
1238      treeOutput->SetBranchAddress("ft1.",&ftrack1);
1239   }
1240   Int_t entries= treeInput->GetEntries();
1241   for (Int_t i=0; i<entries; i++){
1242     treeInput->GetEntry(i);
1243     treeOutput->Fill();
1244   }
1245 }
1246
1247
1248
1249 void AliTPCcalibCosmic::MakeFitTree(TTree * treeInput, TTreeSRedirector *pcstream, const TObjArray * corrArray, Int_t step, Int_t run){
1250   //
1251   // Make fit tree
1252   // refit the tracks with original points + corrected points for each correction
1253   // Input:
1254   //   treeInput - tree with cosmic tracks
1255   //   pcstream  - debug output
1256
1257   // Algorithm:
1258   // Loop over pair of cosmic tracks:
1259   //   1. Find trigger offset between cosmic event and "physic" trigger
1260   //   2. Refit tracks with current transformation
1261   //   3. Refit tracks using additional "primitive" distortion on top
1262   // Best correction estimated as linear combination of corrections
1263   // minimizing the observed distortions
1264   // Observed distortions - matching in the y,z, snp, theta and 1/pt
1265   //
1266   const Double_t kResetCov=20.;
1267   const Double_t kMaxDelta[5]={1,40,0.03,0.01,0.2};
1268   const Double_t kSigma=2.;
1269   const Double_t kMaxTime=1050;
1270   const Double_t kMaxSnp=0.8;
1271   Int_t ncorr=corrArray->GetEntries();
1272   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1273   AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
1274   AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
1275   Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
1276   transform->SetCurrentRun(run);
1277   transform->SetCurrentTimeStamp(TMath::Nint(time));
1278   Double_t covar[15];
1279   for (Int_t i=0;i<15;i++) covar[i]=0;
1280   covar[0]=kSigma*kSigma;
1281   covar[2]=kSigma*kSigma;
1282   covar[5]=kSigma*kSigma/Float_t(150*150);
1283   covar[9]=kSigma*kSigma/Float_t(150*150);
1284   covar[14]=0.2*0.2;
1285   Double_t *distortions = new Double_t[ncorr+1];
1286
1287   AliESDtrack *track0=new AliESDtrack;
1288   AliESDtrack *track1=new AliESDtrack;
1289   AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
1290   AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
1291   treeInput->SetBranchAddress("t0.",&track0);
1292   treeInput->SetBranchAddress("t1.",&track1);
1293   treeInput->SetBranchAddress("ft0.",&ftrack0);
1294   treeInput->SetBranchAddress("ft1.",&ftrack1);
1295   Int_t entries= treeInput->GetEntries();
1296   for (Int_t i=0; i<entries; i+=step){
1297     treeInput->GetEntry(i);
1298     if (i%20==0) printf("%d\n",i);
1299     if (!ftrack0->GetTPCOut()) continue;
1300     if (!ftrack1->GetTPCOut()) continue;
1301     AliTPCseed *seed0=0;
1302     AliTPCseed *seed1=0;
1303     TObject *calibObject;
1304     for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
1305       if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1306     }
1307     for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
1308       if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1309     }
1310     if (!seed0) continue;
1311     if (!seed1) continue;
1312     if (TMath::Abs(seed0->GetSnp())>kMaxSnp) continue;
1313     if (TMath::Abs(seed1->GetSnp())>kMaxSnp) continue;
1314     //
1315     //
1316     Int_t nclA0=0, nclC0=0;     // number of clusters
1317     Int_t nclA1=0, nclC1=0;     // number of clusters
1318     Int_t ncl0=0,ncl1=0;
1319     Double_t rmin0=300, rmax0=-300;  // variables to estimate the time0 - trigger offset
1320     Double_t rmin1=300, rmax1=-300;
1321     Double_t tmin0=2000, tmax0=-2000;
1322     Double_t tmin1=2000, tmax1=-2000;
1323     //
1324     //
1325     // calculate trigger offset usig "missing clusters"
1326     for (Int_t irow=0; irow<159; irow++){
1327       AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1328       if (cluster0 &&cluster0->GetX()>10){
1329         if (cluster0->GetX()<rmin0) { rmin0=cluster0->GetX(); tmin0=cluster0->GetTimeBin();}
1330         if (cluster0->GetX()>rmax0) { rmax0=cluster0->GetX(); tmax0=cluster0->GetTimeBin();}
1331         ncl0++;
1332         if (cluster0->GetDetector()%36<18) nclA0++;
1333         if (cluster0->GetDetector()%36>=18) nclC0++;
1334       }
1335       AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1336       if (cluster1&&cluster1->GetX()>10){
1337         if (cluster1->GetX()<rmin1) { rmin1=cluster1->GetX();  tmin1=cluster1->GetTimeBin();}
1338         if (cluster1->GetX()>rmax1) { rmax1=cluster1->GetX(); tmax1=cluster1->GetTimeBin();}
1339         ncl1++;
1340         if (cluster1->GetDetector()%36<18) nclA1++;
1341         if (cluster1->GetDetector()%36>=18) nclC1++;
1342       }
1343     }
1344     Int_t cosmicType=0;  // types of cosmic topology
1345     if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side
1346     if ((nclA0<nclC0) && (nclA1<nclC1)) cosmicType=1; // CC side
1347     if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=2; // AC side
1348     if ((nclA0<nclC0) && (nclA1>nclC1)) cosmicType=3; // CA side
1349     //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=6; // AC side out of time
1350     //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=7; // CA side out of time
1351     //
1352     Double_t deltaTime = 0;   // correction for trigger not in time - equalizing the time arival
1353     if ((cosmicType>1)&&TMath::Abs(track1->GetZ()-track0->GetZ())>4){
1354       cosmicType+=4;
1355       deltaTime=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth();
1356       if (nclA0>nclC0) deltaTime*=-1; // if A side track
1357     }
1358     //
1359     TVectorD vectorDT(8);
1360     Int_t crossCounter=0;
1361     Double_t deltaTimeCross = AliTPCcalibCosmic::GetDeltaTime(rmin0, rmax0, rmin1, rmax1, tmin0, tmax0, tmin1, tmax1, TMath::Abs(track0->GetY()),vectorDT);
1362     Bool_t isOKTrigger=kTRUE;
1363     for (Int_t ic=0; ic<6;ic++) {
1364       if (TMath::Abs(vectorDT[ic])>0) {
1365         if (vectorDT[ic]+vectorDT[6]<0) isOKTrigger=kFALSE;
1366         if (vectorDT[ic]+vectorDT[7]>kMaxTime) isOKTrigger=kFALSE;
1367         if (isOKTrigger){
1368           crossCounter++;
1369         }
1370       }
1371     }
1372     Double_t deltaTimeCluster=deltaTime;
1373     if ((cosmicType==0 || cosmicType==1) && crossCounter>0){
1374       deltaTimeCluster=deltaTimeCross;
1375       cosmicType+=8;
1376     }
1377     if (nclA0*nclC0>0 || nclA1*nclC1>0) cosmicType+=16;  // mixed A side C side - bad for visualization
1378     //
1379     // Apply current transformation
1380     //
1381     //
1382     for (Int_t irow=0; irow<159; irow++){
1383       AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1384       if (cluster0 &&cluster0->GetX()>10){
1385         Double_t x0[3]={cluster0->GetRow(),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster};
1386         Int_t index0[1]={cluster0->GetDetector()};
1387         transform->Transform(x0,index0,0,1);
1388         cluster0->SetX(x0[0]);
1389         cluster0->SetY(x0[1]);
1390         cluster0->SetZ(x0[2]);
1391         //
1392       }
1393       AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1394       if (cluster1&&cluster1->GetX()>10){
1395         Double_t x1[3]={cluster1->GetRow(),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster};
1396         Int_t index1[1]={cluster1->GetDetector()};
1397         transform->Transform(x1,index1,0,1);
1398         cluster1->SetX(x1[0]);
1399         cluster1->SetY(x1[1]);
1400         cluster1->SetZ(x1[2]);
1401       }
1402     }
1403     //
1404     //
1405     Double_t alpha=track0->GetAlpha();   // rotation frame
1406     Double_t cos = TMath::Cos(alpha);
1407     Double_t sin = TMath::Sin(alpha);
1408     Double_t mass =  TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
1409     AliExternalTrackParam  btrack0=*(ftrack0->GetTPCOut());
1410     AliExternalTrackParam  btrack1=*(ftrack1->GetTPCOut());
1411     btrack0.Rotate(alpha);
1412     btrack1.Rotate(alpha);
1413     // change the sign for track 1
1414     Double_t* par1=(Double_t*)btrack0.GetParameter();
1415     par1[3]*=-1;
1416     par1[4]*=-1;
1417     btrack0.AddCovariance(covar);
1418     btrack1.AddCovariance(covar);
1419     btrack0.ResetCovariance(kResetCov);
1420     btrack1.ResetCovariance(kResetCov);
1421     Bool_t isOK=kTRUE;
1422     Bool_t isOKT=kTRUE;
1423     TObjArray tracks0(ncorr+1);
1424     TObjArray tracks1(ncorr+1);
1425     //
1426     Double_t dEdx0Tot=seed0->CookdEdxAnalytical(0.02,0.6,kTRUE);
1427     Double_t dEdx1Tot=seed1->CookdEdxAnalytical(0.02,0.6,kTRUE);
1428     Double_t dEdx0Max=seed0->CookdEdxAnalytical(0.02,0.6,kFALSE);
1429     Double_t dEdx1Max=seed1->CookdEdxAnalytical(0.02,0.6,kFALSE);
1430     //if (TMath::Abs((dEdx0Max+1)/(dEdx0Tot+1)-1.)>0.1) isOK=kFALSE;
1431     //if (TMath::Abs((dEdx1Max+1)/(dEdx1Tot+1)-1.)>0.1) isOK=kFALSE;
1432     ncl0=0; ncl1=0;
1433     for (Int_t icorr=-1; icorr<ncorr; icorr++){
1434       AliExternalTrackParam  rtrack0=btrack0;
1435       AliExternalTrackParam  rtrack1=btrack1;
1436       AliTPCCorrection *corr = 0;
1437       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
1438       //
1439       for (Int_t irow=159; irow>0; irow--){
1440         AliTPCclusterMI *cluster=seed0->GetClusterPointer(irow);
1441         if (!cluster) continue;
1442         if (!isOKT) break;
1443         Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1444         transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD);  // transform to global
1445         Float_t  r[3]={rD[0],rD[1],rD[2]};
1446         if (corr){
1447           corr->DistortPoint(r, cluster->GetDetector());
1448         }
1449         //
1450         Double_t cov[3]={0.01,0.,0.01};
1451         Double_t lx =cos*r[0]+sin*r[1];
1452         Double_t ly =-sin*r[0]+cos*r[1];
1453         rD[1]=ly; rD[0]=lx; rD[2]=r[2];  //transform to track local
1454         if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, lx,mass,1.,kFALSE)) isOKT=kFALSE;;
1455         if (!rtrack0.Update(&rD[1],cov)) isOKT =kFALSE;
1456         if (icorr<0) ncl0++;
1457       }
1458       //
1459       for (Int_t irow=159; irow>0; irow--){
1460         AliTPCclusterMI *cluster=seed1->GetClusterPointer(irow);
1461         if (!cluster) continue;
1462         if (!isOKT) break;
1463         Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1464         transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD);
1465         Float_t  r[3]={rD[0],rD[1],rD[2]};
1466         if (corr){
1467           corr->DistortPoint(r, cluster->GetDetector());
1468         }
1469         //
1470         Double_t cov[3]={0.01,0.,0.01};
1471         Double_t lx =cos*r[0]+sin*r[1];
1472         Double_t ly =-sin*r[0]+cos*r[1];
1473         rD[1]=ly; rD[0]=lx; rD[2]=r[2];
1474         if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, lx,mass,1.,kFALSE)) isOKT=kFALSE;
1475         if (!rtrack1.Update(&rD[1],cov)) isOKT=kFALSE;
1476         if (icorr<0) ncl1++;
1477       }
1478       if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,10.,kFALSE)) isOKT=kFALSE;
1479       if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,10.,kFALSE)) isOKT=kFALSE;
1480       if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,1.,kFALSE))  isOKT=kFALSE;
1481       if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,1.,kFALSE))  isOKT=kFALSE;
1482       const Double_t *param0=rtrack0.GetParameter();
1483       const Double_t *param1=rtrack1.GetParameter();
1484       for (Int_t ipar=0; ipar<4; ipar++){
1485         if (TMath::Abs(param1[ipar]-param0[ipar])>kMaxDelta[ipar]) isOK=kFALSE;
1486       }
1487       tracks0.AddAt(rtrack0.Clone(), icorr+1);
1488       tracks1.AddAt(rtrack1.Clone(), icorr+1);
1489       AliExternalTrackParam out0=*(ftrack0->GetTPCOut());
1490       AliExternalTrackParam out1=*(ftrack1->GetTPCOut());
1491       Int_t nentries=TMath::Min(ncl0,ncl1);
1492
1493       if (icorr<0) {
1494         (*pcstream)<<"cosmic"<<
1495           "isOK="<<isOK<<              // correct all propagation update and also residuals
1496           "isOKT="<<isOKT<<            // correct all propagation update
1497           "isOKTrigger="<<isOKTrigger<< // correct? estimate of trigger offset
1498           "id="<<cosmicType<<
1499           //
1500           //
1501           "cross="<<crossCounter<<
1502           "vDT.="<<&vectorDT<<
1503           //
1504           "dTime="<<deltaTime<<        // delta time using the A-c side cross
1505           "dTimeCross="<<deltaTimeCross<< // delta time using missing clusters
1506           //
1507           "dEdx0Max="<<dEdx0Max<<
1508           "dEdx0Tot="<<dEdx0Tot<<
1509           "dEdx1Max="<<dEdx1Max<<
1510           "dEdx1Tot="<<dEdx1Tot<<
1511           //
1512           "track0.="<<track0<<         // original track 0
1513           "track1.="<<track1<<         // original track 1
1514           "out0.="<<&out0<<             // outer track  0
1515           "out1.="<<&out1<<             // outer track  1
1516           "rtrack0.="<<&rtrack0<<      // refitted track with current transform
1517           "rtrack1.="<<&rtrack1<<     //
1518           "ncl0="<<ncl0<<
1519           "ncl1="<<ncl1<<
1520           "entries="<<nentries<<       // number of clusters
1521           "\n";
1522       }
1523     }
1524     //
1525
1526     if (isOK){
1527       Int_t nentries=TMath::Min(ncl0,ncl1);
1528       for (Int_t ipar=0; ipar<5; ipar++){
1529         for (Int_t icorr=-1; icorr<ncorr; icorr++){
1530           AliTPCCorrection *corr = 0;
1531           if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
1532           //
1533           AliExternalTrackParam *param0=(AliExternalTrackParam *) tracks0.At(icorr+1);
1534           AliExternalTrackParam *param1=(AliExternalTrackParam *) tracks1.At(icorr+1);
1535           distortions[icorr+1]=param1->GetParameter()[ipar]-param0->GetParameter()[ipar];
1536           if (icorr>=0){
1537             distortions[icorr+1]-=distortions[0];
1538           }
1539           //
1540           if (icorr<0){
1541             Double_t bz=AliTrackerBase::GetBz();
1542             Double_t gxyz[3];
1543             param0->GetXYZ(gxyz);
1544             Int_t dtype=20;
1545             Double_t theta=param0->GetParameter()[3];
1546             Double_t phi = alpha;
1547             Double_t snp = track0->GetInnerParam()->GetSnp();
1548             Double_t mean= distortions[0];
1549             Int_t index = param0->GetIndex(ipar,ipar);
1550             Double_t rms=TMath::Sqrt(param1->GetCovariance()[index]+param1->GetCovariance()[index]);
1551             if (crossCounter<1) rms*=1;
1552             Double_t sector=9*phi/TMath::Pi();
1553             Double_t dsec   = sector-TMath::Nint(sector+0.5);
1554             Double_t gx=gxyz[0],gy=gxyz[1],gz=gxyz[2];
1555             Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
1556             Double_t dRrec=0;
1557             //      Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
1558             Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
1559
1560             (*pcstream)<<"fit"<<  // dump valus for fit
1561               "run="<<run<<       //run number
1562               "bz="<<bz<<         // magnetic filed used
1563               "dtype="<<dtype<<   // detector match type 20
1564               "ptype="<<ipar<<   // parameter type
1565               "theta="<<theta<<   // theta
1566               "phi="<<phi<<       // phi
1567               "snp="<<snp<<       // snp
1568               "mean="<<mean<<     // mean dist value
1569               "rms="<<rms<<       // rms
1570               "sector="<<sector<<
1571               "dsec="<<dsec<<
1572               //
1573               "refX="<<refX<<      // reference radius
1574               "gx="<<gx<<         // global position
1575               "gy="<<gy<<         // global position
1576               "gz="<<gz<<         // global position
1577               "dRrec="<<dRrec<<      // delta Radius in reconstruction
1578               "pt="<<pt<<            //1/pt
1579               "id="<<cosmicType<<     //type of the cosmic used
1580               "entries="<<nentries;// number of entries in bin
1581             (*pcstream)<<"cosmicDebug"<<
1582               "p0.="<<param0<<    // dump distorted track 0
1583               "p1.="<<param1;    // dump distorted track 1
1584           }
1585           if (icorr>=0){
1586             (*pcstream)<<"fit"<<
1587               Form("%s=",corr->GetName())<<distortions[icorr+1];    // dump correction value
1588             (*pcstream)<<"cosmicDebug"<<
1589               Form("%s=",corr->GetName())<<distortions[icorr+1]<<    // dump correction value
1590               Form("%sp0.=",corr->GetName())<<param0<<    // dump distorted track 0
1591               Form("%sp1.=",corr->GetName())<<param1;    // dump distorted track 1
1592           }
1593         } //loop corrections
1594         (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1595         (*pcstream)<<"cosmicDebug"<<"isOK="<<isOK<<"\n";
1596       } //loop over parameters
1597     } // dump results
1598   }//loop tracks
1599 }
1600
1601
1602
1603 Double_t AliTPCcalibCosmic::GetDeltaTime(Double_t rmin0, Double_t rmax0, Double_t rmin1, Double_t rmax1, Double_t tmin0, Double_t tmax0, Double_t tmin1, Double_t tmax1, Double_t dcaR, TVectorD &vectorDT)
1604 {
1605   //
1606   // Estimate trigger offset between random cosmic event and "physics" trigger
1607   // Efficiency about 50 % of cases:
1608   // Cases:
1609   // 0. Tracks crossing A side and C side - no match in z - 30 % of cases
1610   // 1. Track only on one side and  dissapear at small or at high radiuses - 50 % of cases
1611   //    1.a) rmax<Rc    && tmax<Tcmax && tmax>tmin    => deltaT=Tcmax-tmax
1612   //    1.b) rmin>Rcmin && tmin<Tcmax && tmin>tmax    => deltaT=Tcmax-tmin
1613   //                      // case the z matching gives proper time
1614   //    1.c) rmax<Rc    && tmax>Tcmin && tmax<tmin    => deltaT=-tmax
1615   //
1616   // check algorithm:
1617   //    TCut cutStop = "(min(rmax0,rmax1)<235||abs(rmin0-rmin1)>10)"; // tracks not registered for full time
1618
1619   // Combinations:
1620   // 0-1 - forbidden
1621   // 0-2 - forbidden
1622   // 0-3 - occur - wrong correlation
1623   // 1-2 - occur - wrong correlation
1624   // 1-3 - forbidden
1625   // 2-3 - occur - small number of outlyers -20%
1626   // Frequency:
1627   //   0 - 106
1628   //   1 - 265
1629   //   2 - 206
1630   //   3 - 367
1631   //
1632   const Double_t kMaxRCut=235;  // max radius
1633   const Double_t kMinRCut=TMath::Max(dcaR,90.);  // min radius
1634   const Double_t kMaxDCut=30;   // max distance for minimal radius
1635   const Double_t kMinTime=110;
1636   const Double_t kMaxTime=950;
1637   Double_t deltaT=0;
1638   Int_t counter=0;
1639   vectorDT[6]=TMath::Min(TMath::Min(tmin0,tmin1),TMath::Min(tmax0,tmax1));
1640   vectorDT[7]=TMath::Max(TMath::Max(tmin0,tmin1),TMath::Max(tmax0,tmax1));
1641   if (TMath::Min(rmax0,rmax1)<kMaxRCut){
1642     // max cross - deltaT>0
1643     if (rmax0<kMaxRCut && tmax0 <kMaxTime && tmax0>tmin0) vectorDT[0]=kMaxTime-tmax0; // disapear at CE
1644     if (rmax1<kMaxRCut && tmax1 <kMaxTime && tmax1>tmin1) vectorDT[1]=kMaxTime-tmax1; // disapear at CE
1645     // min cross - deltaT<0 - OK they are correlated
1646     if (rmax0<kMaxRCut && tmax0 >kMinTime && tmax0<tmin0) vectorDT[2]=-tmax0;         // disapear at ROC
1647     if (rmax1<kMaxRCut && tmax1 >kMinTime && tmax1<tmin1) vectorDT[3]=-tmax1;         // disapear at ROC
1648   }
1649   if (rmin0> kMinRCut+kMaxDCut && tmin0 <kMaxTime && tmin0>tmax0) vectorDT[4]=kMaxTime-tmin0;
1650   if (rmin1> kMinRCut+kMaxDCut && tmin1 <kMaxTime && tmin1>tmax1) vectorDT[5]=kMaxTime-tmin1;
1651   Bool_t isOK=kTRUE;
1652   for (Int_t i=0; i<6;i++) {
1653     if (TMath::Abs(vectorDT[i])>0) {
1654       counter++;
1655       if (vectorDT[i]+vectorDT[6]<0) isOK=kFALSE;
1656       if (vectorDT[i]+vectorDT[7]>kMaxTime) isOK=kFALSE;
1657       if (isOK) deltaT=vectorDT[i];
1658     }
1659   }
1660   return deltaT;
1661 }