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