[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTime.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
    Comments to be written here:
    1. What do we calibrate.

Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.

  AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);


    2. How to interpret results
    3. Simple example

a.) determine the required time range:

AliXRDPROOFtoolkit tool;
TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
chain->Draw("GetTimeStamp()")

b) analyse calibration object on Proof in calibration train 

AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift, IntegrationTimeDeDx);

c) plot results

TFile f("CalibObjects.root");
AliTPCcalibTime *cal = (AliTPCcalibTime *)f->Get("TPCCalib")->FindObject("cosmicTime");
cal->GetHistVdrift()->Projection(1,0)->Draw()

    4. Analysis using debug streamers.    

    gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
    gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
    AliXRDPROOFtoolkit tool;
    TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
    chainTime->Lookup();


*/


#include "Riostream.h"
#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "THnSparse.h"
#include "TList.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TF1.h"
#include "TVectorD.h"
#include "TProfile.h"
#include "TGraphErrors.h"
#include "TCanvas.h"

#include "AliTPCclusterMI.h"
#include "AliTPCseed.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliESDfriend.h"
#include "AliESDInputHandler.h"
#include "AliAnalysisManager.h"

#include "AliTracker.h"
#include "AliMagFMaps.h"
#include "AliTPCCalROC.h"

#include "AliLog.h"

#include "AliTPCcalibTime.h"

#include "TTreeStream.h"
#include "AliTPCTracklet.h"
#include "TTimeStamp.h"
#include "AliTPCcalibDB.h"
#include "AliTPCcalibLaser.h"

ClassImp(AliTPCcalibTime)


AliTPCcalibTime::AliTPCcalibTime() 
  :AliTPCcalibBase(), 
   fTriggerMask(0),
   fHistDeDxTgl(0),
   fHistDeDx(0),
   fHistVdrift(0),
   fIntegrationTimeDeDx(0),
   fIntegrationTimeVdrift(0),  
   fLaser(0),       // pointer to laser calibration
   fDz(0),          // current delta z
   fdEdx(0),        // current dEdx
   fdEdxRatio(0),   // current dEdx ratio
   fTl(0),          // current tan(lambda)
   fCutMaxD(5),        // maximal distance in rfi ditection
   fCutMaxDz(20),        // maximal distance in rfi ditection
   fCutTheta(0.03),    // maximal distan theta
   fCutMinDir(-0.99)   // direction vector products

{  
  AliInfo("Default Constructor");  
}


AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, ULong64_t TriggerMask, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeDeDx, Int_t deltaIntegrationTimeVdrift)
  :AliTPCcalibBase(),
   fTriggerMask(0),
   fHistDeDxTgl(0),
   fHistDeDx(0),
   fHistVdrift(0),
   fIntegrationTimeDeDx(0),
   fIntegrationTimeVdrift(0),  
   fLaser(0),       // pointer to laser calibration
   fDz(0),          // current delta z
   fdEdx(0),        // current dEdx
   fdEdxRatio(0),   // current dEdx ratio
   fTl(0),          // current tan(lambda)
   fCutMaxD(5),        // maximal distance in rfi ditection
   fCutMaxDz(20),        // maximal distance in rfi ditection
   fCutTheta(0.03),    // maximal distan theta
   fCutMinDir(-0.99)   // direction vector products
{
  
  SetName(name);
  SetTitle(title);

  AliInfo("Non Default Constructor");

  fTriggerMask = TriggerMask;

  fIntegrationTimeDeDx = deltaIntegrationTimeDeDx;
  fIntegrationTimeVdrift = deltaIntegrationTimeVdrift;

  Double_t deltaTime = EndTime - StartTime;
  
  Int_t binsVdrift[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeVdrift), 100};
  Double_t xminVdrift[2] = {StartTime, -20};
  Double_t xmaxVdrift[2] = {EndTime, 20};
  fHistVdrift = new THnSparseF("HistVdrift","vDrift; time;#Delta z",2,binsVdrift,xminVdrift,xmaxVdrift);

  Int_t binsDeDxTgl[3] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),30,100};
  Double_t xminDeDxTgl[3] = {StartTime,-1,0.7};
  Double_t xmaxDeDxTgl[3] = {EndTime,1,1.3};
  fHistDeDxTgl = new THnSparseF("HistDeDxTgl","dEdx vs tgl;time;tgl;dEdxUp/dEdxLow",3,binsDeDxTgl,xminDeDxTgl,xmaxDeDxTgl); 

  Int_t binsDeDx[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),100};
  Double_t xminDeDx[2] = {StartTime,1};
  Double_t xmaxDeDx[2] = {EndTime,100};
  fHistDeDx = new THnSparseF("HistDeDx","dEdx l;time;dEdx",2,binsDeDx,xminDeDx,xmaxDeDx);

}


AliTPCcalibTime::~AliTPCcalibTime(){
  //
  //
  //
}
void AliTPCcalibTime::ResetCurrent(){
  //
  // reset current values
  //
  fDz=0;          // current delta z
  fdEdx=0;        // current dEdx
  fdEdxRatio=0;   // current dEdx ratio
  fTl=0;          // current tan(lambda)

}


void AliTPCcalibTime::Process(AliESDEvent *event) {
  //
  //
  //
  Int_t ntracks=event->GetNumberOfTracks();
  if (ntracks<2) return;
  ResetCurrent();
  //
  ProcessCosmic(event);
  if (fTrigger==16){
    if (!fLaser) fLaser =  new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
    fLaser->Process(event);
  }
  //  
  // fill debug streamer
  if (fStreamLevel>0 && fDz!=0){
    TTreeSRedirector *cstream = GetDebugStreamer();
    if (cstream){
      TTimeStamp tstamp(fTime);
      Float_t valuePressure0  = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
      Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
      Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
      Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
      Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
      Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
      TVectorD vdriftA, vdriftC,vdriftAC;
      if (fLaser && fTrigger==16) {
	if (fLaser->fFitAside)  vdriftA=*(fLaser->fFitAside);
	if (fLaser->fFitCside)  vdriftC=*(fLaser->fFitCside);
	if (fLaser->fFitACside) vdriftAC=*(fLaser->fFitACside);
      }
      (*cstream)<<"timeInfo"<<
	"run="<<fRun<<              //  run number
	"event="<<fEvent<<          //  event number
	"time="<<fTime<<            //  time stamp of event
	"trigger="<<fTrigger<<      //  trigger
	"mag="<<fMagF<<             //  magnetic field
	// Environment values
	"press0="<<valuePressure0<<
	"press1="<<valuePressure1<<
	"pt0="<<ptrelative0<<
	"pt1="<<ptrelative1<<
	"temp0="<<temp0<<
	"temp1="<<temp1<<
	//
	// accumulated values
	//
	"fDz="<<fDz<<          //! current delta z
	"fdEdx="<<fdEdx<<        //! current dEdx
	"fdEdxRatio="<<fdEdxRatio<<   //! current dEdx ratio
	"fTl="<<fTl<<          //! current tan(lambda)
	//
	//laser
	//
	"laserA.="<<&vdriftA<<
	"laserC.="<<&vdriftC<<
	"laserAC.="<<&vdriftAC<<
	"\n";
    }
  }
}


void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event) {

  if (!event) {
    Printf("ERROR: ESD not available");
    return;
  }  
  if (event->GetTimeStamp() == 0 ) {
    Printf("no time stamp!");
    return;
  }

  if (fTriggerMask != 0 && event->GetTriggerMask() != fTriggerMask) return;

  UInt_t time = event->GetTimeStamp();

  //
  // Find cosmic pairs
  // 
  // Track0 is choosen in upper TPC part
  // Track1 is choosen in lower TPC part
  //
  Int_t ntracks=event->GetNumberOfTracks();
  if (ntracks==0) return;
  if (ntracks > 10) return; // temporary debug to remove LASER events


  if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
  
  TObjArray  tpcSeeds(ntracks);
  Double_t vtxx[3]={0,0,0};
  Double_t svtxx[3]={0.000001,0.000001,100.};
  AliESDVertex vtx(vtxx,svtxx);
  //
  // track loop
  //
  for (Int_t i=0;i<ntracks;++i) {
   AliESDtrack *track = event->GetTrack(i);

   const AliExternalTrackParam * trackIn = track->GetInnerParam();
   const AliExternalTrackParam * trackOut = track->GetOuterParam();
   if (!trackIn) continue;
   if (!trackOut) continue;
   
   AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
   TObject *calibObject;
   AliTPCseed *seed = 0;
   for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
     if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
   }
   if (seed) { 
     tpcSeeds.AddAt(seed,i);
     if (track->GetTPCNcls() > 50) {
       Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
       Double_t TPCsignal = seed->CookdEdxNorm(0.0,0.6,1,0,159,0x0,kTRUE,kTRUE);
       Double_t vecDeDx[2] = {time, TPCsignal};
       if (meanP > 12) {
	 fdEdx = TPCsignal;
	 fHistDeDx->Fill(vecDeDx);
       }
     }
   }
  }
  
  if (ntracks<2) return;
  //
  // Find pairs
  //
  for (Int_t i=0;i<ntracks;++i) {
    AliESDtrack *track0 = event->GetTrack(i);     
    // track0 - choosen upper part
    if (!track0) continue;
    if (!track0->GetOuterParam()) continue;
    if (track0->GetOuterParam()->GetAlpha()<0) continue;
    Double_t d1[3];
    track0->GetDirection(d1);    
    for (Int_t j=0;j<ntracks;++j) {
      if (i==j) continue;
       AliESDtrack *track1 = event->GetTrack(j);   
       //track 1 lower part
       if (!track1) continue;
       if (!track1->GetOuterParam()) continue;
       if (track1->GetOuterParam()->GetAlpha()>0) continue;
       //
       Double_t d2[3];
       track1->GetDirection(d2);
      
       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
       if (! seed0) continue;
       if (! seed1) continue;
       Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
       Float_t d0  = track0->GetLinearD(0,0);
       Float_t d1  = track1->GetLinearD(0,0);
       //
       // conservative cuts - convergence to be guarantied
       // applying before track propagation
       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
       if (dir>fCutMinDir) continue;               // direction vector product
       Float_t bz = AliTracker::GetBz();
       Float_t dvertex0[2];   //distance to 0,0
       Float_t dvertex1[2];   //distance to 0,0 
       track0->GetDZ(0,0,0,bz,dvertex0);
       track1->GetDZ(0,0,0,bz,dvertex1);
       if (TMath::Abs(dvertex0[1])>250) continue;
       if (TMath::Abs(dvertex1[1])>250) continue;
       //
       //
       //
       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
       AliExternalTrackParam param0(*track0);
       AliExternalTrackParam param1(*track1);
       //
       // Propagate using Magnetic field and correct fo material budget
       //
       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
       //
       // Propagate rest to the 0,0 DCA - z should be ignored
       //
       //Bool_t b0 = ;
       param0.PropagateToDCA(&vtx,bz,1000);
       //Bool_t b1 = 
       param1.PropagateToDCA(&vtx,bz,1000);
       //      
       param0.GetDZ(0,0,0,bz,dvertex0);
       param1.GetDZ(0,0,0,bz,dvertex1);
       //
       Double_t xyz0[3];//,pxyz0[3];
       Double_t xyz1[3];//,pxyz1[3];
       param0.GetXYZ(xyz0);
       param1.GetXYZ(xyz1);
       Bool_t isPair = IsPair(&param0,&param1);

       Double_t z0 = track0->GetOuterParam()->GetZ();
       Double_t z1 = track1->GetOuterParam()->GetZ();

       Double_t z0inner = track0->GetInnerParam()->GetZ();
       Double_t z1inner = track1->GetInnerParam()->GetZ();

       if (isPair && z0>0 && z0inner>0 && z1<0 && z1inner<0) {
	 Double_t vecVdrift[2] = {time, param0.GetZ() - param1.GetZ()};
	 
	 if (track0->GetTPCNcls() > 80) {
	   fHistVdrift->Fill(vecVdrift);
	   fDz = param0.GetZ() - param1.GetZ();
	 }
       }
       if (isPair && z0 > 0 && z1 > 0) {
	 if (track1->GetTPCNcls()> 110 && track0->GetTPCNcls()> 110 && seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE) > 0) {
	   Float_t dedxratio = seed0->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)/seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE);
	   Double_t vecDeDxTgl[3] = {time, track0->GetTgl(), seed0->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)/seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)};
	   fHistDeDxTgl->Fill(vecDeDxTgl);
	   fdEdxRatio = dedxratio;
	   fTl = track0->GetTgl() ;
	 }
       }
       
    } // end 2nd order loop        
  } // end 1st order loop

}


void AliTPCcalibTime::Analyze() {
  //
  //
  //
  TH2D * hVdrift = GetHistVdrift()->Projection(1,0);
  hVdrift->Draw();
}


Long64_t AliTPCcalibTime::Merge(TCollection *li) {

  TIterator* iter = li->MakeIterator();
  AliTPCcalibTime* cal = 0;

  while ((cal = (AliTPCcalibTime*)iter->Next())) {
    if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
      Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
      return -1;
    }

    // add histograms here...
    fHistDeDxTgl->Add(cal->GetHistDeDxVsTgl());
    fHistVdrift->Add(cal->GetHistVdrift());
    fHistDeDx->Add(cal->GetHistDeDx());

  }
  
  return 0;
  
}


Bool_t  AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
  //
  //
  /*
  // 0. Same direction - OPOSITE  - cutDir +cutT    
  TCut cutDir("cutDir","dir<-0.99")
  // 1. 
  TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
  //
  // 2. The same rphi 
  TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
  //
  //
  //
  TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");  
  // 1/Pt diff cut
  */
  const Double_t *p0 = tr0->GetParameter();
  const Double_t *p1 = tr1->GetParameter();
  if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
  if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
  if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz)  return kFALSE;
  Double_t d0[3], d1[3];
  tr0->GetDirection(d0);    
  tr1->GetDirection(d1);       
  if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
  //
  return kTRUE;  
}
Commit	Line	Data
c74c5f6c	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/*
	17	Comments to be written here:
	18	1. What do we calibrate.
	19
	20	Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
	21
	22	AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
	23
	24
	25	2. How to interpret results
	26	3. Simple example
	27
	28	a.) determine the required time range:
	29
	30	AliXRDPROOFtoolkit tool;
	31	TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
	32	chain->Draw("GetTimeStamp()")
	33
	34	b) analyse calibration object on Proof in calibration train
	35
	36	AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift, IntegrationTimeDeDx);
	37
	38	c) plot results
	39
	40	TFile f("CalibObjects.root");
828383be	41	AliTPCcalibTime cal = (AliTPCcalibTime )f->Get("TPCCalib")->FindObject("cosmicTime");
c74c5f6c	42	cal->GetHistVdrift()->Projection(1,0)->Draw()
	43
	44	4. Analysis using debug streamers.
	45
da6c0bc9	46	gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
	47	gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
	48	AliXRDPROOFtoolkit tool;
	49	TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
	50	chainTime->Lookup();
	51
	52
c74c5f6c	53	*/
	54
	55
	56	#include "Riostream.h"
	57	#include "TChain.h"
	58	#include "TTree.h"
	59	#include "TH1F.h"
	60	#include "TH2F.h"
	61	#include "TH3F.h"
	62	#include "THnSparse.h"
	63	#include "TList.h"
	64	#include "TMath.h"
	65	#include "TCanvas.h"
	66	#include "TFile.h"
	67	#include "TF1.h"
	68	#include "TVectorD.h"
	69	#include "TProfile.h"
	70	#include "TGraphErrors.h"
	71	#include "TCanvas.h"
	72
	73	#include "AliTPCclusterMI.h"
	74	#include "AliTPCseed.h"
	75	#include "AliESDVertex.h"
	76	#include "AliESDEvent.h"
	77	#include "AliESDfriend.h"
	78	#include "AliESDInputHandler.h"
	79	#include "AliAnalysisManager.h"
	80
	81	#include "AliTracker.h"
	82	#include "AliMagFMaps.h"
	83	#include "AliTPCCalROC.h"
	84
	85	#include "AliLog.h"
	86
	87	#include "AliTPCcalibTime.h"
	88
	89	#include "TTreeStream.h"
	90	#include "AliTPCTracklet.h"
da6c0bc9	91	#include "TTimeStamp.h"
	92	#include "AliTPCcalibDB.h"
	93	#include "AliTPCcalibLaser.h"
c74c5f6c	94
	95	ClassImp(AliTPCcalibTime)
	96
	97
	98	AliTPCcalibTime::AliTPCcalibTime()
da6c0bc9	99	:AliTPCcalibBase(),
da6c0bc9	100	fTriggerMask(0),
c74c5f6c	101	fHistDeDxTgl(0),
	102	fHistDeDx(0),
	103	fHistVdrift(0),
	104	fIntegrationTimeDeDx(0),
da6c0bc9	105	fIntegrationTimeVdrift(0),
	106	fLaser(0), // pointer to laser calibration
	107	fDz(0), // current delta z
	108	fdEdx(0), // current dEdx
	109	fdEdxRatio(0), // current dEdx ratio
	110	fTl(0), // current tan(lambda)
c74c5f6c	111	fCutMaxD(5), // maximal distance in rfi ditection
da6c0bc9	112	fCutMaxDz(20), // maximal distance in rfi ditection
c74c5f6c	113	fCutTheta(0.03), // maximal distan theta
	114	fCutMinDir(-0.99) // direction vector products
	115
	116	{
	117	AliInfo("Default Constructor");
	118	}
	119
	120
	121	AliTPCcalibTime::AliTPCcalibTime(const Text_t name, const Text_t title, ULong64_t TriggerMask, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeDeDx, Int_t deltaIntegrationTimeVdrift)
	122	:AliTPCcalibBase(),
	123	fTriggerMask(0),
	124	fHistDeDxTgl(0),
	125	fHistDeDx(0),
	126	fHistVdrift(0),
	127	fIntegrationTimeDeDx(0),
da6c0bc9	128	fIntegrationTimeVdrift(0),
	129	fLaser(0), // pointer to laser calibration
	130	fDz(0), // current delta z
	131	fdEdx(0), // current dEdx
	132	fdEdxRatio(0), // current dEdx ratio
	133	fTl(0), // current tan(lambda)
c74c5f6c	134	fCutMaxD(5), // maximal distance in rfi ditection
da6c0bc9	135	fCutMaxDz(20), // maximal distance in rfi ditection
c74c5f6c	136	fCutTheta(0.03), // maximal distan theta
	137	fCutMinDir(-0.99) // direction vector products
	138	{
	139
	140	SetName(name);
	141	SetTitle(title);
	142
	143	AliInfo("Non Default Constructor");
	144
	145	fTriggerMask = TriggerMask;
	146
	147	fIntegrationTimeDeDx = deltaIntegrationTimeDeDx;
	148	fIntegrationTimeVdrift = deltaIntegrationTimeVdrift;
	149
	150	Double_t deltaTime = EndTime - StartTime;
	151
828383be	152	Int_t binsVdrift[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeVdrift), 100};
c74c5f6c	153	Double_t xminVdrift[2] = {StartTime, -20};
	154	Double_t xmaxVdrift[2] = {EndTime, 20};
	155	fHistVdrift = new THnSparseF("HistVdrift","vDrift; time;#Delta z",2,binsVdrift,xminVdrift,xmaxVdrift);
	156
828383be	157	Int_t binsDeDxTgl[3] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),30,100};
c74c5f6c	158	Double_t xminDeDxTgl[3] = {StartTime,-1,0.7};
	159	Double_t xmaxDeDxTgl[3] = {EndTime,1,1.3};
	160	fHistDeDxTgl = new THnSparseF("HistDeDxTgl","dEdx vs tgl;time;tgl;dEdxUp/dEdxLow",3,binsDeDxTgl,xminDeDxTgl,xmaxDeDxTgl);
	161
0ffd2ae1	162	Int_t binsDeDx[2] = {TMath::Nint(deltaTime/deltaIntegrationTimeDeDx),100};
c74c5f6c	163	Double_t xminDeDx[2] = {StartTime,1};
	164	Double_t xmaxDeDx[2] = {EndTime,100};
	165	fHistDeDx = new THnSparseF("HistDeDx","dEdx l;time;dEdx",2,binsDeDx,xminDeDx,xmaxDeDx);
	166
	167	}
	168
	169
	170
	171	AliTPCcalibTime::~AliTPCcalibTime(){
	172	//
	173	//
	174	//
	175	}
da6c0bc9	176	void AliTPCcalibTime::ResetCurrent(){
	177	//
	178	// reset current values
	179	//
	180	fDz=0; // current delta z
	181	fdEdx=0; // current dEdx
	182	fdEdxRatio=0; // current dEdx ratio
	183	fTl=0; // current tan(lambda)
	184
	185	}
	186
c74c5f6c	187
	188	void AliTPCcalibTime::Process(AliESDEvent *event) {
	189	//
	190	//
	191	//
da6c0bc9	192	Int_t ntracks=event->GetNumberOfTracks();
	193	if (ntracks<2) return;
	194	ResetCurrent();
	195	//
c74c5f6c	196	ProcessCosmic(event);
da6c0bc9	197	if (fTrigger==16){
	198	if (!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
	199	fLaser->Process(event);
	200	}
	201	//
	202	// fill debug streamer
	203	if (fStreamLevel>0 && fDz!=0){
	204	TTreeSRedirector *cstream = GetDebugStreamer();
	205	if (cstream){
	206	TTimeStamp tstamp(fTime);
	207	Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
	208	Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
	209	Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
	210	Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
	211	Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
	212	Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
	213	TVectorD vdriftA, vdriftC,vdriftAC;
	214	if (fLaser && fTrigger==16) {
	215	if (fLaser->fFitAside) vdriftA=*(fLaser->fFitAside);
	216	if (fLaser->fFitCside) vdriftC=*(fLaser->fFitCside);
	217	if (fLaser->fFitACside) vdriftAC=*(fLaser->fFitACside);
	218	}
	219	(*cstream)<<"timeInfo"<<
	220	"run="<<fRun<< // run number
	221	"event="<<fEvent<< // event number
	222	"time="<<fTime<< // time stamp of event
	223	"trigger="<<fTrigger<< // trigger
	224	"mag="<<fMagF<< // magnetic field
	225	// Environment values
	226	"press0="<<valuePressure0<<
	227	"press1="<<valuePressure1<<
	228	"pt0="<<ptrelative0<<
	229	"pt1="<<ptrelative1<<
	230	"temp0="<<temp0<<
	231	"temp1="<<temp1<<
	232	//
	233	// accumulated values
	234	//
	235	"fDz="<<fDz<< //! current delta z
	236	"fdEdx="<<fdEdx<< //! current dEdx
	237	"fdEdxRatio="<<fdEdxRatio<< //! current dEdx ratio
	238	"fTl="<<fTl<< //! current tan(lambda)
	239	//
	240	//laser
	241	//
	242	"laserA.="<<&vdriftA<<
	243	"laserC.="<<&vdriftC<<
	244	"laserAC.="<<&vdriftAC<<
	245	"\n";
	246	}
	247	}
c74c5f6c	248	}
	249
	250
	251
	252	void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event) {
	253
	254	if (!event) {
	255	Printf("ERROR: ESD not available");
	256	return;
	257	}
	258	if (event->GetTimeStamp() == 0 ) {
	259	Printf("no time stamp!");
	260	return;
	261	}
	262
	263	if (fTriggerMask != 0 && event->GetTriggerMask() != fTriggerMask) return;
	264
	265	UInt_t time = event->GetTimeStamp();
	266
	267	//
	268	// Find cosmic pairs
	269	//
	270	// Track0 is choosen in upper TPC part
	271	// Track1 is choosen in lower TPC part
	272	//
	273	Int_t ntracks=event->GetNumberOfTracks();
	274	if (ntracks==0) return;
	275	if (ntracks > 10) return; // temporary debug to remove LASER events
	276
	277
	278	if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
	279	AliESDfriend ESDfriend=static_cast<AliESDfriend>(event->FindListObject("AliESDfriend"));
	280
	281	TObjArray tpcSeeds(ntracks);
	282	Double_t vtxx[3]={0,0,0};
	283	Double_t svtxx[3]={0.000001,0.000001,100.};
	284	AliESDVertex vtx(vtxx,svtxx);
	285	//
	286	// track loop
	287	//
	288	for (Int_t i=0;i<ntracks;++i) {
	289	AliESDtrack *track = event->GetTrack(i);
	290
	291	const AliExternalTrackParam * trackIn = track->GetInnerParam();
	292	const AliExternalTrackParam * trackOut = track->GetOuterParam();
	293	if (!trackIn) continue;
	294	if (!trackOut) continue;
	295
	296	AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
	297	TObject *calibObject;
	298	AliTPCseed *seed = 0;
	299	for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
	300	if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
	301	}
	302	if (seed) {
	303	tpcSeeds.AddAt(seed,i);
	304	if (track->GetTPCNcls() > 50) {
	305	Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
	306	Double_t TPCsignal = seed->CookdEdxNorm(0.0,0.6,1,0,159,0x0,kTRUE,kTRUE);
	307	Double_t vecDeDx[2] = {time, TPCsignal};
da6c0bc9	308	if (meanP > 12) {
	309	fdEdx = TPCsignal;
	310	fHistDeDx->Fill(vecDeDx);
	311	}
c74c5f6c	312	}
	313	}
	314	}
	315
	316	if (ntracks<2) return;
	317	//
	318	// Find pairs
	319	//
	320	for (Int_t i=0;i<ntracks;++i) {
	321	AliESDtrack *track0 = event->GetTrack(i);
	322	// track0 - choosen upper part
	323	if (!track0) continue;
	324	if (!track0->GetOuterParam()) continue;
	325	if (track0->GetOuterParam()->GetAlpha()<0) continue;
	326	Double_t d1[3];
	327	track0->GetDirection(d1);
	328	for (Int_t j=0;j<ntracks;++j) {
	329	if (i==j) continue;
	330	AliESDtrack *track1 = event->GetTrack(j);
	331	//track 1 lower part
	332	if (!track1) continue;
	333	if (!track1->GetOuterParam()) continue;
	334	if (track1->GetOuterParam()->GetAlpha()>0) continue;
	335	//
	336	Double_t d2[3];
	337	track1->GetDirection(d2);
	338
	339	AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
	340	AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
	341	if (! seed0) continue;
	342	if (! seed1) continue;
	343	Float_t dir = (d1[0]d2[0] + d1[1]d2[1] + d1[2]*d2[2]);
	344	Float_t d0 = track0->GetLinearD(0,0);
	345	Float_t d1 = track1->GetLinearD(0,0);
	346	//
	347	// conservative cuts - convergence to be guarantied
	348	// applying before track propagation
	349	if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
	350	if (dir>fCutMinDir) continue; // direction vector product
	351	Float_t bz = AliTracker::GetBz();
	352	Float_t dvertex0[2]; //distance to 0,0
	353	Float_t dvertex1[2]; //distance to 0,0
	354	track0->GetDZ(0,0,0,bz,dvertex0);
	355	track1->GetDZ(0,0,0,bz,dvertex1);
	356	if (TMath::Abs(dvertex0[1])>250) continue;
	357	if (TMath::Abs(dvertex1[1])>250) continue;
	358	//
	359	//
	360	//
	361	Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
	362	AliExternalTrackParam param0(*track0);
	363	AliExternalTrackParam param1(*track1);
	364	//
	365	// Propagate using Magnetic field and correct fo material budget
	366	//
	367	AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
	368	AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
	369	//
	370	// Propagate rest to the 0,0 DCA - z should be ignored
	371	//
828383be	372	//Bool_t b0 = ;
	373	param0.PropagateToDCA(&vtx,bz,1000);
	374	//Bool_t b1 =
	375	param1.PropagateToDCA(&vtx,bz,1000);
c74c5f6c	376	//
	377	param0.GetDZ(0,0,0,bz,dvertex0);
	378	param1.GetDZ(0,0,0,bz,dvertex1);
	379	//
	380	Double_t xyz0[3];//,pxyz0[3];
	381	Double_t xyz1[3];//,pxyz1[3];
	382	param0.GetXYZ(xyz0);
	383	param1.GetXYZ(xyz1);
	384	Bool_t isPair = IsPair(&param0,&param1);
	385
	386	Double_t z0 = track0->GetOuterParam()->GetZ();
	387	Double_t z1 = track1->GetOuterParam()->GetZ();
	388
	389	Double_t z0inner = track0->GetInnerParam()->GetZ();
	390	Double_t z1inner = track1->GetInnerParam()->GetZ();
	391
	392	if (isPair && z0>0 && z0inner>0 && z1<0 && z1inner<0) {
	393	Double_t vecVdrift[2] = {time, param0.GetZ() - param1.GetZ()};
	394
da6c0bc9	395	if (track0->GetTPCNcls() > 80) {
	396	fHistVdrift->Fill(vecVdrift);
	397	fDz = param0.GetZ() - param1.GetZ();
	398	}
c74c5f6c	399	}
	400	if (isPair && z0 > 0 && z1 > 0) {
	401	if (track1->GetTPCNcls()> 110 && track0->GetTPCNcls()> 110 && seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE) > 0) {
da6c0bc9	402	Float_t dedxratio = seed0->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)/seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE);
c74c5f6c	403	Double_t vecDeDxTgl[3] = {time, track0->GetTgl(), seed0->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)/seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)};
c74c5f6c	404	fHistDeDxTgl->Fill(vecDeDxTgl);
da6c0bc9	405	fdEdxRatio = dedxratio;
da6c0bc9	406	fTl = track0->GetTgl() ;
c74c5f6c	407	}
	408	}
	409
	410	} // end 2nd order loop
	411	} // end 1st order loop
	412
	413	}
	414
	415
	416	void AliTPCcalibTime::Analyze() {
828383be	417	//
	418	//
	419	//
c74c5f6c	420	TH2D * hVdrift = GetHistVdrift()->Projection(1,0);
0ffd2ae1	421	hVdrift->Draw();
c74c5f6c	422	}
	423
	424
	425	Long64_t AliTPCcalibTime::Merge(TCollection *li) {
	426
	427	TIterator* iter = li->MakeIterator();
	428	AliTPCcalibTime* cal = 0;
	429
	430	while ((cal = (AliTPCcalibTime*)iter->Next())) {
	431	if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
	432	Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
	433	return -1;
	434	}
	435
	436	// add histograms here...
	437	fHistDeDxTgl->Add(cal->GetHistDeDxVsTgl());
	438	fHistVdrift->Add(cal->GetHistVdrift());
	439	fHistDeDx->Add(cal->GetHistDeDx());
	440
	441	}
	442
	443	return 0;
	444
	445	}
	446
	447
	448
	449	Bool_t AliTPCcalibTime::IsPair(AliExternalTrackParam tr0, AliExternalTrackParam tr1){
	450	//
	451	//
	452	/*
	453	// 0. Same direction - OPOSITE - cutDir +cutT
	454	TCut cutDir("cutDir","dir<-0.99")
	455	// 1.
	456	TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
	457	//
	458	// 2. The same rphi
	459	TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
	460	//
	461	//
	462	//
	463	TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
	464	// 1/Pt diff cut
	465	*/
	466	const Double_t *p0 = tr0->GetParameter();
	467	const Double_t *p1 = tr1->GetParameter();
	468	if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
	469	if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
da6c0bc9	470	if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
c74c5f6c	471	Double_t d0[3], d1[3];
	472	tr0->GetDirection(d0);
	473	tr1->GetDirection(d1);
	474	if (d0[0]d1[0] + d0[1]d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
	475	//
	476	return kTRUE;
	477	}