[u/mrichter/AliRoot.git] / STEER / AliTPCComparisonMI.C

/**************************************************************************
 * 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.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  Time Projection Chamber                                                  //
//  Comparison macro for TPC                                                 //
//  responsible: 
//  marian.ivanov@cern.ch                                                    //
//
//
/*
//


.L $ALICE_ROOT/STEER/AliGenInfo.C+
.L $ALICE_ROOT/STEER/AliTPCComparisonMI.C+


TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",-1,1,0);
t2->Exec();


TCut cprim("cprim","MC.fVDist[3]<1");
TCut csec("cprim","MC.fVDist[3]>1");
TCut crec("crec","fReconstructed==1");
TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
AliTPCComparisonDraw comp; 
comp.SetIO();


//example
comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&fPdg==-211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
comp.fRes->Draw();
comp.fMean->Draw();  
comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1,"fReconstructed>0",10,0,1.5)
comp.fRes->Draw();
*/


#if !defined(__CINT__) || defined(__MAKECINT__)
#include <stdio.h>
#include <string.h>
//ROOT includes
#include "Rtypes.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TCut.h"
#include "TString.h"
#include "TBenchmark.h"
#include "TStopwatch.h"
#include "TParticle.h"
#include "TSystem.h"
#include "TTimer.h"
#include "TVector3.h"
#include "TPad.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TF1.h"
#include "TText.h"
#include "Getline.h"
#include "TStyle.h"

//ALIROOT includes
#include "AliRun.h"
#include "AliStack.h"
#include "AliTPCtrack.h"
#include "AliSimDigits.h"
#include "AliTPCParam.h"
#include "AliTPC.h"
#include "AliTPCLoader.h"
#include "AliDetector.h"
#include "AliTrackReference.h"
#include "AliRun.h"
#include "AliTPCParamSR.h"
#include "AliTracker.h"
#include "AliComplexCluster.h"
#include "AliMagF.h"
#endif
#include "AliGenInfo.h"
#include "AliTPCComparisonMI.h"

//
//
  

void AliTPCRecInfo::Update(AliMCInfo* info,AliTPCParam * param, Bool_t reconstructed, Int_t direction)
{
  //
  //
  //calculates derived variables
  //
  
  if (reconstructed==kFALSE){
    fReconstructed = kFALSE;
    return;
  }
  fReconstructed = kTRUE;
  // Find the nearest track reference
  //
  Double_t radius = TMath::Sqrt(fTrack.GetX()*fTrack.GetX()+fTrack.GetY()*fTrack.GetY());
  TClonesArray * tpcreferences = info->fTPCReferences;
  Int_t nentries = tpcreferences->GetEntriesFast();
  AliTrackReference * ref = 0;
  Double_t dr = 1000000;
  for (Int_t i=0;i<nentries;i++){
    AliTrackReference * refi = (AliTrackReference*)tpcreferences->UncheckedAt(i);
    if (!refi) continue;
    //if  ( (direction<0) && (refi->R()<radius-10)) continue;
    //if  ( (direction>0) && (refi->R()>radius+10)) continue;    
    if  (TMath::Abs(refi->R()-radius)<dr){
      dr = TMath::Abs(refi->R()-radius);
      ref = refi;
    }
  }  
  if (ref==0) {
    fReconstructed = kFALSE;
    return;
  }
  //
  //
  fdEdx = fTrack.GetdEdx();
  //
  // propagate track to nearest reference in direction
  TVector3 local = TPCCmpTr::TR2Local(ref,param);
  local.GetXYZ(fTRLocalCoord);
  //
  // rotate if neccessary
  Float_t ymax = local.X()*TMath::Tan(0.5*param->GetInnerAngle());  
  Float_t y = fTrack.GetYat(local.X());
  if (y > ymax) {
    fTrack.Rotate(param->GetInnerAngle());
  } else if (y <-ymax) {
    fTrack.Rotate(-param->GetInnerAngle());
  }   
  Double_t xtrack = fTrack.GetX();
  Double_t delta = (local.X()-xtrack)*0.1;
  for (Int_t i=0;i<9;i++){
    fTrack.PropagateTo(xtrack+delta*float(i));
  }
  fTrack.PropagateTo(local.X());

  Double_t par[5], cov[15], localX;
  fTrack.GetExternalParameters(localX,par);
  fTrack.GetExternalCovariance(cov);
  //
  //
  fRecPhi=TMath::ASin(par[2]) + fTrack.GetAlpha();
  Double_t phi2 =TMath::ATan2(ref->Py(),ref->Px());
  if (phi2<0) phi2+=2*TMath::Pi();
  fGenPhi =phi2;
  //
  if (fRecPhi<0) fRecPhi+=2*TMath::Pi();
  if (fRecPhi>=2*TMath::Pi()) fRecPhi-=2*TMath::Pi();
  fLambda = TMath::ATan(par[3]);
  fRecPt_1 = TMath::Abs(par[4]);  
  //TPC
  //
  //
  Double_t phi=TMath::ATan2(par[0],localX) + fTrack.GetAlpha();
  Double_t r=TMath::Sqrt(localX*localX + par[0]*par[0]);
  fTPCinR1[0]= r*TMath::Cos(phi); 
  fTPCinR1[1]= r*TMath::Sin(phi); 
  fTPCinR1[2]= par[1];
  fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
  fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);	
  //
  fTPCinP1[0] = fTrack.Px();
  fTPCinP1[1] = fTrack.Py();
  fTPCinP1[2] = fTrack.Pz();
  fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
  //
  fTPCinR0[0] = ref->X();	
  fTPCinR0[1] = ref->Y();	
  fTPCinR0[2] = ref->Z();
  fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
  fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
  //
  fTPCinP0[0]=ref->Px();
  fTPCinP0[1]=ref->Py();
  fTPCinP0[2]=ref->Pz();
  fTPCinP0[3]=ref->Pt();
  //
  //
  if (fTPCinP1[3]>0.0000001){
    fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
    fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);  
    //
    fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
    fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
  }    
  //
  fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3];  //delta rfi
  fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
  fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]);              //delta z
  fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
  //
  fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
  fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
  //
  fTPCDelta[3] = TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]);
  fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
  fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
  Double_t sign = (fTrack.GetC()>0) ? 1:-1;
  fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]);      
  
}


////////////////////////////////////////////////////////////////////////
TPCCmpTr::TPCCmpTr()
{
  Reset();
}

////////////////////////////////////////////////////////////////////////
TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
		   const char* fnCmp,
		   const char* fnGalice, Int_t direction,
		   Int_t nEvents, Int_t firstEvent)
{
  Reset();
  //  fFnGenTracks = fnGenTracks;
  //  fFnCmp = fnCmp;
  sprintf(fFnGenTracks,"%s",fnGenTracks);
  sprintf(fFnCmp,"%s",fnCmp);

  fFirstEventNr = firstEvent;
  fEventNr = firstEvent;
  fNEvents = nEvents;
  fDirection = direction;
  //
  fLoader = AliRunLoader::Open(fnGalice);
  if (gAlice){
    //delete gAlice->GetRunLoader();
    delete gAlice;
    gAlice = 0x0;
  }
  if (fLoader->LoadgAlice()){
    cerr<<"Error occured while l"<<endl;
  }
  Int_t nall = fLoader->GetNumberOfEvents();
  if (nEvents==0) {
    nEvents =nall;
    fNEvents=nall;
    fFirstEventNr=0;
  }    

  if (nall<=0){
    cerr<<"no events available"<<endl;
    fEventNr = 0;
    return;
  }
  if (firstEvent+nEvents>nall) {
    fEventNr = nall-firstEvent;
    cerr<<"restricted number of events availaible"<<endl;
  }
  AliMagF * magf = gAlice->Field();
  AliTracker::SetFieldMap(magf);

}


////////////////////////////////////////////////////////////////////////
TPCCmpTr::~TPCCmpTr()
{
  //if (fLoader) {
  //  delete fLoader;
  //}
}

//////////////////////////////////////////////////////////////
Int_t TPCCmpTr::SetIO()
{
  //
  // 
  CreateTreeCmp();
  if (!fTreeCmp) return 1;
  fParamTPC = GetTPCParam();
  //
  if (!ConnectGenTree()) {
    cerr<<"Cannot connect tree with generated tracks"<<endl;
    return 1;
  }
  return 0;
}

//////////////////////////////////////////////////////////////

Int_t TPCCmpTr::SetIO(Int_t eventNr)
{
  //
  // 
  // SET INPUT
  //gAlice->GetEvent(eventNr);
  //  fLoader->SetEventNumber(eventNr);  
  //fLoader = AliRunLoader::Open("galice.root");
  fLoader->GetEvent(eventNr);  
  //
  AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
  tpcl->LoadTracks();
  fTreeRecTracks = tpcl->TreeT();
  
  //
  return 0;
}


////////////////////////////////////////////////////////////////////////
void TPCCmpTr::Reset()
{
  fEventNr = 0;
  fNEvents = 0;
  fTreeCmp = 0;
  //
  fFileGenTracks = 0;
  fDebug = 0;
  //
  fParamTPC = 0;
  fTreeRecTracks = 0;
  fTreePoints =0;

  fTPCTrack = 0; 
  fTracks   = 0;
  fTrackPoints =0;
}

////////////////////////////////////////////////////////////////////////
Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
{
  fNEvents = nEvents;
  fFirstEventNr = firstEventNr;
  return Exec();
}

////////////////////////////////////////////////////////////////////////
Int_t TPCCmpTr::Exec()
{
  TStopwatch timer;
  timer.Start();

  if (SetIO()==1) 
    return 1;
   
  cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
  for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
       eventNr++) {
    fEventNr=eventNr;
    SetIO(fEventNr);
    fNParticles = gAlice->GetEvent(fEventNr);    

    fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
    fFakeRecTracks = new Int_t[fNParticles];
    fMultiRecTracks = new Int_t[fNParticles];
    for (Int_t i = 0; i<fNParticles; i++) {
      fIndexRecTracks[4*i] = -1;
      fIndexRecTracks[4*i+1] = -1;
      fIndexRecTracks[4*i+2] = -1;
      fIndexRecTracks[4*i+3] = -1;

      fFakeRecTracks[i] = 0;
      fMultiRecTracks[i] = 0;
    }
  
    cout<<"Start to process event "<<fEventNr<<endl;
    cout<<"\tfNParticles = "<<fNParticles<<endl;
    if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
    if (TreeTLoop(eventNr)>0) return 1;

    if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
    if (TreeGenLoop(eventNr)>0) return 1;
    if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;

    delete [] fIndexRecTracks;
    delete [] fFakeRecTracks;
    delete [] fMultiRecTracks;
  }

  CloseOutputFile();

  cerr<<"Exec finished"<<endl;
  timer.Stop();
  timer.Print();
  return 0;

}
////////////////////////////////////////////////////////////////////////
Bool_t TPCCmpTr::ConnectGenTree()
{
//
// connect all branches from the genTracksTree
// use the same variables as for the new cmp tree, it may work
//
  fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
  if (!fFileGenTracks) {
    cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
    return kFALSE;
  }
  fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
  if (!fTreeGenTracks) {
    cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
	<<fFnGenTracks<<endl;
    return kFALSE;
  }
  //
  fMCInfo = new AliMCInfo;
  fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);


  if (fDebug > 1) {
    cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
  }
  return kTRUE;
}


////////////////////////////////////////////////////////////////////////
void TPCCmpTr::CreateTreeCmp() 
{
  fFileCmp = TFile::Open(fFnCmp,"RECREATE");
  if (!fFileCmp) {
    cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
    return;
  }


  fTreeCmp    = new TTree("TPCcmpTracks","TPCcmpTracks");
  //
  fMCInfo = new AliMCInfo;
  fRecInfo = new AliTPCRecInfo;
  //
  fTPCTrack = new AliTPCtrack;
   //
  fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
  fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
  fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
  //
  fTreeCmp->AutoSave(); 

}
////////////////////////////////////////////////////////////////////////
void TPCCmpTr::CloseOutputFile()  
{
  if (!fFileCmp) {
    cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
    return;
  }
  fFileCmp->cd();
  fTreeCmp->Write();    
  delete fTreeCmp;
  
  fFileCmp->Close();
  delete fFileCmp;
  return;
}
////////////////////////////////////////////////////////////////////////

TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
			    AliTPCParam *paramTPC) {

  Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
  Int_t index[4];
  paramTPC->Transform0to1(x,index);
  paramTPC->Transform1to2(x,index);
  return TVector3(x);
}
////////////////////////////////////////////////////////////////////////

Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
{
  //
  // loop over all TPC reconstructed tracks and store info in memory
  //
  TStopwatch  timer;
  timer.Start();
  
  if (!fTreeRecTracks) {
    cerr<<"Can't get a tree with TPC rec. tracks  "<<endl;
    return 1;
  }
  //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
  
  Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
  if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
		      <<nEntries<<endl;
  TBranch * br= fTreeRecTracks->GetBranch("tracks");
  br->SetAddress(&fTPCTrack);
  TBranch *brp = 0;
  if (fTreePoints) brp = fTreePoints->GetBranch("debug");

  if (fTracks){
    fTracks->Delete();    
    delete fTracks;
  }
  if (fTrackPoints){
    fTrackPoints->Delete();
    delete fTrackPoints;
    fTrackPoints =0;
  }
  fTracks      = new TObjArray(nEntries);
  if (brp){
    fTrackPoints = new TObjArray(nEntries);
  }
  else fTrackPoints = 0;

  //
  //load tracks to the memory
  for (Int_t i=0; i<nEntries;i++){
    AliTPCtrack * track = new AliTPCtrack;
    br->SetAddress(&track);
    br->GetEntry(i);
    fTracks->AddAt(track,i);
  }
  //
  //load track points to the memory
  if (brp) for (Int_t i=0; i<nEntries;i++){
    TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
    brp->SetAddress(&arr);
    brp->GetEntry(i);
    if (arr!=0)
      for (Int_t j=0;j<arr->GetEntriesFast();j++){
	AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
	if (point && point->fID>=0){
	  fTrackPoints->AddAt(arr,point->fID);
	  break;
	}
      }    
  }
  //

  //
  for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
    //br->GetEntry(iEntry);
    fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
    //
    Int_t label = fTPCTrack->GetLabel();
    Int_t absLabel = abs(label);
    if (absLabel < fNParticles) {
      //      fIndexRecTracks[absLabel] =  iEntry;
      if (label < 0) fFakeRecTracks[absLabel]++;
      
      if (fMultiRecTracks[absLabel]>0){
	if (fMultiRecTracks[absLabel]<4)
	  fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry; 	
      }
      else      
	fIndexRecTracks[absLabel*4] =  iEntry;
      fMultiRecTracks[absLabel]++;
    }
  }  
  printf("Time spended in TreeTLoop\n");
  timer.Print();
  
  if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;

  return 0;
}

////////////////////////////////////////////////////////////////////////
Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
{
//
// loop over all entries for a given event, find corresponding 
// rec. track and store in the fTreeCmp
//
  TStopwatch timer;
  timer.Start();
  Int_t entry = 0;
  Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
  TBranch * branch = fTreeCmp->GetBranch("RC");
  branch->SetAddress(&fRecInfo); // set all pointers

  while (entry < nParticlesTR) {
    fTreeGenTracks->GetEntry(entry);
    entry++;
    if (eventNr < fMCInfo->fEventNr) continue;
    if (eventNr > fMCInfo->fEventNr) continue;
    //
    if (fDebug > 2 && fMCInfo->fLabel < 10) {
      cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
    }
    fRecInfo->Reset();
    //
    if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
      fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
      
      //      if (nBytes > 0) {
      if (fTPCTrack) {
	//
	//
	//if multifound find track with biggest pt - returning tracks - loos energy
	//
	if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
	  Double_t pt = TMath::Abs(1/fTPCTrack->Get1Pt());
	  for (Int_t i=1;i<4;i++){
	    if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
	      AliTPCtrack * track = (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
	      Double_t pt2 = TMath::Abs(1/track->Get1Pt());	      
	      if  (pt2>pt)
		fTPCTrack = track;		   
	    }
	  }
	}
	//
	fRecInfo->fTP=0;
	fRecInfo->fTrack =*fTPCTrack; 
	fRecInfo->fReconstructed = 1;
	fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
	fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
	fRecInfo->Update(fMCInfo, fParamTPC, kTRUE, 1);
      }      
    } 
    else{
      fRecInfo->Update(fMCInfo, fParamTPC, kFALSE, 1);      
    }    
    fTreeCmp->Fill();
  }
  fTreeCmp->AutoSave();
  fTracks->Delete();
  fTPCTrack =0;
  if (fTrackPoints){
    fTrackPoints->Delete();
    delete fTrackPoints;
    fTrackPoints =0;
  } 
  printf("Time spended in TreeGenLoop\n");
  timer.Print();
  if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;  
  return 0;
}


////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////

void AliTPCComparisonDraw::SetIO(const char *fname)
{
  TFile *file = TFile::Open(fname);
  if (!file){
    printf("Comparison file couldn't be generated\n"); 
    return;
   }
  //
  fTree = (TTree*) file->Get("TPCcmpTracks");
  if (!fTree) {
    printf("no track comparison tree found\n");
    file->Close();
    delete file;
  }

}

void AliTPCComparisonDraw::Eff()
{
  
}

void AliTPCComparisonDraw::ResPt()
{
}
Commit	Line	Data
ae7d73d2	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	///////////////////////////////////////////////////////////////////////////////
	18	// //
	19	// Time Projection Chamber //
	20	// Comparison macro for TPC //
	21	// responsible:
	22	// marian.ivanov@cern.ch //
	23	//
	24	//
	25	/*
	26	//
	27
	28
	29	.L $ALICE_ROOT/STEER/AliGenInfo.C+
	30	.L $ALICE_ROOT/STEER/AliTPCComparisonMI.C+
	31
	32
	33	TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",-1,1,0);
	34	t2->Exec();
	35
	36
	37	TCut cprim("cprim","MC.fVDist[3]<1");
	38	TCut csec("cprim","MC.fVDist[3]>1");
	39	TCut crec("crec","fReconstructed==1");
	40	TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
	41	TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVxMC.fParticle.fVx+MC.fParticle.fVyMC.fParticle.fVy))<1");
	42	TCut csens("csens","abs(sqrt(fVDist[0]2+fVDist[1]2)-170)<50");
	43	TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
	44	AliTPCComparisonDraw comp;
	45	comp.SetIO();
	46
	47
	48	//example
	49	comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&fPdg==-211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
	50	comp.fRes->Draw();
	51	comp.fMean->Draw();
	52	comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1,"fReconstructed>0",10,0,1.5)
	53	comp.fRes->Draw();
	54	*/
	55
	56
	57
	58
	59
	60	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	61	#include <stdio.h>
	62	#include <string.h>
	63	//ROOT includes
	64	#include "Rtypes.h"
65	#include "TFile.h"
66	#include "TTree.h"
67	#include "TChain.h"
68	#include "TCut.h"
69	#include "TString.h"
70	#include "TBenchmark.h"
71	#include "TStopwatch.h"
72	#include "TParticle.h"
73	#include "TSystem.h"
74	#include "TTimer.h"
75	#include "TVector3.h"
76	#include "TPad.h"
77	#include "TCanvas.h"
78	#include "TH1F.h"
79	#include "TH2F.h"
80	#include "TF1.h"
81	#include "TText.h"
82	#include "Getline.h"
83	#include "TStyle.h"
84
85	//ALIROOT includes
86	#include "AliRun.h"
87	#include "AliStack.h"
88	#include "AliTPCtrack.h"
89	#include "AliSimDigits.h"
90	#include "AliTPCParam.h"
91	#include "AliTPC.h"
92	#include "AliTPCLoader.h"
93	#include "AliDetector.h"
94	#include "AliTrackReference.h"
95	#include "AliRun.h"
96	#include "AliTPCParamSR.h"
97	#include "AliTracker.h"
98	#include "AliComplexCluster.h"
99	#include "AliMagF.h"
100	#endif
101	#include "AliGenInfo.h"
102	#include "AliTPCComparisonMI.h"
103
104	//
105	//
106
107
108
109	void AliTPCRecInfo::Update(AliMCInfo* info,AliTPCParam * param, Bool_t reconstructed, Int_t direction)
110	{
111	//
112	//
113	//calculates derived variables
114	//
115
116	if (reconstructed==kFALSE){
117	fReconstructed = kFALSE;
118	return;
119	}
120	fReconstructed = kTRUE;
121	// Find the nearest track reference
122	//
123	Double_t radius = TMath::Sqrt(fTrack.GetX()fTrack.GetX()+fTrack.GetY()fTrack.GetY());
124	TClonesArray * tpcreferences = info->fTPCReferences;
125	Int_t nentries = tpcreferences->GetEntriesFast();
126	AliTrackReference * ref = 0;
127	Double_t dr = 1000000;
128	for (Int_t i=0;i<nentries;i++){
129	AliTrackReference * refi = (AliTrackReference*)tpcreferences->UncheckedAt(i);
130	if (!refi) continue;
131	//if ( (direction<0) && (refi->R()<radius-10)) continue;
132	//if ( (direction>0) && (refi->R()>radius+10)) continue;
133	if (TMath::Abs(refi->R()-radius)<dr){
134	dr = TMath::Abs(refi->R()-radius);
135	ref = refi;
136	}
137	}
138	if (ref==0) {
139	fReconstructed = kFALSE;
140	return;
141	}
142	//
143	//
144	fdEdx = fTrack.GetdEdx();
145	//
146	// propagate track to nearest reference in direction
147	TVector3 local = TPCCmpTr::TR2Local(ref,param);
148	local.GetXYZ(fTRLocalCoord);
149	//
150	// rotate if neccessary
151	Float_t ymax = local.X()TMath::Tan(0.5param->GetInnerAngle());
152	Float_t y = fTrack.GetYat(local.X());
153	if (y > ymax) {
154	fTrack.Rotate(param->GetInnerAngle());
155	} else if (y <-ymax) {
156	fTrack.Rotate(-param->GetInnerAngle());
157	}
158	Double_t xtrack = fTrack.GetX();
159	Double_t delta = (local.X()-xtrack)*0.1;
160	for (Int_t i=0;i<9;i++){
161	fTrack.PropagateTo(xtrack+delta*float(i));
162	}
163	fTrack.PropagateTo(local.X());
164
165	Double_t par[5], cov[15], localX;
166	fTrack.GetExternalParameters(localX,par);
167	fTrack.GetExternalCovariance(cov);
168	//
169	//
170	fRecPhi=TMath::ASin(par[2]) + fTrack.GetAlpha();
171	Double_t phi2 =TMath::ATan2(ref->Py(),ref->Px());
172	if (phi2<0) phi2+=2*TMath::Pi();
173	fGenPhi =phi2;
174	//
175	if (fRecPhi<0) fRecPhi+=2*TMath::Pi();
176	if (fRecPhi>=2TMath::Pi()) fRecPhi-=2TMath::Pi();
177	fLambda = TMath::ATan(par[3]);
178	fRecPt_1 = TMath::Abs(par[4]);
179	//TPC
180	//
181	//
182	Double_t phi=TMath::ATan2(par[0],localX) + fTrack.GetAlpha();
183	Double_t r=TMath::Sqrt(localXlocalX + par[0]par[0]);
184	fTPCinR1[0]= r*TMath::Cos(phi);
185	fTPCinR1[1]= r*TMath::Sin(phi);
186	fTPCinR1[2]= par[1];
187	fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]fTPCinR1[0]+fTPCinR1[1]fTPCinR1[1]);
188	fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);
189	//
190	fTPCinP1[0] = fTrack.Px();
191	fTPCinP1[1] = fTrack.Py();
192	fTPCinP1[2] = fTrack.Pz();
193	fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]fTPCinP1[0]+fTPCinP1[1]fTPCinP1[1]);
194	//
195	fTPCinR0[0] = ref->X();
196	fTPCinR0[1] = ref->Y();
197	fTPCinR0[2] = ref->Z();
198	fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]fTPCinR0[0]+fTPCinR0[1]fTPCinR0[1]);
199	fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
200	//
201	fTPCinP0[0]=ref->Px();
202	fTPCinP0[1]=ref->Py();
203	fTPCinP0[2]=ref->Pz();
204	fTPCinP0[3]=ref->Pt();
205	//
206	//
207	if (fTPCinP1[3]>0.0000001){
208	fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
209	fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
210	//
211	fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
212	fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
213	}
214	//
215	fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3]; //delta rfi
216	fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
217	fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]); //delta z
218	fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
219	//
220	fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
221	fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
222	//
223	fTPCDelta[3] = TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]);
224	fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
225	fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
226	Double_t sign = (fTrack.GetC()>0) ? 1:-1;
227	fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]);
228
229	}
230
231
232
233
234	////////////////////////////////////////////////////////////////////////
235	TPCCmpTr::TPCCmpTr()
236	{
237	Reset();
238	}
239
240	////////////////////////////////////////////////////////////////////////
241	TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
242	const char* fnCmp,
243	const char* fnGalice, Int_t direction,
244	Int_t nEvents, Int_t firstEvent)
245	{
246	Reset();
247	// fFnGenTracks = fnGenTracks;
248	// fFnCmp = fnCmp;
249	sprintf(fFnGenTracks,"%s",fnGenTracks);
250	sprintf(fFnCmp,"%s",fnCmp);
251
252	fFirstEventNr = firstEvent;
253	fEventNr = firstEvent;
254	fNEvents = nEvents;
255	fDirection = direction;
256	//
257	fLoader = AliRunLoader::Open(fnGalice);
258	if (gAlice){
259	//delete gAlice->GetRunLoader();
260	delete gAlice;
261	gAlice = 0x0;
262	}
263	if (fLoader->LoadgAlice()){
264	cerr<<"Error occured while l"<<endl;
265	}
266	Int_t nall = fLoader->GetNumberOfEvents();
267	if (nEvents==0) {
268	nEvents =nall;
269	fNEvents=nall;
270	fFirstEventNr=0;
271	}
272
273	if (nall<=0){
274	cerr<<"no events available"<<endl;
275	fEventNr = 0;
276	return;
277	}
278	if (firstEvent+nEvents>nall) {
279	fEventNr = nall-firstEvent;
280	cerr<<"restricted number of events availaible"<<endl;
281	}
282	AliMagF * magf = gAlice->Field();
283	AliTracker::SetFieldMap(magf);
284
285	}
286
287
288	////////////////////////////////////////////////////////////////////////
289	TPCCmpTr::~TPCCmpTr()
290	{
291	//if (fLoader) {
292	// delete fLoader;
293	//}
294	}
295
296	//////////////////////////////////////////////////////////////
297	Int_t TPCCmpTr::SetIO()
298	{
299	//
300	//
301	CreateTreeCmp();
302	if (!fTreeCmp) return 1;
303	fParamTPC = GetTPCParam();
304	//
305	if (!ConnectGenTree()) {
306	cerr<<"Cannot connect tree with generated tracks"<<endl;
307	return 1;
308	}
309	return 0;
310	}
311
312	//////////////////////////////////////////////////////////////
313
314	Int_t TPCCmpTr::SetIO(Int_t eventNr)
315	{
316	//
317	//
318	// SET INPUT
319	//gAlice->GetEvent(eventNr);
320	// fLoader->SetEventNumber(eventNr);
321	//fLoader = AliRunLoader::Open("galice.root");
322	fLoader->GetEvent(eventNr);
323	//
324	AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
325	tpcl->LoadTracks();
326	fTreeRecTracks = tpcl->TreeT();
327
328	//
329	return 0;
330	}
331
332
333
334	////////////////////////////////////////////////////////////////////////
335	void TPCCmpTr::Reset()
336	{
337	fEventNr = 0;
338	fNEvents = 0;
339	fTreeCmp = 0;
340	//
341	fFileGenTracks = 0;
342	fDebug = 0;
343	//
344	fParamTPC = 0;
345	fTreeRecTracks = 0;
346	fTreePoints =0;
347
348	fTPCTrack = 0;
349	fTracks = 0;
350	fTrackPoints =0;
351	}
352
353	////////////////////////////////////////////////////////////////////////
354	Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
355	{
356	fNEvents = nEvents;
357	fFirstEventNr = firstEventNr;
358	return Exec();
359	}
360
361	////////////////////////////////////////////////////////////////////////
362	Int_t TPCCmpTr::Exec()
363	{
364	TStopwatch timer;
365	timer.Start();
366
367	if (SetIO()==1)
368	return 1;
369
370	cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
371	for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
372	eventNr++) {
373	fEventNr=eventNr;
374	SetIO(fEventNr);
375	fNParticles = gAlice->GetEvent(fEventNr);
376
377	fIndexRecTracks = new Int_t[fNParticles*4]; //write at maximum 4 tracks corresponding to particle
378	fFakeRecTracks = new Int_t[fNParticles];
379	fMultiRecTracks = new Int_t[fNParticles];
380	for (Int_t i = 0; i<fNParticles; i++) {
381	fIndexRecTracks[4*i] = -1;
382	fIndexRecTracks[4*i+1] = -1;
383	fIndexRecTracks[4*i+2] = -1;
384	fIndexRecTracks[4*i+3] = -1;
385
386	fFakeRecTracks[i] = 0;
387	fMultiRecTracks[i] = 0;
388	}
389
390	cout<<"Start to process event "<<fEventNr<<endl;
391	cout<<"\tfNParticles = "<<fNParticles<<endl;
392	if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
393	if (TreeTLoop(eventNr)>0) return 1;
394
395	if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
396	if (TreeGenLoop(eventNr)>0) return 1;
397	if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
398
399	delete [] fIndexRecTracks;
400	delete [] fFakeRecTracks;
401	delete [] fMultiRecTracks;
402	}
403
404	CloseOutputFile();
405
406	cerr<<"Exec finished"<<endl;
407	timer.Stop();
408	timer.Print();
409	return 0;
410
411	}
412	////////////////////////////////////////////////////////////////////////
413	Bool_t TPCCmpTr::ConnectGenTree()
414	{
415	//
416	// connect all branches from the genTracksTree
417	// use the same variables as for the new cmp tree, it may work
418	//
419	fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
420	if (!fFileGenTracks) {
421	cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
422	return kFALSE;
423	}
424	fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
425	if (!fTreeGenTracks) {
426	cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
427	<<fFnGenTracks<<endl;
428	return kFALSE;
429	}
430	//
431	fMCInfo = new AliMCInfo;
432	fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
433
434
435	if (fDebug > 1) {
436	cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
437	}
438	return kTRUE;
439	}
440
441
442	////////////////////////////////////////////////////////////////////////
443	void TPCCmpTr::CreateTreeCmp()
444	{
445	fFileCmp = TFile::Open(fFnCmp,"RECREATE");
446	if (!fFileCmp) {
447	cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
448	return;
449	}
450
451
452	fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
453	//
454	fMCInfo = new AliMCInfo;
455	fRecInfo = new AliTPCRecInfo;
456	//
457	fTPCTrack = new AliTPCtrack;
458	//
459	fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
460	fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
461	fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
462	//
463	fTreeCmp->AutoSave();
464
465	}
466	////////////////////////////////////////////////////////////////////////
467	void TPCCmpTr::CloseOutputFile()
468	{
469	if (!fFileCmp) {
470	cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
471	return;
472	}
473	fFileCmp->cd();
474	fTreeCmp->Write();
475	delete fTreeCmp;
476
477	fFileCmp->Close();
478	delete fFileCmp;
479	return;
480	}
481	////////////////////////////////////////////////////////////////////////
482
483	TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
484	AliTPCParam *paramTPC) {
485
486	Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
487	Int_t index[4];
488	paramTPC->Transform0to1(x,index);
489	paramTPC->Transform1to2(x,index);
490	return TVector3(x);
491	}
492	////////////////////////////////////////////////////////////////////////
493
494	Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
495	{
496	//
497	// loop over all TPC reconstructed tracks and store info in memory
498	//
499	TStopwatch timer;
500	timer.Start();
501
502	if (!fTreeRecTracks) {
503	cerr<<"Can't get a tree with TPC rec. tracks "<<endl;
504	return 1;
505	}
506	//fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
507
508	Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
509	if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
510	<<nEntries<<endl;
511	TBranch * br= fTreeRecTracks->GetBranch("tracks");
512	br->SetAddress(&fTPCTrack);
513	TBranch *brp = 0;
514	if (fTreePoints) brp = fTreePoints->GetBranch("debug");
515
516	if (fTracks){
517	fTracks->Delete();
518	delete fTracks;
519	}
520	if (fTrackPoints){
521	fTrackPoints->Delete();
522	delete fTrackPoints;
523	fTrackPoints =0;
524	}
525	fTracks = new TObjArray(nEntries);
526	if (brp){
527	fTrackPoints = new TObjArray(nEntries);
528	}
529	else fTrackPoints = 0;
530
531	//
532	//load tracks to the memory
533	for (Int_t i=0; i<nEntries;i++){
534	AliTPCtrack * track = new AliTPCtrack;
535	br->SetAddress(&track);
536	br->GetEntry(i);
537	fTracks->AddAt(track,i);
538	}
539	//
540	//load track points to the memory
541	if (brp) for (Int_t i=0; i<nEntries;i++){
542	TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
543	brp->SetAddress(&arr);
544	brp->GetEntry(i);
545	if (arr!=0)
546	for (Int_t j=0;j<arr->GetEntriesFast();j++){
547	AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
548	if (point && point->fID>=0){
549	fTrackPoints->AddAt(arr,point->fID);
550	break;
551	}
552	}
553	}
554	//
555
556	//
557	for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
558	//br->GetEntry(iEntry);
559	fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
560	//
561	Int_t label = fTPCTrack->GetLabel();
562	Int_t absLabel = abs(label);
563	if (absLabel < fNParticles) {
564	// fIndexRecTracks[absLabel] = iEntry;
565	if (label < 0) fFakeRecTracks[absLabel]++;
566
567	if (fMultiRecTracks[absLabel]>0){
568	if (fMultiRecTracks[absLabel]<4)
569	fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] = iEntry;
570	}
571	else
572	fIndexRecTracks[absLabel*4] = iEntry;
573	fMultiRecTracks[absLabel]++;
574	}
575	}
576	printf("Time spended in TreeTLoop\n");
577	timer.Print();
578
579	if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
580
581	return 0;
582	}
583
584	////////////////////////////////////////////////////////////////////////
585	Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
586	{
587	//
588	// loop over all entries for a given event, find corresponding
589	// rec. track and store in the fTreeCmp
590	//
591	TStopwatch timer;
592	timer.Start();
593	Int_t entry = 0;
594	Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
595	TBranch * branch = fTreeCmp->GetBranch("RC");
596	branch->SetAddress(&fRecInfo); // set all pointers
597
598	while (entry < nParticlesTR) {
599	fTreeGenTracks->GetEntry(entry);
600	entry++;
601	if (eventNr < fMCInfo->fEventNr) continue;
602	if (eventNr > fMCInfo->fEventNr) continue;
603	//
604	if (fDebug > 2 && fMCInfo->fLabel < 10) {
605	cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
606	}
607	fRecInfo->Reset();
608	//
609	if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
610	fTPCTrack= (AliTPCtrack)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel4]);
611
612	// if (nBytes > 0) {
613	if (fTPCTrack) {
614	//
615	//
616	//if multifound find track with biggest pt - returning tracks - loos energy
617	//
618	if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
619	Double_t pt = TMath::Abs(1/fTPCTrack->Get1Pt());
620	for (Int_t i=1;i<4;i++){
621	if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
622	AliTPCtrack * track = (AliTPCtrack)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel4+i]);
623	Double_t pt2 = TMath::Abs(1/track->Get1Pt());
624	if (pt2>pt)
625	fTPCTrack = track;
626	}
627	}
628	}
629	//
630	fRecInfo->fTP=0;
631	fRecInfo->fTrack =*fTPCTrack;
632	fRecInfo->fReconstructed = 1;
633	fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
634	fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
635	fRecInfo->Update(fMCInfo, fParamTPC, kTRUE, 1);
636	}
637	}
638	else{
639	fRecInfo->Update(fMCInfo, fParamTPC, kFALSE, 1);
640	}
641	fTreeCmp->Fill();
642	}
643	fTreeCmp->AutoSave();
644	fTracks->Delete();
645	fTPCTrack =0;
646	if (fTrackPoints){
647	fTrackPoints->Delete();
648	delete fTrackPoints;
649	fTrackPoints =0;
650	}
651	printf("Time spended in TreeGenLoop\n");
652	timer.Print();
653	if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
654	return 0;
655	}
656
657
658	////////////////////////////////////////////////////////////////////////
659	////////////////////////////////////////////////////////////////////////
660
661	void AliTPCComparisonDraw::SetIO(const char *fname)
662	{
663	TFile *file = TFile::Open(fname);
664	if (!file){
665	printf("Comparison file couldn't be generated\n");
666	return;
667	}
668	//
669	fTree = (TTree*) file->Get("TPCcmpTracks");
670	if (!fTree) {
671	printf("no track comparison tree found\n");
672	file->Close();
673	delete file;
674	}
675
676	}
677
678	void AliTPCComparisonDraw::Eff()
679	{
680
681	}
682
683	void AliTPCComparisonDraw::ResPt()
684	{
685	}