TPC/AliTPCcalibTime.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 Comments to be written here:
  18
  19 1. What do we calibrate.
  20
  21   Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
  22
  23   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
  24
  25 2. How to interpret results
  26
  27 3. Simple example
  28
  29   a) determine the required time range:
  30
  31   AliXRDPROOFtoolkit tool;
  32   TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
  33   chain->Draw("GetTimeStamp()")
  34
  35   b) analyse calibration object on Proof in calibration train
  36
  37   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift);
  38
  39   c) plot results
  40   .x ~/NimStyle.C
  41   gSystem->Load("libANALYSIS");
  42   gSystem->Load("libTPCcalib");
  43
  44   TFile f("CalibObjectsTrain1.root");
  45   AliTPCcalibTime *calib = (AliTPCcalibTime *)f->Get("calibTime");
  46   calib->GetHistoDrift("all")->Projection(2,0)->Draw()
  47   calib->GetFitDrift("all")->Draw("lp")
  48
  49 4. Analysis using debug streamers.
  50
  51   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  52   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  53   AliXRDPROOFtoolkit tool;
  54   AliXRDPROOFtoolkit::FilterList("timeitstpc.txt","* itstpc",1)
  55   AliXRDPROOFtoolkit::FilterList("timetoftpc.txt","* pointMatch",1)
  56   AliXRDPROOFtoolkit::FilterList("time.txt","* trackInfo",1)
  57   AliXRDPROOFtoolkit::FilterList("timelaser.txt","* laserInfo",1)
  58
  59   TChain * chainTPCITS = tool.MakeChainRandom("timeitstpc.txt.Good","itstpc",0,10000);
  60   TChain * chainTPCTOF = tool.MakeChainRandom("timetoftpc.txt.Good","pointMatch",0,500);
  61   TChain * chainTime = tool.MakeChainRandom("time.txt.Good","trackInfo",0,10000);
  62   TChain * chainLaser = tool.MakeChainRandom("timelaser.txt.Good","laserInfo",0,10000);
  63   chainTime->Lookup();
  64   chainLaser->Lookup();
  65 */
  66
  67 #include "Riostream.h"
  68 #include "TChain.h"
  69 #include "TTree.h"
  70 #include "TH1F.h"
  71 #include "TH2F.h"
  72 #include "TH3F.h"
  73 #include "THnSparse.h"
  74 #include "TList.h"
  75 #include "TMath.h"
  76 #include "TCanvas.h"
  77 #include "TFile.h"
  78 #include "TF1.h"
  79 #include "TVectorD.h"
  80 #include "TProfile.h"
  81 #include "TGraphErrors.h"
  82 #include "TCanvas.h"
  83 #include "AliTPCclusterMI.h"
  84 #include "AliTPCseed.h"
  85 #include "AliESDVertex.h"
  86 #include "AliESDEvent.h"
  87 #include "AliESDfriend.h"
  88 #include "AliESDInputHandler.h"
  89 #include "AliAnalysisManager.h"
  90
  91 #include "AliTracker.h"
  92 #include "AliMagF.h"
  93 #include "AliTPCCalROC.h"
  94
  95 #include "AliLog.h"
  96
  97 #include "AliTPCcalibTime.h"
  98 #include "AliRelAlignerKalman.h"
  99
 100 #include "TTreeStream.h"
 101 #include "AliTPCTracklet.h"
 102 #include "TTimeStamp.h"
 103 #include "AliTPCcalibDB.h"
 104 #include "AliTPCcalibLaser.h"
 105 #include "AliDCSSensorArray.h"
 106 #include "AliDCSSensor.h"
 107
 108 #include "TDatabasePDG.h"
 109 #include "AliTrackPointArray.h"
 110
 111 ClassImp(AliTPCcalibTime)
 112
 113
 114 AliTPCcalibTime::AliTPCcalibTime()
 115   :AliTPCcalibBase(),
 116    fLaser(0),       // pointer to laser calibration
 117    fDz(0),          // current delta z
 118    fCutMaxD(3),        // maximal distance in rfi ditection
 119    fCutMaxDz(25),      // maximal distance in rfi ditection
 120    fCutTheta(0.03),    // maximal distan theta
 121    fCutMinDir(-0.99),  // direction vector products
 122    fCutTracks(10),
 123    fArrayDz(0),          //NEW! Tmap of V drifts for different triggers
 124    fAlignITSTPC(0),      //alignemnt array ITS TPC match
 125    fAlignTRDTPC(0),      //alignemnt array TRD TPC match
 126    fAlignTOFTPC(0),      //alignemnt array TOF TPC match
 127    fTimeBins(0),
 128    fTimeStart(0),
 129    fTimeEnd(0),
 130    fPtBins(0),
 131    fPtStart(0),
 132    fPtEnd(0),
 133    fVdriftBins(0),
 134    fVdriftStart(0),
 135    fVdriftEnd(0),
 136    fRunBins(0),
 137    fRunStart(0),
 138    fRunEnd(0)
 139 //   fBinsVdrift(fTimeBins,fPtBins,fVdriftBins),
 140 //   fXminVdrift(fTimeStart,fPtStart,fVdriftStart),
 141 //   fXmaxVdrift(fTimeEnd,fPtEnd,fVdriftEnd)
 142 {
 143   AliInfo("Default Constructor");
 144   for (Int_t i=0;i<3;i++) {
 145     fHistVdriftLaserA[i]=0;
 146     fHistVdriftLaserC[i]=0;
 147   }
 148   for (Int_t i=0;i<10;i++) {
 149     fCosmiMatchingHisto[i]=0;
 150   }
 151 }
 152
 153 AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeVdrift)
 154   :AliTPCcalibBase(),
 155    fLaser(0),            // pointer to laser calibration
 156    fDz(0),               // current delta z
 157    fCutMaxD(5*0.5356),   // maximal distance in rfi ditection
 158    fCutMaxDz(40),   // maximal distance in rfi ditection
 159    fCutTheta(5*0.004644),// maximal distan theta
 160    fCutMinDir(-0.99),    // direction vector products
 161    fCutTracks(10),
 162    fArrayDz(0),            //Tmap of V drifts for different triggers
 163    fAlignITSTPC(0),      //alignemnt array ITS TPC match
 164    fAlignTRDTPC(0),      //alignemnt array TRD TPC match
 165    fAlignTOFTPC(0),      //alignemnt array TOF TPC match
 166    fTimeBins(0),
 167    fTimeStart(0),
 168    fTimeEnd(0),
 169    fPtBins(0),
 170    fPtStart(0),
 171    fPtEnd(0),
 172    fVdriftBins(0),
 173    fVdriftStart(0),
 174    fVdriftEnd(0),
 175    fRunBins(0),
 176    fRunStart(0),
 177    fRunEnd(0)
 178 {
 179   SetName(name);
 180   SetTitle(title);
 181   for (Int_t i=0;i<3;i++) {
 182     fHistVdriftLaserA[i]=0;
 183     fHistVdriftLaserC[i]=0;
 184   }
 185
 186   AliInfo("Non Default Constructor");
 187   fTimeBins   =(EndTime-StartTime)/deltaIntegrationTimeVdrift;
 188   fTimeStart  =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi();
 189   fTimeEnd    =EndTime;   //(((TObjString*)(mapGRP->GetValue("fAliceStopTime")))->GetString()).Atoi();
 190   fPtBins     = 400;
 191   fPtStart    = -0.04;
 192   fPtEnd      =  0.04;
 193   fVdriftBins = 500;
 194   fVdriftStart= -0.1;
 195   fVdriftEnd  =  0.1;
 196   fRunBins    = 100001;
 197   fRunStart   = -1.5;
 198   fRunEnd     = 99999.5;
 199
 200   Int_t    binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins };
 201   Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart  , fRunStart};
 202   Double_t xmaxVdriftLaser[4] = {fTimeEnd  , fPtEnd  , fVdriftEnd    , fRunEnd  };
 203   TString axisTitle[4]={
 204     "T",
 205     "#delta_{P/T}",
 206     "value",
 207     "run"
 208   };
 209   TString histoName[3]={
 210     "Loffset",
 211     "Lcorr",
 212     "Lgy"
 213   };
 214
 215
 216   for (Int_t i=0;i<3;i++) {
 217     fHistVdriftLaserA[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
 218     fHistVdriftLaserC[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
 219     fHistVdriftLaserA[i]->SetName(histoName[i]);
 220     fHistVdriftLaserC[i]->SetName(histoName[i]);
 221     for (Int_t iaxis=0; iaxis<4;iaxis++){
 222       fHistVdriftLaserA[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
 223       fHistVdriftLaserC[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
 224     }
 225   }
 226   fBinsVdrift[0] = fTimeBins;
 227   fBinsVdrift[1] = fPtBins;
 228   fBinsVdrift[2] = fVdriftBins;
 229   fBinsVdrift[3] = fRunBins;
 230   fXminVdrift[0] = fTimeStart;
 231   fXminVdrift[1] = fPtStart;
 232   fXminVdrift[2] = fVdriftStart;
 233   fXminVdrift[3] = fRunStart;
 234   fXmaxVdrift[0] = fTimeEnd;
 235   fXmaxVdrift[1] = fPtEnd;
 236   fXmaxVdrift[2] = fVdriftEnd;
 237   fXmaxVdrift[3] = fRunEnd;
 238
 239   fArrayDz=new TObjArray();
 240   fAlignITSTPC = new TObjArray;      //alignemnt array ITS TPC match
 241   fAlignTRDTPC = new TObjArray;      //alignemnt array ITS TPC match
 242   fAlignTOFTPC = new TObjArray;      //alignemnt array ITS TPC match
 243
 244   // fArrayDz->AddLast(fHistVdriftLaserA[0]);
 245 //   fArrayDz->AddLast(fHistVdriftLaserA[1]);
 246 //   fArrayDz->AddLast(fHistVdriftLaserA[2]);
 247 //   fArrayDz->AddLast(fHistVdriftLaserC[0]);
 248 //   fArrayDz->AddLast(fHistVdriftLaserC[1]);
 249 //   fArrayDz->AddLast(fHistVdriftLaserC[2]);
 250
 251   fCosmiMatchingHisto[0]=new TH1F("Cosmics matching","p0-all"   ,100,-10*0.5356  ,10*0.5356  );
 252   fCosmiMatchingHisto[1]=new TH1F("Cosmics matching","p1-all"   ,100,-10*4.541   ,10*4.541   );
 253   fCosmiMatchingHisto[2]=new TH1F("Cosmics matching","p2-all"   ,100,-10*0.01134 ,10*0.01134 );
 254   fCosmiMatchingHisto[3]=new TH1F("Cosmics matching","p3-all"   ,100,-10*0.004644,10*0.004644);
 255   fCosmiMatchingHisto[4]=new TH1F("Cosmics matching","p4-all"   ,100,-10*0.03773 ,10*0.03773 );
 256   fCosmiMatchingHisto[5]=new TH1F("Cosmics matching","p0-isPair",100,-10*0.5356  ,10*0.5356  );
 257   fCosmiMatchingHisto[6]=new TH1F("Cosmics matching","p1-isPair",100,-10*4.541   ,10*4.541   );
 258   fCosmiMatchingHisto[7]=new TH1F("Cosmics matching","p2-isPair",100,-10*0.01134 ,10*0.01134 );
 259   fCosmiMatchingHisto[8]=new TH1F("Cosmics matching","p3-isPair",100,-10*0.004644,10*0.004644);
 260   fCosmiMatchingHisto[9]=new TH1F("Cosmics matching","p4-isPair",100,-10*0.03773 ,10*0.03773 );
 261 //  Char_t nameHisto[3]={'p','0','\n'};
 262 //  for (Int_t i=0;i<10;i++){
 263 //    fCosmiMatchingHisto[i]=new TH1F("Cosmics matching",nameHisto,8192,0,0);
 264 //    nameHisto[1]++;
 265 //    if(i==4) nameHisto[1]='0';
 266 //  }
 267 }
 268
 269 AliTPCcalibTime::~AliTPCcalibTime(){
 270   //
 271   // Destructor
 272   //
 273   for(Int_t i=0;i<3;i++){
 274     if(fHistVdriftLaserA[i]){
 275       delete fHistVdriftLaserA[i];
 276       fHistVdriftLaserA[i]=NULL;
 277     }
 278     if(fHistVdriftLaserC[i]){
 279       delete fHistVdriftLaserC[i];
 280       fHistVdriftLaserC[i]=NULL;
 281     }
 282   }
 283   if(fArrayDz){
 284     fArrayDz->SetOwner();
 285     fArrayDz->Delete();
 286     delete fArrayDz;
 287     fArrayDz=NULL;
 288   }
 289   for(Int_t i=0;i<5;i++){
 290     if(fCosmiMatchingHisto[i]){
 291       delete fCosmiMatchingHisto[i];
 292       fCosmiMatchingHisto[i]=NULL;
 293     }
 294   }
 295   fAlignITSTPC->Delete();
 296   fAlignTRDTPC->Delete();
 297   fAlignTOFTPC->Delete();
 298 }
 299
 300 Bool_t AliTPCcalibTime::IsLaser(AliESDEvent */*event*/){
 301   return kTRUE; //More accurate creteria to be added
 302 }
 303 Bool_t AliTPCcalibTime::IsCosmics(AliESDEvent */*event*/){
 304   return kTRUE; //More accurate creteria to be added
 305 }
 306 Bool_t AliTPCcalibTime::IsBeam(AliESDEvent */*event*/){
 307   return kTRUE; //More accurate creteria to be added
 308 }
 309 void AliTPCcalibTime::ResetCurrent(){
 310   fDz=0; //Reset current dz
 311 }
 312 void AliTPCcalibTime::Process(AliESDEvent *event){
 313   if(!event) return;
 314   if (event->GetNumberOfTracks()<2) return;
 315   ResetCurrent();
 316   if(IsLaser  (event)) ProcessLaser (event);
 317   if(IsCosmics(event)) ProcessCosmic(event);
 318   if(IsBeam   (event)) ProcessBeam  (event);
 319 }
 320
 321 void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){
 322   //
 323   // Fit drift velocity using laser
 324   //
 325   // 0. cuts
 326   const Int_t    kMinTracks     = 40;    // minimal number of laser tracks
 327   const Int_t    kMinTracksSide = 20;    // minimal number of tracks per side
 328   const Float_t  kMaxDeltaZ     = 30.;   // maximal trigger delay
 329   const Float_t  kMaxDeltaV     = 0.05;  // maximal deltaV
 330   const Float_t  kMaxRMS        = 0.1;   // maximal RMS of tracks
 331   //
 332   /*
 333     TCut cutRMS("sqrt(laserA.fElements[4])<0.1&&sqrt(laserC.fElements[4])<0.1");
 334     TCut cutZ("abs(laserA.fElements[0]-laserC.fElements[0])<3");
 335     TCut cutV("abs(laserA.fElements[1]-laserC.fElements[1])<0.01");
 336     TCut cutY("abs(laserA.fElements[2]-laserC.fElements[2])<2");
 337     TCut cutAll = cutRMS+cutZ+cutV+cutY;
 338   */
 339   if (event->GetNumberOfTracks()<kMinTracks) return;
 340   //
 341   if(!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
 342   fLaser->Process(event);
 343   if (fLaser->GetNtracks()<kMinTracks) return;   // small amount of tracks cut
 344   if (fLaser->fFitAside->GetNrows()==0  && fLaser->fFitCside->GetNrows()==0) return;  // no fit neither a or C side
 345   //
 346   // debug streamer  - activate stream level
 347   // Use it for tuning of the cuts
 348   //
 349   // cuts to be applied
 350   //
 351   Int_t isReject[2]={0,0};
 352   //
 353   // not enough tracks
 354   if (TMath::Abs((*fLaser->fFitAside)[3]) < kMinTracksSide) isReject[0]|=1;
 355   if (TMath::Abs((*fLaser->fFitCside)[3]) < kMinTracksSide) isReject[1]|=1;
 356   // unreasonable z offset
 357   if (TMath::Abs((*fLaser->fFitAside)[0])>kMaxDeltaZ)  isReject[0]|=2;
 358   if (TMath::Abs((*fLaser->fFitCside)[0])>kMaxDeltaZ)  isReject[1]|=2;
 359   // unreasonable drift velocity
 360   if (TMath::Abs((*fLaser->fFitAside)[1]-1)>kMaxDeltaV)  isReject[0]|=4;
 361   if (TMath::Abs((*fLaser->fFitCside)[1]-1)>kMaxDeltaV)  isReject[1]|=4;
 362   // big chi2
 363   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitAside)[4]))>kMaxRMS ) isReject[0]|=8;
 364   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitCside)[4]))>kMaxRMS ) isReject[1]|=8;
 365
 366
 367
 368
 369   if (fStreamLevel>0){
 370     printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
 371
 372     TTreeSRedirector *cstream = GetDebugStreamer();
 373     if (cstream){
 374       TTimeStamp tstamp(fTime);
 375       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 376       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 377       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 378       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 379       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 380       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 381       TVectorD vecGoofie(20);
 382       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 383       if (goofieArray){
 384         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 385           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 386           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 387         }
 388       }
 389       (*cstream)<<"laserInfo"<<
 390         "run="<<fRun<<              //  run number
 391         "event="<<fEvent<<          //  event number
 392         "time="<<fTime<<            //  time stamp of event
 393         "trigger="<<fTrigger<<      //  trigger
 394         "mag="<<fMagF<<             //  magnetic field
 395         // Environment values
 396         "press0="<<valuePressure0<<
 397         "press1="<<valuePressure1<<
 398         "pt0="<<ptrelative0<<
 399         "pt1="<<ptrelative1<<
 400         "temp0="<<temp0<<
 401         "temp1="<<temp1<<
 402         "vecGoofie.="<<&vecGoofie<<
 403         //laser
 404         "rejectA="<<isReject[0]<<
 405         "rejectC="<<isReject[1]<<
 406         "laserA.="<<fLaser->fFitAside<<
 407         "laserC.="<<fLaser->fFitCside<<
 408         "laserAC.="<<fLaser->fFitACside<<
 409         "trigger="<<event->GetFiredTriggerClasses()<<
 410         "\n";
 411     }
 412   }
 413   //
 414   // fill histos
 415   //
 416   TVectorD vdriftA(5), vdriftC(5),vdriftAC(5);
 417   vdriftA=*(fLaser->fFitAside);
 418   vdriftC=*(fLaser->fFitCside);
 419   vdriftAC=*(fLaser->fFitACside);
 420   Int_t npointsA=0, npointsC=0;
 421   Float_t chi2A=0, chi2C=0;
 422   npointsA= TMath::Nint(vdriftA[3]);
 423   chi2A= vdriftA[4];
 424   npointsC= TMath::Nint(vdriftC[3]);
 425   chi2C= vdriftC[4];
 426
 427   TTimeStamp tstamp(fTime);
 428   Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 429   Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 430   Double_t driftA=0, driftC=0;
 431   if (vdriftA[1]>1.-kMaxDeltaV) driftA = 1./vdriftA[1]-1.;
 432   if (vdriftC[1]>1.-kMaxDeltaV) driftC = 1./vdriftC[1]-1.;
 433   //
 434   Double_t vecDriftLaserA[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftA,event->GetRunNumber()};
 435   Double_t vecDriftLaserC[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftC,event->GetRunNumber()};
 436   //  Double_t vecDrift[4]      ={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitACside))[1])-1,event->GetRunNumber()};
 437
 438   for (Int_t icalib=0;icalib<3;icalib++){
 439     if (icalib==0){ //z0 shift
 440       vecDriftLaserA[2]=vdriftA[0]/250.;
 441       vecDriftLaserC[2]=vdriftC[0]/250.;
 442     }
 443     if (icalib==1){ //vdrel shift
 444       vecDriftLaserA[2]=driftA;
 445       vecDriftLaserC[2]=driftC;
 446     }
 447     if (icalib==2){ //gy shift - full gy - full drift
 448       vecDriftLaserA[2]=vdriftA[2]/250.;
 449       vecDriftLaserC[2]=vdriftC[2]/250.;
 450     }
 451     if (isReject[0]==0) fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
 452     if (isReject[1]==0) fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
 453   }
 454
 455 //   THnSparse* curHist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 456 //   TString shortName=curHist->ClassName();
 457 //   shortName+="_MEAN_DRIFT_LASER_";
 458 //   delete curHist;
 459 //   curHist=NULL;
 460 //   TString name="";
 461
 462 //   name=shortName;
 463 //   name+=event->GetFiredTriggerClasses();
 464 //   name.ToUpper();
 465 //   curHist=(THnSparseF*)fArrayDz->FindObject(name);
 466 //   if(!curHist){
 467 //     curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 468 //     fArrayDz->AddLast(curHist);
 469 //   }
 470 //   curHist->Fill(vecDrift);
 471
 472 //   name=shortName;
 473 //   name+="ALL";
 474 //   name.ToUpper();
 475 //   curHist=(THnSparseF*)fArrayDz->FindObject(name);
 476 //   if(!curHist){
 477 //     curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 478 //     fArrayDz->AddLast(curHist);
 479 //   }
 480 //   curHist->Fill(vecDrift);
 481 }
 482
 483 void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){
 484   if (!event) {
 485     Printf("ERROR: ESD not available");
 486     return;
 487   }
 488   if (event->GetTimeStamp() == 0 ) {
 489     Printf("no time stamp!");
 490     return;
 491   }
 492
 493   //fd
 494   // Find cosmic pairs
 495   //
 496   // Track0 is choosen in upper TPC part
 497   // Track1 is choosen in lower TPC part
 498   //
 499   Int_t ntracks=event->GetNumberOfTracks();
 500   if (ntracks==0) return;
 501   if (ntracks > fCutTracks) return;
 502
 503   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
 504   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 505
 506   TObjArray  tpcSeeds(ntracks);
 507   Double_t vtxx[3]={0,0,0};
 508   Double_t svtxx[3]={0.000001,0.000001,100.};
 509   AliESDVertex vtx(vtxx,svtxx);
 510   //
 511   // track loop
 512   //
 513   for (Int_t i=0;i<ntracks;++i) {
 514     AliESDtrack *track = event->GetTrack(i);
 515
 516     const AliExternalTrackParam * trackIn = track->GetInnerParam();
 517     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 518     if (!trackIn) continue;
 519     if (!trackOut) continue;
 520
 521     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 522     if (friendTrack) ProcessAlignITS(track,friendTrack);
 523     if (friendTrack) ProcessAlignTRD(track,friendTrack);
 524     if (friendTrack) ProcessAlignTOF(track,friendTrack);
 525     TObject *calibObject;
 526     AliTPCseed *seed = 0;
 527     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 528     if (seed) tpcSeeds.AddAt(seed,i);
 529   }
 530   if (ntracks<2) return;
 531   //
 532   // Find pairs
 533   //
 534
 535   for (Int_t i=0;i<ntracks;++i) {
 536     AliESDtrack *track0 = event->GetTrack(i);
 537     // track0 - choosen upper part
 538     if (!track0) continue;
 539     if (!track0->GetOuterParam()) continue;
 540     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 541     Double_t d1[3];
 542     track0->GetDirection(d1);
 543     for (Int_t j=0;j<ntracks;++j) {
 544       if (i==j) continue;
 545       AliESDtrack *track1 = event->GetTrack(j);
 546       //track 1 lower part
 547       if (!track1) continue;
 548       if (!track1->GetOuterParam()) continue;
 549       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 550       //
 551       Double_t d2[3];
 552       track1->GetDirection(d2);
 553
 554       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 555       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 556       if (! seed0) continue;
 557       if (! seed1) continue;
 558       Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
 559       Float_t dist0  = track0->GetLinearD(0,0);
 560       Float_t dist1  = track1->GetLinearD(0,0);
 561       //
 562       // conservative cuts - convergence to be guarantied
 563       // applying before track propagation
 564       if (TMath::Abs(dist0+dist1)>fCutMaxD) continue;   // distance to the 0,0
 565       if (dir>fCutMinDir) continue;               // direction vector product
 566       Float_t bz = AliTracker::GetBz();
 567       Float_t dvertex0[2];   //distance to 0,0
 568       Float_t dvertex1[2];   //distance to 0,0
 569       track0->GetDZ(0,0,0,bz,dvertex0);
 570       track1->GetDZ(0,0,0,bz,dvertex1);
 571       if (TMath::Abs(dvertex0[1])>250) continue;
 572       if (TMath::Abs(dvertex1[1])>250) continue;
 573       //
 574       //
 575       //
 576       Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
 577       AliExternalTrackParam param0(*track0);
 578       AliExternalTrackParam param1(*track1);
 579       //
 580       // Propagate using Magnetic field and correct fo material budget
 581       //
 582       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
 583       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
 584       //
 585       // Propagate rest to the 0,0 DCA - z should be ignored
 586       //
 587       //Bool_t b0 = ;
 588       param0.PropagateToDCA(&vtx,bz,1000);
 589       //Bool_t b1 =
 590       param1.PropagateToDCA(&vtx,bz,1000);
 591       param0.GetDZ(0,0,0,bz,dvertex0);
 592       param1.GetDZ(0,0,0,bz,dvertex1);
 593       Double_t xyz0[3];
 594       Double_t xyz1[3];
 595       param0.GetXYZ(xyz0);
 596       param1.GetXYZ(xyz1);
 597       Bool_t isPair = IsPair(&param0,&param1);
 598       Bool_t isCross = IsCross(track0, track1);
 599       Bool_t isSame = IsSame(track0, track1);
 600
 601       THnSparse* hist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 602       TString shortName=hist->ClassName();
 603       shortName+="_MEAN_VDRIFT_COSMICS_";
 604       delete hist;
 605       hist=NULL;
 606
 607       if(isSame || (isCross && isPair)){
 608         if (track0->GetTPCNcls() > 80) {
 609           fDz = param0.GetZ() - param1.GetZ();
 610           if(track0->GetOuterParam()->GetZ()<0) fDz=-fDz;
 611           TTimeStamp tstamp(fTime);
 612           Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 613           Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 614           Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
 615           THnSparse* curHist=NULL;
 616           TString name="";
 617
 618           name=shortName;
 619           name+=event->GetFiredTriggerClasses();
 620           name.ToUpper();
 621           curHist=(THnSparseF*)fArrayDz->FindObject(name);
 622           if(!curHist){
 623             curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 624             fArrayDz->AddLast(curHist);
 625           }
 626 //        curHist=(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
 627 //        if(!curHist){
 628 //          curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 629 //          fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
 630 //        }
 631           curHist->Fill(vecDrift);
 632
 633           name=shortName;
 634           name+="ALL";
 635           name.ToUpper();
 636           curHist=(THnSparseF*)fArrayDz->FindObject(name);
 637           if(!curHist){
 638             curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 639             fArrayDz->AddLast(curHist);
 640           }
 641 //        curHist=(THnSparseF*)(fMapDz->GetValue("all"));
 642 //        if(!curHist){
 643 //          curHist=new THnSparseF("all","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 644 //          fMapDz->Add(new TObjString("all"),curHist);
 645 //        }
 646           curHist->Fill(vecDrift);
 647         }
 648       }
 649       TTreeSRedirector *cstream = GetDebugStreamer();
 650       if (fStreamLevel>0){
 651         if (cstream){
 652         (*cstream)<<"trackInfo"<<
 653           "tr0.="<<track0<<
 654           "tr1.="<<track1<<
 655           "p0.="<<&param0<<
 656           "p1.="<<&param1<<
 657           "isPair="<<isPair<<
 658           "isCross="<<isCross<<
 659           "isSame="<<isSame<<
 660           "fDz="<<fDz<<
 661           "fRun="<<fRun<<
 662           "fTime="<<fTime<<
 663           "\n";
 664         }
 665       }
 666     } // end 2nd order loop
 667   } // end 1st order loop
 668
 669   if (fStreamLevel>0){
 670     TTreeSRedirector *cstream = GetDebugStreamer();
 671     if (cstream){
 672       TTimeStamp tstamp(fTime);
 673       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 674       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 675       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 676       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 677       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 678       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 679       TVectorD vecGoofie(20);
 680       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 681       if (goofieArray){
 682         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 683           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 684           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 685         }
 686       }
 687       (*cstream)<<"timeInfo"<<
 688         "run="<<fRun<<              //  run number
 689         "event="<<fEvent<<          //  event number
 690         "time="<<fTime<<            //  time stamp of event
 691         "trigger="<<fTrigger<<      //  trigger
 692         "mag="<<fMagF<<             //  magnetic field
 693         // Environment values
 694         "press0="<<valuePressure0<<
 695         "press1="<<valuePressure1<<
 696         "pt0="<<ptrelative0<<
 697         "pt1="<<ptrelative1<<
 698         "temp0="<<temp0<<
 699         "temp1="<<temp1<<
 700         "vecGoofie.=<<"<<&vecGoofie<<
 701         //
 702         // accumulated values
 703         //
 704         "fDz="<<fDz<<          //! current delta z
 705         "trigger="<<event->GetFiredTriggerClasses()<<
 706         "\n";
 707     }
 708   }
 709   printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
 710 }
 711
 712 void AliTPCcalibTime::ProcessBeam(AliESDEvent */*event*/){
 713 }
 714
 715 void AliTPCcalibTime::Analyze(){}
 716
 717 THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name){
 718   TIterator* iterator = fArrayDz->MakeIterator();
 719   iterator->Reset();
 720   TString newName=name;
 721   newName.ToUpper();
 722   THnSparse* newHist=new THnSparseF(newName,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 723   THnSparse* addHist=NULL;
 724   while((addHist=(THnSparseF*)iterator->Next())){
 725   if(!addHist) continue;
 726     TString histName=addHist->GetName();
 727     if(!histName.Contains(newName)) continue;
 728     addHist->Print();
 729     newHist->Add(addHist);
 730   }
 731   return newHist;
 732 }
 733
 734 TObjArray* AliTPCcalibTime::GetHistoDrift(){
 735   return fArrayDz;
 736 }
 737
 738 TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){
 739   THnSparse* histoDrift=GetHistoDrift(name);
 740   TGraphErrors* graphDrift=NULL;
 741   if(histoDrift){
 742     graphDrift=FitSlices(histoDrift,2,0,400,100,0.05,0.95, kTRUE);
 743     TString end=histoDrift->GetName();
 744     Int_t pos=end.Index("_");
 745     end=end(pos,end.Capacity()-pos);
 746     TString graphName=graphDrift->ClassName();
 747     graphName+=end;
 748     graphName.ToUpper();
 749     graphDrift->SetName(graphName);
 750   }
 751   return graphDrift;
 752 }
 753
 754 TObjArray* AliTPCcalibTime::GetGraphDrift(){
 755   TObjArray* arrayGraphDrift=new TObjArray();
 756   TIterator* iterator=fArrayDz->MakeIterator();
 757   iterator->Reset();
 758   THnSparse* addHist=NULL;
 759   while((addHist=(THnSparseF*)iterator->Next())) arrayGraphDrift->AddLast(GetGraphDrift(addHist->GetName()));
 760   return arrayGraphDrift;
 761 }
 762
 763 AliSplineFit* AliTPCcalibTime::GetFitDrift(const char* name){
 764   TGraph* graphDrift=GetGraphDrift(name);
 765   AliSplineFit* fitDrift=NULL;
 766   if(graphDrift && graphDrift->GetN()){
 767     fitDrift=new AliSplineFit();
 768     fitDrift->SetGraph(graphDrift);
 769     fitDrift->SetMinPoints(graphDrift->GetN()+1);
 770     fitDrift->InitKnots(graphDrift,2,0,0.001);
 771     fitDrift->SplineFit(0);
 772     TString end=graphDrift->GetName();
 773     Int_t pos=end.Index("_");
 774     end=end(pos,end.Capacity()-pos);
 775     TString fitName=fitDrift->ClassName();
 776     fitName+=end;
 777     fitName.ToUpper();
 778     //fitDrift->SetName(fitName);
 779     delete graphDrift;
 780     graphDrift=NULL;
 781   }
 782   return fitDrift;
 783 }
 784
 785 //TObjArray* AliTPCcalibTime::GetFitDrift(){
 786 //  TObjArray* arrayFitDrift=new TObjArray();
 787 //  TIterator* iterator = fArrayDz->MakeIterator();
 788 //  iterator->Reset();
 789 //  THnSparse* addHist=NULL;
 790 //  while((addHist=(THnSparseF*)iterator->Next())) arrayFitDrift->AddLast(GetFitDrift(addHist->GetName()));
 791 //  return arrayFitDrift;
 792 //}
 793
 794 Long64_t AliTPCcalibTime::Merge(TCollection *li) {
 795   TIterator* iter = li->MakeIterator();
 796   AliTPCcalibTime* cal = 0;
 797
 798   while ((cal = (AliTPCcalibTime*)iter->Next())) {
 799     if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
 800       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 801       return -1;
 802     }
 803     for (Int_t imeas=0; imeas<3; imeas++){
 804       if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
 805         fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
 806         fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
 807       }
 808     }
 809     TObjArray* addArray=cal->GetHistoDrift();
 810     if(!addArray) return 0;
 811     TIterator* iterator = addArray->MakeIterator();
 812     iterator->Reset();
 813     THnSparse* addHist=NULL;
 814     while((addHist=(THnSparseF*)iterator->Next())){
 815       if(!addHist) continue;
 816       addHist->Print();
 817       THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName());
 818       if(!localHist){
 819         localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 820         fArrayDz->AddLast(localHist);
 821       }
 822       localHist->Add(addHist);
 823     }
 824 //    TMap * addMap=cal->GetHistoDrift();
 825 //    if(!addMap) return 0;
 826 //    TIterator* iterator = addMap->MakeIterator();
 827 //    iterator->Reset();
 828 //    TPair* addPair=0;
 829 //    while((addPair=(TPair *)(addMap->FindObject(iterator->Next())))){
 830 //      THnSparse* addHist=dynamic_cast<THnSparseF*>(addPair->Value());
 831 //      if (!addHist) continue;
 832 //      addHist->Print();
 833 //      THnSparse* localHist=dynamic_cast<THnSparseF*>(fMapDz->GetValue(addHist->GetName()));
 834 //      if(!localHist){
 835 //        localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
 836 //        fMapDz->Add(new TObjString(addHist->GetName()),localHist);
 837 //      }
 838 //      localHist->Add(addHist);
 839 //    }
 840     for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
 841   }
 842   return 0;
 843 }
 844
 845 Bool_t  AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 846   /*
 847   // 0. Same direction - OPOSITE  - cutDir +cutT
 848   TCut cutDir("cutDir","dir<-0.99")
 849   // 1.
 850   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 851   //
 852   // 2. The same rphi
 853   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 854   //
 855   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 856   // 1/Pt diff cut
 857   */
 858   const Double_t *p0 = tr0->GetParameter();
 859   const Double_t *p1 = tr1->GetParameter();
 860   fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
 861   fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
 862   fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
 863   fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
 864   fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
 865
 866   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 867   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 868   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz)  return kFALSE;
 869   Double_t d0[3], d1[3];
 870   tr0->GetDirection(d0);
 871   tr1->GetDirection(d1);
 872   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 873
 874   fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
 875   fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
 876   fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
 877   fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
 878   fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
 879
 880   return kTRUE;
 881 }
 882 Bool_t AliTPCcalibTime::IsCross(AliESDtrack *tr0, AliESDtrack *tr1){
 883   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;
 884 }
 885
 886 Bool_t AliTPCcalibTime::IsSame(AliESDtrack */*tr0*/, AliESDtrack */*tr1*/){
 887   // To be implemented
 888   return kFALSE;
 889 }
 890
 891 /*
 892 chainDrift->Draw("p0.fP[0]+p1.fP[0]","isPair");
 893   mean ~-0.02  ~-0.03913
 894   RMS  ~ 0.5   ~ 0.5356    --> 3    (fCutMaxD)
 895
 896 chainDrift->Draw("p0.fP[1]-p1.fP[1]","isPair");
 897   mean         ~ 0.1855
 898   RMS          ~ 4.541     -->25    (fCutMaxDz)
 899
 900 chainDrift->Draw("p1.fAlpha-p0.fAlpha+pi","isPair");
 901 //chainDrift->Draw("p1.fAlpha+p0.fAlpha","isPair");
 902 //chainDrift->Draw("p1.fP[2]-p0.fP[2]+pi","isPair");
 903 //chainDrift->Draw("p1.fP[2]+p0.fP[2]","isPair");
 904   mean ~ 0     ~ 0.001898
 905   RMS  ~ 0.009 ~ 0.01134   --> 0.06
 906
 907 chainDrift->Draw("p0.fP[3]+p1.fP[3]","isPair");
 908   mean ~ 0.0013 ~ 0.001539
 909   RMS  ~ 0.003  ~ 0.004644 --> 0.03  (fCutTheta)
 910
 911 chainDrift->Draw("p1.fP[4]+p0.fP[4]>>his(100,-0.2,0.2)","isPair")
 912   mean ~ 0.012  ~-0.0009729
 913   RMS  ~ 0.036  ~ 0.03773  --> 0.2
 914 */
 915
 916
 917 void  AliTPCcalibTime::ProcessAlignITS(AliESDtrack* track, AliESDfriendTrack *friendTrack){
 918   //
 919   // Process track
 920   // Update TPC-ITS alignment
 921   //
 922   const Int_t    kMinTPC  = 80;
 923   const Int_t    kMinITS  = 3;
 924   const Double_t kMinZ    = 10;
 925   const Double_t kMaxDy   = 2;
 926   const Double_t kMaxAngle= 0.02;
 927   //
 928   Int_t dummycl[1000];
 929   if (track->GetITSclusters(dummycl)<kMinITS) return;  // minimal amount of clusters
 930   if (track->GetTPCNcls()<kMinTPC) return;  // minimal amount of clusters cut
 931   //
 932   if (!friendTrack->GetITSOut()) return;
 933   if (!track->GetInnerParam())   return;
 934   if (!track->GetOuterParam())   return;
 935   // exclude crossing track
 936   if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0)   return;
 937   if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ)   return;
 938   //
 939   AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetInnerParam()));
 940   AliExternalTrackParam pITS(*(friendTrack->GetITSOut()));
 941   //
 942   //
 943   //
 944   Int_t htime = fTime/3600; //time in hours
 945   if (fAlignITSTPC->GetEntries()<htime){
 946     fAlignITSTPC->Expand(htime*2+20);
 947   }
 948   AliRelAlignerKalman* align =  (AliRelAlignerKalman*)fAlignITSTPC->At(htime);
 949   if (!align){
 950     align=new AliRelAlignerKalman();
 951     align->SetOutRejSigma(2.);
 952     //align->SetRejectOutliers(kFALSE);
 953     align->SetMagField(fMagF);
 954     fAlignITSTPC->AddAt(align,htime);
 955   }
 956   pITS.Rotate(pTPC.GetAlpha());
 957   pITS.PropagateTo(pTPC.GetX(),fMagF);
 958   if (TMath::Abs(pITS.GetY()-pTPC.GetY())>kMaxDy) return;
 959   if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
 960   if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
 961   align->AddTrackParams(&pITS,&pTPC);
 962   align->SetTimeStamp(fTime);
 963   //  align->SetRunNumber(fRun );
 964   static Int_t entry=0;
 965   entry++;
 966   //  Int_t nupdates=align->GetNUpdates();
 967   Int_t nupdates=entry;
 968   align->SetOutRejSigma(1.+1./Double_t(nupdates));
 969   TTreeSRedirector *cstream = GetDebugStreamer();
 970   if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
 971     TTimeStamp tstamp(fTime);
 972     Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 973     Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 974     Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 975     Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 976     Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 977     Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 978     TVectorD vecGoofie(20);
 979     AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 980     if (goofieArray){
 981       for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 982         AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 983         if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 984       }
 985     }
 986     TVectorD gpTPC(3), gdTPC(3);
 987     TVectorD gpITS(3), gdITS(3);
 988     pTPC.GetXYZ(gpTPC.GetMatrixArray());
 989     pTPC.GetDirection(gdTPC.GetMatrixArray());
 990     pITS.GetXYZ(gpITS.GetMatrixArray());
 991     pITS.GetDirection(gdITS.GetMatrixArray());
 992     (*cstream)<<"itstpc"<<
 993       "run="<<fRun<<              //  run number
 994       "event="<<fEvent<<          //  event number
 995       "time="<<fTime<<            //  time stamp of event
 996       "trigger="<<fTrigger<<      //  trigger
 997       "mag="<<fMagF<<             //  magnetic field
 998       // Environment values
 999       "press0="<<valuePressure0<<
1000       "press1="<<valuePressure1<<
1001       "pt0="<<ptrelative0<<
1002       "pt1="<<ptrelative1<<
1003       "temp0="<<temp0<<
1004       "temp1="<<temp1<<
1005       "vecGoofie.="<<&vecGoofie<<
1006       "entry="<<entry<<  // current entry
1007       //
1008       "a.="<<align<<     // current alignment
1009       "pITS.="<<&pITS<<  // track param ITS
1010       "pTPC.="<<&pTPC<<  // track param TPC
1011       "gpTPC.="<<&gpTPC<<
1012       "gdTPC.="<<&gdTPC<<
1013       "gpITS.="<<&gpITS<<
1014       "gdITS.="<<&gdITS<<
1015       "\n";
1016   }
1017
1018 }
1019 void  AliTPCcalibTime::ProcessAlignTRD(AliESDtrack* track, AliESDfriendTrack *friendTrack){
1020   //
1021   // Process track
1022   // Update TPC-TRD alignment
1023   //
1024   const Int_t    kMinTPC  = 80;
1025   const Int_t    kMinTRD  = 60;
1026   const Double_t kMinZ    = 10;
1027   const Double_t kMaxDy   = 2;
1028   const Double_t kMaxAngle= 0.02;
1029   //
1030   Int_t dummycl[1000];
1031   if (track->GetTRDclusters(dummycl)<kMinTRD) return;  // minimal amount of clusters
1032   if (track->GetTPCNcls()<kMinTPC) return;  // minimal amount of clusters cut
1033   //
1034   if (!friendTrack->GetTRDIn()) return;
1035   if (!track->GetInnerParam())   return;
1036   if (!track->GetOuterParam())   return;
1037   // exclude crossing track
1038   if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0)   return;
1039   if (TMath::Abs(track->GetOuterParam()->GetZ())<kMinZ)   return;
1040   //
1041   AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetOuterParam()));
1042   AliExternalTrackParam pTRD(*(friendTrack->GetTRDIn()));
1043   //
1044   //
1045   //
1046   Int_t htime = fTime/3600; //time in hours
1047   if (fAlignTRDTPC->GetEntries()<htime){
1048     fAlignTRDTPC->Expand(htime*2+20);
1049   }
1050   AliRelAlignerKalman* align =  (AliRelAlignerKalman*)fAlignTRDTPC->At(htime);
1051   if (!align){
1052     align=new AliRelAlignerKalman();
1053     align->SetOutRejSigma(2.);
1054     //align->SetRejectOutliers(kFALSE);
1055     align->SetMagField(fMagF);
1056     fAlignTRDTPC->AddAt(align,htime);
1057   }
1058   pTRD.Rotate(pTPC.GetAlpha());
1059   pTRD.PropagateTo(pTPC.GetX(),fMagF);
1060   ((Double_t*)pTRD.GetCovariance())[2]+=3.*3.;
1061   ((Double_t*)pTRD.GetCovariance())[9]+=0.1*0.1;
1062
1063   if (TMath::Abs(pTRD.GetY()-pTPC.GetY())>kMaxDy) return;
1064   if (TMath::Abs(pTRD.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
1065   if (TMath::Abs(pTRD.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
1066   align->AddTrackParams(&pTRD,&pTPC);
1067   align->SetTimeStamp(fTime);
1068   //  align->SetRunNumber(fRun );
1069   static Int_t entry=0;
1070   entry++;
1071   //  Int_t nupdates=align->GetNUpdates();
1072   Int_t nupdates=entry;
1073   align->SetOutRejSigma(1.+1./Double_t(nupdates));
1074   TTreeSRedirector *cstream = GetDebugStreamer();
1075   if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
1076     TTimeStamp tstamp(fTime);
1077     Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1078     Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1079     Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1080     Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1081     Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1082     Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1083     TVectorD vecGoofie(20);
1084     AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1085     if (goofieArray){
1086       for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1087         AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1088         if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1089       }
1090     }
1091     TVectorD gpTPC(3), gdTPC(3);
1092     TVectorD gpTRD(3), gdTRD(3);
1093     pTPC.GetXYZ(gpTPC.GetMatrixArray());
1094     pTPC.GetDirection(gdTPC.GetMatrixArray());
1095     pTRD.GetXYZ(gpTRD.GetMatrixArray());
1096     pTRD.GetDirection(gdTRD.GetMatrixArray());
1097     (*cstream)<<"itstpc"<<
1098       "run="<<fRun<<              //  run number
1099       "event="<<fEvent<<          //  event number
1100       "time="<<fTime<<            //  time stamp of event
1101       "trigger="<<fTrigger<<      //  trigger
1102       "mag="<<fMagF<<             //  magnetic field
1103       // Environment values
1104       "press0="<<valuePressure0<<
1105       "press1="<<valuePressure1<<
1106       "pt0="<<ptrelative0<<
1107       "pt1="<<ptrelative1<<
1108       "temp0="<<temp0<<
1109       "temp1="<<temp1<<
1110       "vecGoofie.="<<&vecGoofie<<
1111       "entry="<<entry<<  // current entry
1112       //
1113       "a.="<<align<<     // current alignment
1114       "pTRD.="<<&pTRD<<  // track param TRD
1115       "pTPC.="<<&pTPC<<  // track param TPC
1116       "gpTPC.="<<&gpTPC<<
1117       "gdTPC.="<<&gdTPC<<
1118       "gpTRD.="<<&gpTRD<<
1119       "gdTRD.="<<&gdTRD<<
1120       "\n";
1121   }
1122
1123 }
1124
1125
1126 void  AliTPCcalibTime::ProcessAlignTOF(AliESDtrack* track, AliESDfriendTrack *friendTrack){
1127   //
1128   //process TOF-TPC  alignment
1129   //
1130   Int_t kminNcl=80;
1131   Float_t kMaxDy=6;
1132   Float_t kMaxDz=10;
1133   if (track->GetTPCNcls()<kminNcl) return;
1134   if (track->GetOuterParam()==0) return;
1135   if (track->GetInnerParam()==0) return;
1136   if (track->GetTOFsignal()<=0)  return;
1137   //
1138   AliExternalTrackParam *paramOut = new AliExternalTrackParam(*(track->GetOuterParam()));
1139   AliExternalTrackParam *param=0;
1140   const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
1141   if (!points) return;
1142   Int_t npoints = points->GetNPoints();
1143   AliTrackPoint point;
1144   //Double_t alpha=
1145   Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
1146   TTreeSRedirector * cstream =  GetDebugStreamer();
1147   //
1148   //
1149   //
1150   for (Int_t ipoint=0;ipoint<npoints;ipoint++){
1151     //
1152     points->GetPoint(point,ipoint);
1153     Float_t xyz[3];
1154     point.GetXYZ(xyz);
1155     Double_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1156     if (r<300)  continue;
1157     if (r>400)  continue;
1158     param=paramOut;
1159     if (!param) continue;
1160     AliTracker::PropagateTrackToBxByBz(param,r,mass,2.,kTRUE);
1161     AliTracker::PropagateTrackToBxByBz(param,r,mass,0.1,kTRUE);
1162     AliTrackPoint lpoint = point.Rotate(param->GetAlpha());
1163     param->PropagateTo(lpoint.GetX(),fMagF);
1164     //
1165     //
1166     // this is ugly - we need AliCluster constructor
1167     //
1168     AliExternalTrackParam &pTPC=*param;
1169     AliExternalTrackParam pTOF(*param);
1170     ((Double_t*)pTOF.GetParameter())[0] =lpoint.GetY();
1171     ((Double_t*)pTOF.GetParameter())[1] =lpoint.GetZ();
1172     pTOF.ResetCovariance(20);
1173     ((Double_t*)pTOF.GetCovariance())[0]+=3.*3.;
1174     ((Double_t*)pTOF.GetCovariance())[2]+=3.*3.;
1175     ((Double_t*)pTOF.GetCovariance())[5]+=0.1*0.1;
1176     ((Double_t*)pTOF.GetCovariance())[9]+=0.1*0.1;
1177     if (TMath::Abs(pTOF.GetY()-pTPC.GetY())>kMaxDy) continue;
1178     if (TMath::Abs(pTOF.GetZ()-pTPC.GetZ())>kMaxDz) continue;
1179     //
1180     Int_t htime = fTime/3600; //time in hours
1181
1182     if (fAlignTOFTPC->GetEntries()<htime){
1183       fAlignTOFTPC->Expand(htime*2+20);
1184     }
1185     AliRelAlignerKalman* align =  (AliRelAlignerKalman*)fAlignTOFTPC->At(htime);
1186     if (!align){
1187       align=new AliRelAlignerKalman();
1188       align->SetOutRejSigma(2.);
1189       //align->SetRejectOutliers(kFALSE);
1190       align->SetMagField(fMagF);
1191       fAlignTOFTPC->AddAt(align,htime);
1192     }
1193     pTOF.Rotate(pTPC.GetAlpha());
1194     pTOF.PropagateTo(pTPC.GetX(),fMagF);
1195     align->AddTrackParams(&pTOF,&pTPC);
1196     align->SetTimeStamp(fTime);
1197     static Int_t entry=0;
1198     entry++;
1199     //    Int_t nupdates=align->GetNUpdates();
1200     Int_t nupdates=entry;
1201     align->SetOutRejSigma(1.+1./Double_t(nupdates));
1202
1203     //
1204     //
1205     if (cstream) {
1206       (*cstream) << "pointMatch" <<
1207         "run="<<fRun<<              //  run number
1208         "event="<<fEvent<<          //  event number
1209         "time="<<fTime<<            //  time stamp of event
1210         "trigger="<<fTrigger<<      //  trigger
1211         "mag="<<fMagF<<             //  magnetic field
1212         //
1213         "a.="<<align<<     // current alignment
1214         "p.="<<&point<<
1215         "lp.="<<&lpoint<<
1216         "pTPC.="<<&pTPC<<
1217         "pTOF.="<<&pTOF<<
1218         "\n";
1219     }
1220
1221
1222
1223   }
1224   delete paramOut;
1225 }