[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCReconstruction.cxx


#include <TDatabasePDG.h>
#include <TString.h>
#include <TSystem.h>
#include <TROOT.h>
#include <TFile.h>
#include <TPRegexp.h>

#include <AliExternalTrackParam.h>
#include <AliTPCcalibDB.h>
#include <AliTPCclusterMI.h>
#include <AliTPCSpaceCharge3D.h>
#include <AliTrackerBase.h>
#include <AliTrackPointArray.h>
#include <AliLog.h>
#include <AliTPCParam.h>
#include <AliTPCROC.h>
#include <TTreeStream.h>

#include "AliToyMCTrack.h"
#include "AliToyMCEvent.h"

#include "AliToyMCReconstruction.h"


//____________________________________________________________________________________
AliToyMCReconstruction::AliToyMCReconstruction() : TObject()
, fSeedingRow(140)
, fSeedingDist(10)
, fClusterType(0)
, fCorrectionType(kNoCorrection)
, fDoTrackFit(kTRUE)
, fUseMaterial(kFALSE)
, fIdealTracking(kFALSE)
, fTime0(-1)
, fCreateT0seed(kFALSE)
, fStreamer(0x0)
, fTree(0x0)
, fEvent(0x0)
, fTPCParam(0x0)
, fSpaceCharge(0x0)
{
  //
  //  ctor
  //
  fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
}

//____________________________________________________________________________________
AliToyMCReconstruction::~AliToyMCReconstruction()
{
  //
  //  dtor
  //

  delete fStreamer;
//   delete fTree;
}

//____________________________________________________________________________________
void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
{
  //
  //
  //

  TFile f(file);
  if (!f.IsOpen() || f.IsZombie()) {
    printf("ERROR: couldn't open the file '%s'\n", file);
    return;
  }
  
 fTree=(TTree*)f.Get("toyMCtree");
  if (!fTree) {
    printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
    return;
  }

  fEvent=0x0;
  fTree->SetBranchAddress("event",&fEvent);
  
  // read spacecharge from the Userinfo ot the tree
  InitSpaceCharge();
  
  TString debugName=file;
  debugName.ReplaceAll(".root","");
  debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
                        fUseMaterial,fIdealTracking,fClusterType,
                        Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
  debugName.Append(".debug.root");
  
  gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
  if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
  
  gROOT->cd();

  static AliExternalTrackParam dummySeedT0;
  static AliExternalTrackParam dummySeed;
  static AliExternalTrackParam dummyTrack;

  AliExternalTrackParam t0seed;
  AliExternalTrackParam seed;
  AliExternalTrackParam track;
  AliExternalTrackParam tOrig;

  AliExternalTrackParam *dummy;
  
  Int_t maxev=fTree->GetEntries();
  if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
  
  for (Int_t iev=0; iev<maxev; ++iev){
    printf("==============  Processing Event %6d =================\n",iev);
    fTree->GetEvent(iev);
    for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
      printf(" > ======  Processing Track %6d ========  \n",itr);
      const AliToyMCTrack *tr=fEvent->GetTrack(itr);
      tOrig = *tr;

      
      // set dummy 
      t0seed    = dummySeedT0;
      seed      = dummySeed;
      track     = dummyTrack;
      
      Float_t z0=fEvent->GetZ();
      Float_t t0=fEvent->GetT0();

      Float_t vdrift=GetVDrift();
      Float_t zLength=GetZLength(0);

      // crate time0 seed, steered by fCreateT0seed
      printf("t0 seed\n");
      fTime0=-1.;
      fCreateT0seed=kTRUE;
      dummy = GetSeedFromTrack(tr);
      
      if (dummy) {
        t0seed = *dummy;
        delete dummy;

        // crate real seed using the time 0 from the first seed
        // set fCreateT0seed now to false to get the seed in z coordinates
        fTime0 = t0seed.GetZ()-zLength/vdrift;
        fCreateT0seed = kFALSE;
        printf("seed (%.2g)\n",fTime0);
        dummy  = GetSeedFromTrack(tr);
        if (dummy) {
          seed = *dummy;
          delete dummy;

          // create fitted track
          if (fDoTrackFit){
            printf("track\n");
            dummy = GetFittedTrackFromSeed(tr, &seed);
            track = *dummy;
            delete dummy;
          }

          // propagate seed to 0
          const Double_t kMaxSnp = 0.85;
          const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
//           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
          
        }
      }

      Int_t ctype(fCorrectionType);
      
      if (fStreamer) {
        (*fStreamer) << "Tracks" <<
        "iev="         << iev             <<
        "z0="          << z0              <<
        "t0="          << t0              <<
        "fTime0="      << fTime0          <<
        "itr="         << itr             <<
        "clsType="     << fClusterType    <<
        "corrType="    << ctype           <<
        "seedRow="     << fSeedingRow     <<
        "seedDist="    << fSeedingDist    <<
        "vdrift="      << vdrift          <<
        "zLength="     << zLength         <<
        "t0seed.="     << &t0seed         <<
        "seed.="       << &seed           <<
        "track.="      << &track          <<
        "tOrig.="      << &tOrig          <<
        "\n";
      }
      
      
    }
  }

  delete fStreamer;
  fStreamer=0x0;

  delete fEvent;
  fEvent = 0x0;
  
  delete fTree;
  fTree=0x0;
  f.Close();
}

//____________________________________________________________________________________
AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr)
{
  //
  // if we don't have a valid time0 informaion (fTime0) available yet
  // assume we create a seed for the time0 estimate
  //

  // seed point informaion
  AliTrackPoint    seedPoint[3];
  const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
  
  // number of clusters to loop over
  const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
  
  UChar_t nextSeedRow=fSeedingRow;
  Int_t   nseeds=0;
  
  //assumes sorted clusters
  for (Int_t icl=0;icl<ncls;++icl) {
    const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
    if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
    if (!cl) continue;
    // use row in sector
    const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
    // skip clusters without proper pad row
    if (row>200) continue;
    
    //check seeding row
    // if we are in the last row and still miss a seed we use the last row
    //   even if the row spacing will not be equal
    if (row>=nextSeedRow || icl==ncls-1){
      seedCluster[nseeds]=cl;
      SetTrackPointFromCluster(cl, seedPoint[nseeds]);
      ++nseeds;
      nextSeedRow+=fSeedingDist;
      
      if (nseeds==3) break;
    }
  }
  
  // check we really have 3 seeds
  if (nseeds!=3) {
    AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
    return 0x0;
  }
  
  // do cluster correction for fCorrectionType:
  //   0 - no correction
  //   1 - TPC center
  //   2 - average eta
  //   3 - ideal
  // assign the cluster abs time as z component to all seeds
  for (Int_t iseed=0; iseed<3; ++iseed) {
    Float_t xyz[3]={0,0,0};
    seedPoint[iseed].GetXYZ(xyz);
    const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
    
    const Int_t sector=seedCluster[iseed]->GetDetector();
    const Int_t sign=1-2*((sector/18)%2);
    
    if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
      printf("correction type: %d\n",(Int_t)fCorrectionType);

      // the settings below are for the T0 seed
      // for known T0 the z position is already calculated in SetTrackPointFromCluster
      if ( fCreateT0seed ){
        if ( fCorrectionType == kTPCCenter  ) xyz[2] = 125.*sign;
        //!!! TODO: is this the correct association?
        if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
      }
      
      if ( fCorrectionType == kIdeal      ) xyz[2] = seedCluster[iseed]->GetZ();
      
      //!!! TODO: to be replaced with the proper correction
      fSpaceCharge->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
    }

    // after the correction set the time bin as z-Position in case of a T0 seed
    if ( fCreateT0seed )
      xyz[2]=seedCluster[iseed]->GetTimeBin();
    
    seedPoint[iseed].SetXYZ(xyz);
  }
  
  const Double_t kMaxSnp = 0.85;
  const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
  
  AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
  seed->ResetCovariance(10);

  if (fCreateT0seed){
    // if fTime0 < 0 we assume that we create a seed for the T0 estimate
    AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
    if (TMath::Abs(seed->GetX())>3) {
      printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
    }
  }
  
  return seed;
  
}

//____________________________________________________________________________________
void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
{
  //
  // make AliTrackPoint out of AliTPCclusterMI
  //
  
  if (!cl) return;
    Float_t xyz[3]={0.,0.,0.};
  //   ClusterToSpacePoint(cl,xyz);
  //   cl->GetGlobalCov(cov);
  //TODO: what to do with the covariance matrix???
  //TODO: the problem is that it is used in GetAngle in AliTrackPoint
  //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
  //TODO: for the moment simply assign 1 permill squared
  // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
  //   Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
  //                   xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
  //   cl->GetGlobalXYZ(xyz);
  //   cl->GetGlobalCov(cov);
  // voluem ID to add later ....
  //   p.SetXYZ(xyz);
  //   p.SetCov(cov);
  AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint();
  p=*tp;
  delete tp;
  //   cl->Print();
  //   p.Print();
  p.SetVolumeID(cl->GetDetector());
  
  
  if ( !fCreateT0seed && !fIdealTracking ) {
    p.GetXYZ(xyz);
    const Int_t sector=cl->GetDetector();
    const Int_t sign=1-2*((sector/18)%2);
    const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
    printf(" z:  %.2f  %.2f\n",xyz[2],zT0);
    xyz[2]=zT0;
    p.SetXYZ(xyz);
  }
  
  
  //   p.Rotate(p.GetAngle()).Print();
}

//____________________________________________________________________________________
void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
{
  //
  // convert the cluster to a space point in global coordinates
  //
  if (!cl) return;
  xyz[0]=cl->GetRow();
  xyz[1]=cl->GetPad();
  xyz[2]=cl->GetTimeBin(); // this will not be correct at all
  Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
  //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
  fTPCParam->Transform8to4(xyz,i);
  //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
  fTPCParam->Transform4to3(xyz,i);
  //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
  fTPCParam->Transform2to1(xyz,i);
  //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
}

//____________________________________________________________________________________
AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
{
  //
  //
  //

  // create track
  AliExternalTrackParam *track = new AliExternalTrackParam(*seed);

  Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();

  const AliTPCROC * roc = AliTPCROC::Instance();
  
  const Double_t kRTPC0  = roc->GetPadRowRadii(0,0);
  const Double_t kRTPC1  = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
  const Double_t kMaxSnp = 0.85;
  const Double_t kMaxR   = 500.;
  const Double_t kMaxZ   = 500.;
  
  //   const Double_t kMaxZ0=220;
//   const Double_t kZcut=3;
  
  const Double_t refX = tr->GetX();
  
  const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
  
  // loop over all other points and add to the track
  for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
    AliTrackPoint pIn;
    const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
    if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
    SetTrackPointFromCluster(cl, pIn);
    if (fCorrectionType != kNoCorrection){
      Float_t xyz[3]={0,0,0};
      pIn.GetXYZ(xyz);
//       if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
      fSpaceCharge->CorrectPoint(xyz, cl->GetDetector());
      pIn.SetXYZ(xyz);
    }
    // rotate the cluster to the local detector frame
    track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
    AliTrackPoint prot = pIn.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
    if (TMath::Abs(prot.GetX())<kRTPC0) continue;
    if (TMath::Abs(prot.GetX())>kRTPC1) continue;
    //
    Bool_t res=kTRUE;
    if (fUseMaterial) res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp);
    else res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,kFALSE);

    if (!res) break;
    
    if (TMath::Abs(track->GetZ())>kMaxZ) break;
    if (TMath::Abs(track->GetX())>kMaxR) break;
//     if (TMath::Abs(track->GetZ())<kZcut)continue;
    //
    Double_t pointPos[2]={0,0};
    Double_t pointCov[3]={0,0,0};
    pointPos[0]=prot.GetY();//local y
    pointPos[1]=prot.GetZ();//local z
    pointCov[0]=prot.GetCov()[3];//simay^2
    pointCov[1]=prot.GetCov()[4];//sigmayz
    pointCov[2]=prot.GetCov()[5];//sigmaz^2
    
    if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
  }

  if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp);
  else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,kFALSE);

  // rotate fittet track to the frame of the original track and propagate to same reference
  track->Rotate(tr->GetAlpha());
  
  if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp);
  else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
  
  return track;
}

//____________________________________________________________________________________
void AliToyMCReconstruction::InitSpaceCharge()
{
  //
  // Init the space charge map
  //

  TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
  if (fTree) {
    TList *l=fTree->GetUserInfo();
    for (Int_t i=0; i<l->GetEntries(); ++i) {
      TString s(l->At(i)->GetName());
      if (s.Contains("SC_")) {
        filename=s;
        break;
      }
    }
  }

  printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
  TFile f(filename.Data());
  fSpaceCharge=(AliTPCSpaceCharge3D*)f.Get("map");
  
  //   fSpaceCharge = new AliTPCSpaceCharge3D();
  //   fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
  //   fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
  // //   fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
  //   fSpaceCharge->InitSpaceCharge3DDistortion();
  
}

//____________________________________________________________________________________
Double_t AliToyMCReconstruction::GetVDrift() const
{
  //
  //
  //
  return fTPCParam->GetDriftV();
}

//____________________________________________________________________________________
Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
{
  //
  //
  //
  if (roc<0 || roc>71) return -1;
  return fTPCParam->GetZLength(roc);
}

//____________________________________________________________________________________
TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
  TString s=gSystem->GetFromPipe(Form("ls %s",files));

  TTree *tFirst=0x0;
  TObjArray *arrFiles=s.Tokenize("\n");
  
  for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
    TString name(arrFiles->At(ifile)->GetName());
    
    TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
    TObjArray *arrMatch=0x0;
    arrMatch=reg.MatchS(name);
    
    if (!tFirst) {
      TFile *f=TFile::Open(name.Data());
      if (!f) continue;
      TTree *t=(TTree*)f->Get("Tracks");
      if (!t) {
        delete f;
        continue;
      }
      
      t->SetName(arrMatch->At(1)->GetName());
      tFirst=t;
    } else {
      tFirst->AddFriend(Form("t%s=Tracks",arrMatch->At(1)->GetName()), name.Data());
//       tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
    }
  }

  tFirst->GetListOfFriends()->Print();
  return tFirst;
}
Commit	Line	Data
d1cf83f5	1
	2	#include <TDatabasePDG.h>
	3	#include <TString.h>
4a777885	4	#include <TSystem.h>
	5	#include <TROOT.h>
	6	#include <TFile.h>
9e98dea8	7	#include <TPRegexp.h>
d1cf83f5	8
	9	#include <AliExternalTrackParam.h>
	10	#include <AliTPCcalibDB.h>
	11	#include <AliTPCclusterMI.h>
	12	#include <AliTPCSpaceCharge3D.h>
	13	#include <AliTrackerBase.h>
	14	#include <AliTrackPointArray.h>
	15	#include <AliLog.h>
	16	#include <AliTPCParam.h>
	17	#include <AliTPCROC.h>
	18	#include <TTreeStream.h>
	19
	20	#include "AliToyMCTrack.h"
4a777885	21	#include "AliToyMCEvent.h"
d1cf83f5	22
	23	#include "AliToyMCReconstruction.h"
	24
	25
	26	//____________________________________________________________________________________
	27	AliToyMCReconstruction::AliToyMCReconstruction() : TObject()
	28	, fSeedingRow(140)
	29	, fSeedingDist(10)
	30	, fClusterType(0)
	31	, fCorrectionType(kNoCorrection)
	32	, fDoTrackFit(kTRUE)
	33	, fUseMaterial(kFALSE)
600eaa0d	34	, fIdealTracking(kFALSE)
d1cf83f5	35	, fTime0(-1)
600eaa0d	36	, fCreateT0seed(kFALSE)
d1cf83f5	37	, fStreamer(0x0)
d1cf83f5	38	, fTree(0x0)
4a777885	39	, fEvent(0x0)
d1cf83f5	40	, fTPCParam(0x0)
	41	, fSpaceCharge(0x0)
	42	{
	43	//
	44	// ctor
	45	//
	46	fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
	47	}
	48
	49	//____________________________________________________________________________________
	50	AliToyMCReconstruction::~AliToyMCReconstruction()
	51	{
	52	//
	53	// dtor
	54	//
	55
	56	delete fStreamer;
4a777885	57	// delete fTree;
	58	}
	59
	60	//____________________________________________________________________________________
	61	void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
	62	{
	63	//
	64	//
	65	//
	66
	67	TFile f(file);
	68	if (!f.IsOpen() \|\| f.IsZombie()) {
	69	printf("ERROR: couldn't open the file '%s'\n", file);
	70	return;
	71	}
	72
	73	fTree=(TTree*)f.Get("toyMCtree");
	74	if (!fTree) {
	75	printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
	76	return;
	77	}
	78
	79	fEvent=0x0;
	80	fTree->SetBranchAddress("event",&fEvent);
	81
	82	// read spacecharge from the Userinfo ot the tree
	83	InitSpaceCharge();
	84
	85	TString debugName=file;
	86	debugName.ReplaceAll(".root","");
9e98dea8	87	debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
600eaa0d	88	fUseMaterial,fIdealTracking,fClusterType,
600eaa0d	89	Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
4a777885	90	debugName.Append(".debug.root");
	91
	92	gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
	93	if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
	94
	95	gROOT->cd();
	96
	97	static AliExternalTrackParam dummySeedT0;
	98	static AliExternalTrackParam dummySeed;
	99	static AliExternalTrackParam dummyTrack;
	100
	101	AliExternalTrackParam t0seed;
	102	AliExternalTrackParam seed;
	103	AliExternalTrackParam track;
	104	AliExternalTrackParam tOrig;
	105
	106	AliExternalTrackParam *dummy;
	107
	108	Int_t maxev=fTree->GetEntries();
	109	if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
	110
	111	for (Int_t iev=0; iev<maxev; ++iev){
	112	printf("============== Processing Event %6d =================\n",iev);
	113	fTree->GetEvent(iev);
	114	for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
	115	printf(" > ====== Processing Track %6d ======== \n",itr);
	116	const AliToyMCTrack *tr=fEvent->GetTrack(itr);
	117	tOrig = *tr;
	118
	119
	120	// set dummy
	121	t0seed = dummySeedT0;
	122	seed = dummySeed;
	123	track = dummyTrack;
	124
	125	Float_t z0=fEvent->GetZ();
	126	Float_t t0=fEvent->GetT0();
	127
	128	Float_t vdrift=GetVDrift();
	129	Float_t zLength=GetZLength(0);
	130
600eaa0d	131	// crate time0 seed, steered by fCreateT0seed
4a777885	132	printf("t0 seed\n");
4a777885	133	fTime0=-1.;
600eaa0d	134	fCreateT0seed=kTRUE;
4a777885	135	dummy = GetSeedFromTrack(tr);
	136
	137	if (dummy) {
	138	t0seed = *dummy;
	139	delete dummy;
	140
	141	// crate real seed using the time 0 from the first seed
600eaa0d	142	// set fCreateT0seed now to false to get the seed in z coordinates
4a777885	143	fTime0 = t0seed.GetZ()-zLength/vdrift;
600eaa0d	144	fCreateT0seed = kFALSE;
600eaa0d	145	printf("seed (%.2g)\n",fTime0);
4a777885	146	dummy = GetSeedFromTrack(tr);
	147	if (dummy) {
	148	seed = *dummy;
	149	delete dummy;
	150
	151	// create fitted track
	152	if (fDoTrackFit){
	153	printf("track\n");
	154	dummy = GetFittedTrackFromSeed(tr, &seed);
	155	track = *dummy;
	156	delete dummy;
	157	}
600eaa0d	158
	159	// propagate seed to 0
	160	const Double_t kMaxSnp = 0.85;
	161	const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
	162	// AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
	163
4a777885	164	}
	165	}
	166
	167	Int_t ctype(fCorrectionType);
	168
	169	if (fStreamer) {
	170	(*fStreamer) << "Tracks" <<
	171	"iev=" << iev <<
	172	"z0=" << z0 <<
	173	"t0=" << t0 <<
600eaa0d	174	"fTime0=" << fTime0 <<
4a777885	175	"itr=" << itr <<
	176	"clsType=" << fClusterType <<
	177	"corrType=" << ctype <<
	178	"seedRow=" << fSeedingRow <<
	179	"seedDist=" << fSeedingDist <<
	180	"vdrift=" << vdrift <<
	181	"zLength=" << zLength <<
	182	"t0seed.=" << &t0seed <<
	183	"seed.=" << &seed <<
	184	"track.=" << &track <<
	185	"tOrig.=" << &tOrig <<
	186	"\n";
	187	}
	188
	189
	190	}
	191	}
	192
	193	delete fStreamer;
	194	fStreamer=0x0;
600eaa0d	195
	196	delete fEvent;
	197	fEvent = 0x0;
4a777885	198
600eaa0d	199	delete fTree;
	200	fTree=0x0;
	201	f.Close();
d1cf83f5	202	}
	203
	204	//____________________________________________________________________________________
	205	AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr)
	206	{
	207	//
	208	// if we don't have a valid time0 informaion (fTime0) available yet
	209	// assume we create a seed for the time0 estimate
	210	//
	211
	212	// seed point informaion
	213	AliTrackPoint seedPoint[3];
	214	const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
	215
	216	// number of clusters to loop over
	217	const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
	218
	219	UChar_t nextSeedRow=fSeedingRow;
	220	Int_t nseeds=0;
	221
	222	//assumes sorted clusters
	223	for (Int_t icl=0;icl<ncls;++icl) {
	224	const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
	225	if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
	226	if (!cl) continue;
	227	// use row in sector
	228	const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
	229	// skip clusters without proper pad row
	230	if (row>200) continue;
	231
	232	//check seeding row
	233	// if we are in the last row and still miss a seed we use the last row
	234	// even if the row spacing will not be equal
	235	if (row>=nextSeedRow \|\| icl==ncls-1){
	236	seedCluster[nseeds]=cl;
	237	SetTrackPointFromCluster(cl, seedPoint[nseeds]);
	238	++nseeds;
	239	nextSeedRow+=fSeedingDist;
	240
	241	if (nseeds==3) break;
	242	}
	243	}
	244
	245	// check we really have 3 seeds
	246	if (nseeds!=3) {
	247	AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
	248	return 0x0;
	249	}
	250
	251	// do cluster correction for fCorrectionType:
	252	// 0 - no correction
	253	// 1 - TPC center
	254	// 2 - average eta
	255	// 3 - ideal
	256	// assign the cluster abs time as z component to all seeds
	257	for (Int_t iseed=0; iseed<3; ++iseed) {
	258	Float_t xyz[3]={0,0,0};
	259	seedPoint[iseed].GetXYZ(xyz);
600eaa0d	260	const Float_t r=TMath::Sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
d1cf83f5	261
600eaa0d	262	const Int_t sector=seedCluster[iseed]->GetDetector();
600eaa0d	263	const Int_t sign=1-2*((sector/18)%2);
d1cf83f5	264
d1cf83f5	265	if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
4a777885	266	printf("correction type: %d\n",(Int_t)fCorrectionType);
600eaa0d	267
	268	// the settings below are for the T0 seed
	269	// for known T0 the z position is already calculated in SetTrackPointFromCluster
	270	if ( fCreateT0seed ){
	271	if ( fCorrectionType == kTPCCenter ) xyz[2] = 125.*sign;
	272	//!!! TODO: is this the correct association?
	273	if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.TMath::DegToRad())r*sign;
	274	}
4a777885	275
d1cf83f5	276	if ( fCorrectionType == kIdeal ) xyz[2] = seedCluster[iseed]->GetZ();
600eaa0d	277
d1cf83f5	278	//!!! TODO: to be replaced with the proper correction
	279	fSpaceCharge->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
	280	}
	281
600eaa0d	282	// after the correction set the time bin as z-Position in case of a T0 seed
	283	if ( fCreateT0seed )
	284	xyz[2]=seedCluster[iseed]->GetTimeBin();
d1cf83f5	285
	286	seedPoint[iseed].SetXYZ(xyz);
	287	}
	288
	289	const Double_t kMaxSnp = 0.85;
	290	const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
	291
	292	AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
	293	seed->ResetCovariance(10);
	294
600eaa0d	295	if (fCreateT0seed){
d1cf83f5	296	// if fTime0 < 0 we assume that we create a seed for the T0 estimate
	297	AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
	298	if (TMath::Abs(seed->GetX())>3) {
	299	printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
	300	}
	301	}
	302
	303	return seed;
	304
	305	}
	306
	307	//____________________________________________________________________________________
	308	void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
	309	{
	310	//
	311	// make AliTrackPoint out of AliTPCclusterMI
	312	//
	313
	314	if (!cl) return;
4a777885	315	Float_t xyz[3]={0.,0.,0.};
d1cf83f5	316	// ClusterToSpacePoint(cl,xyz);
	317	// cl->GetGlobalCov(cov);
	318	//TODO: what to do with the covariance matrix???
	319	//TODO: the problem is that it is used in GetAngle in AliTrackPoint
	320	//TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
	321	//TODO: for the moment simply assign 1 permill squared
	322	// in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
	323	// Float_t cov[6]={xyz[0]xyz[0]1e-6,xyz[0]xyz[1]1e-6,xyz[0]xyz[2]1e-6,
	324	// xyz[1]xyz[1]1e-6,xyz[1]xyz[2]1e-6,xyz[2]xyz[2]1e-6};
	325	// cl->GetGlobalXYZ(xyz);
	326	// cl->GetGlobalCov(cov);
	327	// voluem ID to add later ....
	328	// p.SetXYZ(xyz);
	329	// p.SetCov(cov);
	330	AliTrackPoint tp=const_cast<AliTPCclusterMI>(cl)->MakePoint();
	331	p=*tp;
	332	delete tp;
	333	// cl->Print();
	334	// p.Print();
	335	p.SetVolumeID(cl->GetDetector());
600eaa0d	336
	337
	338	if ( !fCreateT0seed && !fIdealTracking ) {
	339	p.GetXYZ(xyz);
	340	const Int_t sector=cl->GetDetector();
	341	const Int_t sign=1-2*((sector/18)%2);
	342	const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)GetVDrift() )sign;
	343	printf(" z: %.2f %.2f\n",xyz[2],zT0);
	344	xyz[2]=zT0;
	345	p.SetXYZ(xyz);
4a777885	346	}
600eaa0d	347
600eaa0d	348
d1cf83f5	349	// p.Rotate(p.GetAngle()).Print();
	350	}
	351
	352	//____________________________________________________________________________________
	353	void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
	354	{
	355	//
	356	// convert the cluster to a space point in global coordinates
	357	//
	358	if (!cl) return;
	359	xyz[0]=cl->GetRow();
	360	xyz[1]=cl->GetPad();
	361	xyz[2]=cl->GetTimeBin(); // this will not be correct at all
	362	Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
	363	// printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
	364	fTPCParam->Transform8to4(xyz,i);
	365	// printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
	366	fTPCParam->Transform4to3(xyz,i);
	367	// printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
	368	fTPCParam->Transform2to1(xyz,i);
	369	// printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
	370	}
	371
	372	//____________________________________________________________________________________
	373	AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack tr, const AliExternalTrackParam seed)
	374	{
	375	//
	376	//
	377	//
	378
	379	// create track
	380	AliExternalTrackParam track = new AliExternalTrackParam(seed);
	381
	382	Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
	383
	384	const AliTPCROC * roc = AliTPCROC::Instance();
	385
	386	const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
	387	const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
	388	const Double_t kMaxSnp = 0.85;
	389	const Double_t kMaxR = 500.;
	390	const Double_t kMaxZ = 500.;
	391
	392	// const Double_t kMaxZ0=220;
	393	// const Double_t kZcut=3;
	394
	395	const Double_t refX = tr->GetX();
	396
	397	const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
	398
	399	// loop over all other points and add to the track
4a777885	400	for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
d1cf83f5	401	AliTrackPoint pIn;
	402	const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
	403	if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
	404	SetTrackPointFromCluster(cl, pIn);
	405	if (fCorrectionType != kNoCorrection){
	406	Float_t xyz[3]={0,0,0};
	407	pIn.GetXYZ(xyz);
600eaa0d	408	// if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
d1cf83f5	409	fSpaceCharge->CorrectPoint(xyz, cl->GetDetector());
	410	pIn.SetXYZ(xyz);
	411	}
	412	// rotate the cluster to the local detector frame
	413	track->Rotate(((cl->GetDetector()%18)20+10)TMath::DegToRad());
	414	AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
	415	if (TMath::Abs(prot.GetX())<kRTPC0) continue;
	416	if (TMath::Abs(prot.GetX())>kRTPC1) continue;
	417	//
	418	Bool_t res=kTRUE;
	419	if (fUseMaterial) res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp);
	420	else res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
	421
	422	if (!res) break;
	423
	424	if (TMath::Abs(track->GetZ())>kMaxZ) break;
	425	if (TMath::Abs(track->GetX())>kMaxR) break;
	426	// if (TMath::Abs(track->GetZ())<kZcut)continue;
	427	//
	428	Double_t pointPos[2]={0,0};
	429	Double_t pointCov[3]={0,0,0};
	430	pointPos[0]=prot.GetY();//local y
	431	pointPos[1]=prot.GetZ();//local z
	432	pointCov[0]=prot.GetCov()[3];//simay^2
	433	pointCov[1]=prot.GetCov()[4];//sigmayz
	434	pointCov[2]=prot.GetCov()[5];//sigmaz^2
	435
	436	if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
	437	}
	438
	439	if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp);
	440	else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
	441
	442	// rotate fittet track to the frame of the original track and propagate to same reference
	443	track->Rotate(tr->GetAlpha());
	444
	445	if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp);
	446	else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
	447
	448	return track;
	449	}
	450
	451	//____________________________________________________________________________________
	452	void AliToyMCReconstruction::InitSpaceCharge()
	453	{
	454	//
	455	// Init the space charge map
	456	//
	457
	458	TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
	459	if (fTree) {
	460	TList *l=fTree->GetUserInfo();
	461	for (Int_t i=0; i<l->GetEntries(); ++i) {
	462	TString s(l->At(i)->GetName());
	463	if (s.Contains("SC_")) {
	464	filename=s;
	465	break;
	466	}
	467	}
	468	}
	469
	470	printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
	471	TFile f(filename.Data());
	472	fSpaceCharge=(AliTPCSpaceCharge3D*)f.Get("map");
473
474	// fSpaceCharge = new AliTPCSpaceCharge3D();
475	// fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
476	// fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
477	// // fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
478	// fSpaceCharge->InitSpaceCharge3DDistortion();
479
480	}
481
482	//____________________________________________________________________________________
483	Double_t AliToyMCReconstruction::GetVDrift() const
484	{
485	//
486	//
487	//
488	return fTPCParam->GetDriftV();
489	}
490
491	//____________________________________________________________________________________
492	Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
493	{
494	//
495	//
496	//
497	if (roc<0 \|\| roc>71) return -1;
498	return fTPCParam->GetZLength(roc);
499	}
500
9e98dea8	501	//____________________________________________________________________________________
	502	TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
	503	TString s=gSystem->GetFromPipe(Form("ls %s",files));
	504
	505	TTree *tFirst=0x0;
	506	TObjArray *arrFiles=s.Tokenize("\n");
	507
	508	for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
	509	TString name(arrFiles->At(ifile)->GetName());
	510
	511	TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
	512	TObjArray *arrMatch=0x0;
	513	arrMatch=reg.MatchS(name);
	514
	515	if (!tFirst) {
	516	TFile *f=TFile::Open(name.Data());
	517	if (!f) continue;
	518	TTree t=(TTree)f->Get("Tracks");
	519	if (!t) {
	520	delete f;
	521	continue;
	522	}
	523
	524	t->SetName(arrMatch->At(1)->GetName());
	525	tFirst=t;
	526	} else {
	527	tFirst->AddFriend(Form("t%s=Tracks",arrMatch->At(1)->GetName()), name.Data());
	528	// tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
	529	}
	530	}
	531
376089b6	532	tFirst->GetListOfFriends()->Print();
9e98dea8	533	return tFirst;
9e98dea8	534	}
d1cf83f5	535