TPC/AliTPCcalibTime.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /*
  18 Comments to be written here:
  19
  20 1. What do we calibrate.
  21
  22   Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
  23
  24   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
  25
  26 2. How to interpret results
  27
  28 3. Simple example
  29
  30   a) determine the required time range:
  31
  32   AliXRDPROOFtoolkit tool;
  33   TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
  34   chain->Draw("GetTimeStamp()")
  35
  36   b) analyse calibration object on Proof in calibration train
  37
  38   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift);
  39
  40   c) plot results
  41   .x ~/NimStyle.C
  42   gSystem->Load("libANALYSIS");
  43   gSystem->Load("libTPCcalib");
  44
  45   TFile f("CalibObjectsTrain1.root");
  46   AliTPCcalibTime *calib = (AliTPCcalibTime *)f->Get("calibTime");
  47   calib->GetHistoDrift("all")->Projection(2,0)->Draw()
  48   calib->GetFitDrift("all")->Draw("lp")
  49
  50 4. Analysis using debug streamers.
  51
  52   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  53   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  54   AliXRDPROOFtoolkit tool;
  55   TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
  56   chainTime->Lookup();
  57   TChain * chainLaser = tool.MakeChain("time.txt","timeLaser",0,10200);
  58   chainLaser->Lookup();
  59 */
  60
  61 #include "Riostream.h"
  62 #include "TChain.h"
  63 #include "TTree.h"
  64 #include "TH1F.h"
  65 #include "TH2F.h"
  66 #include "TH3F.h"
  67 #include "THnSparse.h"
  68 #include "TList.h"
  69 #include "TMath.h"
  70 #include "TCanvas.h"
  71 #include "TFile.h"
  72 #include "TF1.h"
  73 #include "TVectorD.h"
  74 #include "TProfile.h"
  75 #include "TGraphErrors.h"
  76 #include "TCanvas.h"
  77
  78 #include "AliTPCclusterMI.h"
  79 #include "AliTPCseed.h"
  80 #include "AliESDVertex.h"
  81 #include "AliESDEvent.h"
  82 #include "AliESDfriend.h"
  83 #include "AliESDInputHandler.h"
  84 #include "AliAnalysisManager.h"
  85
  86 #include "AliTracker.h"
  87 #include "AliMagF.h"
  88 #include "AliTPCCalROC.h"
  89
  90 #include "AliLog.h"
  91
  92 #include "AliTPCcalibTime.h"
  93
  94 #include "TTreeStream.h"
  95 #include "AliTPCTracklet.h"
  96 #include "TTimeStamp.h"
  97 #include "AliTPCcalibDB.h"
  98 #include "AliTPCcalibLaser.h"
  99 #include "AliDCSSensorArray.h"
 100 #include "AliDCSSensor.h"
 101
 102 ClassImp(AliTPCcalibTime)
 103
 104
 105 AliTPCcalibTime::AliTPCcalibTime()
 106   :AliTPCcalibBase(),
 107    fLaser(0),       // pointer to laser calibration
 108    fDz(0),          // current delta z
 109    fCutMaxD(3),        // maximal distance in rfi ditection
 110    fCutMaxDz(25),      // maximal distance in rfi ditection
 111    fCutTheta(0.03),    // maximal distan theta
 112    fCutMinDir(-0.99),  // direction vector products
 113    fCutTracks(10),
 114    fArrayDz(0),          //NEW! Tmap of V drifts for different triggers
 115    fTimeBins(0),
 116    fTimeStart(0),
 117    fTimeEnd(0),
 118    fPtBins(0),
 119    fPtStart(0),
 120    fPtEnd(0),
 121    fVdriftBins(0),
 122    fVdriftStart(0),
 123    fVdriftEnd(0),
 124    fRunBins(0),
 125    fRunStart(0),
 126    fRunEnd(0)
 127 //   fBinsVdrift(fTimeBins,fPtBins,fVdriftBins),
 128 //   fXminVdrift(fTimeStart,fPtStart,fVdriftStart),
 129 //   fXmaxVdrift(fTimeEnd,fPtEnd,fVdriftEnd)
 130 {
 131   AliInfo("Default Constructor");
 132   for (Int_t i=0;i<3;i++) {
 133     fHistVdriftLaserA[i]=0;
 134     fHistVdriftLaserC[i]=0;
 135   }
 136   for (Int_t i=0;i<10;i++) {
 137     fCosmiMatchingHisto[i]=0;
 138   }
 139 }
 140
 141 AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeVdrift)
 142   :AliTPCcalibBase(),
 143    fLaser(0),            // pointer to laser calibration
 144    fDz(0),               // current delta z
 145    fCutMaxD(5*0.5356),   // maximal distance in rfi ditection
 146    fCutMaxDz(40),   // maximal distance in rfi ditection
 147    fCutTheta(5*0.004644),// maximal distan theta
 148    fCutMinDir(-0.99),    // direction vector products
 149    fCutTracks(10),
 150    fArrayDz(0),            //Tmap of V drifts for different triggers
 151    fTimeBins(0),
 152    fTimeStart(0),
 153    fTimeEnd(0),
 154    fPtBins(0),
 155    fPtStart(0),
 156    fPtEnd(0),
 157    fVdriftBins(0),
 158    fVdriftStart(0),
 159    fVdriftEnd(0),
 160    fRunBins(0),
 161    fRunStart(0),
 162    fRunEnd(0)
 163 {
 164   SetName(name);
 165   SetTitle(title);
 166   for (Int_t i=0;i<3;i++) {
 167     fHistVdriftLaserA[i]=0;
 168     fHistVdriftLaserC[i]=0;
 169   }
 170
 171   AliInfo("Non Default Constructor");
 172   fTimeBins   =(EndTime-StartTime)/deltaIntegrationTimeVdrift;
 173   fTimeStart  =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi();
 174   fTimeEnd    =EndTime;   //(((TObjString*)(mapGRP->GetValue("fAliceStopTime")))->GetString()).Atoi();
 175   fPtBins     = 400;
 176   fPtStart    = -0.04;
 177   fPtEnd      =  0.04;
 178   fVdriftBins = 500;
 179   fVdriftStart= -0.1;
 180   fVdriftEnd  =  0.1;
 181   fRunBins    = 100001;
 182   fRunStart   = -1.5;
 183   fRunEnd     = 99999.5;
 184
 185   Int_t    binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins };
 186   Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart  , fRunStart};
 187   Double_t xmaxVdriftLaser[4] = {fTimeEnd  , fPtEnd  , fVdriftEnd    , fRunEnd  };
 188   TString axisTitle[4]={
 189     "T",
 190     "#delta_{P/T}",
 191     "value",
 192     "run"
 193   };
 194   TString histoName[3]={
 195     "Loffset",
 196     "Lcorr",
 197     "Lgy"
 198   };
 199
 200
 201   for (Int_t i=0;i<3;i++) {
 202     fHistVdriftLaserA[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
 203     fHistVdriftLaserC[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
 204     fHistVdriftLaserA[i]->SetName(histoName[i]);
 205     fHistVdriftLaserC[i]->SetName(histoName[i]);
 206     for (Int_t iaxis=0; iaxis<4;iaxis++){
 207       fHistVdriftLaserA[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
 208       fHistVdriftLaserC[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
 209     }
 210   }
 211   fBinsVdrift[0] = fTimeBins;
 212   fBinsVdrift[1] = fPtBins;
 213   fBinsVdrift[2] = fVdriftBins;
 214   fBinsVdrift[3] = fRunBins;
 215   fXminVdrift[0] = fTimeStart;
 216   fXminVdrift[1] = fPtStart;
 217   fXminVdrift[2] = fVdriftStart;
 218   fXminVdrift[3] = fRunStart;
 219   fXmaxVdrift[0] = fTimeEnd;
 220   fXmaxVdrift[1] = fPtEnd;
 221   fXmaxVdrift[2] = fVdriftEnd;
 222   fXmaxVdrift[3] = fRunEnd;
 223
 224   fArrayDz=new TObjArray();
 225   fArrayDz->AddLast(fHistVdriftLaserA[0]);
 226   fArrayDz->AddLast(fHistVdriftLaserA[1]);
 227   fArrayDz->AddLast(fHistVdriftLaserA[2]);
 228   fArrayDz->AddLast(fHistVdriftLaserC[0]);
 229   fArrayDz->AddLast(fHistVdriftLaserC[1]);
 230   fArrayDz->AddLast(fHistVdriftLaserC[2]);
 231
 232   fCosmiMatchingHisto[0]=new TH1F("Cosmics matching","p0-all"   ,100,-10*0.5356  ,10*0.5356  );
 233   fCosmiMatchingHisto[1]=new TH1F("Cosmics matching","p1-all"   ,100,-10*4.541   ,10*4.541   );
 234   fCosmiMatchingHisto[2]=new TH1F("Cosmics matching","p2-all"   ,100,-10*0.01134 ,10*0.01134 );
 235   fCosmiMatchingHisto[3]=new TH1F("Cosmics matching","p3-all"   ,100,-10*0.004644,10*0.004644);
 236   fCosmiMatchingHisto[4]=new TH1F("Cosmics matching","p4-all"   ,100,-10*0.03773 ,10*0.03773 );
 237   fCosmiMatchingHisto[5]=new TH1F("Cosmics matching","p0-isPair",100,-10*0.5356  ,10*0.5356  );
 238   fCosmiMatchingHisto[6]=new TH1F("Cosmics matching","p1-isPair",100,-10*4.541   ,10*4.541   );
 239   fCosmiMatchingHisto[7]=new TH1F("Cosmics matching","p2-isPair",100,-10*0.01134 ,10*0.01134 );
 240   fCosmiMatchingHisto[8]=new TH1F("Cosmics matching","p3-isPair",100,-10*0.004644,10*0.004644);
 241   fCosmiMatchingHisto[9]=new TH1F("Cosmics matching","p4-isPair",100,-10*0.03773 ,10*0.03773 );
 242 //  Char_t nameHisto[3]={'p','0','\n'};
 243 //  for (Int_t i=0;i<10;i++){
 244 //    fCosmiMatchingHisto[i]=new TH1F("Cosmics matching",nameHisto,8192,0,0);
 245 //    nameHisto[1]++;
 246 //    if(i==4) nameHisto[1]='0';
 247 //  }
 248 }
 249
 250 AliTPCcalibTime::~AliTPCcalibTime(){
 251   //
 252   // Destructor
 253   //
 254   for(Int_t i=0;i<3;i++){
 255     if(fHistVdriftLaserA[i]){
 256       delete fHistVdriftLaserA[i];
 257       fHistVdriftLaserA[i]=NULL;
 258     }
 259     if(fHistVdriftLaserC[i]){
 260       delete fHistVdriftLaserC[i];
 261       fHistVdriftLaserC[i]=NULL;
 262     }
 263   }
 264   if(fArrayDz){
 265     fArrayDz->SetOwner();
 266     fArrayDz->Delete();
 267     delete fArrayDz;
 268     fArrayDz=NULL;
 269   }
 270   for(Int_t i=0;i<5;i++){
 271     if(fCosmiMatchingHisto[i]){
 272       delete fCosmiMatchingHisto[i];
 273       fCosmiMatchingHisto[i]=NULL;
 274     }
 275   }
 276 }
 277
 278 Bool_t AliTPCcalibTime::IsLaser(AliESDEvent */*event*/){
 279   return kTRUE; //More accurate creteria to be added
 280 }
 281 Bool_t AliTPCcalibTime::IsCosmics(AliESDEvent */*event*/){
 282   return kTRUE; //More accurate creteria to be added
 283 }
 284 Bool_t AliTPCcalibTime::IsBeam(AliESDEvent */*event*/){
 285   return kTRUE; //More accurate creteria to be added
 286 }
 287 void AliTPCcalibTime::ResetCurrent(){
 288   fDz=0; //Reset current dz
 289 }
 290 void AliTPCcalibTime::Process(AliESDEvent *event){
 291   if(!event) return;
 292   if (event->GetNumberOfTracks()<2) return;
 293   ResetCurrent();
 294   if(IsLaser  (event)) ProcessLaser (event);
 295   if(IsCosmics(event)) ProcessCosmic(event);
 296   if(IsBeam   (event)) ProcessBeam  (event);
 297 }
 298
 299 void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){
 300   //
 301   // Fit drift velocity using laser
 302   //
 303   // 0. cuts
 304   const Int_t    kMinTracks     = 40;    // minimal number of laser tracks
 305   const Int_t    kMinTracksSide = 20;    // minimal number of tracks per side
 306   const Float_t  kMaxDeltaZ     = 30.;   // maximal trigger delay
 307   const Float_t  kMaxDeltaV     = 0.05;  // maximal deltaV
 308   const Float_t  kMaxRMS        = 0.1;   // maximal RMS of tracks
 309   //
 310   /*
 311     TCut cutRMS("sqrt(laserA.fElements[4])<0.1&&sqrt(laserC.fElements[4])<0.1");
 312     TCut cutZ("abs(laserA.fElements[0]-laserC.fElements[0])<3");
 313     TCut cutV("abs(laserA.fElements[1]-laserC.fElements[1])<0.01");
 314     TCut cutY("abs(laserA.fElements[2]-laserC.fElements[2])<2");
 315     TCut cutAll = cutRMS+cutZ+cutV+cutY;
 316   */
 317   if (event->GetNumberOfTracks()<kMinTracks) return;
 318   //
 319   if(!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
 320   fLaser->Process(event);
 321   if (fLaser->GetNtracks()<kMinTracks) return;   // small amount of tracks cut
 322   if (fLaser->fFitAside->GetNrows()==0  && fLaser->fFitCside->GetNrows()==0) return;  // no fit neither a or C side
 323   //
 324   // debug streamer  - activate stream level
 325   // Use it for tuning of the cuts
 326   //
 327   // cuts to be applied
 328   //
 329   Int_t isReject[2]={0,0};
 330   //
 331   // not enough tracks
 332   if (TMath::Abs((*fLaser->fFitAside)[3]) < kMinTracksSide) isReject[0]|=1;
 333   if (TMath::Abs((*fLaser->fFitCside)[3]) < kMinTracksSide) isReject[1]|=1;
 334   // unreasonable z offset
 335   if (TMath::Abs((*fLaser->fFitAside)[0])>kMaxDeltaZ)  isReject[0]|=2;
 336   if (TMath::Abs((*fLaser->fFitCside)[0])>kMaxDeltaZ)  isReject[1]|=2;
 337   // unreasonable drift velocity
 338   if (TMath::Abs((*fLaser->fFitAside)[1]-1)>kMaxDeltaV)  isReject[0]|=4;
 339   if (TMath::Abs((*fLaser->fFitCside)[1]-1)>kMaxDeltaV)  isReject[1]|=4;
 340   // big chi2
 341   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitAside)[4]))>kMaxRMS ) isReject[0]|=8;
 342   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitCside)[4]))>kMaxRMS ) isReject[1]|=8;
 343
 344
 345
 346
 347   if (fStreamLevel>0){
 348     printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
 349
 350     TTreeSRedirector *cstream = GetDebugStreamer();
 351     if (cstream){
 352       TTimeStamp tstamp(fTime);
 353       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 354       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 355       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 356       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 357       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 358       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 359       TVectorD vecGoofie(20);
 360       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 361       if (goofieArray){
 362         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 363           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 364           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 365         }
 366       }
 367       (*cstream)<<"laserInfo"<<
 368         "run="<<fRun<<              //  run number
 369         "event="<<fEvent<<          //  event number
 370         "time="<<fTime<<            //  time stamp of event
 371         "trigger="<<fTrigger<<      //  trigger
 372         "mag="<<fMagF<<             //  magnetic field
 373         // Environment values
 374         "press0="<<valuePressure0<<
 375         "press1="<<valuePressure1<<
 376         "pt0="<<ptrelative0<<
 377         "pt1="<<ptrelative1<<
 378         "temp0="<<temp0<<
 379         "temp1="<<temp1<<
 380         "vecGoofie.=<<"<<&vecGoofie<<
 381         //laser
 382         "rejectA="<<isReject[0]<<
 383         "rejectC="<<isReject[1]<<
 384         "laserA.="<<fLaser->fFitAside<<
 385         "laserC.="<<fLaser->fFitCside<<
 386         "laserAC.="<<fLaser->fFitACside<<
 387         "trigger="<<event->GetFiredTriggerClasses()<<
 388         "\n";
 389     }
 390   }
 391   //
 392   // fill histos
 393   //
 394   TVectorD vdriftA(5), vdriftC(5),vdriftAC(5);
 395   vdriftA=*(fLaser->fFitAside);
 396   vdriftC=*(fLaser->fFitCside);
 397   vdriftAC=*(fLaser->fFitACside);
 398   Int_t npointsA=0, npointsC=0;
 399   Float_t chi2A=0, chi2C=0;
 400   npointsA= TMath::Nint(vdriftA[3]);
 401   chi2A= vdriftA[4];
 402   npointsC= TMath::Nint(vdriftC[3]);
 403   chi2C= vdriftC[4];
 404
 405   TTimeStamp tstamp(fTime);
 406   Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 407   Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 408   Double_t driftA=0, driftC=0;
 409   if (vdriftA[1]>1.-kMaxDeltaV) driftA = 1./vdriftA[1]-1.;
 410   if (vdriftC[1]>1.-kMaxDeltaV) driftC = 1./vdriftC[1]-1.;
 411   //
 412   Double_t vecDriftLaserA[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftA,event->GetRunNumber()};
 413   Double_t vecDriftLaserC[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftC,event->GetRunNumber()};
 414   //  Double_t vecDrift[4]      ={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitACside))[1])-1,event->GetRunNumber()};
 415
 416   for (Int_t icalib=0;icalib<3;icalib++){
 417     if (icalib==0){ //z0 shift
 418       vecDriftLaserA[2]=vdriftA[0]/250.;
 419       vecDriftLaserC[2]=vdriftC[0]/250.;
 420     }
 421     if (icalib==1){ //vdrel shift
 422       vecDriftLaserA[2]=driftA;
 423       vecDriftLaserC[2]=driftC;
 424     }
 425     if (icalib==2){ //gy shift - full gy - full drift
 426       vecDriftLaserA[2]=vdriftA[2]/250.;
 427       vecDriftLaserC[2]=vdriftC[2]/250.;
 428     }
 429     if (npointsA>kMinTracks) fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
 430     if (npointsC>kMinTracks) fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
 431   }
 432
 433 //   THnSparse* curHist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 434 //   TString shortName=curHist->ClassName();
 435 //   shortName+="_MEAN_DRIFT_LASER_";
 436 //   delete curHist;
 437 //   curHist=NULL;
 438 //   TString name="";
 439
 440 //   name=shortName;
 441 //   name+=event->GetFiredTriggerClasses();
 442 //   name.ToUpper();
 443 //   curHist=(THnSparseF*)fArrayDz->FindObject(name);
 444 //   if(!curHist){
 445 //     curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 446 //     fArrayDz->AddLast(curHist);
 447 //   }
 448 //   curHist->Fill(vecDrift);
 449
 450 //   name=shortName;
 451 //   name+="ALL";
 452 //   name.ToUpper();
 453 //   curHist=(THnSparseF*)fArrayDz->FindObject(name);
 454 //   if(!curHist){
 455 //     curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 456 //     fArrayDz->AddLast(curHist);
 457 //   }
 458 //   curHist->Fill(vecDrift);
 459 }
 460
 461 void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){
 462   if (!event) {
 463     Printf("ERROR: ESD not available");
 464     return;
 465   }
 466   if (event->GetTimeStamp() == 0 ) {
 467     Printf("no time stamp!");
 468     return;
 469   }
 470
 471   //fd
 472   // Find cosmic pairs
 473   //
 474   // Track0 is choosen in upper TPC part
 475   // Track1 is choosen in lower TPC part
 476   //
 477   Int_t ntracks=event->GetNumberOfTracks();
 478   if (ntracks==0) return;
 479   if (ntracks > fCutTracks) return;
 480
 481   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 482   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 483
 484   TObjArray  tpcSeeds(ntracks);
 485   Double_t vtxx[3]={0,0,0};
 486   Double_t svtxx[3]={0.000001,0.000001,100.};
 487   AliESDVertex vtx(vtxx,svtxx);
 488   //
 489   // track loop
 490   //
 491   for (Int_t i=0;i<ntracks;++i) {
 492     AliESDtrack *track = event->GetTrack(i);
 493
 494     const AliExternalTrackParam * trackIn = track->GetInnerParam();
 495     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 496     if (!trackIn) continue;
 497     if (!trackOut) continue;
 498
 499     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 500     TObject *calibObject;
 501     AliTPCseed *seed = 0;
 502     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 503     if (seed) tpcSeeds.AddAt(seed,i);
 504   }
 505   if (ntracks<2) return;
 506   //
 507   // Find pairs
 508   //
 509
 510   for (Int_t i=0;i<ntracks;++i) {
 511     AliESDtrack *track0 = event->GetTrack(i);
 512     // track0 - choosen upper part
 513     if (!track0) continue;
 514     if (!track0->GetOuterParam()) continue;
 515     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 516     Double_t d1[3];
 517     track0->GetDirection(d1);
 518     for (Int_t j=0;j<ntracks;++j) {
 519       if (i==j) continue;
 520       AliESDtrack *track1 = event->GetTrack(j);
 521       //track 1 lower part
 522       if (!track1) continue;
 523       if (!track1->GetOuterParam()) continue;
 524       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 525       //
 526       Double_t d2[3];
 527       track1->GetDirection(d2);
 528
 529       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 530       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 531       if (! seed0) continue;
 532       if (! seed1) continue;
 533       Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
 534       Float_t dist0  = track0->GetLinearD(0,0);
 535       Float_t dist1  = track1->GetLinearD(0,0);
 536       //
 537       // conservative cuts - convergence to be guarantied
 538       // applying before track propagation
 539       if (TMath::Abs(dist0+dist1)>fCutMaxD) continue;   // distance to the 0,0
 540       if (dir>fCutMinDir) continue;               // direction vector product
 541       Float_t bz = AliTracker::GetBz();
 542       Float_t dvertex0[2];   //distance to 0,0
 543       Float_t dvertex1[2];   //distance to 0,0
 544       track0->GetDZ(0,0,0,bz,dvertex0);
 545       track1->GetDZ(0,0,0,bz,dvertex1);
 546       if (TMath::Abs(dvertex0[1])>250) continue;
 547       if (TMath::Abs(dvertex1[1])>250) continue;
 548       //
 549       //
 550       //
 551       Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
 552       AliExternalTrackParam param0(*track0);
 553       AliExternalTrackParam param1(*track1);
 554       //
 555       // Propagate using Magnetic field and correct fo material budget
 556       //
 557       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 558       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 559       //
 560       // Propagate rest to the 0,0 DCA - z should be ignored
 561       //
 562       //Bool_t b0 = ;
 563       param0.PropagateToDCA(&vtx,bz,1000);
 564       //Bool_t b1 =
 565       param1.PropagateToDCA(&vtx,bz,1000);
 566       param0.GetDZ(0,0,0,bz,dvertex0);
 567       param1.GetDZ(0,0,0,bz,dvertex1);
 568       Double_t xyz0[3];
 569       Double_t xyz1[3];
 570       param0.GetXYZ(xyz0);
 571       param1.GetXYZ(xyz1);
 572       Bool_t isPair = IsPair(&param0,&param1);
 573       Bool_t isCross = IsCross(track0, track1);
 574       Bool_t isSame = IsSame(track0, track1);
 575
 576       THnSparse* hist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 577       TString shortName=hist->ClassName();
 578       shortName+="_MEAN_VDRIFT_COSMICS_";
 579       delete hist;
 580       hist=NULL;
 581
 582       if(isSame || (isCross && isPair)){
 583         if (track0->GetTPCNcls() > 80) {
 584           fDz = param0.GetZ() - param1.GetZ();
 585           if(track0->GetOuterParam()->GetZ()<0) fDz=-fDz;
 586           TTimeStamp tstamp(fTime);
 587           Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 588           Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 589           Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
 590           THnSparse* curHist=NULL;
 591           TString name="";
 592
 593           name=shortName;
 594           name+=event->GetFiredTriggerClasses();
 595           name.ToUpper();
 596           curHist=(THnSparseF*)fArrayDz->FindObject(name);
 597           if(!curHist){
 598             curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 599             fArrayDz->AddLast(curHist);
 600           }
 601 //        curHist=(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
 602 //        if(!curHist){
 603 //          curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 604 //          fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
 605 //        }
 606           curHist->Fill(vecDrift);
 607
 608           name=shortName;
 609           name+="ALL";
 610           name.ToUpper();
 611           curHist=(THnSparseF*)fArrayDz->FindObject(name);
 612           if(!curHist){
 613             curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 614             fArrayDz->AddLast(curHist);
 615           }
 616 //        curHist=(THnSparseF*)(fMapDz->GetValue("all"));
 617 //        if(!curHist){
 618 //          curHist=new THnSparseF("all","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 619 //          fMapDz->Add(new TObjString("all"),curHist);
 620 //        }
 621           curHist->Fill(vecDrift);
 622         }
 623       }
 624       TTreeSRedirector *cstream = GetDebugStreamer();
 625       if (fStreamLevel>0){
 626         if (cstream){
 627         (*cstream)<<"trackInfo"<<
 628           "tr0.="<<track0<<
 629           "tr1.="<<track1<<
 630           "p0.="<<&param0<<
 631           "p1.="<<&param1<<
 632           "isPair="<<isPair<<
 633           "isCross="<<isCross<<
 634           "isSame="<<isSame<<
 635           "fDz="<<fDz<<
 636           "fRun="<<fRun<<
 637           "fTime="<<fTime<<
 638           "\n";
 639         }
 640       }
 641     } // end 2nd order loop
 642   } // end 1st order loop
 643
 644   if (fStreamLevel>0){
 645     TTreeSRedirector *cstream = GetDebugStreamer();
 646     if (cstream){
 647       TTimeStamp tstamp(fTime);
 648       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 649       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 650       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 651       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 652       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 653       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 654       TVectorD vecGoofie(20);
 655       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 656       if (goofieArray){
 657         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 658           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 659           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 660         }
 661       }
 662       (*cstream)<<"timeInfo"<<
 663         "run="<<fRun<<              //  run number
 664         "event="<<fEvent<<          //  event number
 665         "time="<<fTime<<            //  time stamp of event
 666         "trigger="<<fTrigger<<      //  trigger
 667         "mag="<<fMagF<<             //  magnetic field
 668         // Environment values
 669         "press0="<<valuePressure0<<
 670         "press1="<<valuePressure1<<
 671         "pt0="<<ptrelative0<<
 672         "pt1="<<ptrelative1<<
 673         "temp0="<<temp0<<
 674         "temp1="<<temp1<<
 675         "vecGoofie.=<<"<<&vecGoofie<<
 676         //
 677         // accumulated values
 678         //
 679         "fDz="<<fDz<<          //! current delta z
 680         "trigger="<<event->GetFiredTriggerClasses()<<
 681         "\n";
 682     }
 683   }
 684   printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
 685 }
 686
 687 void AliTPCcalibTime::ProcessBeam(AliESDEvent */*event*/){
 688 }
 689
 690 void AliTPCcalibTime::Analyze(){}
 691
 692 THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name){
 693   TIterator* iterator = fArrayDz->MakeIterator();
 694   iterator->Reset();
 695   TString newName=name;
 696   newName.ToUpper();
 697   THnSparse* newHist=new THnSparseF(newName,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 698   THnSparse* addHist=NULL;
 699   while((addHist=(THnSparseF*)iterator->Next())){
 700   if(!addHist) continue;
 701     TString histName=addHist->GetName();
 702     if(!histName.Contains(newName)) continue;
 703     addHist->Print();
 704     newHist->Add(addHist);
 705   }
 706   return newHist;
 707 }
 708
 709 TObjArray* AliTPCcalibTime::GetHistoDrift(){
 710   return fArrayDz;
 711 }
 712
 713 TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){
 714   THnSparse* histoDrift=GetHistoDrift(name);
 715   TGraphErrors* graphDrift=NULL;
 716   if(histoDrift){
 717     graphDrift=FitSlices(histoDrift,2,0,400,100,0.05,0.95, kTRUE);
 718     TString end=histoDrift->GetName();
 719     Int_t pos=end.Index("_");
 720     end=end(pos,end.Capacity()-pos);
 721     TString graphName=graphDrift->ClassName();
 722     graphName+=end;
 723     graphName.ToUpper();
 724     graphDrift->SetName(graphName);
 725   }
 726   return graphDrift;
 727 }
 728
 729 TObjArray* AliTPCcalibTime::GetGraphDrift(){
 730   TObjArray* arrayGraphDrift=new TObjArray();
 731   TIterator* iterator=fArrayDz->MakeIterator();
 732   iterator->Reset();
 733   THnSparse* addHist=NULL;
 734   while((addHist=(THnSparseF*)iterator->Next())) arrayGraphDrift->AddLast(GetGraphDrift(addHist->GetName()));
 735   return arrayGraphDrift;
 736 }
 737
 738 AliSplineFit* AliTPCcalibTime::GetFitDrift(const char* name){
 739   TGraph* graphDrift=GetGraphDrift(name);
 740   AliSplineFit* fitDrift=NULL;
 741   if(graphDrift && graphDrift->GetN()){
 742     fitDrift=new AliSplineFit();
 743     fitDrift->SetGraph(graphDrift);
 744     fitDrift->SetMinPoints(graphDrift->GetN()+1);
 745     fitDrift->InitKnots(graphDrift,2,0,0.001);
 746     fitDrift->SplineFit(0);
 747     TString end=graphDrift->GetName();
 748     Int_t pos=end.Index("_");
 749     end=end(pos,end.Capacity()-pos);
 750     TString fitName=fitDrift->ClassName();
 751     fitName+=end;
 752     fitName.ToUpper();
 753     //fitDrift->SetName(fitName);
 754     delete graphDrift;
 755     graphDrift=NULL;
 756   }
 757   return fitDrift;
 758 }
 759
 760 //TObjArray* AliTPCcalibTime::GetFitDrift(){
 761 //  TObjArray* arrayFitDrift=new TObjArray();
 762 //  TIterator* iterator = fArrayDz->MakeIterator();
 763 //  iterator->Reset();
 764 //  THnSparse* addHist=NULL;
 765 //  while((addHist=(THnSparseF*)iterator->Next())) arrayFitDrift->AddLast(GetFitDrift(addHist->GetName()));
 766 //  return arrayFitDrift;
 767 //}
 768
 769 Long64_t AliTPCcalibTime::Merge(TCollection *li) {
 770   TIterator* iter = li->MakeIterator();
 771   AliTPCcalibTime* cal = 0;
 772
 773   while ((cal = (AliTPCcalibTime*)iter->Next())) {
 774     if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
 775       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 776       return -1;
 777     }
 778     for (Int_t imeas=0; imeas<3; imeas++){
 779       if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
 780         fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
 781         fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
 782       }
 783     }
 784     TObjArray* addArray=cal->GetHistoDrift();
 785     if(!addArray) return 0;
 786     TIterator* iterator = addArray->MakeIterator();
 787     iterator->Reset();
 788     THnSparse* addHist=NULL;
 789     while((addHist=(THnSparseF*)iterator->Next())){
 790       if(!addHist) continue;
 791       addHist->Print();
 792       THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName());
 793       if(!localHist){
 794         localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 795         fArrayDz->AddLast(localHist);
 796       }
 797       localHist->Add(addHist);
 798     }
 799 //    TMap * addMap=cal->GetHistoDrift();
 800 //    if(!addMap) return 0;
 801 //    TIterator* iterator = addMap->MakeIterator();
 802 //    iterator->Reset();
 803 //    TPair* addPair=0;
 804 //    while((addPair=(TPair *)(addMap->FindObject(iterator->Next())))){
 805 //      THnSparse* addHist=dynamic_cast<THnSparseF*>(addPair->Value());
 806 //      if (!addHist) continue;
 807 //      addHist->Print();
 808 //      THnSparse* localHist=dynamic_cast<THnSparseF*>(fMapDz->GetValue(addHist->GetName()));
 809 //      if(!localHist){
 810 //        localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 811 //        fMapDz->Add(new TObjString(addHist->GetName()),localHist);
 812 //      }
 813 //      localHist->Add(addHist);
 814 //    }
 815     for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
 816   }
 817   return 0;
 818 }
 819
 820 Bool_t  AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 821   /*
 822   // 0. Same direction - OPOSITE  - cutDir +cutT
 823   TCut cutDir("cutDir","dir<-0.99")
 824   // 1.
 825   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 826   //
 827   // 2. The same rphi
 828   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 829   //
 830   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 831   // 1/Pt diff cut
 832   */
 833   const Double_t *p0 = tr0->GetParameter();
 834   const Double_t *p1 = tr1->GetParameter();
 835   fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
 836   fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
 837   fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
 838   fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
 839   fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
 840
 841   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 842   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 843   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz)  return kFALSE;
 844   Double_t d0[3], d1[3];
 845   tr0->GetDirection(d0);
 846   tr1->GetDirection(d1);
 847   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 848
 849   fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
 850   fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
 851   fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
 852   fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
 853   fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
 854
 855   return kTRUE;
 856 }
 857 Bool_t AliTPCcalibTime::IsCross(AliESDtrack *tr0, AliESDtrack *tr1){
 858   return  tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0 && tr0->GetInnerParam()->GetZ()*tr1->GetInnerParam()->GetZ()<0 && tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0 && tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0;
 859 }
 860
 861 Bool_t AliTPCcalibTime::IsSame(AliESDtrack */*tr0*/, AliESDtrack */*tr1*/){
 862   // To be implemented
 863   return kFALSE;
 864 }
 865
 866 /*
 867 chainDrift->Draw("p0.fP[0]+p1.fP[0]","isPair");
 868   mean ~-0.02  ~-0.03913
 869   RMS  ~ 0.5   ~ 0.5356    --> 3    (fCutMaxD)
 870
 871 chainDrift->Draw("p0.fP[1]-p1.fP[1]","isPair");
 872   mean         ~ 0.1855
 873   RMS          ~ 4.541     -->25    (fCutMaxDz)
 874
 875 chainDrift->Draw("p1.fAlpha-p0.fAlpha+pi","isPair");
 876 //chainDrift->Draw("p1.fAlpha+p0.fAlpha","isPair");
 877 //chainDrift->Draw("p1.fP[2]-p0.fP[2]+pi","isPair");
 878 //chainDrift->Draw("p1.fP[2]+p0.fP[2]","isPair");
 879   mean ~ 0     ~ 0.001898
 880   RMS  ~ 0.009 ~ 0.01134   --> 0.06
 881
 882 chainDrift->Draw("p0.fP[3]+p1.fP[3]","isPair");
 883   mean ~ 0.0013 ~ 0.001539
 884   RMS  ~ 0.003  ~ 0.004644 --> 0.03  (fCutTheta)
 885
 886 chainDrift->Draw("p1.fP[4]+p0.fP[4]>>his(100,-0.2,0.2)","isPair")
 887   mean ~ 0.012  ~-0.0009729
 888   RMS  ~ 0.036  ~ 0.03773  --> 0.2
 889 */
 890