From c92725b70caa0ab534fb3f66b38f75a98fc567cd Mon Sep 17 00:00:00 2001 From: hristov Date: Fri, 18 May 2007 16:45:01 +0000 Subject: [PATCH] First version of some classes for comparison and performance studies (Marian) --- PWG1/AliESDComparisonMI.C | 1699 +++++++++++++++++++++++++++++++++++++ PWG1/AliESDComparisonMI.h | 233 +++++ PWG1/AliGenInfo.cxx | 1360 +++++++++++++++++++++++++++++ PWG1/AliGenInfo.h | 237 ++++++ PWG1/AliTreeDraw.cxx | 387 +++++++++ PWG1/AliTreeDraw.h | 62 ++ PWG1/PWG1LinkDef.h | 16 + PWG1/libPWG1.pkg | 9 + PWG1/rootlogon.C | 20 + 9 files changed, 4023 insertions(+) create mode 100644 PWG1/AliESDComparisonMI.C create mode 100644 PWG1/AliESDComparisonMI.h create mode 100644 PWG1/AliGenInfo.cxx create mode 100644 PWG1/AliGenInfo.h create mode 100644 PWG1/AliTreeDraw.cxx create mode 100644 PWG1/AliTreeDraw.h create mode 100644 PWG1/PWG1LinkDef.h create mode 100644 PWG1/libPWG1.pkg create mode 100644 PWG1/rootlogon.C diff --git a/PWG1/AliESDComparisonMI.C b/PWG1/AliESDComparisonMI.C new file mode 100644 index 00000000000..39c5b6cdb98 --- /dev/null +++ b/PWG1/AliESDComparisonMI.C @@ -0,0 +1,1699 @@ +/************************************************************************** + * 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 ESD // +// responsible: +// marian.ivanov@cern.ch // +// +// + +/* +marian.ivanov@cern.ch +Usage: + + +.L $ALICE_ROOT/STEER/AliGenInfo.C+ +//be sure you created genTracks file before +.L $ALICE_ROOT/STEER/AliESDComparisonMI.C+ +// +ESDCmpTr *t2 = new ESDCmpTr("genTracks.root","cmpESDTracks.root","galice.root",-1,0,0); +t2->Exec(); + +// +//some cuts definition +TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01") +//TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5") +//TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9"); +TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5"); +TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5"); + + +TCut crec("crec","fReconstructed==1"); +TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25"); +TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1"); + +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)"); +TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.") + +AliESDComparisonDraw comp; +comp.SetIO(); +TFile f("genHits.root"); +TTree * treel = (TTree*)f.Get("HitLines"); +if (treel) comp->fTree->AddFriend(treel,"L"); + +// +//example +comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)"); +comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy"); +comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy"); +comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()"); +comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)"); +comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); +comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); + + +TH1F his("his","his",100,0,20); +TH1F hpools("hpools","hpools",100,-7,7); +TH1F hfake("hfake","hfake",1000,0,150); +TProfile profp0("profp0","profp0",20,-0.4,0.9) + +comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06) +comp.fRes->Draw(); +comp.fMean->Draw(); + +comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06) +comp.fRes->Draw(); + +comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5) +comp.fRes->Draw(); + +comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio<0.1",10,0.2,1.5) +comp.fRes->Draw(); +comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio>0.1",10,0.2,1.5) +comp.fRes->Draw(); + +comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fESDTrack.fITSFakeRatio==0") + +TH1F his("his","his",100,0,20); +TH1F hpools("hpools","hpools",100,-7,7); + +TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1); +TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4); +TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3); +TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2); + +comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim) +comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim) +comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim) +comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim) + + +comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1") +comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1") +comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1") +comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1") + +comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx0","fTPCOn&&abs(fPdg)==211&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); +comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx1","fTPCOn&&abs(fPdg)==2212&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); +comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx2","fTPCOn&&abs(fPdg)==321&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); +comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx3","fTPCOn&&abs(fPdg)==11&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); + +hdedx3->SetXTitle("P(GeV/c)"); +hdedx3->SetYTitle("dEdx(unit)"); +hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same"); + +comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8) + +TProfile prof("prof","prof",10,0.5,5); + + + + +*/ + + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +//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 "AliESDtrack.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" +#include "AliESD.h" +#include "AliESDfriend.h" +#include "AliESDtrack.h" +#include "AliITStrackMI.h" +#include "AliTRDtrack.h" +#include "AliHelix.h" +#include "AliESDVertex.h" +#include "AliExternalTrackParam.h" +#include "AliESDkink.h" +#include "AliESDv0.h" +#include "AliV0.h" + +#endif +#include "AliGenInfo.h" +#include "AliESDComparisonMI.h" + + + + + +void MakeAliases(AliESDComparisonDraw&comp) +{ + // + // aliases definition + // + comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)"); + comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy"); + comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy"); + comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()"); + comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)"); + comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); + comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); + + comp.fTree->SetAlias("trddedx","(RC.fESDTrack.fTRDsignals[0]+RC.fESDTrack.fTRDsignals[1]+RC.fESDTrack.fTRDsignals[2]+RC.fESDTrack.fTRDsignals[3]+RC.fESDTrack.fTRDsignals[4]+RC.fESDTrack.fTRDsignals[5])/6."); + + comp.fTree->SetAlias("dtofmc2","fESDTrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)"); + comp.fTree->SetAlias("dtofrc2","(fESDTrack.fTrackTime[2]-fESDTrack.fTOFsignal)"); + + comp.fTree->SetAlias("psum","fESDTrack.fTOFr[4]+fESDTrack.fTOFr[3]+fESDTrack.fTOFr[2]+fESDTrack.fTOFr[1]+fESDTrack.fTOFr[0]"); + comp.fTree->SetAlias("P0","fESDTrack.fTOFr[0]/psum"); + comp.fTree->SetAlias("P1","fESDTrack.fTOFr[1]/psum"); + comp.fTree->SetAlias("P2","fESDTrack.fTOFr[2]/psum"); + comp.fTree->SetAlias("P3","fESDTrack.fTOFr[3]/psum"); + comp.fTree->SetAlias("P4","fESDTrack.fTOFr[4]/psum"); + comp.fTree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)"); +} + + + + +void AliESDRecInfo::UpdatePoints(AliESDtrack*track) +{ + // + // + Int_t iclusters[200]; + Float_t density[160]; + for (Int_t i=0;i<160;i++) density[i]=-1.; + fTPCPoints[0]= 160; + fTPCPoints[1] = -1; + // + if (fTPCPoints[0]GetTPCclusters(iclusters); + + Int_t ngood=0; + Int_t undeff=0; + Int_t nall =0; + Int_t range=20; + for (Int_t i=0;i<160;i++){ + Int_t last = i-range; + if (nall=0){ + if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--; + if (iclusters[last]==-1) undeff--; + } + if (iclusters[i]>0&& (iclusters[i]&0x8000)==0) ngood++; + if (iclusters[i]==-1) undeff++; + if (nall==range &&undeffmaxdens){ + maxdens=density[i]; + indexmax=i; + } + } + // + //max dens point + fTPCPoints[3] = maxdens; + fTPCPoints[1] = indexmax; + // + // last point + for (Int_t i=indexmax;i<160;i++){ + if (density[i]<0) continue; + if (density[i]0;i--){ + if (density[i]<0) continue; + if (density[i]0;i--){ + if (iclusters[i]==-1) continue; //dead zone + nall++; + if (iclusters[i]>0) ngood++; + if (nall>20) break; + } + fTPCPoints[4] = Float_t(ngood)/Float_t(nall); + // + if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1; + // + // + // check TRDPoints + /* + nclusters=track->GetTRDclusters(iclusters); + for (Int_t i=nclusters;i>0;i--){ + + } + */ + + +} + +// +// +void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed, AliESD */*event*/) +{ + // + // + //calculates derived variables + // + // + UpdatePoints(&fESDTrack); + fBestTOFmatch=1000; + AliTrackReference * ref = &(info->fTrackRef); + fTPCinR0[0] = info->fTrackRef.X(); + fTPCinR0[1] = info->fTrackRef.Y(); + fTPCinR0[2] = info->fTrackRef.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(); + fTPCinP0[4] = ref->P(); + fDeltaP = (ref->P()-info->fParticle.P())/info->fParticle.P(); + // + // + if (fTPCinP0[3]>0.0000001){ + // + fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]); + fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]); + } + // + // + fITSinP0[0]=info->fParticle.Px(); + fITSinP0[1]=info->fParticle.Py(); + fITSinP0[2]=info->fParticle.Pz(); + fITSinP0[3]=info->fParticle.Pt(); + // + fITSinR0[0]=info->fParticle.Vx(); + fITSinR0[1]=info->fParticle.Vy(); + fITSinR0[2]=info->fParticle.Vz(); + fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]); + fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]); + // + // + if (fITSinP0[3]>0.0000001){ + fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]); + fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]); + } + // + for (Int_t i=0;i<4;i++) fStatus[i] =0; + fReconstructed = kFALSE; + fTPCOn = kFALSE; + fITSOn = kFALSE; + fTRDOn = kFALSE; + if (reconstructed==kFALSE) return; + + fLabels[0] = info->fLabel; + fLabels[1] = info->fPrimPart; + fReconstructed = kTRUE; + fTPCOn = ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE; + fITSOn = ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE; + fTRDOn = ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE; + // + // + if ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0){ + fStatus[1] =3; + } + else{ + if ((fESDTrack.GetStatus()&AliESDtrack::kTPCout)>0){ + fStatus[1] =2; + } + else{ + if ((fESDTrack.GetStatus()&AliESDtrack::kTPCin)>0) + fStatus[1]=1; + } + } + // + if ((fESDTrack.GetStatus()&AliESDtrack::kITSout)>0){ + fStatus[0] =2; + } + else{ + if ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0){ + fStatus[0] =1; + } + else{ + fStatus[0]=0; + } + } + + // + // + if ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0){ + fStatus[2] =2; + } + else{ + if ((fESDTrack.GetStatus()&AliESDtrack::kTRDout)>0){ + fStatus[2] =1; + } + } + if ((fESDTrack.GetStatus()&AliESDtrack::kTRDStop)>0){ + fStatus[2] =10; + } + + // + //TOF + // + if (((fESDTrack.GetStatus()&AliESDtrack::kTOFout)>0)){ + // + // best tof match + Double_t times[5]; + fESDTrack.GetIntegratedTimes(times); + for (Int_t i=0;i<5;i++){ + if ( TMath::Abs(fESDTrack.GetTOFsignal()-times[i]) fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE; // decay product of original particle + fStatus[3]=1; + } + if (toffake) fStatus[3] =3; //total fake + if (tofdaughter) fStatus[3]=2; //fake because of decay + }else{ + fStatus[3]=0; + } + + + if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){ + //TPC + fESDTrack.GetInnerXYZ(fTPCinR1); + fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]); + fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]); + fESDTrack.GetInnerPxPyPz(fTPCinP1); + fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]); + fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]); + // + // + if (fTPCinP1[3]>0.000000000000001){ + fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]); + fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]); + } + Double_t cov[15], param[5],x, alpha; + fESDTrack.GetInnerExternalCovariance(cov); + fESDTrack.GetInnerExternalParameters(alpha, x,param); + if (x<50) return ; + // + 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 = (param[4]>0)? 1.:-1; + fSign =sign; + fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(TMath::Abs(cov[14])); + } + if (fITSOn){ + // ITS + Double_t param[5],x; + fESDTrack.GetExternalParameters(x,param); + // fESDTrack.GetConstrainedExternalParameters(x,param); + Double_t cov[15]; + fESDTrack.GetExternalCovariance(cov); + //fESDTrack.GetConstrainedExternalCovariance(cov); + if (TMath::Abs(param[4])<0.0000000001) return; + + fESDTrack.GetXYZ(fITSinR1); + fESDTrack.GetPxPyPz(fITSinP1); + fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]); + // + fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]); + fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]); + // + // + if (fITSinP1[3]>0.0000001){ + fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]); + fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]); + } + // + // + fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3]; //delta rfi + fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]); + fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]); //delta z + fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]); + fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]); + fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]); + fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1])); + fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]); + fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]); + Double_t sign = (param[4]>0) ? 1:-1; + fSign = sign; + fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]); + } + +} + + +void AliESDRecV0Info::Update(Float_t vertex[3]) +{ + + if ( (fT1.fStatus[1]>0)&& (fT2.fStatus[1]>0)){ + Float_t distance1,distance2; + Double_t xx[3],pp[3]; + // + Double_t xd[3],pd[3],signd; + Double_t xm[3],pm[3],signm; + // + // + if (fT1.fITSOn&&fT2.fITSOn){ + for (Int_t i=0;i<3;i++){ + xd[i] = fT2.fITSinR1[i]; + pd[i] = fT2.fITSinP1[i]; + xm[i] = fT1.fITSinR1[i]; + pm[i] = fT1.fITSinP1[i]; + } + } + else{ + + for (Int_t i=0;i<3;i++){ + xd[i] = fT2.fTPCinR1[i]; + pd[i] = fT2.fTPCinP1[i]; + xm[i] = fT1.fTPCinR1[i]; + pm[i] = fT1.fTPCinP1[i]; + } + } + // + // + signd = fT2.fSign<0 ? -1:1; + signm = fT1.fSign<0 ? -1:1; + + AliHelix dhelix1(xd,pd,signd); + dhelix1.GetMomentum(0,pp,0); + dhelix1.Evaluate(0,xx); + // + // Double_t x2[3],p2[3]; + // + AliHelix mhelix(xm,pm,signm); + // + //find intersection linear + // + Double_t phase[2][2],radius[2]; + Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200); + Double_t delta1=10000,delta2=10000; + + if (points==1){ + fRs[0] = TMath::Sqrt(radius[0]); + fRs[1] = TMath::Sqrt(radius[0]); + } + if (points==2){ + fRs[0] =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1])); + fRs[1] =TMath::Max(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1])); + } + + if (points>0){ + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + if (points==2){ + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + if (points==1){ + fRs[0] = TMath::Sqrt(radius[0]); + fRs[1] = TMath::Sqrt(radius[0]); + fDistMinR = delta1; + } + if (points==2){ + if (radius[0]0){ + dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + if (points==2){ + dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + + distance2 = TMath::Min(delta1,delta2); + if (distance2>100) fDist2 =100; + return; + if (delta10){ + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + if (points==2){ + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + if (points==1){ + fMinR = TMath::Sqrt(radius[0]); + fDistMinR = delta1; + } + if (points==2){ + if (radius[0]0){ + dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + if (points==2){ + dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + + distance2 = TMath::Min(delta1,delta2); + if (delta1GetRunLoader(); + delete gAlice; + gAlice = 0x0; + } + if (fLoader->LoadgAlice()){ + cerr<<"Error occured while l"<GetNumberOfEvents(); + if (nEvents==0) { + nEvents =nall; + fNEvents=nall; + fFirstEventNr=0; + } + + if (nall<=0){ + cerr<<"no events available"<nall) { + fEventNr = nall-firstEvent; + cerr<<"restricted number of events availaible"<Field(); + AliTracker::SetFieldMap(magf,0); + +} + + +//////////////////////////////////////////////////////////////////////// +ESDCmpTr::~ESDCmpTr() +{ + if (fLoader) { + delete fLoader; + } +} + +////////////////////////////////////////////////////////////// +Int_t ESDCmpTr::SetIO() +{ + // + // + CreateTreeCmp(); + if (!fTreeCmp) return 1; + fParamTPC = GetTPCParam(); + // + if (!ConnectGenTree()) { + cerr<<"Cannot connect tree with generated tracks"<SetBranchStatus("*",1); + tree->SetBranchAddress("ESD", &fEvent); + tree->SetBranchAddress("ESDfriend.",&fESDfriend); + tree->GetEntry(eventNr); + fEvent->SetESDfriend(fESDfriend); + } + + + /* + Char_t ename[100]; + sprintf(ename,"%d",eventNr); + fEvent = (AliESD*)f.Get(ename); + if (!fEvent){ + sprintf(ename,"ESD%d",eventNr); + fEvent = (AliESD*)f.Get(ename); + } + + TTree* tree = (TTree*) f.Get("esdTree"); + if (!tree) { + Error("CheckESD", "no ESD tree found"); + return kFALSE; + } + tree->SetBranchAddress("ESD", &fEvent); + tree->GetEntry(eventNr); + */ + + if (!fEvent) return 1; + + return 0; +} + + + +//////////////////////////////////////////////////////////////////////// +void ESDCmpTr::Reset() +{ + fEventNr = 0; + fNEvents = 0; + fTreeCmp = 0; + fTreeCmpKinks =0; + fTreeCmpV0 =0; + // fFnCmp = "cmpTracks.root"; + fFileGenTracks = 0; + fDebug = 0; + // + fParamTPC = 0; + fEvent =0; +} + +//////////////////////////////////////////////////////////////////////// +Int_t ESDCmpTr::Exec(Int_t nEvents, Int_t firstEventNr) +{ + fNEvents = nEvents; + fFirstEventNr = firstEventNr; + return Exec(); +} + +//////////////////////////////////////////////////////////////////////// +Int_t ESDCmpTr::Exec() +{ + TStopwatch timer; + timer.Start(); + + if (SetIO()==1) + return 1; + + fNextTreeGenEntryToRead = 0; + fNextKinkToRead = 0; + fNextV0ToRead =0; + cerr<<"fFirstEventNr, fNEvents: "<GetEvent(fEventNr); + + fIndexRecTracks = new Short_t[fNParticles*20]; //write at maximum 4 tracks corresponding to particle + fIndexRecKinks = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle + fIndexRecV0 = new Short_t[fNParticles*20]; //write at maximum 20 tracks corresponding to particle + + fFakeRecTracks = new Short_t[fNParticles]; + fMultiRecTracks = new Short_t[fNParticles]; + fMultiRecKinks = new Short_t[fNParticles]; + fMultiRecV0 = new Short_t[fNParticles]; + + for (Int_t i = 0; i0) return 1; + + if (fDebug>2) cout<<"\tStart loop over tree genTracks"<0) return 1; + BuildKinkInfo0(eventNr); + BuildV0Info(eventNr); + fRecArray->Delete(); + + if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<Get("genTracksTree"); + if (!fTreeGenTracks) { + cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " + <SetBranchAddress("MC",&fMCInfo); + // + // + fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree"); + if (!fTreeGenKinks) { + cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " + <SetBranchAddress("MC",&fGenKinkInfo); + } + + fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree"); + if (!fTreeGenV0) { + cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " + <SetBranchAddress("MC",&fGenV0Info); + } + // + if (fDebug > 1) { + cout<<"Number of gen. tracks with TR: "<GetEntries()<Branch("MC","AliMCInfo",&fMCInfo,256000); + fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo,256000); + // fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack); + delete esdTrack; + // + // + fTreeCmpKinks = new TTree("ESDcmpKinks","ESDcmpKinks"); + fGenKinkInfo = new AliGenKinkInfo; + fRecKinkInfo = new AliESDRecKinkInfo; + fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo,256000); + fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo,256000); + // + // + fTreeCmpV0 = new TTree("ESDcmpV0","ESDcmpV0"); + fGenV0Info = new AliGenV0Info; + fRecV0Info = new AliESDRecV0Info; + fTreeCmpV0->Branch("MC.","AliGenV0Info", &fGenV0Info,256000); + fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info,256000); + // + fTreeCmp->AutoSave(); + fTreeCmpKinks->AutoSave(); + fTreeCmpV0->AutoSave(); +} +//////////////////////////////////////////////////////////////////////// +void ESDCmpTr::CloseOutputFile() +{ + if (!fFileCmp) { + cerr<<"File "<cd(); + fTreeCmp->Write(); + delete fTreeCmp; + + fFileCmp->Close(); + delete fFileCmp; + return; +} +//////////////////////////////////////////////////////////////////////// + +TVector3 ESDCmpTr::TR2Local(AliTrackReference *trackRef, + AliTPCParam *paramTPC) { + + Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()}; + Int_t index[4]; + paramTPC->Transform0to1(x,index); + paramTPC->Transform1to2Ideal(x,index); + return TVector3(x); +} +//////////////////////////////////////////////////////////////////////// + +Int_t ESDCmpTr::TreeTLoop() +{ + // + // loop over all ESD reconstructed tracks and store info in memory + // + // + loop over all reconstructed kinks + TStopwatch timer; + timer.Start(); + // + Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks(); + Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks(); + Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0s(); + fSignedKinks = new Short_t[nKinks]; + fSignedV0 = new Short_t[nV0MIs]; + // + // load kinks to the memory + for (Int_t i=0; iGetKink(i); + fSignedKinks[i]=0; + if (kink->fStatus<0) continue; + } + // + for (Int_t i=0; iGetV0(i); + fSignedV0[i]=0; + if (v0MI->fStatus<0) continue; + } + + // + // + AliESDtrack * track=0; + for (Int_t iEntry=0; iEntryUncheckedAt(iEntry); + track = (AliESDtrack*)fEvent->GetTrack(iEntry); + // + Int_t label = track->GetLabel(); + Int_t absLabel = abs(label); + if (absLabel < fNParticles) { + // fIndexRecTracks[absLabel] = iEntry; + if (label < 0) fFakeRecTracks[absLabel]++; + if (fMultiRecTracks[absLabel]>0){ + if (fMultiRecTracks[absLabel]<20) + fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] = iEntry; + } + else + fIndexRecTracks[absLabel*20] = iEntry; + fMultiRecTracks[absLabel]++; + } + } + // sort reconstructed kinks + // + AliESDkink * kink=0; + for (Int_t iEntry=0; iEntryGetKink(iEntry); + if (!kink) continue; + // + Int_t label0 = TMath::Abs(kink->GetLabel(0)); + Int_t label1 = TMath::Abs(kink->GetLabel(1)); + Int_t absLabel = TMath::Min(label0,label1); + if (absLabel < fNParticles) { + if (fMultiRecKinks[absLabel]>0){ + if (fMultiRecKinks[absLabel]<20) + fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] = iEntry; + } + else + fIndexRecKinks[absLabel*20] = iEntry; + fMultiRecKinks[absLabel]++; + } + } + // --sort reconstructed V0 + // + AliV0 * v0MI=0; + for (Int_t iEntry=0; iEntryGetV0(iEntry); + if (!v0MI) continue; + // + // Int_t label0 = TMath::Abs(v0MI->fLab[0]); + //Int_t label1 = TMath::Abs(v0MI->fLab[1]); + // + for (Int_t i=0;i<2;i++){ + // Int_t absLabel = TMath::Min(label0,label1); + // Int_t absLabel = TMath::Abs(v0MI->fLab[i]); + Int_t absLabel = TMath::Abs(0); + if (absLabel < fNParticles) { + if (fMultiRecV0[absLabel]>0){ + if (fMultiRecV0[absLabel]<20) + fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] = iEntry; + } + else + fIndexRecV0[absLabel*20] = iEntry; + fMultiRecV0[absLabel]++; + } + } + } + + + printf("Time spended in TreeTLoop\n"); + timer.Print(); + + if (fDebug > 2) cerr<<"end of TreeTLoop"<GetEntriesFast(); + cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<GetBranch("RC"); + TBranch * branchF = fTreeCmp->GetBranch("F"); + + branch->SetAddress(&fRecInfo); // set all pointers + fRecArray = new TObjArray(fNParticles); + AliESDtrack dummytrack; // + AliESDfriendTrack dummytrackF; // + + while (entry < nParticlesTR) { + fTreeGenTracks->GetEntry(entry); + entry++; + if (eventNr < fMCInfo->fEventNr) continue; + if (eventNr > fMCInfo->fEventNr) continue;; + // + fNextTreeGenEntryToRead = entry-1; + if (fDebug > 2 && fMCInfo->fLabel < 10) { + cerr<<"Fill track with a label "<fLabel<fNTPCRef<1) continue; // not TPCref + // + fRecInfo->Reset(); + AliESDtrack * track=0; + fRecInfo->fReconstructed =0; + TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC); + local.GetXYZ(fRecInfo->fTRLocalCoord); + // + if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) { + //track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]); + track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]); + // + // + // find nearest track if multifound + //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge); + // + Int_t status = 0; + if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++; + if ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++; + if ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++; + + // + if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){ + // + Double_t p[3]; + track->GetInnerPxPyPz(p); + Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2]; + // + for (Int_t i=1;i<20;i++){ + if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){ + AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]); + if (!track2) continue; + //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; // + //if (sign2<0) continue; + track2->GetInnerPxPyPz(p); + Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2]; + /* + if (sign<0){ + sign = sign2; + track = track2; + maxp = mom; + continue; + } + */ + // + Int_t status2 = 0; + if ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++; + if ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++; + if ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++; + if (status2fESDTrack)) AliESDtrack(*track); + if (track->GetFriendTrack()){ + new( &(fRecInfo->fTrackF)) AliESDfriendTrack(*(track->GetFriendTrack())); + } + }else{ + fRecInfo->fESDTrack = dummytrack; + fRecInfo->fTrackF = dummytrackF; + } + + + + if (track->GetITStrack()) + //fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack()); + new (&(fRecInfo->fITStrack)) AliITStrackMI(*((AliITStrackMI*)track->GetITStrack())); + else{ + fRecInfo->fITStrack = *track; + } + if (track->GetTRDtrack()){ + fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack()); + } + else{ + fRecInfo->fTRDtrack.SetdEdx(-1); + } + fRecInfo->fReconstructed = 1; + fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel]; + fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel]; + // + fRecInfo->Update(fMCInfo,fParamTPC,kTRUE, fEvent); + } + else{ + fRecInfo->fESDTrack = dummytrack; + fRecInfo->Update(fMCInfo,fParamTPC,kFALSE, fEvent); + } + fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel); + fTreeCmp->Fill(); + } + fTreeCmp->AutoSave(); + //fTracks->Delete(); + printf("Time spended in TreeGenLoop\n"); + timer.Print(); + if (fDebug > 2) cerr<<"end of TreeGenLoop"<GetEntriesFast(); + cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<GetBranch("RC."); + branch->SetAddress(&fRecKinkInfo); // set all pointers + + // + while (entry < nParticlesTR) { + fTreeGenKinks->GetEntry(entry); + entry++; + if (eventNr < fGenKinkInfo->fMCm.fEventNr) continue; + if (eventNr > fGenKinkInfo->fMCm.fEventNr) continue;; + // + fNextKinkToRead = entry-1; + // + // + AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCm.fLabel); + AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCd.fLabel); + fRecKinkInfo->fT1 = (*fRecInfo1); + fRecKinkInfo->fT2 = (*fRecInfo2); + fRecKinkInfo->fStatus =0; + if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1; + if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2; + if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1; + + if (fRecKinkInfo->fStatus==3){ + fRecKinkInfo->Update(); + } + Int_t label = TMath::Min(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel); + Int_t label2 = TMath::Max(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel); + + AliESDkink *kink=0; + fRecKinkInfo->fRecStatus =0; + fRecKinkInfo->fMultiple = fMultiRecKinks[label]; + fRecKinkInfo->fKinkMultiple=0; + // + if (fMultiRecKinks[label]>0){ + + // for (Int_t j=0;j=0;j--){ + Int_t index = fIndexRecKinks[label*20+j]; + //AliESDkink *kink2 = (AliESDkink*)fKinks->At(index); + AliESDkink *kink2 = (AliESDkink*)fEvent->GetKink(index); + if (TMath::Abs(kink2->fLab[0])==label &&TMath::Abs(kink2->fLab[1])==label2) { + fRecKinkInfo->fKinkMultiple++; + fSignedKinks[index]=1; + Int_t c0=0; + if (kink){ + // if (kink->fTRDOn) c0++; + //if (kink->fITSOn) c0++; + if (kink->GetStatus(2)>0) c0++; + if (kink->GetStatus(0)>0) c0++; + } + Int_t c2=0; + // if (kink2->fTRDOn) c2++; + //if (kink2->fITSOn) c2++; + if (kink2->GetStatus(2)>0) c2++; + if (kink2->GetStatus(0)>0) c2++; + + if (c2fLab[1])==label &&TMath::Abs(kink2->fLab[0])==label2) { + fRecKinkInfo->fKinkMultiple++; + fSignedKinks[index]=1; + Int_t c0=0; + if (kink){ + //if (kink->fTRDOn) c0++; + //if (kink->fITSOn) c0++; + if (kink->GetStatus(2)>0) c0++; + if (kink->GetStatus(0)>0) c0++; + + } + Int_t c2=0; + // if (kink2->fTRDOn) c2++; + //if (kink2->fITSOn) c2++; + if (kink2->GetStatus(2)>0) c2++; + if (kink2->GetStatus(0)>0) c2++; + + if (c2fKink = *kink; + fRecKinkInfo->fRecStatus=1; + } + fTreeCmpKinks->Fill(); + } + // Int_t nkinks = fKinks->GetEntriesFast(); + Int_t nkinks = fEvent->GetNumberOfKinks(); + for (Int_t i=0;iAt(i); + AliESDkink *kink = (AliESDkink*)fEvent->GetKink(i); + if (!kink) continue; + // + fRecKinkInfo->fKink = *kink; + fRecKinkInfo->fRecStatus =-2; + // + AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[0])); + AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[1])); + if (fRecInfo1 && fRecInfo2){ + fRecKinkInfo->fT1 = (*fRecInfo1); + fRecKinkInfo->fT2 = (*fRecInfo2); + fRecKinkInfo->fRecStatus =-1; + } + fTreeCmpKinks->Fill(); + } + } + + + fTreeCmpKinks->AutoSave(); + printf("Time spended in BuilKinkInfo Loop\n"); + timer.Print(); + if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<GetEntriesFast(); + cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<GetBranch("RC."); + branch->SetAddress(&fRecV0Info); // set all pointers + const AliESDVertex * esdvertex = fEvent->GetVertex(); + Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()}; + + // + while (entry < nParticlesTR) { + fTreeGenV0->GetEntry(entry); + entry++; + if (eventNr < fGenV0Info->fMCm.fEventNr) continue; + if (eventNr > fGenV0Info->fMCm.fEventNr) continue;; + // + fNextV0ToRead = entry-1; + // + // + AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCm.fLabel); + AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCd.fLabel); + if (fGenV0Info->fMCm.fCharge*fGenV0Info->fMCd.fCharge>0) continue; // interactions + if (!fRecInfo1 || !fRecInfo2) continue; + fRecV0Info->fT1 = (*fRecInfo1); + fRecV0Info->fT2 = (*fRecInfo2); + fRecV0Info->fV0Status =0; + if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1; + if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2; + + if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1; + + + if (abs(fRecV0Info->fV0Status)==3){ + fRecV0Info->Update(vertex); + { + // + // TPC V0 Info + Double_t x,alpha, param[5],cov[15]; + if ( fRecV0Info->fT1.fESDTrack.GetInnerParam() && fRecV0Info->fT2.fESDTrack.GetInnerParam()){ + fRecV0Info->fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param); + fRecV0Info->fT1.fESDTrack.GetInnerExternalCovariance(cov); + AliExternalTrackParam paramP(x,alpha,param,cov); + // + fRecV0Info->fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param); + fRecV0Info->fT2.fESDTrack.GetInnerExternalCovariance(cov); + AliExternalTrackParam paramM(x,alpha,param,cov); + // + fRecV0Info->fV0tpc.SetParamN(paramM); + fRecV0Info->fV0tpc.SetParamP(paramP); + Double_t pid1[5],pid2[5]; + fRecV0Info->fT1.fESDTrack.GetESDpid(pid1); + fRecV0Info->fT1.fESDTrack.GetESDpid(pid2); + // + //fRecV0Info->fV0tpc.UpdatePID(pid1,pid2); + fRecV0Info->fV0tpc.Update(vertex); + + // + // + fRecV0Info->fT1.fESDTrack.GetExternalParameters(x,param); + fRecV0Info->fT1.fESDTrack.GetExternalCovariance(cov); + alpha = fRecV0Info->fT1.fESDTrack.GetAlpha(); + new (¶mP) AliExternalTrackParam(x,alpha,param,cov); + // + fRecV0Info->fT2.fESDTrack.GetExternalParameters(x,param); + fRecV0Info->fT2.fESDTrack.GetExternalCovariance(cov); + alpha = fRecV0Info->fT2.fESDTrack.GetAlpha(); + new (¶mM) AliExternalTrackParam(x,alpha,param,cov); + // + fRecV0Info->fV0its.SetParamN(paramM); + fRecV0Info->fV0its.SetParamP(paramP); + // fRecV0Info->fV0its.UpdatePID(pid1,pid2); + fRecV0Info->fV0its.Update(vertex); + } + } + if (TMath::Abs(fGenV0Info->fMCm.fPdg)==11 &&TMath::Abs(fGenV0Info->fMCd.fPdg)==11){ + if (fRecV0Info->fDist2>10){ + fRecV0Info->Update(vertex); + } + if (fRecV0Info->fDist2>10){ + fRecV0Info->Update(vertex); + } + } + } + // + // take the V0 from reconstruction + + Int_t label = TMath::Min(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel); + Int_t label2 = TMath::Max(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel); + AliV0 *v0MI=0; + fRecV0Info->fRecStatus =0; + fRecV0Info->fMultiple = fMultiRecV0[label]; + fRecV0Info->fV0Multiple=0; + // + if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){ + + // for (Int_t j=0;j=0;j--){ + Int_t index = fIndexRecV0[label*20+j]; + if (index<0) continue; + AliV0 *v0MI2 = (AliV0*)fEvent->GetV0(index); +// if (TMath::Abs(v0MI2->fLab[0])==label &&TMath::Abs(v0MI2->fLab[1])==label2) { +// v0MI =v0MI2; +// fRecV0Info->fV0Multiple++; +// fSignedV0[index]=1; +// } +// if (TMath::Abs(v0MI2->fLab[1])==label &&TMath::Abs(v0MI2->fLab[0])==label2) { +// v0MI =v0MI2; +// fRecV0Info->fV0Multiple++; +// fSignedV0[index]=1; +// } + } + } + if (v0MI){ + fRecV0Info->fV0rec = *v0MI; + fRecV0Info->fRecStatus=1; + } + + fTreeCmpV0->Fill(); + } + // + // write fake v0s + + Int_t nV0MIs = fEvent->GetNumberOfV0s(); + for (Int_t i=0;iGetV0(i); + if (!v0MI) continue; + // + fRecV0Info->fV0rec = *v0MI; + fRecV0Info->fV0Status =-10; + fRecV0Info->fRecStatus =-2; + // + // AliESDRecInfo* fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[0])); +// AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[1])); +// if (fRecInfo1 && fRecInfo2){ +// fRecV0Info->fT1 = (*fRecInfo1); +// fRecV0Info->fT2 = (*fRecInfo2); +// fRecV0Info->fRecStatus =-1; +// } + fRecV0Info->Update(vertex); + fTreeCmpV0->Fill(); + } + } + + + + fTreeCmpV0->AutoSave(); + printf("Time spended in BuilV0Info Loop\n"); + timer.Print(); + if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<Get("ESDcmpTracks"); + if (!fTree) { + printf("no track comparison tree found\n"); + file->Close(); + delete file; + } +} + + diff --git a/PWG1/AliESDComparisonMI.h b/PWG1/AliESDComparisonMI.h new file mode 100644 index 00000000000..bf8f7971b4d --- /dev/null +++ b/PWG1/AliESDComparisonMI.h @@ -0,0 +1,233 @@ + +class AliESDComparisonDraw: public AliComparisonDraw{ +public: + AliESDComparisonDraw(){fTree = 0;} + void SetIO(const char *fname = "cmpESDTracks.root"); + ClassDef(AliESDComparisonDraw,1) +}; +ClassImp(AliESDComparisonDraw) + + +///////////////////////////////////////////////////////////////////////// +class AliESDRecInfo: public TObject { + +public: + AliESDRecInfo(){} + ~AliESDRecInfo(){} + void UpdatePoints(AliESDtrack* track); + void Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed, AliESD *event); + void Reset(); + Float_t fTPCPoints[10]; //start , biggest end points,max density .. density at the last 30 pad-rows + Double_t fTPCinR0[5]; //generated position of the track at inner tpc - radius [3] and fi [4] + Double_t fTPCinR1[5]; //reconstructed postion of the track - radius [3] and fi [ + Double_t fTPCinP0[5]; //generated position of the track at inner tpc + Double_t fTPCinP1[5]; //reconstructed postion of the track + Double_t fTPCAngle0[2]; // generated angle + Double_t fTPCAngle1[2]; //refconstructed angle + Double_t fTPCDelta[5]; // deltas + Double_t fTPCPools[5]; // pools + Double_t fITSinR0[5]; //generated position of the track at inner tpc + Double_t fITSinR1[5]; //reconstructed postion of the track + Double_t fITSinP0[5]; //generated position of the track at inner tpc + Double_t fITSinP1[5]; //reconstructed postion of the track + Double_t fITSAngle0[2]; // generated angle + Double_t fITSAngle1[2]; //refconstructed angle + Double_t fITSDelta[5]; // deltas + Double_t fITSPools[5]; // pools + AliESDtrack fESDTrack; // tpc track + AliESDfriendTrack fTrackF; // friend track + AliITStrackMI fITStrack; //its track + AliTRDtrack fTRDtrack; //its track + Float_t fBestTOFmatch; //best matching between times + Float_t fTRLocalCoord[3]; //local coordinates of the track ref. + Int_t fReconstructed; //flag if track was reconstructed + Int_t fFake; // fake track + Int_t fMultiple; // number of reconstructions + Bool_t fTPCOn; // TPC refitted inward + Int_t fStatus[4]; // status -0 not found - 1 -only in - 2 -in-out -3 -in -out-refit + Bool_t fITSOn; // ITS refitted inward + Bool_t fTRDOn; // ITS refitted inward + Float_t fDeltaP; //delta of momenta + Double_t fSign; // sign + Int_t fLabels[2]; // labels + ClassDef(AliESDRecInfo,2) // container for +}; +ClassImp(AliESDRecInfo) + +void AliESDRecInfo::Reset() +{ + fMultiple =0; + fFake =0; + fReconstructed=0; +} + + +///////////////////////////////////////////////////////// +///////////////////////////////////////////////////////// +///////////////////////////////////////////////////////// + + +class AliESDRecV0Info: public TObject { +public: + void Update(Float_t vertex[3]); + AliESDRecInfo fT1; //track1 + AliESDRecInfo fT2; //track2 + Double_t fDist1; //info about closest distance according closest MC - linear DCA + Double_t fDist2; //info about closest distance parabolic DCA + Double_t fInvMass; //reconstructed invariant mass - + // + Double_t fPdr[3]; //momentum at vertex daughter - according approx at DCA + Double_t fXr[3]; //rec. position according helix + // + Double_t fRs[2]; // minimum radius in rphi intersection + Double_t fDistMinR; // distance at minimal radius + Double_t fPm[3]; //momentum at the vertex mother + Double_t fAngle[3]; //three angles + Double_t fRr; // rec position of the vertex + Int_t fLab[2]; //MC label of the partecle + Float_t fPointAngleFi; //point angle fi + Float_t fPointAngleTh; //point angle theta + Float_t fPointAngle; //point angle full + Int_t fV0Status; // status of the kink + AliV0 fV0tpc; // Vo information from reconsturction according TPC + AliV0 fV0its; // Vo information from reconsturction according ITS + AliV0 fV0rec; // V0 information form the reconstruction + Int_t fMultiple; + Int_t fV0Multiple; + Int_t fRecStatus; // status form the reconstuction + ClassDef(AliESDRecV0Info,2) // container for +}; + +ClassImp(AliESDRecV0Info) + + +class AliESDRecKinkInfo: public TObject { +public: + void Update(); + AliESDRecInfo fT1; //track1 + AliESDRecInfo fT2; //track2 + AliESDkink fKink; //kink + Double_t fDist1; //info about closest distance according closest MC - linear DCA + Double_t fDist2; //info about closest distance parabolic DCA + Double_t fInvMass; //reconstructed invariant mass - + // + Double_t fPdr[3]; //momentum at vertex daughter - according approx at DCA + Double_t fXr[3]; //rec. position according helix + // + Double_t fPm[3]; //momentum at the vertex mother + Double_t fAngle[3]; //three angles + Double_t fRr; // rec position of the vertex + Double_t fMinR; // minimum radius in rphi intersection + Double_t fDistMinR; // distance at minimal radius + Int_t fLab[2]; //MC label of the partecle + Float_t fPointAngleFi; //point angle fi + Float_t fPointAngleTh; //point angle theta + Float_t fPointAngle; //point angle full + Int_t fStatus; //status -tracks + Int_t fRecStatus; //kink -status- 0 - not found 1-good - fake + Int_t fMultiple; + Int_t fKinkMultiple; + ClassDef(AliESDRecKinkInfo,1) // container for +}; + +ClassImp(AliESDRecKinkInfo) + + + +//////////////////////////////////////////////////////////////////////// +// +// Start of implementation of the class ESDCmpTr +// +//////////////////////////////////////////////////////////////////////// + +class ESDCmpTr { + +public: + ESDCmpTr(); + ESDCmpTr(const char* fnGenTracks, + const char* fnCmpRes ="cmpTracks.root", + const char* fnGalice ="galice.root", Int_t direction=0, + Int_t nEvents=1, Int_t firstEvent=0); + virtual ~ESDCmpTr(); + void Reset(); + Int_t Exec(); + Int_t Exec(Int_t nEvents, Int_t firstEventNr); + Int_t SetIO(); + Int_t SetIO(Int_t eventNr ); + void CreateTreeCmp(); + void CloseOutputFile(); + Bool_t ConnectGenTree(); + Int_t TreeGenLoop(Int_t eventNr); + Int_t TreeTLoop(); + Int_t BuildKinkInfo0(Int_t eventNr); // build kink info 0 + Int_t BuildV0Info(Int_t eventNr); // build kink info 0 + void SetFirstEventNr(Int_t i) {fFirstEventNr = i;} + void SetNEvents(Int_t i) {fNEvents = i;} + void SetDebug(Int_t level) {fDebug = level;} + +// tmp method, should go to TrackReferenceESD + static TVector3 TR2Local(AliTrackReference *trackRef, + AliTPCParam *paramTPC); + +private: + + Int_t fEventNr; //! current event number + Int_t fNEvents; //! number of events to process + Int_t fFirstEventNr; //! first event to process + // + char fFnCmp[1000]; //! output file name with cmp tracks + TFile *fFileCmp; //! output file with cmp tracks + TTree *fTreeCmp; //! output tree with cmp tracks + TTree *fTreeCmpKinks; //! output tree with cmp Kinks + TTree *fTreeCmpV0; //! output tree with cmp V0 + // + char fFnGenTracks[1000]; //! input file name with gen tracks + TFile *fFileGenTracks; + TTree *fTreeGenTracks; + TTree *fTreeGenKinks; // tree with gen kinks + TTree *fTreeGenV0; // tree with gen V0 + // + // + Int_t fDirection; + // + AliRunLoader * fLoader; //! pointer to the run loader + //TTree *fTreeRecTracks; //! tree with reconstructed tracks + // + Short_t *fIndexRecTracks; //! index of particle label in the TreeT_ESD + Short_t *fFakeRecTracks; //! number of fake tracks + Short_t *fMultiRecTracks; //! number of multiple reconstructions + // + Short_t *fIndexRecKinks; //! index of particle label in treeesd + Short_t *fMultiRecKinks; //! number of multiple reconstructions + Short_t *fSignedKinks; //! indicator that kink was not fake + // + Short_t *fIndexRecV0; //! index of particle label in treeesd + Short_t *fMultiRecV0; //! number of multiple reconstructions + Short_t *fSignedV0; //! indicator that kink was not fake + // + TObjArray *fRecArray; // container with rec infos + AliESD *fEvent; //!event + AliESDfriend *fESDfriend; //!event friend + // + AliTPCParam* fParamTPC; //! AliTPCParam + Int_t fNParticles; //! number of particles in the input tree genTracks + Int_t fDebug; //! debug flag + Int_t fNextTreeGenEntryToRead; //! last entry already read from genTracks tree + Int_t fNextKinkToRead; //! last entry already read from genKinks tree + Int_t fNextV0ToRead; //! last entry already read from genV0 tree + // + AliMCInfo* fMCInfo; //! MC information writen per particle + AliGenKinkInfo* fGenKinkInfo; //! MC information writen per Kink + AliGenV0Info* fGenV0Info; //! MC information writen per Kink + AliESDRecInfo* fRecInfo; //! Rec. information writen per particle + AliESDfriendTrack* fFriend; //! friend track + AliESDRecKinkInfo* fRecKinkInfo; //! reconstructed kink info + AliESDRecV0Info* fRecV0Info; //! reconstructed kink info + // + + ClassDef(ESDCmpTr,1) // class which creates and fills tree with ESDGenTrack objects +}; +ClassImp(ESDCmpTr) + + + diff --git a/PWG1/AliGenInfo.cxx b/PWG1/AliGenInfo.cxx new file mode 100644 index 00000000000..0214cf07db9 --- /dev/null +++ b/PWG1/AliGenInfo.cxx @@ -0,0 +1,1360 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + + +/////////////////////////////////////////////////////////////////////////// +/* + +Origin: marian.ivanov@cern.ch + +Generate complex MC information - used for Comparison later on +How to use it? + +.L +AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0) +t->Exec(); + +*/ + +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +//ROOT includes +#include "TROOT.h" +#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 "TH1F.h" +#include "TH2F.h" +#include "TCanvas.h" +#include "TPad.h" +#include "TF1.h" +#include "TView.h" +#include "TView3D.h" +#include "TGeometry.h" +#include "TPolyLine3D.h" + +//ALIROOT includes +#include "AliRun.h" +#include "AliStack.h" +#include "AliSimDigits.h" +#include "AliTPCParam.h" +#include "AliTPC.h" +#include "AliTPCLoader.h" +#include "AliDetector.h" +#include "AliTrackReference.h" +#include "AliTPCParamSR.h" +#include "AliTracker.h" +#include "AliMagF.h" +#include "AliHelix.h" +#include "AliPoints.h" +#include "AliTrackPointArray.h" + +#endif +#include "AliGenInfo.h" +// +// + +ClassImp(AliTPCdigitRow); +ClassImp(AliMCInfo); +ClassImp(AliGenV0Info) +ClassImp(AliGenKinkInfo) +ClassImp(AliGenInfoMaker) + + + +AliTPCParam * GetTPCParam(){ + AliTPCParamSR * par = new AliTPCParamSR; + par->Update(); + return par; +} + + +//_____________________________________________________________________________ +Float_t TPCBetheBloch(Float_t bg) +{ + // + // Bethe-Bloch energy loss formula + // + const Double_t kp1=0.76176e-1; + const Double_t kp2=10.632; + const Double_t kp3=0.13279e-4; + const Double_t kp4=1.8631; + const Double_t kp5=1.9479; + + Double_t dbg = (Double_t) bg; + + Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg); + + Double_t aa = TMath::Power(beta,kp4); + Double_t bb = TMath::Power(1./dbg,kp5); + + bb=TMath::Log(kp3+bb); + + return ((Float_t)((kp2-aa-bb)*kp1/aa)); +} + + + + + +//////////////////////////////////////////////////////////////////////// +AliMCInfo::AliMCInfo() +{ + fTPCReferences = new TClonesArray("AliTrackReference",10); + fITSReferences = new TClonesArray("AliTrackReference",10); + fTRDReferences = new TClonesArray("AliTrackReference",10); + fTOFReferences = new TClonesArray("AliTrackReference",10); + fTRdecay.SetTrack(-1); + fCharge = 0; +} + +AliMCInfo::~AliMCInfo() +{ + if (fTPCReferences) { + delete fTPCReferences; + } + if (fITSReferences){ + delete fITSReferences; + } + if (fTRDReferences){ + delete fTRDReferences; + } + if (fTOFReferences){ + delete fTOFReferences; + } + +} + + + +void AliMCInfo::Update() +{ + // + // + fMCtracks =1; + if (!fTPCReferences) { + fNTPCRef =0; + return; + } + Float_t direction=1; + //Float_t rlast=0; + fNTPCRef = fTPCReferences->GetEntriesFast(); + fNITSRef = fITSReferences->GetEntriesFast(); + fNTRDRef = fTRDReferences->GetEntriesFast(); + fNTOFRef = fTOFReferences->GetEntriesFast(); + + for (Int_t iref =0;irefGetEntriesFast();iref++){ + AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref); + //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y()); + Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside + if (iref==0) direction = newdirection; + if ( newdirection*direction<0){ + //changed direction + direction = newdirection; + fMCtracks+=1; + } + //rlast=r; + } + // + // decay info + fTPCdecay=kFALSE; + if (fTRdecay.GetTrack()>0){ + fDecayCoord[0] = fTRdecay.X(); + fDecayCoord[1] = fTRdecay.Y(); + fDecayCoord[2] = fTRdecay.Z(); + if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){ + fTPCdecay=kTRUE; + } + else{ + fDecayCoord[0] = 0; + fDecayCoord[1] = 0; + fDecayCoord[2] = 0; + } + } + // + // + //digits information update + fRowsWithDigits = fTPCRow.RowsOn(); + fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows + fRowsTrackLength = fTPCRow.Last() - fTPCRow.First(); + // + // + // calculate primary ionization per cm + if (fParticle.GetPDG()){ + fMass = fParticle.GetMass(); + fCharge = fParticle.GetPDG()->Charge(); + if (fTPCReferences->GetEntriesFast()>0){ + fTrackRef = *((AliTrackReference*)fTPCReferences->At(0)); + } + if (fMass>0){ + Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+ + fTrackRef.Py()*fTrackRef.Py()+ + fTrackRef.Pz()*fTrackRef.Pz()); + if (p>0.001){ + Float_t betagama = p /fMass; + fPrim = TPCBetheBloch(betagama); + }else fPrim=0; + } + }else{ + fMass =0; + fPrim =0; + } +} + +///////////////////////////////////////////////////////////////////////////////// +/* +void AliGenV0Info::Update() +{ + fMCPd[0] = fMCd.fParticle.Px(); + fMCPd[1] = fMCd.fParticle.Py(); + fMCPd[2] = fMCd.fParticle.Pz(); + fMCPd[3] = fMCd.fParticle.P(); + // + fMCX[0] = fMCd.fParticle.Vx(); + fMCX[1] = fMCd.fParticle.Vy(); + fMCX[2] = fMCd.fParticle.Vz(); + fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]); + // + fPdg[0] = fMCd.fParticle.GetPdgCode(); + fPdg[1] = fMCm.fParticle.GetPdgCode(); + // + fLab[0] = fMCd.fParticle.GetUniqueID(); + fLab[1] = fMCm.fParticle.GetUniqueID(); + // +} +*/ + +void AliGenV0Info::Update(Float_t vertex[3]) +{ + fMCPd[0] = fMCd.fParticle.Px(); + fMCPd[1] = fMCd.fParticle.Py(); + fMCPd[2] = fMCd.fParticle.Pz(); + fMCPd[3] = fMCd.fParticle.P(); + // + fMCX[0] = fMCd.fParticle.Vx(); + fMCX[1] = fMCd.fParticle.Vy(); + fMCX[2] = fMCd.fParticle.Vz(); + fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]); + // + fPdg[0] = fMCd.fParticle.GetPdgCode(); + fPdg[1] = fMCm.fParticle.GetPdgCode(); + // + fLab[0] = fMCd.fParticle.GetUniqueID(); + fLab[1] = fMCm.fParticle.GetUniqueID(); + // + // + // + Double_t x1[3],p1[3]; + TParticle & pdaughter = fMCd.fParticle; + x1[0] = pdaughter.Vx(); + x1[1] = pdaughter.Vy(); + x1[2] = pdaughter.Vz(); + p1[0] = pdaughter.Px(); + p1[1] = pdaughter.Py(); + p1[2] = pdaughter.Pz(); + Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1; + AliHelix dhelix1(x1,p1,sign); + // + // + Double_t x2[3],p2[3]; + // + TParticle & pmother = fMCm.fParticle; + x2[0] = pmother.Vx(); + x2[1] = pmother.Vy(); + x2[2] = pmother.Vz(); + p2[0] = pmother.Px(); + p2[1] = pmother.Py(); + p2[2] = pmother.Pz(); + // + // + sign = (pmother.GetPDG()->Charge()>0) ? -1:1; + AliHelix mhelix(x2,p2,sign); + + // + // + // + //find intersection linear + // + Double_t distance1, distance2; + Double_t phase[2][2],radius[2]; + Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius); + Double_t delta1=10000,delta2=10000; + if (points>0){ + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + else{ + fInvMass=-1; + return; + } + if (points==2){ + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + distance1 = TMath::Min(delta1,delta2); + // + //find intersection parabolic + // + points = dhelix1.GetRPHIintersections(mhelix, phase, radius); + delta1=10000,delta2=10000; + + if (points>0){ + dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + if (points==2){ + dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + + distance2 = TMath::Min(delta1,delta2); + // + if (delta10) fInvMass = TMath::Sqrt(fInvMass); + + +} + + + +///////////////////////////////////////////////////////////////////////////////// +void AliGenKinkInfo::Update() +{ + fMCPd[0] = fMCd.fParticle.Px(); + fMCPd[1] = fMCd.fParticle.Py(); + fMCPd[2] = fMCd.fParticle.Pz(); + fMCPd[3] = fMCd.fParticle.P(); + // + fMCX[0] = fMCd.fParticle.Vx(); + fMCX[1] = fMCd.fParticle.Vy(); + fMCX[2] = fMCd.fParticle.Vz(); + fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]); + // + fPdg[0] = fMCd.fParticle.GetPdgCode(); + fPdg[1] = fMCm.fParticle.GetPdgCode(); + // + fLab[0] = fMCd.fParticle.GetUniqueID(); + fLab[1] = fMCm.fParticle.GetUniqueID(); + // + // + // + Double_t x1[3],p1[3]; + TParticle & pdaughter = fMCd.fParticle; + x1[0] = pdaughter.Vx(); + x1[1] = pdaughter.Vy(); + x1[2] = pdaughter.Vz(); + p1[0] = pdaughter.Px(); + p1[1] = pdaughter.Py(); + p1[2] = pdaughter.Pz(); + Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1; + AliHelix dhelix1(x1,p1,sign); + // + // + Double_t x2[3],p2[3]; + // + TParticle & pmother = fMCm.fParticle; + x2[0] = pmother.Vx(); + x2[1] = pmother.Vy(); + x2[2] = pmother.Vz(); + p2[0] = pmother.Px(); + p2[1] = pmother.Py(); + p2[2] = pmother.Pz(); + // + AliTrackReference & pdecay = fMCm.fTRdecay; + x2[0] = pdecay.X(); + x2[1] = pdecay.Y(); + x2[2] = pdecay.Z(); + p2[0] = pdecay.Px(); + p2[1] = pdecay.Py(); + p2[2] = pdecay.Pz(); + // + sign = (pmother.GetPDG()->Charge()>0) ? -1:1; + AliHelix mhelix(x2,p2,sign); + + // + // + // + //find intersection linear + // + Double_t distance1, distance2; + Double_t phase[2][2],radius[2]; + Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius); + Double_t delta1=10000,delta2=10000; + if (points>0){ + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + if (points==2){ + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + distance1 = TMath::Min(delta1,delta2); + // + //find intersection parabolic + // + points = dhelix1.GetRPHIintersections(mhelix, phase, radius); + delta1=10000,delta2=10000; + + if (points>0){ + dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + if (points==2){ + dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + + distance2 = TMath::Min(delta1,delta2); + // + if (delta1= 8*kgRowBytes) { + cerr<<"AliTPCdigitRow::SetRow: index "< upto) return total; + if (TESTBIT(fDig[i],j)) total++; + } + } + return total; +} +//////////////////////////////////////////////////////////////////////// +void AliTPCdigitRow::Reset() +{ +// +// resets all rows to zero +// + for (Int_t i = 0; i=0; i--) { + for (Int_t j = 7; j >= 0; j--) { + if TESTBIT(fDig[i],j) return i*8+j; + } + } + return -1; +} +//////////////////////////////////////////////////////////////////////// +Int_t AliTPCdigitRow::First() +{ +// +// returns the first row number with a digit +// returns -1 if now digits +// + for (Int_t i = 0; iGetRunLoader(); + delete gAlice; + gAlice = 0x0; + } + if (fLoader->LoadgAlice()){ + cerr<<"Error occured while l"<GetNumberOfEvents(); + if (nEvents==0) { + nEvents =nall; + fNEvents=nall; + fFirstEventNr=0; + } + + if (nall<=0){ + cerr<<"no events available"<nall) { + fEventNr = nall-firstEvent; + cerr<<"restricted number of events availaible"<Field(); + AliTracker::SetFieldMap(magf,0); +} + + +AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i) +{ + // + if (iUnloadgAlice(); + gAlice = 0; + delete fLoader; + } +} + +Int_t AliGenInfoMaker::SetIO() +{ + // + // + CreateTreeGenTracks(); + if (!fTreeGenTracks) return 1; + // AliTracker::SetFieldFactor(); + + fParamTPC = GetTPCParam(); + // + return 0; +} + +//////////////////////////////////////////////////////////////////////// +Int_t AliGenInfoMaker::SetIO(Int_t eventNr) +{ + // + // + // SET INPUT + fLoader->SetEventNumber(eventNr); + // + fLoader->LoadHeader(); + fLoader->LoadKinematics(); + fStack = fLoader->Stack(); + // + fLoader->LoadTrackRefs(); + fLoader->LoadHits(); + fTreeTR = fLoader->TreeTR(); + // + AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader"); + tpcl->LoadDigits(); + fTreeD = tpcl->TreeD(); + return 0; +} + +Int_t AliGenInfoMaker::CloseIOEvent() +{ + fLoader->UnloadHeader(); + fLoader->UnloadKinematics(); + fLoader->UnloadTrackRefs(); + AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader"); + tpcl->UnloadDigits(); + return 0; +} + +Int_t AliGenInfoMaker::CloseIO() +{ + fLoader->UnloadgAlice(); + return 0; +} + + + +//////////////////////////////////////////////////////////////////////// +Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr) +{ + fNEvents = nEvents; + fFirstEventNr = firstEventNr; + return Exec(); +} + +//////////////////////////////////////////////////////////////////////// +Int_t AliGenInfoMaker::Exec() +{ + TStopwatch timer; + timer.Start(); + Int_t status =SetIO(); + if (status>0) return status; + // + + for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents; + fEventNr++) { + SetIO(fEventNr); + fNParticles = fStack->GetNtrack(); + // + fGenInfo = new AliMCInfo*[fNParticles]; + for (UInt_t i = 0; i0) return 1; + // + if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<0) return 1; + // + if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<0) return 1; + if (fDebug>2) cout<<"\tEnd loop over TreeK"<0) return 1; + if (BuildV0Info()>0) return 1; + //if (BuildHitLines()>0) return 1; + if (fDebug>2) cout<<"\tEnd loop over TreeK"<Branch("MC","AliMCInfo",&info); + delete info; + // + AliGenKinkInfo *kinkinfo = new AliGenKinkInfo; + fTreeKinks = new TTree("genKinksTree","genKinksTree"); + fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo); + delete kinkinfo; + // + AliGenV0Info *v0info = new AliGenV0Info; + fTreeV0 = new TTree("genV0Tree","genV0Tree"); + fTreeV0->Branch("MC","AliGenV0Info",&v0info); + delete v0info; + // + // + AliTrackPointArray * points0 = new AliTrackPointArray; + AliTrackPointArray * points1 = new AliTrackPointArray; + AliTrackPointArray * points2 = new AliTrackPointArray; + fTreeHitLines = new TTree("HitLines","HitLines"); + fTreeHitLines->Branch("TPC.","AliTrackPointArray",&points0,32000,10); + fTreeHitLines->Branch("TRD.","AliTrackPointArray",&points1,32000,10); + fTreeHitLines->Branch("ITS.","AliTrackPointArray",&points2,32000,10); + Int_t label=0; + fTreeHitLines->Branch("Label",&label,"label/I"); + // + fTreeGenTracks->AutoSave(); + fTreeKinks->AutoSave(); + fTreeV0->AutoSave(); + fTreeHitLines->AutoSave(); +} +//////////////////////////////////////////////////////////////////////// +void AliGenInfoMaker::CloseOutputFile() +{ + if (!fFileGenTracks) { + cerr<<"File "<cd(); + fTreeGenTracks->Write(); + delete fTreeGenTracks; + fTreeKinks->Write(); + delete fTreeKinks; + fTreeV0->Write(); + delete fTreeV0; + + fFileGenTracks->Close(); + delete fFileGenTracks; + return; +} + +//////////////////////////////////////////////////////////////////////// +Int_t AliGenInfoMaker::TreeKLoop() +{ +// +// open the file with treeK +// loop over all entries there and save information about some tracks +// + + AliStack * stack = fStack; + if (!stack) {cerr<<"Stack was not found!\n"; return 1;} + + if (fDebug > 0) { + cout<<"There are "<GetBranch("MC"); + TParticle *particle = stack->ParticleFromTreeK(0); + fVPrim[0] = particle->Vx(); + fVPrim[1] = particle->Vy(); + fVPrim[2] = particle->Vz(); + for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) { + // load only particles with TR + AliMCInfo * info = GetInfo(iParticle); + if (!info) continue; + ////////////////////////////////////////////////////////////////////// + info->fLabel = iParticle; + // + info->fParticle = *(stack->Particle(iParticle)); + info->fVDist[0] = info->fParticle.Vx()-fVPrim[0]; + info->fVDist[1] = info->fParticle.Vy()-fVPrim[1]; + info->fVDist[2] = info->fParticle.Vz()-fVPrim[2]; + info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+ + info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]); + // + // + ipdg = info->fParticle.GetPdgCode(); + info->fPdg = ipdg; + ppdg = pdg.GetParticle(ipdg); + info->fEventNr = fEventNr; + info->Update(); + ////////////////////////////////////////////////////////////////////// + br->SetAddress(&info); + fTreeGenTracks->Fill(); + } + fTreeGenTracks->AutoSave(); + if (fDebug > 2) cerr<<"end of TreeKLoop"< 0) { + cout<<"There are "<GetBranch("MC"); + // + AliGenKinkInfo * kinkinfo = new AliGenKinkInfo; + // + for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) { + // load only particles with TR + AliMCInfo * info = GetInfo(iParticle); + if (!info) continue; + if (info->fCharge==0) continue; + if (info->fTRdecay.P()<0.13) continue; //momenta cut + if (info->fTRdecay.R()>500) continue; //R cut - decay outside barrel + TParticle & particle = info->fParticle; + Int_t first = particle.GetDaughter(0); + if (first ==0) continue; + + Int_t last = particle.GetDaughter(1); + if (last ==0) last=first; + AliMCInfo * info2 = 0; + AliMCInfo * dinfo = 0; + + for (Int_t id2=first;id2<=last;id2++){ + info2 = GetInfo(id2); + if (!info2) continue; + TParticle &p2 = info2->fParticle; + Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy()); + if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue; + if (!(p2.GetPDG())) continue; + dinfo =info2; + + kinkinfo->fMCm = (*info); + kinkinfo->fMCd = (*dinfo); + if (kinkinfo->fMCm.fParticle.GetPDG()==0 || kinkinfo->fMCd.fParticle.GetPDG()==0) continue; + kinkinfo->Update(); + br->SetAddress(&kinkinfo); + fTreeKinks->Fill(); + } + /* + if (dinfo){ + kinkinfo->fMCm = (*info); + kinkinfo->fMCd = (*dinfo); + kinkinfo->Update(); + br->SetAddress(&kinkinfo); + fTreeKinks->Fill(); + } + */ + } + fTreeGenTracks->AutoSave(); + if (fDebug > 2) cerr<<"end of Kink Loop"< 0) { + cout<<"There are "<GetBranch("MC"); + // + AliGenV0Info * v0info = new AliGenV0Info; + // + // + for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) { + // load only particles with TR + AliMCInfo * info = GetInfo(iParticle); + if (!info) continue; + if (info->fCharge==0) continue; + // + // + TParticle & particle = info->fParticle; + Int_t mother = particle.GetMother(0); + if (mother <=0) continue; + // + TParticle * motherparticle = stack->Particle(mother); + if (!motherparticle) continue; + // + Int_t last = motherparticle->GetDaughter(1); + if (last==(int)iParticle) continue; + AliMCInfo * info2 = info; + AliMCInfo * dinfo = GetInfo(last); + if (!dinfo) continue; + if (!dinfo->fParticle.GetPDG()) continue; + if (!info2->fParticle.GetPDG()) continue; + if (dinfo){ + v0info->fMCm = (*info); + v0info->fMCd = (*dinfo); + v0info->fMotherP = (*motherparticle); + v0info->Update(fVPrim); + br->SetAddress(&v0info); + fTreeV0->Fill(); + } + } + fTreeV0->AutoSave(); + if (fDebug > 2) cerr<<"end of V0 Loop"< 0) { + cout<<"There are "<GetBranch("TPC."); +// TBranch * brtrd = fTreeHitLines->GetBranch("TRD."); +// TBranch * brits = fTreeHitLines->GetBranch("ITS."); +// TBranch * brlabel = fTreeHitLines->GetBranch("Label"); +// brlabel->SetAddress(&label); +// brtpc->SetAddress(&tpcp); +// brtrd->SetAddress(&trdp); +// brits->SetAddress(&itsp); +// // +// AliDetector *dtpc = gAlice->GetDetector("TPC"); +// AliDetector *dtrd = gAlice->GetDetector("TRD"); +// AliDetector *dits = gAlice->GetDetector("ITS"); + +// for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) { +// // load only particles with TR +// AliMCInfo * info = GetInfo(iParticle); +// if (!info) continue; +// Int_t primpart = info->fPrimPart; +// ipdg = info->fParticle.GetPdgCode(); +// label = iParticle; +// // +// gAlice->ResetHits(); +// tpcp->Reset(); +// itsp->Reset(); +// trdp->Reset(); +// tpcp->fLabel1 = ipdg; +// trdp->fLabel1 = ipdg; +// itsp->fLabel1 = ipdg; +// if (dtpc->TreeH()->GetEvent(primpart)){ +// dtpc->LoadPoints(primpart); +// tpcp->Reset(dtpc,iParticle); +// } +// if (dtrd->TreeH()->GetEvent(primpart)){ +// dtrd->LoadPoints(primpart); +// trdp->Reset(dtrd,iParticle); +// } +// if (dits->TreeH()->GetEvent(primpart)){ +// dits->LoadPoints(primpart); +// itsp->Reset(dits,iParticle); +// } +// // +// fTreeHitLines->Fill(); +// dtpc->ResetPoints(); +// dtrd->ResetPoints(); +// dits->ResetPoints(); +// } +// delete tpcp; +// delete trdp; +// delete itsp; +// fTreeHitLines->AutoSave(); +// if (fDebug > 2) cerr<<"end of TreeKLoop"<GetNInnerSector(); + Int_t rowShift = 0; + Int_t zero=fParamTPC->GetZeroSup()+6; + // + // + AliSimDigits digitsAddress, *digits=&digitsAddress; + fTreeD->GetBranch("Segment")->SetAddress(&digits); + + Int_t sectorsByRows=(Int_t)fTreeD->GetEntries(); + if (fDebug > 1) cout<<"\tsectorsByRows = "<GetEvent(i)) continue; + Int_t sec,row; + fParamTPC->AdjustSectorRow(digits->GetID(),sec,row); + if (fDebug > 1) cout< nInnerSector) rowShift = fParamTPC->GetNRowLow(); + // + digits->ExpandTrackBuffer(); + digits->First(); + do { + Int_t iRow=digits->CurrentRow(); + Int_t iColumn=digits->CurrentColumn(); + Short_t digitValue = digits->CurrentDigit(); + if (digitValue >= zero) { + Int_t label; + for (Int_t j = 0; j<3; j++) { + // label = digits->GetTrackID(iRow,iColumn,j); + label = digits->GetTrackIDFast(iRow,iColumn,j)-2; + if (label >= (Int_t)fNParticles) { //don't label from bakground event + continue; + } + if (label >= 0 && label <= (Int_t)fNParticles) { + if (fDebug > 6 ) { + cout<<"Inner loop: sector, iRow, iColumn, label, value, row " + <fTPCRow.SetRow(row+rowShift); + } + } + } + } + } while (digits->Next()); + } + + if (fDebug > 2) cerr<<"end of TreeDLoop"<GetEntries(); + if (fDebug > 1) cout<<"There are "<GetBranch("TPC")) treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR); + if (treeTR->GetBranch("ITS")) treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR); + if (treeTR->GetBranch("TRD")) treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR); + if (treeTR->GetBranch("TOF")) treeTR->GetBranch("TOF")->SetAddress(&TOFArrayTR); + if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR); + // + // + // + for (Int_t iPrimPart = 0; iPrimPartGetEntry(iPrimPart); + // + // Loop over TPC references + // + for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) { + AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef); + // + if (trackRef->TestBit(BIT(2))){ + //if decay + if (trackRef->P()GetTrack(); + AliMCInfo * info = GetInfo(label); + if (!info) info = MakeInfo(label); + info->fTRdecay = *trackRef; + } + // + if (trackRef->P()GetTrack(); + AliMCInfo * info = GetInfo(label); + if (!info) info = MakeInfo(label); + if (!info) continue; + info->fPrimPart = iPrimPart; + TClonesArray & arr = *(info->fTPCReferences); + new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef); + } + // + // Loop over ITS references + // + for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) { + AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef); + // + // + if (trackRef->P()GetTrack(); + AliMCInfo * info = GetInfo(label); + if ( (!info) && trackRef->Pt()fPrimPart = iPrimPart; + TClonesArray & arr = *(info->fITSReferences); + new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef); + } + // + // Loop over TRD references + // + for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) { + AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef); + // + if (trackRef->P()GetTrack(); + AliMCInfo * info = GetInfo(label); + if ( (!info) && trackRef->Pt()fPrimPart = iPrimPart; + TClonesArray & arr = *(info->fTRDReferences); + new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef); + } + // + // Loop over TOF references + // + for (Int_t iTrackRef = 0; iTrackRef < TOFArrayTR->GetEntriesFast(); iTrackRef++) { + AliTrackReference *trackRef = (AliTrackReference*)TOFArrayTR->At(iTrackRef); + Int_t label = trackRef->GetTrack(); + AliMCInfo * info = GetInfo(label); + if (!info){ + if (trackRef->Pt()Pt()fPrimPart = iPrimPart; + TClonesArray & arr = *(info->fTOFReferences); + new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef); + } + // + // get dacay position + for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) { + AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef); + // + Int_t label = trackRef->GetTrack(); + AliMCInfo * info = GetInfo(label); + if (!info) continue; + if (!trackRef->TestBit(BIT(2))) continue; //if not decay + // if (TMath::Abs(trackRef.X()); + info->fTRdecay = *trackRef; + } + } + // + TPCArrayTR->Delete(); + delete TPCArrayTR; + TRDArrayTR->Delete(); + delete TRDArrayTR; + TOFArrayTR->Delete(); + delete TOFArrayTR; + + ITSArrayTR->Delete(); + delete ITSArrayTR; + RunArrayTR->Delete(); + delete RunArrayTR; + // + return 0; +} + +//////////////////////////////////////////////////////////////////////// +Float_t AliGenInfoMaker::TR2LocalX(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 x[0]; +} +//////////////////////////////////////////////////////////////////////// + + + diff --git a/PWG1/AliGenInfo.h b/PWG1/AliGenInfo.h new file mode 100644 index 00000000000..70d9dbcdc9b --- /dev/null +++ b/PWG1/AliGenInfo.h @@ -0,0 +1,237 @@ +#ifndef ALIGENINFO_H +#define ALIGENINFO_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + + + +////////////////////////////////////////////////////////////////////////////// +// Class AliGenInfo // +// collect together MC info for comparison purposes - effieciency studies and so on// // +// marian.ivanov@cern.ch // +////////////////////////////////////////////////////////////////////////////// + + + +//////////////////////////////////////////////////////////////////////// +// +// Start of implementation of the class AliTPCdigitRow +// +//////////////////////////////////////////////////////////////////////// + +#include +#include "AliTrackReference.h" + +class TFile; +class AliRunLoader; +class AliStack; +class AliTPCParam; + +const Int_t kgRowBytes = 32; + +class AliTPCdigitRow: public TObject { +public: + AliTPCdigitRow(); + virtual ~AliTPCdigitRow(){;} + void SetRow(Int_t row); + Bool_t TestRow(Int_t row); + AliTPCdigitRow & operator=(const AliTPCdigitRow &digOld); + Int_t RowsOn(Int_t upto=8*kgRowBytes); + Int_t Last(); + Int_t First(); + void Reset(); + +//private: + UChar_t fDig[kgRowBytes]; + ClassDef(AliTPCdigitRow,1) // container for digit pattern +}; + + +//////////////////////////////////////////////////////////////////////// +// +// Start of implementation of the class AliMCInfo +// +//////////////////////////////////////////////////////////////////////// + +class AliMCInfo: public TObject { + +public: + AliMCInfo(); + ~AliMCInfo(); + void Update(); + AliTrackReference fTrackRef; // track reference saved in the output tree + AliTrackReference fTrackRefOut; // decay track reference saved in the output tree + AliTrackReference fTRdecay; // track reference at decay point + // + Int_t fPrimPart; // index of primary particle in TreeH + TParticle fParticle; // generated particle + Float_t fMass; // mass of the particle + Float_t fCharge; // + Int_t fLabel; // track label + Int_t fEventNr; // event number + Int_t fMCtracks; // indication of how many times the track is retuturned back + Int_t fPdg; //pdg code + Float_t fDecayCoord[3]; // position of particle decay + Double_t fVDist[4]; //distance of the particle vertex from primary vertex + Bool_t fTPCdecay; //indicates decay in TPC + Int_t fRowsWithDigitsInn; // number of rows with digits in the inner sectors + Int_t fRowsWithDigits; // number of rows with digits in the outer sectors + Int_t fRowsTrackLength; // last - first row with digit + Float_t fPrim; // theoretical dedx in tpc according particle momenta and mass + AliTPCdigitRow fTPCRow; // information about digits row pattern + Int_t fNTPCRef; // tpc references counter + Int_t fNITSRef; // ITS references counter + Int_t fNTRDRef; // TRD references counter + Int_t fNTOFRef; // TOF references counter + TClonesArray * fTPCReferences; //containner with all track references -in the TPC + TClonesArray * fITSReferences; //container with ITS references + TClonesArray * fTRDReferences; //container with TRD references + TClonesArray * fTOFReferences; //container with TRD references + // + ClassDef(AliMCInfo,1) // container for +}; + + + +class AliGenV0Info: public TObject { +public: + AliMCInfo fMCd; //info about daughter particle - second particle for V0 + AliMCInfo fMCm; //info about mother particle - first particle for V0 + TParticle fMotherP; //particle info about mother particle + void Update(Float_t vertex[3]); // put some derived info to special field + Double_t fMCDist1; //info about closest distance according closest MC - linear DCA + Double_t fMCDist2; //info about closest distance parabolic DCA + // + Double_t fMCPdr[3]; //momentum at vertex daughter - according approx at DCA + Double_t fMCPd[4]; //exact momentum from MC info + Double_t fMCX[3]; //exact position of the vertex + Double_t fMCXr[3]; //rec. position according helix + // + Double_t fMCPm[3]; //momentum at the vertex mother + Double_t fMCAngle[3]; //three angels + Double_t fMCRr; // rec position of the vertex + Double_t fMCR; //exact r position of the vertex + Int_t fPdg[2]; //pdg code of mother and daugter particles + Int_t fLab[2]; //MC label of the partecle + // + Double_t fInvMass; //reconstructed invariant mass - + Float_t fPointAngleFi; //point angle fi + Float_t fPointAngleTh; //point angle theta + Float_t fPointAngle; //point angle full + // + ClassDef(AliGenV0Info,1) // container for +}; + + + +class AliGenKinkInfo: public TObject { +public: + AliMCInfo fMCd; //info about daughter particle - second particle for V0 + AliMCInfo fMCm; //info about mother particle - first particle for V0 + void Update(); // put some derived info to special field + Float_t GetQt(); // + Double_t fMCDist1; //info about closest distance according closest MC - linear DCA + Double_t fMCDist2; //info about closest distance parabolic DCA + // + Double_t fMCPdr[3]; //momentum at vertex daughter - according approx at DCA + Double_t fMCPd[4]; //exact momentum from MC info + Double_t fMCX[3]; //exact position of the vertex + Double_t fMCXr[3]; //rec. position according helix + // + Double_t fMCPm[3]; //momentum at the vertex mother + Double_t fMCAngle[3]; //three angels + Double_t fMCRr; // rec position of the vertex + Double_t fMCR; //exact r position of the vertex + Int_t fPdg[2]; //pdg code of mother and daugter particles + Int_t fLab[2]; //MC label of the partecle + ClassDef(AliGenKinkInfo,1) // container for +}; + + + +//////////////////////////////////////////////////////////////////////// +// +// Start of implementation of the class AliGenInfoMaker +// +//////////////////////////////////////////////////////////////////////// + +class AliGenInfoMaker { + +public: + AliGenInfoMaker(); + AliGenInfoMaker(const char * fnGalice, const char* fnRes ="genTracks.root", + Int_t nEvents=1, Int_t firstEvent=0); + virtual ~AliGenInfoMaker(); + void Reset(); + Int_t Exec(); + Int_t Exec(Int_t nEvents, Int_t firstEventNr); + void CreateTreeGenTracks(); + void CloseOutputFile(); + Int_t TreeKLoop(); + Int_t TreeTRLoop(); + Int_t TreeDLoop(); + Int_t BuildKinkInfo(); // build information about MC kinks + Int_t BuildV0Info(); // build information about MC kinks + Int_t BuildHitLines(); // build information about MC kinks + //void FillInfo(Int_t iParticle); + void SetFirstEventNr(Int_t i) {fFirstEventNr = i;} + void SetNEvents(Int_t i) {fNEvents = i;} + void SetDebug(Int_t level) {fDebug = level;} + Int_t SetIO(Int_t eventNr); + Int_t CloseIOEvent(); + Int_t CloseIO(); + Int_t SetIO(); + Float_t TR2LocalX(AliTrackReference *trackRef, + AliTPCParam *paramTPC); + AliMCInfo * GetInfo(UInt_t i){return (i +#include +//ROOT includes +#include "TROOT.h" +#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 "TH1F.h" +#include "TH2F.h" +#include "TCanvas.h" +#include "TPad.h" +#include "TF1.h" +#include "TView.h" +#include "TView3D.h" +#include "TGeometry.h" +#include "TPolyLine3D.h" +#include "TPolyMarker3D.h" + +//ALIROOT includes +#include "AliTrackPointArray.h" +#include "AliTreeDraw.h" + +#endif + +// +// Class for visualization and some statistacal analysis using tree +// To be used in comparisons +// and calib viewers based on tree + + +ClassImp(AliTreeDraw) + + +AliTreeDraw::AliTreeDraw(): + fTree(0), + fRes(0), + fMean(0), + fPoints(0){ + // + // default constructor + // +} + +void AliTreeDraw::ClearHisto(){ + // + // + delete fRes; + delete fMean; + fRes=0; + fMean=0; +} + + + +TH1F * AliTreeDraw::DrawXY(const char * chx, const char *chy, const char* selection, + const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes) +{ + // + Double_t* bins = CreateLogBins(nbins, minx, maxx); + TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy); + char cut[1000]; + sprintf(cut,"%s&&%s",selection,quality); + char expression[1000]; + sprintf(expression,"%s:%s>>hRes2",chy,chx); + fTree->Draw(expression, cut, "groff"); + TH1F* hMean=0; + TH1F* hRes = CreateResHisto(hRes2, &hMean); + AliLabelAxes(hRes, chx, chy); + // + delete hRes2; + delete[] bins; + ClearHisto(); + fRes = hRes; + fMean = hMean; + return hRes; +} + + + +TH1F * AliTreeDraw::DrawLogXY(const char * chx, const char *chy, const char* selection, + const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes) +{ + // + // + // + Double_t* bins = CreateLogBins(nbins, minx, maxx); + TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy); + char cut[1000]; + sprintf(cut,"%s&&%s",selection,quality); + char expression[1000]; + sprintf(expression,"%s:%s>>hRes2",chy,chx); + fTree->Draw(expression, cut, "groff"); + TH1F* hMean=0; + TH1F* hRes = CreateResHisto(hRes2, &hMean); + AliLabelAxes(hRes, chx, chy); + // + delete hRes2; + delete[] bins; + ClearHisto(); + fRes = hRes; + fMean = hMean; + return hRes; +} + +/////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////// +TH1F * AliTreeDraw::Eff(const char *variable, const char* selection, const char * quality, + Int_t nbins, Float_t min, Float_t max) +{ + // + // + TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max); + TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max); + char inputGen[1000]; + sprintf(inputGen,"%s>>hGen", variable); + fTree->Draw(inputGen, selection, "groff"); + char selectionRec[256]; + sprintf(selectionRec, "%s && %s", selection, quality); + char inputRec[1000]; + sprintf(inputRec,"%s>>hRec", variable); + fTree->Draw(inputRec, selectionRec, "groff"); + // + TH1F* hEff = CreateEffHisto(hGen, hRec); + AliLabelAxes(hEff, variable, "#epsilon [%]"); + fRes = hEff; + delete hRec; + delete hGen; + return hEff; +} + + + +/////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////// +TH1F * AliTreeDraw::EffLog(const char *variable, const char* selection, const char * quality, + Int_t nbins, Float_t min, Float_t max) +{ + // + // + Double_t* bins = CreateLogBins(nbins, min, max); + TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins); + TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins); + char inputGen[1000]; + sprintf(inputGen,"%s>>hGen", variable); + fTree->Draw(inputGen, selection, "groff"); + char selectionRec[256]; + sprintf(selectionRec, "%s && %s", selection, quality); + char inputRec[1000]; + sprintf(inputRec,"%s>>hRec", variable); + fTree->Draw(inputRec, selectionRec, "groff"); + // + TH1F* hEff = CreateEffHisto(hGen, hRec); + AliLabelAxes(hEff, variable, "#epsilon [%]"); + fRes = hEff; + delete hRec; + delete hGen; + delete[] bins; + return hEff; +} + + +/////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////// + +Double_t* AliTreeDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax) +{ + Double_t* bins = new Double_t[nBins+1]; + bins[0] = xMin; + Double_t factor = pow(xMax/xMin, 1./nBins); + for (Int_t i = 1; i <= nBins; i++) + bins[i] = factor * bins[i-1]; + return bins; +} + + + + +void AliTreeDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle) +{ + // + histo->GetXaxis()->SetTitle(xAxisTitle); + histo->GetYaxis()->SetTitle(yAxisTitle); +} + + +TH1F* AliTreeDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec) +{ + // + Int_t nBins = hGen->GetNbinsX(); + TH1F* hEff = (TH1F*) hGen->Clone("hEff"); + hEff->SetTitle(""); + hEff->SetStats(kFALSE); + hEff->SetMinimum(0.); + hEff->SetMaximum(110.); + // + for (Int_t iBin = 0; iBin <= nBins; iBin++) { + Double_t nGen = hGen->GetBinContent(iBin); + Double_t nRec = hRec->GetBinContent(iBin); + if (nGen > 0) { + Double_t eff = nRec/nGen; + hEff->SetBinContent(iBin, 100. * eff); + Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen); + // if (error == 0) error = sqrt(0.1/nGen); + // + if (error == 0) error = 0.0001; + hEff->SetBinError(iBin, 100. * error); + } else { + hEff->SetBinContent(iBin, 100. * 0.5); + hEff->SetBinError(iBin, 100. * 0.5); + } + } + return hEff; +} + + + +TH1F* AliTreeDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean, Bool_t drawBinFits, + Bool_t overflowBinFits) +{ + TVirtualPad* currentPad = gPad; + TAxis* axis = hRes2->GetXaxis(); + Int_t nBins = axis->GetNbins(); + TH1F* hRes, *hMean; + if (axis->GetXbins()->GetSize()){ + hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray()); + hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray()); + } + else{ + hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax()); + hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax()); + + } + hRes->SetStats(false); + hRes->SetOption("E"); + hRes->SetMinimum(0.); + // + hMean->SetStats(false); + hMean->SetOption("E"); + + // create the fit function + TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3); + + fitFunc->SetLineWidth(2); + fitFunc->SetFillStyle(0); + // create canvas for fits + TCanvas* canBinFits = NULL; + Int_t nPads = (overflowBinFits) ? nBins+2 : nBins; + Int_t nx = Int_t(sqrt(nPads-1.));// + 1; + Int_t ny = (nPads-1) / nx + 1; + if (drawBinFits) { + canBinFits = (TCanvas*)gROOT->FindObject("canBinFits"); + if (canBinFits) delete canBinFits; + canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700); + canBinFits->Divide(nx, ny); + } + + // loop over x bins and fit projection + Int_t dBin = ((overflowBinFits) ? 1 : 0); + for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) { + if (drawBinFits) canBinFits->cd(bin + dBin); + TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin); + // + if (hBin->GetEntries() > 5) { + fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS()); + hBin->Fit(fitFunc,"s"); + Double_t sigma = TMath::Abs(fitFunc->GetParameter(2)); + + if (sigma > 0.){ + hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2))); + hMean->SetBinContent(bin, fitFunc->GetParameter(1)); + } + else{ + hRes->SetBinContent(bin, 0.); + hMean->SetBinContent(bin,0); + } + hRes->SetBinError(bin, fitFunc->GetParError(2)); + hMean->SetBinError(bin, fitFunc->GetParError(1)); + + // + // + + } else { + hRes->SetBinContent(bin, 0.); + hRes->SetBinError(bin, 0.); + hMean->SetBinContent(bin, 0.); + hMean->SetBinError(bin, 0.); + } + + + if (drawBinFits) { + char name[256]; + if (bin == 0) { + sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin)); + } else if (bin == nBins+1) { + sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle()); + } else { + sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin), + axis->GetTitle(), axis->GetBinUpEdge(bin)); + } + canBinFits->cd(bin + dBin); + hBin->SetTitle(name); + hBin->SetStats(kTRUE); + hBin->DrawCopy("E"); + canBinFits->Update(); + canBinFits->Modified(); + canBinFits->Update(); + } + + delete hBin; + } + + delete fitFunc; + currentPad->cd(); + *phMean = hMean; + return hRes; +} + + + + +void AliTreeDraw::GetPoints3D(const char * label, const char * chpoints, + const char* selection, TTree * tpoints, Int_t color,Float_t rmin){ + // + // load selected points from tree + // + if (!fPoints) fPoints = new TObjArray; + if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label"); + TBranch * br = tpoints->GetBranch(chpoints); + if (!br) return; + AliTrackPointArray * points = new AliTrackPointArray; + br->SetAddress(&points); + + Int_t npoints = fTree->Draw(label,selection); + Float_t xyz[30000]; + rmin*=rmin; + for (Int_t i=0;iGetV1()[i]; + tpoints->GetEntryWithIndex(index,index); + if (points->GetNPoints()<2) continue; + Int_t goodpoints=0; + for (Int_t i=0;iGetNPoints(); i++){ + xyz[goodpoints*3] = points->GetX()[i]; + xyz[goodpoints*3+1] = points->GetY()[i]; + xyz[goodpoints*3+2] = points->GetZ()[i]; + if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++; + } + TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz); + marker->SetMarkerColor(color); + marker->SetMarkerStyle(1); + fPoints->AddLast(marker); + } +} + diff --git a/PWG1/AliTreeDraw.h b/PWG1/AliTreeDraw.h new file mode 100644 index 00000000000..9ea3b81a173 --- /dev/null +++ b/PWG1/AliTreeDraw.h @@ -0,0 +1,62 @@ +#ifndef ALITREEDRAW_H +#define ALITREEDRAW_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + + + +////////////////////////////////////////////////////////////////////////////// +// Class AliTreeDraw // +// // +// marian.ivanov@cern.ch // +////////////////////////////////////////////////////////////////////////////// + + + +#include +#include + +class TH1; +class TH1F; +class TH2F; +class TTree; + +class AliTreeDraw: public TObject{ +public: + AliTreeDraw(); + ~AliTreeDraw(){;} + void SetTree(TTree *tree){fTree=tree;} + const TH1 * GetRes() const{ return (TH1*)fRes;} + const TH1 * GetMean() const{ return (TH1*)fMean;} + const TObjArray* GetPoints() const {return fPoints;} + void ClearHisto(); + void ClearPoints(){if (fPoints) fPoints->Clear();} + // + TH1F * DrawXY(const char * chx, const char *chy, const char* selection, + const char * quality,Int_t nbins, Float_t minx, Float_t maxx, + Float_t miny, Float_t maxy, Int_t nBinsRes=30); + TH1F * DrawLogXY(const char * chx, const char *chy, const char* selection, + const char * quality, Int_t nbins,Float_t minx, Float_t maxx, + Float_t miny, Float_t maxy, Int_t nBinsRes=30); + TH1F * Eff(const char *variable, const char* selection, const char * quality, + Int_t nbins,Float_t min, Float_t max); + TH1F * EffLog(const char *variable, const char* selection, const char * quality, + Int_t nbins,Float_t min, Float_t max); + // + void GetPoints3D(const char * label, const char * chpoints, const char* selection, TTree * tpoints, Int_t color=6, Float_t rmin=4.); + + private: + static void AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle); + static Double_t* CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax); + static TH1F* CreateEffHisto(TH1F* hGen, TH1F* hRec); + static TH1F* CreateResHisto(TH2F* hRes2, TH1F **phMean, + Bool_t drawBinFits = kTRUE,Bool_t overflowBinFits = kFALSE); + + TTree * fTree; //the tree for visualization - NOT OWNER + TH1F * fRes; //temporary histogram - OWNER + TH1F * fMean; //temporary histogram - OWNER + TObjArray *fPoints;// - OWNER + ClassDef(AliTreeDraw,0) +}; + +#endif diff --git a/PWG1/PWG1LinkDef.h b/PWG1/PWG1LinkDef.h new file mode 100644 index 00000000000..212156bbb58 --- /dev/null +++ b/PWG1/PWG1LinkDef.h @@ -0,0 +1,16 @@ +#ifdef __CINT__ + +#pragma link off all glols; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class AliTreeDraw+; +#pragma link C++ class AliTPCdigitRow+; +#pragma link C++ class AliGenV0Info+; +#pragma link C++ class AliGenKinkInfo+; +#pragma link C++ class AliGenInfoMaker+; +#pragma link C++ class AliMCInfo+; + + + +#endif diff --git a/PWG1/libPWG1.pkg b/PWG1/libPWG1.pkg new file mode 100644 index 00000000000..062658667b7 --- /dev/null +++ b/PWG1/libPWG1.pkg @@ -0,0 +1,9 @@ + +SRCS:= AliTreeDraw.cxx \ + AliGenInfo.cxx + +HDRS:= $(SRCS:.cxx=.h) + +DHDR:= PWG1LinkDef.h + +EINCLUDE:= TPC diff --git a/PWG1/rootlogon.C b/PWG1/rootlogon.C new file mode 100644 index 00000000000..21171f1ff66 --- /dev/null +++ b/PWG1/rootlogon.C @@ -0,0 +1,20 @@ + +{ + printf("\nWELCOME to ALICE\n\n"); + gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER\ + -I$ALICE_ROOT/TPC -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TRD -I$ALICE_ROOT/TOF -I$ALICE_ROOT/RAW -I$ALICE_ROOT/kdtree/include"); + //load library neccesary for TPC alone + gSystem->Load("$(ROOTSYS)/lib/libEG"); + gSystem->Load("$(ROOTSYS)/lib/libGeom"); + gSystem->Load("$(ROOTSYS)/lib/libVMC"); + gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libSTEER"); + gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTPCrec"); + gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTPCsim"); + gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTRDrec"); + gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTRDsim"); + //gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTRD"); + //gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libTRACKING"); + + +} + -- 2.39.3