--- /dev/null
+/**************************************************************************
+ * 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,1,0);
+t2->Exec();
+
+//
+//some cuts definition
+TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.00005&&abs(MC.fVDist[2])<0.0005")
+//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 csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
+
+
+TCut crec("crec","fReconstructed==1");
+TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
+TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
+TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
+TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
+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();
+
+//
+//example
+comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
+
+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.DrawLogXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&abs(fPdg)==211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
+comp.fRes->Draw();
+comp.fMean->Draw();
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120&&abs(fPdg)==211"+cteta1+cpos1+cprim,"fTPCOn",20,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,250); hdedx0->SetMarkerColor(2);
+TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,250); hdedx1->SetMarkerColor(3);
+TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,250); hdedx2->SetMarkerColor(4);
+TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,250); hdedx3->SetMarkerColor(5);
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==211&&fITStrack.fN==6"+cprim)
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim)
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==321&&fITStrack.fN==6"+cprim)
+comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==11&&fITStrack.fN==6"+cprim)
+hdedx0->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx3->Draw("same");
+
+comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
+
+*/
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <string.h>
+//ROOT includes
+#include "Rtypes.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "TString.h"
+#include "TBenchmark.h"
+#include "TStopwatch.h"
+#include "TParticle.h"
+#include "TSystem.h"
+#include "TTimer.h"
+#include "TVector3.h"
+#include "TPad.h"
+#include "TCanvas.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TText.h"
+#include "Getline.h"
+#include "TStyle.h"
+
+//ALIROOT includes
+#include "AliRun.h"
+#include "AliStack.h"
+#include "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 "AliESDtrack.h"
+#include "AliITStrackV2.h"
+#include "AliTRDtrack.h"
+#endif
+#include "AliGenInfo.h"
+#include "AliESDComparisonMI.h"
+
+
+
+//
+//
+void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed)
+{
+ //
+ //
+ //calculates derived variables
+ //
+ //
+ 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]);
+ }
+ //
+ if (reconstructed==kFALSE){
+ fReconstructed = kFALSE;
+ fTPCOn = kFALSE;
+ fITSOn = kFALSE;
+ fTRDOn = kFALSE;
+ return;
+ }
+ 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 (fTPCOn){
+ //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.0000001){
+ fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
+ fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
+ }
+ Double_t cov[15], param[5],x;
+ fESDTrack.GetInnerExternalCovariance(cov);
+ fESDTrack.GetInnerExternalParameters(x,param);
+ //
+ 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;
+ fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(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;
+ fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);
+ }
+
+}
+
+
+void AliESDRecV0Info::Update(Float_t vertex[3])
+{
+ Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
+ Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
+
+ Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
+ Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
+ vnorm2 = TMath::Sqrt(vnorm2);
+ Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
+ Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
+ pnorm2 = TMath::Sqrt(pnorm2);
+
+ fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
+ fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
+ fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+ESDCmpTr::ESDCmpTr()
+{
+ Reset();
+}
+
+////////////////////////////////////////////////////////////////////////
+ESDCmpTr::ESDCmpTr(const char* fnGenTracks,
+ const char* fnCmp,
+ const char* fnGalice, Int_t direction,
+ Int_t nEvents, Int_t firstEvent)
+{
+ Reset();
+ // fFnGenTracks = fnGenTracks;
+ // fFnCmp = fnCmp;
+ sprintf(fFnGenTracks,"%s",fnGenTracks);
+ sprintf(fFnCmp,"%s",fnCmp);
+
+ fFirstEventNr = firstEvent;
+ fEventNr = firstEvent;
+ fNEvents = nEvents;
+ fDirection = direction;
+ //
+ fLoader = AliRunLoader::Open(fnGalice);
+ if (gAlice){
+ //delete gAlice->GetRunLoader();
+ delete gAlice;
+ gAlice = 0x0;
+ }
+ if (fLoader->LoadgAlice()){
+ cerr<<"Error occured while l"<<endl;
+ }
+ Int_t nall = fLoader->GetNumberOfEvents();
+ if (nEvents==0) {
+ nEvents =nall;
+ fNEvents=nall;
+ fFirstEventNr=0;
+ }
+
+ if (nall<=0){
+ cerr<<"no events available"<<endl;
+ fEventNr = 0;
+ return;
+ }
+ if (firstEvent+nEvents>nall) {
+ fEventNr = nall-firstEvent;
+ cerr<<"restricted number of events availaible"<<endl;
+ }
+ AliMagF * magf = gAlice->Field();
+ AliTracker::SetFieldMap(magf);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+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"<<endl;
+ return 1;
+ }
+ return 0;
+}
+
+//////////////////////////////////////////////////////////////
+
+Int_t ESDCmpTr::SetIO(Int_t eventNr)
+{
+ //
+ //
+ // SET INPUT
+ //
+ TFile f("AliESDs.root");
+ //
+
+ TTree* tree = (TTree*) f.Get("esdTree");
+ if (!tree) {
+ 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);
+ }
+ }
+ else{
+ tree->SetBranchAddress("ESD", &fEvent);
+ tree->GetEntry(eventNr);
+ }
+
+
+ /*
+ 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;
+ // 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;
+ cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
+ for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
+ eventNr++) {
+ fEventNr = eventNr;
+ SetIO(fEventNr);
+ fNParticles = gAlice->GetEvent(fEventNr);
+
+ fIndexRecTracks = new Int_t[fNParticles*4]; //write at maximum 4 tracks corresponding to particle
+ fFakeRecTracks = new Int_t[fNParticles];
+ fMultiRecTracks = new Int_t[fNParticles];
+ for (Int_t i = 0; i<fNParticles; i++) {
+ fIndexRecTracks[4*i] = -1;
+ fIndexRecTracks[4*i+1] = -1;
+ fIndexRecTracks[4*i+2] = -1;
+ fIndexRecTracks[4*i+3] = -1;
+
+ fFakeRecTracks[i] = 0;
+ fMultiRecTracks[i] = 0;
+ }
+
+ cout<<"Start to process event "<<fEventNr<<endl;
+ cout<<"\tfNParticles = "<<fNParticles<<endl;
+ if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
+ if (TreeTLoop()>0) return 1;
+
+ if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
+ if (TreeGenLoop(eventNr)>0) return 1;
+ if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
+
+ delete [] fIndexRecTracks;
+ delete [] fFakeRecTracks;
+ delete [] fMultiRecTracks;
+ }
+
+ CloseOutputFile();
+
+ cerr<<"Exec finished"<<endl;
+ timer.Stop();
+ timer.Print();
+ return 0;
+
+}
+////////////////////////////////////////////////////////////////////////
+Bool_t ESDCmpTr::ConnectGenTree()
+{
+//
+// connect all branches from the genTracksTree
+// use the same variables as for the new cmp tree, it may work
+//
+ fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
+ if (!fFileGenTracks) {
+ cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
+ return kFALSE;
+ }
+ fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
+ if (!fTreeGenTracks) {
+ cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
+ <<fFnGenTracks<<endl;
+ return kFALSE;
+ }
+ //
+ fMCInfo = new AliMCInfo;
+ fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
+ //
+ if (fDebug > 1) {
+ cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
+ }
+ return kTRUE;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+void ESDCmpTr::CreateTreeCmp()
+{
+ fFileCmp = TFile::Open(fFnCmp,"RECREATE");
+ if (!fFileCmp) {
+ cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
+ return;
+ }
+
+
+ fTreeCmp = new TTree("ESDcmpTracks","ESDcmpTracks");
+ //
+ fMCInfo = new AliMCInfo;
+ fRecInfo = new AliESDRecInfo;
+ //
+ AliESDtrack * esdTrack = new AliESDtrack;
+ AliITStrackV2 * itsTrack = new AliITStrackV2;
+ //
+ fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
+ fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo);
+ fTreeCmp->Branch("fESDTrack","AliESDtrack",&esdTrack);
+ // fTreeCmp->Branch("ITS","AliITStrackV2",&itsTrack);
+ delete esdTrack;
+ //
+ fTreeCmp->AutoSave();
+
+}
+////////////////////////////////////////////////////////////////////////
+void ESDCmpTr::CloseOutputFile()
+{
+ if (!fFileCmp) {
+ cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
+ return;
+ }
+ fFileCmp->cd();
+ fTreeCmp->Write();
+ delete fTreeCmp;
+
+ fFileCmp->Close();
+ delete fFileCmp;
+ return;
+}
+////////////////////////////////////////////////////////////////////////
+
+TVector3 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->Transform1to2(x,index);
+ return TVector3(x);
+}
+////////////////////////////////////////////////////////////////////////
+
+Int_t ESDCmpTr::TreeTLoop()
+{
+ //
+ // loop over all ESD reconstructed tracks and store info in memory
+ //
+ TStopwatch timer;
+ timer.Start();
+ //
+ Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();
+ //
+ fTracks = new TObjArray(nEntries);
+ //
+ //load tracks to the memory
+ for (Int_t i=0; i<nEntries;i++){
+ AliESDtrack * track =fEvent->GetTrack(i);
+ fTracks->AddAt(track,i);
+ }
+ //
+ //
+ AliESDtrack * track=0;
+ for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
+ track = (AliESDtrack*)fTracks->UncheckedAt(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]<4)
+ fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] = iEntry;
+ }
+ else
+ fIndexRecTracks[absLabel*4] = iEntry;
+ fMultiRecTracks[absLabel]++;
+ }
+ }
+ printf("Time spended in TreeTLoop\n");
+ timer.Print();
+
+ if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
+{
+//
+// loop over all entries for a given event, find corresponding
+// rec. track and store in the fTreeCmp
+//
+ TStopwatch timer;
+ timer.Start();
+ Int_t entry = fNextTreeGenEntryToRead;
+ Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
+ cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
+ <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
+ TBranch * branch = fTreeCmp->GetBranch("RC");
+ branch->SetAddress(&fRecInfo); // set all pointers
+
+ while (entry < nParticlesTR) {
+ fTreeGenTracks->GetEntry(entry);
+ entry++;
+ if (eventNr < fMCInfo->fEventNr) continue;
+ if (eventNr > fMCInfo->fEventNr) continue;;
+ //
+ fNextTreeGenEntryToRead = entry-1;
+ if (fDebug > 2 && fMCInfo->fLabel < 10) {
+ cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
+ }
+ //
+ fRecInfo->Reset();
+ AliESDtrack * track=0;
+ fRecInfo->fReconstructed =0;
+ TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
+ local.GetXYZ(fRecInfo->fTRLocalCoord);
+ //
+ if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
+ track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
+ //
+ //
+ // find nearest track if multifound
+ if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
+ Double_t xyz[3];
+ track->GetInnerXYZ(xyz);
+ Float_t dz = TMath::Abs(local.Z()-xyz[2]);
+ for (Int_t i=1;i<4;i++){
+ if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
+ AliESDtrack * track2 = (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
+ track2->GetInnerXYZ(xyz);
+ Float_t dz2 = TMath::Abs(local.Z()-xyz[2]);
+ if (TMath::Abs(dz2)<dz)
+ track = track2;
+ }
+ }
+ }
+ //
+ fRecInfo->fESDTrack =*track;
+ if (track->GetITStrack())
+ fRecInfo->fITStrack = *((AliITStrackV2*)track->GetITStrack());
+ 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);
+ }
+ else{
+ fRecInfo->Update(fMCInfo,fParamTPC,kFALSE);
+ }
+
+ fTreeCmp->Fill();
+ }
+ fTreeCmp->AutoSave();
+ fTracks->Delete();
+ printf("Time spended in TreeGenLoop\n");
+ timer.Print();
+ if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
+
+ return 0;
+}
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+
+void AliESDComparisonDraw::SetIO(const char *fname)
+{
+ //
+ TFile* file = TFile::Open(fname);
+ //
+ fTree = (TTree*) file->Get("ESDcmpTracks");
+ if (!fTree) {
+ printf("no track comparison tree found\n");
+ file->Close();
+ delete file;
+ }
+}
+
+
--- /dev/null
+
+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 Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed);
+ 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
+ AliITStrackV2 fITStrack; //its track
+ AliTRDtrack fTRDtrack; //its track
+ 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
+ Bool_t fITSOn; // ITS refitted inward
+ Bool_t fTRDOn; // ITS refitted inward
+ Float_t fDeltaP; //delta of momenta
+ void Reset();
+ ClassDef(AliESDRecInfo,1) // 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 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
+ ClassDef(AliESDRecV0Info,1) // container for
+};
+
+ClassImp(AliESDRecV0Info)
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// 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();
+ 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
+ //
+ char fFnGenTracks[1000]; //! input file name with gen tracks
+ TFile *fFileGenTracks;
+ TTree *fTreeGenTracks;
+ //
+ //
+ Int_t fDirection;
+ //
+ AliRunLoader * fLoader; //! pointer to the run loader
+ //TTree *fTreeRecTracks; //! tree with reconstructed tracks
+ //
+ Int_t *fIndexRecTracks; //! index of particle label in the TreeT_ESD
+ Int_t *fFakeRecTracks; //! number of fake tracks
+ Int_t *fMultiRecTracks; //! number of multiple reconstructions
+ //
+ TObjArray * fTracks; //!container with tracks
+ AliESD *fEvent; //!event
+
+ //
+ 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
+ //
+ AliMCInfo* fMCInfo; //! MC information writen per particle
+ AliESDRecInfo* fRecInfo; //! Rec. information writen per particle
+ //
+
+ ClassDef(ESDCmpTr,1) // class which creates and fills tree with ESDGenTrack objects
+};
+ClassImp(ESDCmpTr)
+
+
+
--- /dev/null
+/**************************************************************************
+ * 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
+
+Macro to generate comples MC information - used for Comparison later on
+How to use it?
+
+.L $ALICE_ROOT/STEER/AliGenInfo.C+
+AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",1,0)
+t->Exec();
+
+*/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <string.h>
+//ROOT includes
+#include "Rtypes.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "TString.h"
+#include "TBenchmark.h"
+#include "TStopwatch.h"
+#include "TParticle.h"
+#include "TSystem.h"
+#include "TTimer.h"
+#include "TVector3.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TPad.h"
+#include "TF1.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"
+#endif
+#include "AliGenInfo.h"
+//
+//
+
+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);
+ fTRdecay.SetTrack(-1);
+}
+
+AliMCInfo::~AliMCInfo()
+{
+ if (fTPCReferences) {
+ delete fTPCReferences;
+ }
+ if (fITSReferences){
+ delete fITSReferences;
+ }
+ if (fTRDReferences){
+ delete fTRDReferences;
+ }
+}
+
+
+
+void AliMCInfo::Update()
+{
+ //
+ //
+ fMCtracks =1;
+ if (!fTPCReferences) {
+ fNTPCRef =0;
+ return;
+ }
+ Float_t direction=1;
+ //Float_t rlast=0;
+ fNTPCRef = fTPCReferences->GetEntriesFast();
+ fNITSRef = fITSReferences->GetEntriesFast();
+
+ for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();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();
+
+}
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////
+//
+// End of implementation of the class AliMCInfo
+//
+////////////////////////////////////////////////////////////////////////
+
+
+
+////////////////////////////////////////////////////////////////////////
+digitRow::digitRow()
+{
+ Reset();
+}
+////////////////////////////////////////////////////////////////////////
+digitRow & digitRow::operator=(const digitRow &digOld)
+{
+ for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
+ return (*this);
+}
+////////////////////////////////////////////////////////////////////////
+void digitRow::SetRow(Int_t row)
+{
+ if (row >= 8*kgRowBytes) {
+ cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
+ return;
+ }
+ Int_t iC = row/8;
+ Int_t iB = row%8;
+ SETBIT(fDig[iC],iB);
+}
+
+////////////////////////////////////////////////////////////////////////
+Bool_t digitRow::TestRow(Int_t row)
+{
+//
+// return kTRUE if row is on
+//
+ Int_t iC = row/8;
+ Int_t iB = row%8;
+ return TESTBIT(fDig[iC],iB);
+}
+////////////////////////////////////////////////////////////////////////
+Int_t digitRow::RowsOn(Int_t upto)
+{
+//
+// returns number of rows with a digit
+// count only rows less equal row number upto
+//
+ Int_t total = 0;
+ for (Int_t i = 0; i<kgRowBytes; i++) {
+ for (Int_t j = 0; j < 8; j++) {
+ if (i*8+j > upto) return total;
+ if (TESTBIT(fDig[i],j)) total++;
+ }
+ }
+ return total;
+}
+////////////////////////////////////////////////////////////////////////
+void digitRow::Reset()
+{
+//
+// resets all rows to zero
+//
+ for (Int_t i = 0; i<kgRowBytes; i++) {
+ fDig[i] <<= 8;
+ }
+}
+////////////////////////////////////////////////////////////////////////
+Int_t digitRow::Last()
+{
+//
+// returns the last row number with a digit
+// returns -1 if now digits
+//
+ for (Int_t i = kgRowBytes-1; i>=0; i--) {
+ for (Int_t j = 7; j >= 0; j--) {
+ if TESTBIT(fDig[i],j) return i*8+j;
+ }
+ }
+ return -1;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t digitRow::First()
+{
+//
+// returns the first row number with a digit
+// returns -1 if now digits
+//
+ for (Int_t i = 0; i<kgRowBytes; i++) {
+ for (Int_t j = 0; j < 8; j++) {
+ if (TESTBIT(fDig[i],j)) return i*8+j;
+ }
+ }
+ return -1;
+}
+
+////////////////////////////////////////////////////////////////////////
+//
+// end of implementation of a class digitRow
+//
+////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////
+AliGenInfoMaker::AliGenInfoMaker()
+{
+ Reset();
+}
+
+////////////////////////////////////////////////////////////////////////
+AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
+ Int_t nEvents, Int_t firstEvent)
+{
+ Reset();
+ fFirstEventNr = firstEvent;
+ fEventNr = firstEvent;
+ fNEvents = nEvents;
+ // fFnRes = fnRes;
+ sprintf(fFnRes,"%s",fnRes);
+ //
+ fLoader = AliRunLoader::Open(fnGalice);
+ if (gAlice){
+ delete gAlice->GetRunLoader();
+ delete gAlice;
+ gAlice = 0x0;
+ }
+ if (fLoader->LoadgAlice()){
+ cerr<<"Error occured while l"<<endl;
+ }
+ Int_t nall = fLoader->GetNumberOfEvents();
+ if (nEvents==0) {
+ nEvents =nall;
+ fNEvents=nall;
+ fFirstEventNr=0;
+ }
+
+ if (nall<=0){
+ cerr<<"no events available"<<endl;
+ fEventNr = 0;
+ return;
+ }
+ if (firstEvent+nEvents>nall) {
+ fEventNr = nall-firstEvent;
+ cerr<<"restricted number of events availaible"<<endl;
+ }
+ AliMagF * magf = gAlice->Field();
+ AliTracker::SetFieldMap(magf);
+}
+
+
+AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
+{
+ //
+ if (i<fNParticles) {
+ if (fGenInfo[i]) return fGenInfo[i];
+ fGenInfo[i] = new AliMCInfo;
+ fNInfos++;
+ return fGenInfo[i];
+ }
+ else
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+void AliGenInfoMaker::Reset()
+{
+ fEventNr = 0;
+ fNEvents = 0;
+ fTreeGenTracks = 0;
+ fFileGenTracks = 0;
+ fGenInfo = 0;
+ fNInfos = 0;
+ //
+ //
+ fDebug = 0;
+ fVPrim[0] = -1000.;
+ fVPrim[1] = -1000.;
+ fVPrim[2] = -1000.;
+ fParamTPC = 0;
+}
+////////////////////////////////////////////////////////////////////////
+AliGenInfoMaker::~AliGenInfoMaker()
+{
+
+ if (fLoader){
+ fLoader->UnloadgAlice();
+ 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();
+ 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; i<fNParticles; i++) {
+ fGenInfo[i]=0;
+ }
+ //
+ cout<<"Start to process event "<<fEventNr<<endl;
+ cout<<"\tfNParticles = "<<fNParticles<<endl;
+ if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
+ if (TreeTRLoop()>0) return 1;
+ //
+ if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
+ if (TreeDLoop()>0) return 1;
+ //
+ if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
+ if (TreeKLoop()>0) return 1;
+ if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
+ for (UInt_t i = 0; i<fNParticles; i++) {
+ if (fGenInfo[i]) delete fGenInfo[i];
+ }
+ delete []fGenInfo;
+ CloseIOEvent();
+ }
+ //
+ CloseIO();
+ CloseOutputFile();
+
+ cerr<<"Exec finished"<<endl;
+
+ timer.Stop();
+ timer.Print();
+ return 0;
+}
+////////////////////////////////////////////////////////////////////////
+void AliGenInfoMaker::CreateTreeGenTracks()
+{
+ fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
+ if (!fFileGenTracks) {
+ cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
+ return;
+ }
+ fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
+ AliMCInfo * info = new AliMCInfo;
+ //
+ fTreeGenTracks->Branch("MC","AliMCInfo",&info);
+ delete info;
+ fTreeGenTracks->AutoSave();
+}
+////////////////////////////////////////////////////////////////////////
+void AliGenInfoMaker::CloseOutputFile()
+{
+ if (!fFileGenTracks) {
+ cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
+ return;
+ }
+ fFileGenTracks->cd();
+ fTreeGenTracks->Write();
+ delete fTreeGenTracks;
+ 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 "<<fNParticles<<" primary and secondary particles in event "
+ <<fEventNr<<endl;
+ }
+ Int_t ipdg = 0;
+ TParticlePDG *ppdg = 0;
+ // not all generators give primary vertex position. Take the vertex
+ // of the particle 0 as primary vertex.
+ TDatabasePDG pdg; //get pdg table
+ //thank you very much root for this
+ TBranch * br = fTreeGenTracks->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"<<endl;
+ return 0;
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::TreeDLoop()
+{
+ //
+ // open the file with treeD
+ // loop over all entries there and save information about some tracks
+ //
+
+ Int_t nInnerSector = fParamTPC->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 = "<<sectorsByRows<<endl;
+ for (Int_t i=0; i<sectorsByRows; i++) {
+ if (!fTreeD->GetEvent(i)) continue;
+ Int_t sec,row;
+ fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
+ if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
+ // here I expect that upper sectors follow lower sectors
+ if (sec > 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 "
+ <<sec<<" "
+ <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
+ <<" "<<row<<endl;
+ }
+ AliMCInfo * info = GetInfo(label);
+ if (info){
+ info->fTPCRow.SetRow(row+rowShift);
+ }
+ }
+ }
+ }
+ } while (digits->Next());
+ }
+
+ if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
+ return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliGenInfoMaker::TreeTRLoop()
+{
+ //
+ // loop over TrackReferences and store the first one for each track
+ //
+ TTree * treeTR = fTreeTR;
+ Int_t nPrimaries = (Int_t) treeTR->GetEntries();
+ if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
+ //
+ //
+ //track references for TPC
+ TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
+ TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
+ TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
+ TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
+ //
+ if (treeTR->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("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
+ //
+ //
+ //
+ for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
+ treeTR->GetEntry(iPrimPart);
+ //
+ // Loop over TPC references
+ //
+ for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
+ AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
+ //
+ if (trackRef->Pt()<fgTPCPtCut) continue;
+ Int_t label = trackRef->GetTrack();
+ AliMCInfo * info = GetInfo(label);
+ if (!info) info = MakeInfo(label);
+ if (!info) continue;
+ 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->Pt()<fgTPCPtCut) continue;
+ Int_t label = trackRef->GetTrack();
+ AliMCInfo * info = GetInfo(label);
+ if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
+ if (!info) info = MakeInfo(label);
+ if (!info) continue;
+ 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->Pt()<fgTPCPtCut) continue;
+ Int_t label = trackRef->GetTrack();
+ AliMCInfo * info = GetInfo(label);
+ if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
+ if (!info) info = MakeInfo(label);
+ if (!info) continue;
+ TClonesArray & arr = *(info->fTRDReferences);
+ 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;
+ 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];
+}
+////////////////////////////////////////////////////////////////////////
+
+
+
+TH1F * AliComparisonDraw::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)
+{
+ //
+ Double_t* bins = CreateLogBins(nbins, minx, maxx);
+ Int_t nBinsRes = 30;
+ 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;
+ fRes = hRes;
+ fMean = hMean;
+ return hRes;
+}
+
+TH1F * AliComparisonDraw::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)
+{
+ //
+ Double_t* bins = CreateLogBins(nbins, minx, maxx);
+ Int_t nBinsRes = 30;
+ 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;
+ fRes = hRes;
+ fMean = hMean;
+ return hRes;
+}
+
+///////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////
+TH1F * AliComparisonDraw::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 * AliComparisonDraw::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* AliComparisonDraw::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 AliComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
+{
+ //
+ histo->GetXaxis()->SetTitle(xAxisTitle);
+ histo->GetYaxis()->SetTitle(yAxisTitle);
+}
+
+
+TH1F* AliComparisonDraw::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* AliComparisonDraw::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
+ //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
+ // TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
+ 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;
+}
+
+
+
--- /dev/null
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class digitRow
+//
+////////////////////////////////////////////////////////////////////////
+const Int_t kgRowBytes = 32;
+
+class digitRow: public TObject {
+
+public:
+ digitRow();
+ virtual ~digitRow(){;}
+ void SetRow(Int_t row);
+ Bool_t TestRow(Int_t row);
+ digitRow & operator=(const digitRow &digOld);
+ Int_t RowsOn(Int_t upto=8*kgRowBytes);
+ Int_t Last();
+ Int_t First();
+ void Reset();
+
+//private:
+ UChar_t fDig[kgRowBytes];
+ ClassDef(digitRow,1) // container for digit pattern
+};
+ClassImp(digitRow)
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// 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
+ //
+ 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
+ digitRow fTPCRow; // information about digits row pattern
+ Int_t fNTPCRef; // tpc references counter
+ Int_t fNITSRef; // tpc references counter
+ TClonesArray * fTPCReferences; //containner with all track references -in the TPC
+ TClonesArray * fITSReferences; //container with ITS references
+ TClonesArray * fTRDReferences; //container with TRD references
+ //
+ ClassDef(AliMCInfo,3) // container for
+};
+ClassImp(AliMCInfo)
+
+
+
+class AliGenV0Info: public TObject {
+public:
+ AliMCInfo fMCd; //info about daughter particle
+ AliMCInfo fMCm; //info about mother particle
+ void Update(); // 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
+ ClassDef(AliGenV0Info,1) // container for
+};
+ClassImp(AliGenV0Info)
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// 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();
+ //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<fNParticles)? fGenInfo[i]:0;}
+ AliMCInfo * MakeInfo(UInt_t i);
+
+public:
+ Int_t fDebug; //! debug flag
+ Int_t fEventNr; //! current event number
+ Int_t fLabel; //! track label
+ Int_t fNEvents; //! number of events to process
+ Int_t fFirstEventNr; //! first event to process
+ UInt_t fNParticles; //! number of particles in TreeK
+ TTree *fTreeGenTracks; //! output tree with generated tracks
+ char fFnRes[1000]; //! output file name with stored tracks
+ TFile *fFileGenTracks; //! output file with stored fTreeGenTracks
+ //
+ AliRunLoader * fLoader; //! pointer to the run loader
+ TTree * fTreeD; //! current tree with digits
+ TTree * fTreeTR; //! current tree with TR
+ AliStack *fStack; //! current stack
+ //
+ AliMCInfo ** fGenInfo; //! array with pointers to gen info
+ Int_t fNInfos; //! number of tracks with infos
+ //
+ AliTPCParam* fParamTPC; //! AliTPCParam
+ Double_t fVPrim[3]; //! primary vertex position
+ // the fVDist[3] contains size of the 3-vector
+
+private:
+
+// some constants for the original non-pareller tracking (by Y.Belikov)
+ static const Int_t seedRow11 = 158; // nRowUp - 1
+ static const Int_t seedRow12 = 139; // nRowUp - 1 - (Int_t) 0.125*nRowUp
+ static const Int_t seedRow21 = 149; // seedRow11 - shift
+ static const Int_t seedRow22 = 130; // seedRow12 - shift
+ static const Double_t kRaddeg = 180./3.14159265358979312;
+ //
+ static const Double_t fgTPCPtCut = 0.1; // do not store particles with generated pT less than this
+ static const Double_t fgITSPtCut = 0.2; // do not store particles with generated pT less than this
+ static const Double_t fgTRDPtCut = 0.2; // do not store particles with generated pT less than this
+
+ ClassDef(AliGenInfoMaker,1) // class which creates and fills tree with TPCGenTrack objects
+};
+ClassImp(AliGenInfoMaker)
+
+
+
+class AliComparisonDraw: public TObject{
+public:
+ 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);
+ 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);
+ 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);
+ //
+ 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);
+ public:
+ TTree * fTree;
+ TH1F * fRes; //temporary file
+ TH1F * fMean; //temporary file
+ ClassDef(AliComparisonDraw,1)
+};
+ClassImp(AliComparisonDraw)
+
+
+
+AliTPCParam * GetTPCParam();
+Float_t TPCBetheBloch(Float_t bg);
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+
+#include "AliITStrackerMI.h"
+#include "AliITSclusterV2.h"
+
+#include "AliITShit.h"
+#include "AliITS.h"
+#include "AliITSgeom.h"
+//#include "AliITSgeometryDetail.h"
+//#include "AliITSparameter.h"
+#include "alles.h"
+#include "TFile.h"
+#include "TStopwatch.h"
+#include "Rtypes.h"
+#include "TTree.h"
+
+
+#include "AliRunLoader.h"
+#include "AliStack.h"
+#include "TF1.h"
+#include "AliTrackReference.h"
+#endif
+#include "AliITSclusterComparison.h"
+
+
+ClassImp(AliITSCI)
+ClassImp(AliITSClusterErrAnal)
+
+
+
+AliITSClusterErrAnal::AliITSClusterErrAnal(Char_t *chloader )
+{
+ //
+ //SET Input loaders
+ if (gAlice){
+ delete gAlice->GetRunLoader();
+ delete gAlice;
+ gAlice = 0x0;
+ }
+ fRunLoader = AliRunLoader::Open(chloader);
+ if (fRunLoader == 0x0){
+ cerr<<"Can not open session"<<endl;
+ return;
+ }
+ fITSLoader = fRunLoader->GetLoader("ITSLoader");
+ if (fITSLoader == 0x0){
+ cerr<<"Can not get ITS Loader"<<endl;
+ return ;
+ }
+ if (fRunLoader->LoadgAlice()){
+ cerr<<"Error occured while l"<<endl;
+ return;
+ }
+ fRunLoader->CdGAFile();
+ fITS = (AliITS*)gAlice->GetDetector("ITS");
+ if (!fITS) {
+ cerr<<"AliITSclusterComparison.C : Can not find the ITS detector !"<<endl;
+ }
+ AliITSgeom *geom = fITS->GetITSgeom();
+ //An instance of the ITS tracker
+ fTracker = new AliITStrackerMI(geom);
+ //
+ //
+ AliITSCI * clinfo = new AliITSCI();
+ fFileA = new TFile("itsclusteranal.root","recreate");
+ fFileA->cd();
+ fTreeA = new TTree("itscl","itscl");
+ fTreeA->Branch("itscl","AliITSCI",&clinfo);
+
+ AliITSclusterV2 * cl = new AliITSclusterV2;
+ fTreeB = new TTree("Clusters","Clusters");
+ fTreeB->Branch("cl","AliITSclusterV2",&cl);
+
+ fClusters = 0;
+ delete clinfo;
+}
+
+void AliITSClusterErrAnal::SetIO(Int_t event)
+{
+ //
+ //set input output for given event
+ fRunLoader->SetEventNumber(event);
+ fRunLoader->LoadHeader();
+ fRunLoader->LoadKinematics();
+ fRunLoader->LoadTrackRefs();
+ fITSLoader->LoadHits();
+ fITSLoader->LoadRecPoints("read");
+ //
+ fStack = fRunLoader->Stack();
+ fHitTree = fITSLoader->TreeH();
+ fClusterTree = fITSLoader->TreeR();
+ fReferenceTree = fRunLoader->TreeTR();
+ //
+}
+
+void AliITSClusterErrAnal::LoadClusters()
+{
+ //
+ //
+ // Load clusters
+ if (fClusters) {
+ fClusters->Delete();
+ delete fClusters;
+ }
+ Int_t nparticles = fStack->GetNtrack();
+ fClusters = new TObjArray(nparticles);
+ //
+ TObjArray *ClusterArray = new TClonesArray("AliITSclusterV2");
+ TObjArray carray(2000);
+ TBranch *branch=fClusterTree->GetBranch("Clusters");
+ if (!branch) {
+ Error("ReadClusters","Can't get the branch !");
+ return;
+ }
+ branch->SetAddress(&ClusterArray);
+ TBranch *brcl = fTreeB->GetBranch("cl");
+ AliITSclusterV2* cluster= new AliITSclusterV2;
+ brcl->SetAddress(&cluster);
+ Int_t nentries = (Int_t)fClusterTree->GetEntries();
+ for (Int_t i=0;i<nentries;i++){
+ fClusterTree->GetEvent(i);
+ Int_t nCluster = ClusterArray->GetEntriesFast();
+ SignDeltas(ClusterArray,-3.8);
+ for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
+ cluster = (AliITSclusterV2*)ClusterArray->UncheckedAt(iCluster);
+ carray.AddLast(new AliITSclusterV2(*cluster));
+ fTreeB->Fill();
+ }
+ }
+ //
+ Int_t nClusters = carray.GetEntriesFast();
+ printf("Total number of clusters %d \n", nClusters);
+ fTracker->LoadClusters(fClusterTree);
+ //
+ //
+ //SORT clusters
+ //
+ Int_t all=0;
+ Int_t noisecl =0;
+ for (Int_t i=0;i<nClusters;i++){
+ AliITSclusterV2 *cl = (AliITSclusterV2 *) carray.UncheckedAt(i);
+ if (cl->GetLabel(0)<0) noisecl++;
+ //
+ for (Int_t itrack=0; itrack<3;itrack++){
+ Int_t lab = cl->GetLabel(itrack);
+ if (lab>=0){
+ TObjArray * array = (TObjArray*)fClusters->At(lab);
+ if (array==0){
+ array = new TObjArray(20);
+ fClusters->AddAt(array,lab);
+ }
+ array->AddLast(cl);
+ all++;
+ }
+ }
+ }
+ printf("Total number of track clusters %d \n", all);
+ printf("Total number of noise clusters %d \n", noisecl);
+}
+
+void AliITSClusterErrAnal::SignDeltas( TObjArray *ClusterArray, Float_t vz)
+{
+ //
+ //
+ Int_t entries = ClusterArray->GetEntriesFast();
+ if (entries<4) return;
+ AliITSclusterV2* cluster = (AliITSclusterV2*)ClusterArray->At(0);
+ Int_t layer = cluster->GetLayer();
+ if (layer>1) return;
+ Int_t index[10000];
+ Int_t ncandidates=0;
+ Float_t r = (layer>0)? 7:4;
+ //
+ for (Int_t i=0;i<entries;i++){
+ AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(i);
+ Float_t nz = 1+(cl0->GetZ()-vz)/r;
+ if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
+ index[ncandidates] = i; //candidate to belong to delta electron track
+ ncandidates++;
+ if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
+ cl0->SetDeltaProbability(1);
+ }
+ }
+ //
+ //
+ //
+ for (Int_t i=0;i<ncandidates;i++){
+ AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(index[i]);
+ if (cl0->GetDeltaProbability()>0.8) continue;
+ //
+ Int_t ncl = 0;
+ Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
+ sumy=sumz=sumy2=sumyz=sumw=0.0;
+ for (Int_t j=0;j<ncandidates;j++){
+ if (i==j) continue;
+ AliITSclusterV2* cl1 = (AliITSclusterV2*)ClusterArray->At(index[j]);
+ //
+ Float_t dz = cl0->GetZ()-cl1->GetZ();
+ Float_t dy = cl0->GetY()-cl1->GetY();
+ if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
+ Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
+ y[ncl] = cl1->GetY();
+ z[ncl] = cl1->GetZ();
+ sumy+= y[ncl]*weight;
+ sumz+= z[ncl]*weight;
+ sumy2+=y[ncl]*y[ncl]*weight;
+ sumyz+=y[ncl]*z[ncl]*weight;
+ sumw+=weight;
+ ncl++;
+ }
+ }
+ if (ncl<4) continue;
+ Float_t det = sumw*sumy2 - sumy*sumy;
+ Float_t delta=1000;
+ if (TMath::Abs(det)>0.01){
+ Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
+ Float_t k = (sumyz*sumw - sumy*sumz)/det;
+ delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
+ }
+ else{
+ Float_t z0 = sumyz/sumy;
+ delta = TMath::Abs(cl0->GetZ()-z0);
+ }
+ if ( delta<0.05) {
+ cl0->SetDeltaProbability(1-20.*delta);
+ }
+ }
+}
+
+void AliITSClusterErrAnal::LoadParticles()
+{
+ //
+ //
+ // fReferences = new TObjArray;
+}
+
+void AliITSClusterErrAnal::SortReferences()
+{
+ //
+ //
+ //
+ printf("Sorting references\n");
+ fReferences = new TObjArray;
+ Int_t ntracks = fStack->GetNtrack();
+ fReferences->Expand(ntracks);
+ Int_t nentries = (Int_t)fReferenceTree->GetEntries();
+ TClonesArray * arr = new TClonesArray("AliTrackReference");
+ TBranch * br = fReferenceTree->GetBranch("ITS");
+ br->SetAddress(&arr);
+ //
+ TClonesArray *labarr=0;
+ Int_t nreferences=0;
+ Int_t nreftracks=0;
+ for (Int_t iprim=0;iprim<nentries;iprim++){
+ if (br->GetEntry(iprim)){
+ for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){
+ AliTrackReference *ref =(AliTrackReference*)arr->At(iref);
+ if (!ref) continue;
+ Int_t lab = ref->GetTrack();
+ //if ( (lab<0) || (lab>ntracks)) continue;
+ //
+ if (fReferences->At(lab)==0) {
+ labarr = new TClonesArray("AliTrackReference");
+ fReferences->AddAt(labarr,lab);
+ nreftracks++;
+ }
+ TClonesArray &larr = *labarr;
+ new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
+ nreferences++;
+ }
+ }
+ }
+ printf("Total number of references = \t%d\n", nreferences);
+ printf("Total number of tracks with references = \t%d\n", nreftracks);
+ printf("End - Sorting references\n");
+
+}
+
+
+void AliITSClusterErrAnal::GetNTeor(Int_t layer, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz)
+{
+ //
+ //get "mean shape"
+ if (layer==0){
+ ny = 1.+TMath::Abs(phi)*3.2;
+ nz = 1.+TMath::Abs(theta)*0.34;
+ return;
+ }
+ if (layer==1){
+ ny = 1.+TMath::Abs(phi)*3.2;
+ nz = 1.+TMath::Abs(theta)*0.28;
+ return;
+ }
+
+ if (layer>3){
+ ny = 2.02+TMath::Abs(phi)*1.95;
+ nz = 2.02+TMath::Abs(phi)*2.35;
+ return;
+ }
+ ny = 6.6-2.7*abs(phi);
+ nz = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta);
+}
+
+
+Int_t AliITSClusterErrAnal::GetError(Int_t layer, const AliITSclusterV2*cl, Float_t theta, Float_t phi,Float_t expQ, Float_t &erry, Float_t &errz)
+{
+ //calculate cluster position error
+ //
+ Float_t nz,ny;
+ GetNTeor(layer, theta,phi,ny,nz);
+ erry = TMath::Sqrt(cl->GetSigmaY2());
+ errz = TMath::Sqrt(cl->GetSigmaZ2());
+ //
+ // PIXELS
+ if (layer<2){
+
+ if (TMath::Abs(ny-cl->GetNy())>0.6) {
+ if (ny<cl->GetNy()){
+ erry*=0.4+TMath::Abs(ny-cl->GetNy());
+ errz*=0.4+TMath::Abs(ny-cl->GetNy());
+ }else{
+ erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
+ errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
+ }
+ }
+ if (TMath::Abs(nz-cl->GetNz())>1.) {
+ erry*=TMath::Abs(nz-cl->GetNz());
+ errz*=TMath::Abs(nz-cl->GetNz());
+ }
+ erry*=0.85;
+ errz*=0.85;
+ erry= TMath::Min(erry,float(0.005));
+ errz= TMath::Min(errz,float(0.03));
+ return 10;
+ }
+ //STRIPS
+ if (layer>3){
+ if (cl->GetNy()==100||cl->GetNz()==100){
+ erry = 0.004;
+ errz = 0.2;
+ return 100;
+ }
+ if (cl->GetNy()+cl->GetNz()>12){
+ erry = 0.06;
+ errz = 0.57;
+ return 100;
+ }
+ Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
+ Float_t chargematch = TMath::Max(double(normq/expQ),2.);
+ //
+ if (cl->GetType()==1 || cl->GetType()==10 ){
+ if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
+ errz = 0.043;
+ erry = 0.00094;
+ return 101;
+ }
+ if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
+ errz = 0.06;
+ erry =0.0013;
+ return 102;
+ }
+ erry = 0.0027;
+ errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15);
+ return 103;
+ }
+ if (cl->GetType()==2 || cl->GetType()==11 ){
+ erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05);
+ errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5);
+ return 104;
+ }
+
+ if (cl->GetType()>100 ){
+ if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
+ errz = 0.05;
+ erry = 0.00096;
+ return 105;
+ }
+ if (cl->GetNy()+cl->GetNz()-nz-ny<1){
+ errz = 0.10;
+ erry = 0.0025;
+ return 106;
+ }
+
+ errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4);
+ erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05);
+ return 107;
+ }
+ Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
+ if (diff<1) diff=1;
+ if (diff>4) diff=4;
+
+ if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
+ errz = 0.14*diff;
+ erry = 0.003*diff;
+ return 108;
+ }
+ erry = 0.04*diff;
+ errz = 0.06*diff;
+ return 109;
+ }
+
+ //DRIFTS
+ Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
+ Float_t chargematch = normq/expQ;
+ Float_t factorz=1;
+ Int_t cnz = cl->GetNz()%10;
+ //charge match
+ if (cl->GetType()==1){
+ if (chargematch<1.25){
+ erry = 0.0028*(1.+6./cl->GetQ()); // gold clusters
+ }
+ else{
+ erry = 0.003*chargematch;
+ if (cl->GetNz()==3) erry*=1.5;
+ }
+ if (chargematch<1.0){
+ errz = 0.0011*(1.+6./cl->GetQ());
+ }
+ else{
+ errz = 0.002*(1+2*(chargematch-1.));
+ }
+ if (cnz>nz+0.6) {
+ erry*=(cnz-nz+0.5);
+ errz*=1.4*(cnz-nz+0.5);
+ }
+ }
+ if (cl->GetType()>1){
+ if (chargematch<1){
+ erry = 0.00385*(1.+6./cl->GetQ()); // gold clusters
+ errz = 0.0016*(1.+6./cl->GetQ());
+ }
+ else{
+ errz = 0.0014*(1+3*(chargematch-1.));
+ erry = 0.003*(1+3*(chargematch-1.));
+ }
+ if (cnz>nz+0.6) {
+ erry*=(cnz-nz+0.5);
+ errz*=1.4*(cnz-nz+0.5);
+ }
+ }
+
+ if (TMath::Abs(cl->GetY())>2.5){
+ factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
+ }
+ if (TMath::Abs(cl->GetY())<1){
+ factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
+ }
+ factorz= TMath::Min(factorz,float(4.));
+ errz*=factorz;
+
+ erry= TMath::Min(erry,float(0.05));
+ errz= TMath::Min(errz,float(0.05));
+ return 200;
+}
+
+
+AliITSclusterV2 * AliITSClusterErrAnal::FindNearestCluster(AliITSCI * clinfo, Float_t dmax)
+{
+ //
+ //
+ AliTrackReference * ref = &(clinfo->fRef);
+ TObjArray * clusters = (TObjArray*)(fClusters->At(ref->GetTrack()));
+ if (!clusters) return 0;
+ Int_t nclusters = clusters->GetEntriesFast();
+ Float_t maxd = dmax;
+ AliITSclusterV2 * cluster=0;
+ Float_t radius = ref->R();
+ clinfo->fNClusters=0;
+ for (Int_t icluster=0;icluster<nclusters;icluster++){
+ AliITSclusterV2 * cl = (AliITSclusterV2*)clusters->At(icluster);
+ if (!cl) continue;
+ Float_t dz = cl->GetZ()-ref->Z();
+ //
+ if (TMath::Abs(dz)>dmax) continue;
+ Int_t idet=cl->GetDetectorIndex();
+ Double_t r=-1;
+ Float_t dy =0;
+ // for (Int_t i=cl->fLayer;i<=cl->fLayer-1;i++){
+
+ AliITStrackerMI::AliITSdetector & det = fTracker->GetDetector(cl->fLayer,idet);
+ if (TMath::Abs(radius-det.GetR())<1.5) {
+ Bool_t isOK=kFALSE;
+ for (Int_t itrack=0;itrack<3;itrack++){
+ if (cl->fTracks[itrack]== clinfo->fRef.GetTrack()) {
+ clinfo->fNClusters++;
+ isOK= kTRUE;
+ }
+ }
+ if (!isOK) continue;
+ r = det.GetR();
+ Double_t phicl= det.GetPhi();
+ Double_t ca=TMath::Cos(phicl), sa=TMath::Sin(phicl);
+ Double_t x = double(ref->X())*ca +double(ref->Y())*sa;
+ Double_t y = -double(ref->X())*sa +double(ref->Y())*ca;
+ if (TMath::Abs(x-r-0.012)>0.1) continue;
+ //
+ Double_t px = ref->Px(), py = ref->Py();
+ Double_t lpx = double(px)*ca +double(py)*sa;
+ Double_t lpy = -double(px)*sa +double(py)*ca;
+ Double_t lphi = TMath::ATan2(lpy,lpx);
+
+ Double_t theta = ref->Pz()/ref->Pt();
+ Double_t dx = x-r;
+ dy = y -cl->GetY()-dx*lphi;
+ dz = ref->Z()-cl->GetZ()-dx*theta;
+ //dx-=0.012;
+ Double_t d = TMath::Sqrt(dy*dy+dz*dz);
+ if (d<maxd){
+ cluster=cl;
+ maxd=d;
+ clinfo->Update(fTracker);
+ clinfo->fCl = *cluster;
+ clinfo->fLx = x;
+ clinfo->fLy = y;
+ clinfo->fLayer = cl->fLayer;
+ clinfo->fPhiCl = phicl;
+ clinfo->fRCl = det.GetR();
+ clinfo->fPhiL = lphi;
+ clinfo->fDz = dz;
+ clinfo->fDy = dy;
+ cl->SetSigmaZ2(TMath::Abs(cl->GetSigmaZ2()));
+ cl->SetSigmaY2(TMath::Abs(cl->GetSigmaY2()));
+ clinfo->fPoolZ = dz/TMath::Sqrt(cl->GetSigmaZ2());
+ clinfo->fPoolY = dy/TMath::Sqrt(cl->GetSigmaY2());
+ clinfo->fPoolY2 = dy/TMath::Sqrt(cl->GetSigmaY2());
+ clinfo->fErrY = TMath::Sqrt(cl->GetSigmaY2());
+ clinfo->fErrZ = TMath::Sqrt(cl->GetSigmaZ2());
+ Float_t erry,errz,ny,nz;
+ clinfo->fErrType = GetError(cl->fLayer,cl,theta,lphi,clinfo->fExpectedQ,erry,errz);
+ GetNTeor(cl->fLayer, theta,lphi,ny,nz);
+ clinfo->fNormQ = cl->GetQ()/TMath::Sqrt(1+theta*theta+lpy*lpy/(lpx*lpx));
+ clinfo->fTNy = ny;
+ clinfo->fTNz = nz;
+ clinfo->fErrY = erry;
+ clinfo->fErrZ = errz;
+ clinfo->fPoolY2 = dy/clinfo->fErrY;
+ clinfo->fPoolZ2 = dz/clinfo->fErrZ;
+ clinfo->fLabPos =-1;
+ Double_t xyz[3];
+ det.GetGlobalXYZ(cl,xyz);
+ clinfo->fGx = xyz[0];
+ clinfo->fGy = xyz[1];
+ for (Int_t i=0;i<3;i++) if (cl->GetLabel(i)==clinfo->fRef.GetTrack()) clinfo->fLabPos = i;
+ }
+ }
+ }
+ // }
+ return cluster;
+}
+
+
+
+Int_t AliITSClusterErrAnal::Analyze(Int_t trackmax) {
+ //
+ //
+ // dummy cluster to be fill if not cluster info
+ //
+ AliITSCI * clinfo = new AliITSCI;
+ AliITSclusterV2 dummy;
+ TBranch * branch = fTreeA->GetBranch("itscl");
+ branch->SetAddress(&clinfo);
+ Int_t nall =0;
+ Int_t withcl =0;
+ Int_t ntracks = fReferences->GetEntriesFast();
+ ntracks = TMath::Min(trackmax,ntracks);
+ for (Int_t itrack=0;itrack<ntracks;itrack++){
+ TClonesArray * references = (TClonesArray*) fReferences->At(itrack);
+ if (!references) continue;
+ Int_t nreferences = references->GetEntriesFast();
+ if (nreferences<4) continue;
+ TObjArray * clarray = (TObjArray*)fClusters->At(itrack);
+ for (Int_t iref=0;iref<nreferences;iref++){
+ AliTrackReference* trackref = (AliTrackReference*)references->At(iref);
+ if (trackref&& trackref->P()>0.01){
+ TParticle * p = fStack->Particle(trackref->GetTrack());
+ if (p) clinfo->fP = *p;
+ clinfo->fRef = *trackref;
+ clinfo->fStatus = 0;
+ clinfo->Update(fTracker);
+ nall++;
+ AliITSclusterV2 * cluster = FindNearestCluster(clinfo,3.0);
+ if (cluster){
+ //clinfo->fCl = *cluster;
+ clinfo->fStatus=1;
+ withcl++;
+ }else{
+ clinfo->fCl= dummy;
+ }
+ clinfo->Update(fTracker);
+ fTreeA->Fill();
+ }
+ }
+ }
+
+ fFileA->cd();
+ fTreeA->Write();
+ fTreeB->Write();
+ fFileA->Close();
+ return 0;
+}
+
+
+
+void AliITSClusterErrAnal::MakeSeeds(Double_t zv)
+{
+ //
+ //
+ //
+
+ AliITStrackerMI::AliITSlayer & layer3 = fTracker->GetLayer(3);
+ AliITStrackerMI::AliITSlayer & layer2 = fTracker->GetLayer(2);
+ AliITStrackerMI::AliITSlayer & layer1 = fTracker->GetLayer(1);
+ AliITStrackerMI::AliITSlayer & layer0 = fTracker->GetLayer(0);
+ Int_t ncl3 = layer3.GetNumberOfClusters();
+ Int_t ncl2 = layer2.GetNumberOfClusters();
+ Int_t ncl1 = layer1.GetNumberOfClusters();
+ Int_t ncl0 = layer0.GetNumberOfClusters();
+ TH1F *hc = new TH1F("hc","hc",100,-0.01,0.01);
+ TH1F *hrc = new TH1F("hrc","hrc",100,-1,1);
+ TH1F *hl = new TH1F("hl","hl",100,-0.05,0.05);
+ TH1F *hlr = new TH1F("hlr","hlr",1000,-30.,30.);
+ TH2F *hcl = new TH2F("hcl","hcl",1000,-0.01,0.01,1000,-0.03,0.03);
+ TH2F *hcc = new TH2F("hcc","hcc",1000,-0.01,0.01,1000,-0.03,0.03);
+ TH1F *hr = new TH1F("hr","hr",100,0,2000.);
+ //
+ TH1F *hpoolc = new TH1F("hpoolc","hpoolc",100,-5,5);
+ TH1F *hpooll = new TH1F("hpooll","hpooll",100,-5,5);
+ //
+ Double_t ratio23 = layer2.GetR()/layer3.GetR();
+ Double_t ratio12 = layer1.GetR()/layer2.GetR();
+ Double_t referenceR = (layer0.GetR()+layer1.GetR()+layer2.GetR()+layer3.GetR())*0.25;
+ //
+ Double_t *clusterxyzr3 = new Double_t[ncl3*6];
+ Double_t *clusterxyzr2 = new Double_t[ncl2*6];
+ Double_t *clusterxyzr1 = new Double_t[ncl1*6];
+ Double_t *clusterxyzr0 = new Double_t[ncl0*6];
+ //
+ // Get global coordinates 3
+ for (Int_t icl3=0;icl3<ncl3;icl3++){
+ AliITSclusterV2* cl3 = layer3.GetCluster(icl3);
+ Double_t *xyz3 = &clusterxyzr3[6*icl3];
+ AliITStrackerMI::AliITSdetector & det = layer3.GetDetector(cl3->GetDetectorIndex());
+ det.GetGlobalXYZ(cl3,xyz3);
+ xyz3[3] = TMath::Sqrt(xyz3[0]*xyz3[0]+xyz3[1]*xyz3[1]);
+ Double_t fi3 = TMath::ATan2(xyz3[1],xyz3[0]);
+ xyz3[4] = TMath::Cos(fi3);
+ xyz3[5] = TMath::Sin(fi3);
+ }
+ //
+ // Get global coordinates 2
+ for (Int_t icl2=0;icl2<ncl2;icl2++){
+ AliITSclusterV2* cl2 = layer2.GetCluster(icl2);
+ Double_t *xyz2 = &clusterxyzr2[6*icl2];
+ AliITStrackerMI::AliITSdetector & det = layer2.GetDetector(cl2->GetDetectorIndex());
+ det.GetGlobalXYZ(cl2,xyz2);
+ xyz2[3] = TMath::Sqrt(xyz2[0]*xyz2[0]+xyz2[1]*xyz2[1]);
+ Double_t fi2 = TMath::ATan2(xyz2[1],xyz2[0]);
+ xyz2[4] = TMath::Cos(fi2);
+ xyz2[5] = TMath::Sin(fi2);
+ }
+ //
+ // Get global coordinates 1
+ for (Int_t icl1=0;icl1<ncl1;icl1++){
+ AliITSclusterV2* cl1 = layer1.GetCluster(icl1);
+ Double_t *xyz1 = &clusterxyzr1[6*icl1];
+ AliITStrackerMI::AliITSdetector & det = layer1.GetDetector(cl1->GetDetectorIndex());
+ det.GetGlobalXYZ(cl1,xyz1);
+ xyz1[3] = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
+ Double_t fi1 = TMath::ATan2(xyz1[1],xyz1[0]);
+ xyz1[4] = TMath::Cos(fi1);
+ xyz1[5] = TMath::Sin(fi1);
+ }
+ //
+ // Get global coordinates 0
+ for (Int_t icl0=0;icl0<ncl0;icl0++){
+ AliITSclusterV2* cl0 = layer0.GetCluster(icl0);
+ Double_t *xyz0 = &clusterxyzr0[6*icl0];
+ AliITStrackerMI::AliITSdetector & det = layer0.GetDetector(cl0->GetDetectorIndex());
+ det.GetGlobalXYZ(cl0,xyz0);
+ xyz0[3] = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
+ Double_t fi0 = TMath::ATan2(xyz0[1],xyz0[0]);
+ xyz0[4] = TMath::Cos(fi0);
+ xyz0[5] = TMath::Sin(fi0);
+
+ }
+ //
+ //
+ Int_t ntracks=0;
+ Int_t ntracks0=0;
+ Int_t good0=0;
+ Int_t accepted0=0;
+ //
+ Int_t nfoundtracks=0;
+ Int_t good1 =0;;
+ Int_t accepted1=0;
+ Int_t accepted2=0;
+ //
+ Int_t *trackindex = new Int_t[ncl3*4*300];
+ const AliITSclusterV2 **trackclusters = new const AliITSclusterV2*[ncl3*4*300];
+ //
+ for (Int_t icl3=0;icl3<ncl3;icl3++){
+ AliITSclusterV2* cl3 = layer3.GetCluster(icl3);
+ Double_t *xyz3 = &clusterxyzr3[icl3*6];
+ Double_t &r3 = clusterxyzr3[icl3*6+3];
+ ntracks=0; //number of possible tracks corresponding to given cluster
+ //
+ //
+ // find pairs in layer 2
+ Float_t zc = (xyz3[2]-zv)*ratio23+zv;
+ Float_t lphi = TMath::ATan2(xyz3[1],xyz3[0]);
+ layer2.SelectClusters(zc-1,zc+1, layer2.GetR()*lphi-1, layer2.GetR()*lphi+1.);
+ const AliITSclusterV2 * cl2=0;
+ Int_t npairs2=0;
+ Int_t ci2[100];
+ while ( (cl2=layer2.GetNextCluster(ci2[npairs2]))!=0){
+ Int_t index2 = ci2[npairs2];
+ Double_t *xyz2 = &clusterxyzr2[index2*6];
+ Double_t &r2 = clusterxyzr2[index2*6+3];
+ //
+ if (TMath::Abs((xyz3[2]-zv)*r2/r3-(xyz2[2]-zv))>0.3) continue;
+ npairs2++;
+ }
+ //
+ //
+ // find pairs in layer 2
+ for (Int_t jcl2 =0;jcl2<npairs2; jcl2++){
+ Int_t index2 = ci2[jcl2];
+ Double_t *xyz2 = &clusterxyzr2[index2*6];
+ Float_t zc = (xyz2[2]-zv)*ratio12+zv;
+ Float_t lphi = TMath::ATan2(xyz2[1],xyz2[0]);
+ cl2 = layer2.GetCluster(index2);
+ //
+ //
+ layer1.SelectClusters(zc-1,zc+1, layer1.GetR()*lphi-1, layer1.GetR()*lphi+1.);
+ const AliITSclusterV2 * cl1=0;
+ Int_t ci1;
+ //
+ Double_t lr2 = (xyz3[0]-xyz2[0])*(xyz3[0]-xyz2[0])+(xyz3[1]-xyz2[1])*(xyz3[1]-xyz2[1]);
+ Double_t lr = TMath::Sqrt(lr2);
+ Double_t sinphi = (xyz3[1]-xyz2[1])/lr;
+ Double_t cosphi = (xyz3[0]-xyz2[0])/lr;
+ Double_t t1 = (xyz3[2]-xyz2[2])/lr;
+ //
+ while ( (cl1=layer1.GetNextCluster(ci1))!=0){
+ if (cl1->GetQ()==0) continue;
+ Int_t index1 = ci1;
+ Double_t *xyz1 = &clusterxyzr1[index1*6];
+ AliITSclusterV2* cl1 = layer1.GetCluster(index1);
+
+ if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) good0++;
+ Double_t loc1[2];
+ //
+ loc1[0] = (xyz1[0]-xyz2[0])*cosphi+(xyz1[1]-xyz2[1])*sinphi;
+ loc1[1] = -(xyz1[0]-xyz2[0])*sinphi+(xyz1[1]-xyz2[1])*cosphi;
+ //
+ Double_t c1 = 2.*loc1[1]/(loc1[0]*(loc1[0]-lr)); //curvature approximated using 2 derivation
+ Double_t c=c1;
+ Double_t t2 = (xyz2[2]-xyz1[2])/TMath::Sqrt(loc1[0]*loc1[0]+loc1[1]*loc1[1]);
+ if (TMath::Abs(c)>0.008) continue;
+ if (TMath::Abs(t2-t1)>0.02) continue;
+ //
+ trackclusters[ntracks*4] = 0;
+ trackclusters[ntracks*4+1] = cl1;
+ trackclusters[ntracks*4+2] = cl2;
+ trackclusters[ntracks*4+3] = cl3;
+ trackindex[ntracks*4+0] = 0;
+ trackindex[ntracks*4+1] = index1;
+ trackindex[ntracks*4+2] = index2;
+ trackindex[ntracks*4+3] = icl3;
+ if (ntracks>=300*ncl3){
+ break;
+ }
+ if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) {
+ //printf("%f\n",c);
+
+ //hc->Fill(c);
+ //hl->Fill(t2-t1);
+ //hr->Fill(1/c);
+ //hcc->Fill(c,TMath::Abs(c1));
+ //hcl->Fill(c1,t2-t1);
+
+ accepted0++;
+ }
+ ntracks++;
+ }
+ } // end - find pairs in layer 2
+ //
+ ntracks0+=ntracks;
+ //hcc->Draw();
+ Int_t ntracks2=0;
+ Int_t cindex2[10000*4];
+ Float_t pools[10000];
+
+ for (Int_t itrack=0;itrack<ntracks;itrack++){
+ const AliITSclusterV2 *cl0, *cl1,*cl2,*cl3;
+ cl3 = trackclusters[itrack*4+3];
+ cl2 = trackclusters[itrack*4+2];
+ cl1 = trackclusters[itrack*4+1];
+ Int_t index1 = trackindex[itrack*4+1];
+ Int_t index2 = trackindex[itrack*4+2];
+ Int_t index3 = trackindex[itrack*4+3];
+ //
+ Double_t *xyz1 = &clusterxyzr1[index1*6],*xyz2 = &clusterxyzr2[index2*6],*xyz3=&clusterxyzr3[index3*6];
+ Double_t &r1 = clusterxyzr1[index1*6+3],&r2 = clusterxyzr2[index2*6+3], &r3 = clusterxyzr3[index3*6+3];
+
+ Double_t r0 = layer0.GetR();
+ Float_t zc = xyz1[2] + (xyz2[2]-xyz1[2])*(r0-r1)/(r2-r1);
+ Float_t lphi = TMath::ATan2(xyz1[1],xyz1[0]);
+
+ layer0.SelectClusters(zc-1,zc+1, layer0.GetR()*lphi-0.3, layer0.GetR()*lphi+0.3);
+ Int_t ci0;
+ //
+ //
+ Double_t lr2 = (xyz2[0]-xyz1[0])*(xyz2[0]-xyz1[0])+(xyz2[1]-xyz1[1])*(xyz2[1]-xyz1[1]);
+ Double_t lr = TMath::Sqrt(lr2);
+ Double_t sinphi = (xyz2[1]-xyz1[1])/lr;
+ Double_t cosphi = (xyz2[0]-xyz1[0])/lr;
+ //
+ Double_t loc3[2]; // cluster3 in local coordindate
+ loc3[0] = (xyz3[0]-xyz1[0])*cosphi+(xyz3[1]-xyz1[1])*sinphi;
+ loc3[1] = -(xyz3[0]-xyz1[0])*sinphi+(xyz3[1]-xyz1[1])*cosphi;
+
+ Double_t tan1 = (xyz3[2]-xyz2[2])/TMath::Sqrt((xyz3[0]-xyz2[0])*(xyz3[0]-xyz2[0])+(xyz3[1]-xyz2[1])*(xyz3[1]-xyz2[1]));
+ Double_t c3 = 2.*loc3[1]/(loc3[0]*(loc3[0]-lr)); //curvature approximated using 2 derivation
+ //
+ while ( (cl0=layer0.GetNextCluster(ci0))!=0){
+ if (cl0->GetQ()==0) continue;
+ Double_t *xyz0 = &clusterxyzr0[ci0*6];
+ if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0)
+ && cl3->GetLabel(0)==cl1->GetLabel(0)&&cl3->GetLabel(0)==cl0->GetLabel(0)){
+ good1++;
+ }
+ Double_t loc0[2]; // cluster3 in local coordindate
+ loc0[0] = (xyz0[0]-xyz1[0])*cosphi+(xyz0[1]-xyz1[1])*sinphi;
+ loc0[1] = -(xyz0[0]-xyz1[0])*sinphi+(xyz0[1]-xyz1[1])*cosphi;
+ //
+ Double_t c0 = 2.*loc0[1]/(loc0[0]*(loc0[0]-lr)); //curvature approximated using 2 derivation
+ Double_t tan2 = (xyz1[2]-xyz0[2])/TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])+(xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]));
+ //if (TMath::Abs(tan1-tan2)>0.05) continue;
+ //if (TMath::Abs(c3-c0)>0.006) continue;
+ Double_t poolphi = (c3-c0)/0.0011;
+ poolphi*=poolphi;
+ Double_t poolr = ((c3-c0)/(TMath::Abs(c3+c0)+0.001))/0.14;
+ Double_t poolth = (tan1-tan2)/0.0094;
+ poolth*=poolth;
+ Double_t poolthr = ((tan1-tan2)/(TMath::Abs(c3+c0)+0.005))/0.83;
+ if (poolr*poolr+poolthr*poolthr>40) continue;
+ //if ( TMath::Abs((c3-c0)/(c3+c0))>0.6) continue;
+ cindex2[ntracks2*4+0] = ci0;
+ cindex2[ntracks2*4+1] = index1;
+ cindex2[ntracks2*4+2] = index2;
+ cindex2[ntracks2*4+3] = index3;
+ pools[ntracks2] = poolr*poolr+poolthr*poolthr;
+ //
+ ntracks2++;
+ if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0)
+ && cl3->GetLabel(0)==cl1->GetLabel(0)){
+ hl->Fill(tan1-tan2);
+ hrc->Fill(poolr);
+ hc->Fill((c3-c0));
+ hlr->Fill(poolthr);
+ hpoolc->Fill(poolr);
+ hpooll->Fill(poolthr);
+ accepted1++;
+ }
+ }
+ }
+ Int_t sindex[10000];
+ TMath::Sort(ntracks2,pools,sindex,kFALSE);
+ Float_t minpool = pools[sindex[0]];
+ Int_t max = TMath::Min(ntracks2,20);
+ for (Int_t itrack=0;itrack<ntracks2;itrack++){
+ if (itrack>=max) continue;
+ //if (pools[sindex[itrack]]>minpool+5) continue;
+ const AliITSclusterV2 *cl0, *cl1,*cl2,*cl3;
+ cl3 = layer3.GetCluster(cindex2[sindex[itrack]*4+3]);
+ cl2 = layer2.GetCluster(cindex2[sindex[itrack]*4+2]);
+ cl1 = layer1.GetCluster(cindex2[sindex[itrack]*4+1]);
+ cl0 = layer0.GetCluster(cindex2[sindex[itrack]*4+0]);
+ if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0)
+ && cl3->GetLabel(0)==cl1->GetLabel(0)){
+ accepted2++;
+ }
+ nfoundtracks++;
+ }
+ }
+ hc->Draw();
+ printf("First pass \t%d\t%d\t%d\n",ntracks0,good0,accepted0);
+ printf("Found tracks\t%d\t%d\t%d\t%d\n",nfoundtracks,good1, accepted1,accepted2);
+
+}
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "AliTrackReference.h"
+#endif
+
+
+class AliITSCI: public TObject {
+ public :
+ AliITSCI(){;}
+ virtual ~AliITSCI(){;}
+ Float_t GetExpectedQ();
+ //
+ AliITSclusterV2 fCl;
+ AliTrackReference fRef;
+ TParticle fP;
+ Char_t fStatus;
+ //
+ Float_t fBeta2;
+ Float_t fVertex;
+ Int_t fNClusters;
+ Int_t fLabPos;
+ Float_t fExpectedQ;
+ Float_t fNormQ;
+ Int_t fPDG;
+ Int_t fNTracks;
+ Float_t fPt;
+ Float_t fPtotal;
+ Float_t fR;
+ Float_t fCharge;
+ Bool_t fIsPrim;
+ Double_t fDy;
+ Double_t fDz;
+ Double_t fPoolY;
+ Double_t fPoolY2;
+ Double_t fPoolZ;
+ Double_t fPoolZ2;
+ Double_t fErrY;
+ Double_t fErrZ;
+ Int_t fErrType;
+ Double_t fPhiCl;
+ Double_t fRCl;
+ Double_t fLx;
+ Double_t fLy;
+ Double_t fGx;
+ Double_t fGy;
+ Float_t fTNy;
+ Float_t fTNz;
+ Float_t fTheta;
+ Double_t fPhi;
+ Double_t fPhiL;
+ Int_t fLayer;
+ void Update(AliITStrackerMI*tracker);
+ ClassDef(AliITSCI,1)
+};
+
+class AliITSClusterErrAnal: public TObject{
+public:
+ AliITSClusterErrAnal(Char_t *chloader = "galice.root");
+ void SetIO(Int_t event);
+ Int_t Analyze(Int_t trackmax);
+ void LoadClusters();
+ void LoadParticles();
+ void SignDeltas( TObjArray *ClusterArray, Float_t vz);
+ void SortReferences();
+ AliITSclusterV2 * FindNearestCluster(AliITSCI * clinfo, Float_t dmax=0.5);
+ void GetNTeor(Int_t layer, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz);
+ Int_t GetError(Int_t layer, const AliITSclusterV2*cl, Float_t theta, Float_t phi, Float_t expQ, Float_t &erry, Float_t &errz);
+ void MakeSeeds(Double_t zv);
+public:
+ AliRunLoader * fRunLoader;
+ AliLoader * fITSLoader;
+ TTree * fHitTree;
+ TTree * fClusterTree;
+ TTree * fReferenceTree;
+ AliITS * fITS;
+ //
+ TTree * fTreeA;
+ TTree * fTreeB;
+ TFile * fFileA;
+ AliITStrackerMI *fTracker;
+ AliStack *fStack;
+ TObjArray *fClusters; //array of clusters
+ TObjArray fExactPoints;
+ TObjArray *fReferences;
+ TObjArray *fParticles;
+ ClassDef(AliITSClusterErrAnal,1)
+};
+
+
+void AliITSCI::Update(AliITStrackerMI*tracker)
+{
+ //
+ //
+ fPt = fRef.Pt();
+ fR = fRef.R();
+ fCharge = fP.GetPDG()->Charge();
+ //fIsPrim = (fP.GetMother(0)>0)? kFALSE :kTRUE;
+ Double_t vertex = TMath::Sqrt( fP.Vx()*fP.Vx()
+ +fP.Vy()*fP.Vy());
+ fVertex = vertex;
+ fIsPrim = (vertex<1) ? kTRUE:kFALSE;
+ fTheta = fRef.Pz()/fRef.Pt();
+ fPhi = TMath::ATan2(fRef.Py(),fRef.Px());
+ fPDG = fP.GetPdgCode();
+ fExpectedQ = GetExpectedQ();
+}
+
+
+Float_t AliITSCI::GetExpectedQ()
+{
+ //
+ //
+ //
+ //fPtotal = fP.P();
+ fPtotal = fRef.P();
+
+ Double_t p2= fPtotal*fPtotal;
+ Float_t mass = fP.GetPDG()->Mass();
+ Double_t beta2=p2/(p2 + mass*mass);
+ fBeta2 = beta2;
+ Float_t d=0.003;
+ // Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
+ //dE+=0.33e-3/(beta2)*d;
+ //if (beta2/(1-beta2)>3.5*3.5)
+ // dE=0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2)*d;
+ Double_t dE = 0.22*0.153e-3*(39.2-55.6*beta2+28.7*beta2*beta2+27.41/beta2)*d;
+
+ return dE*10000000.; //normalization to ADC
+ //return dE; //normalization to ADC
+}
+
+/*
+
+.L $ALICE_ROOT/COMPARISON/AliITSclusterComparison.C+
+.L $ALICE_ROOT/COMPARISON/AliGenInfo.C+
+
+AliITSClusterErrAnal anal;
+anal.SetIO(0);
+anal.LoadParticles();
+anal.SortReferences();
+anal.LoadClusters();
+
+
+anal.Analyze(30000000);
+
+
+anal.MakeSeeds(-3.63); >seed1.txt
+
+
+TFile f("itsclusteranal.root");
+TTree * tree= (TTree*)f.Get("itscl")
+TTree * tree2= (TTree*)f.Get("Clusters")
+
+AliComparisonDraw comp;
+comp.fTree = tree;
+
+TH1F * hpools = new TH1F("pools","pools",100,-6,6); hpools->SetLineColor(4); hpools->SetFillColor(0);
+TH1F * hpools0 = new TH1F("pools0","pools0",100,-6,6); hpools0->SetLineColor(2); hpools0->SetFillColor(0);
+
+TH1F * hsigma = new TH1F("hsigma","hsigma",1000,0,0.3)
+
+TCut cprim("cprim","fCl.fTracks[0]<110000&&fCl.fTracks[1]<110000&&fCl.fTracks[2]<110000");
+
+TCut cprim("cprim","fCl.fTracks[0]<24000&&fCl.fTracks[1]<24000&&fCl.fTracks[2]<24000");
+TCut cpt("cpt","sqrt(fRef.fPx**2+fRef.fPy**2)>0.05");
+TCut cphi("cphi","abs(fPhiL)<0.9");
+TCut ctheta("ctheta","abs(fTheta)<1.1");
+TCut cangle = cphi+ctheta;
+TCut cpools5("cpools5","abs(fPoolY)<5&&abs(fPoolZ)<5");
+TCut cpools52("cpools52","abs(fPoolY2)<5&&abs(fPoolZ2)<5");
+TCut cpools62("cpools52","sqrt(fPoolY2**2+fPoolZ2**2)<6");
+TCut cgold("cgold","abs(fCl.fNy+fCl.fNz-fTNy-fTNz)<1.5");
+TCut cback("cback","fRef.fX*fRef.fPx+fRef.fPy*fRef.fY<0");
+
+comp.fTree->Draw("fPoolY2>>pools","fStatus==1&&abs(fLayer-0.5)<0.6&&fIsPrim"+cpt+cangle+cpools62,"")
+
+comp.DrawXY("fErrType","fPoolZ2","fStatus==1&&abs(fLayer-4.5)<0.6","1",8,99.9,108,-5,5)
+
+comp.DrawXY("fCl.fType","fPoolY2","fStatus==1&&abs(fLayer-4.5)<0.6&&fIsPrim","1",4,102.9,106.1,-5,5)
+comp.DrawXY("abs(fCl.fChargeRatio)","fPoolY2","fStatus==1&&abs(fLayer-4.5)<0.6&&fIsPrim","1",4,0.0,1,-5,5)
+
+TProfile prof1("prof1","prof1",30,0.1,1); prof1.SetLineColor(2);
+TProfile prof2("prof2","prof2",30,0.1,1); prof2.SetLineColor(3);
+
+comp.fTree->Draw("fNormQ/fExpectedQ:fBeta2>>prof1","fStatus==1&&fCl.fLayer>1&&fNormQ/fExpectedQ<1.5&&fCl.fNz<4&&abs(fPDG)==211&&fIsPrim")
+comp.fTree->Draw("fNormQ/fExpectedQ:fBeta2>>prof2","fStatus==1&&fCl.fLayer>1&&fNormQ/fExpectedQ<1.5&&fCl.fNz<4&&abs(fPDG)!=211&&fIsPrim")
+
+prof1->Draw();prof2->Draw("same");
+
+
+TH1F * hdeltam = new TH1F("hdeltam" ,"hdeltam",100,-3,10); hdeltam->SetLineColor(4); hdeltam->SetFillColor(0);
+TH1F * hdelta0 = new TH1F("hdelta0" ,"hdelta0",100,-3,10); hdelta0->SetLineColor(2); hdelta0->SetFillColor(0);
+
+
+
+*/
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// Time Projection Chamber //
+// Comparison macro for TPC //
+// responsible:
+// marian.ivanov@cern.ch //
+//
+//
+/*
+//
+
+
+.L $ALICE_ROOT/STEER/AliGenInfo.C+
+.L $ALICE_ROOT/STEER/AliTPCComparisonMI.C+
+
+
+TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",-1,1,0);
+t2->Exec();
+
+
+TCut cprim("cprim","MC.fVDist[3]<1");
+TCut csec("cprim","MC.fVDist[3]>1");
+TCut crec("crec","fReconstructed==1");
+TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
+TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
+TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
+TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
+AliTPCComparisonDraw comp;
+comp.SetIO();
+
+
+//example
+comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&fPdg==-211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
+comp.fRes->Draw();
+comp.fMean->Draw();
+comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1,"fReconstructed>0",10,0,1.5)
+comp.fRes->Draw();
+*/
+
+
+
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <string.h>
+//ROOT includes
+#include "Rtypes.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TCut.h"
+#include "TString.h"
+#include "TBenchmark.h"
+#include "TStopwatch.h"
+#include "TParticle.h"
+#include "TSystem.h"
+#include "TTimer.h"
+#include "TVector3.h"
+#include "TPad.h"
+#include "TCanvas.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TText.h"
+#include "Getline.h"
+#include "TStyle.h"
+
+//ALIROOT includes
+#include "AliRun.h"
+#include "AliStack.h"
+#include "AliTPCtrack.h"
+#include "AliSimDigits.h"
+#include "AliTPCParam.h"
+#include "AliTPC.h"
+#include "AliTPCLoader.h"
+#include "AliDetector.h"
+#include "AliTrackReference.h"
+#include "AliRun.h"
+#include "AliTPCParamSR.h"
+#include "AliTracker.h"
+#include "AliComplexCluster.h"
+#include "AliMagF.h"
+#endif
+#include "AliGenInfo.h"
+#include "AliTPCComparisonMI.h"
+
+//
+//
+
+
+
+void AliTPCRecInfo::Update(AliMCInfo* info,AliTPCParam * param, Bool_t reconstructed, Int_t direction)
+{
+ //
+ //
+ //calculates derived variables
+ //
+
+ if (reconstructed==kFALSE){
+ fReconstructed = kFALSE;
+ return;
+ }
+ fReconstructed = kTRUE;
+ // Find the nearest track reference
+ //
+ Double_t radius = TMath::Sqrt(fTrack.GetX()*fTrack.GetX()+fTrack.GetY()*fTrack.GetY());
+ TClonesArray * tpcreferences = info->fTPCReferences;
+ Int_t nentries = tpcreferences->GetEntriesFast();
+ AliTrackReference * ref = 0;
+ Double_t dr = 1000000;
+ for (Int_t i=0;i<nentries;i++){
+ AliTrackReference * refi = (AliTrackReference*)tpcreferences->UncheckedAt(i);
+ if (!refi) continue;
+ //if ( (direction<0) && (refi->R()<radius-10)) continue;
+ //if ( (direction>0) && (refi->R()>radius+10)) continue;
+ if (TMath::Abs(refi->R()-radius)<dr){
+ dr = TMath::Abs(refi->R()-radius);
+ ref = refi;
+ }
+ }
+ if (ref==0) {
+ fReconstructed = kFALSE;
+ return;
+ }
+ //
+ //
+ fdEdx = fTrack.GetdEdx();
+ //
+ // propagate track to nearest reference in direction
+ TVector3 local = TPCCmpTr::TR2Local(ref,param);
+ local.GetXYZ(fTRLocalCoord);
+ //
+ // rotate if neccessary
+ Float_t ymax = local.X()*TMath::Tan(0.5*param->GetInnerAngle());
+ Float_t y = fTrack.GetYat(local.X());
+ if (y > ymax) {
+ fTrack.Rotate(param->GetInnerAngle());
+ } else if (y <-ymax) {
+ fTrack.Rotate(-param->GetInnerAngle());
+ }
+ Double_t xtrack = fTrack.GetX();
+ Double_t delta = (local.X()-xtrack)*0.1;
+ for (Int_t i=0;i<9;i++){
+ fTrack.PropagateTo(xtrack+delta*float(i));
+ }
+ fTrack.PropagateTo(local.X());
+
+ Double_t par[5], cov[15], localX;
+ fTrack.GetExternalParameters(localX,par);
+ fTrack.GetExternalCovariance(cov);
+ //
+ //
+ fRecPhi=TMath::ASin(par[2]) + fTrack.GetAlpha();
+ Double_t phi2 =TMath::ATan2(ref->Py(),ref->Px());
+ if (phi2<0) phi2+=2*TMath::Pi();
+ fGenPhi =phi2;
+ //
+ if (fRecPhi<0) fRecPhi+=2*TMath::Pi();
+ if (fRecPhi>=2*TMath::Pi()) fRecPhi-=2*TMath::Pi();
+ fLambda = TMath::ATan(par[3]);
+ fRecPt_1 = TMath::Abs(par[4]);
+ //TPC
+ //
+ //
+ Double_t phi=TMath::ATan2(par[0],localX) + fTrack.GetAlpha();
+ Double_t r=TMath::Sqrt(localX*localX + par[0]*par[0]);
+ fTPCinR1[0]= r*TMath::Cos(phi);
+ fTPCinR1[1]= r*TMath::Sin(phi);
+ fTPCinR1[2]= par[1];
+ fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
+ fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);
+ //
+ fTPCinP1[0] = fTrack.Px();
+ fTPCinP1[1] = fTrack.Py();
+ fTPCinP1[2] = fTrack.Pz();
+ fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
+ //
+ fTPCinR0[0] = ref->X();
+ fTPCinR0[1] = ref->Y();
+ fTPCinR0[2] = ref->Z();
+ fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
+ fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
+ //
+ fTPCinP0[0]=ref->Px();
+ fTPCinP0[1]=ref->Py();
+ fTPCinP0[2]=ref->Pz();
+ fTPCinP0[3]=ref->Pt();
+ //
+ //
+ if (fTPCinP1[3]>0.0000001){
+ fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
+ fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
+ //
+ fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
+ fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
+ }
+ //
+ fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3]; //delta rfi
+ fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
+ fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]); //delta z
+ fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
+ //
+ fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
+ fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
+ //
+ fTPCDelta[3] = TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]);
+ fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
+ fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
+ Double_t sign = (fTrack.GetC()>0) ? 1:-1;
+ fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]);
+
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+TPCCmpTr::TPCCmpTr()
+{
+ Reset();
+}
+
+////////////////////////////////////////////////////////////////////////
+TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
+ const char* fnCmp,
+ const char* fnGalice, Int_t direction,
+ Int_t nEvents, Int_t firstEvent)
+{
+ Reset();
+ // fFnGenTracks = fnGenTracks;
+ // fFnCmp = fnCmp;
+ sprintf(fFnGenTracks,"%s",fnGenTracks);
+ sprintf(fFnCmp,"%s",fnCmp);
+
+ fFirstEventNr = firstEvent;
+ fEventNr = firstEvent;
+ fNEvents = nEvents;
+ fDirection = direction;
+ //
+ fLoader = AliRunLoader::Open(fnGalice);
+ if (gAlice){
+ //delete gAlice->GetRunLoader();
+ delete gAlice;
+ gAlice = 0x0;
+ }
+ if (fLoader->LoadgAlice()){
+ cerr<<"Error occured while l"<<endl;
+ }
+ Int_t nall = fLoader->GetNumberOfEvents();
+ if (nEvents==0) {
+ nEvents =nall;
+ fNEvents=nall;
+ fFirstEventNr=0;
+ }
+
+ if (nall<=0){
+ cerr<<"no events available"<<endl;
+ fEventNr = 0;
+ return;
+ }
+ if (firstEvent+nEvents>nall) {
+ fEventNr = nall-firstEvent;
+ cerr<<"restricted number of events availaible"<<endl;
+ }
+ AliMagF * magf = gAlice->Field();
+ AliTracker::SetFieldMap(magf);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+TPCCmpTr::~TPCCmpTr()
+{
+ //if (fLoader) {
+ // delete fLoader;
+ //}
+}
+
+//////////////////////////////////////////////////////////////
+Int_t TPCCmpTr::SetIO()
+{
+ //
+ //
+ CreateTreeCmp();
+ if (!fTreeCmp) return 1;
+ fParamTPC = GetTPCParam();
+ //
+ if (!ConnectGenTree()) {
+ cerr<<"Cannot connect tree with generated tracks"<<endl;
+ return 1;
+ }
+ return 0;
+}
+
+//////////////////////////////////////////////////////////////
+
+Int_t TPCCmpTr::SetIO(Int_t eventNr)
+{
+ //
+ //
+ // SET INPUT
+ //gAlice->GetEvent(eventNr);
+ // fLoader->SetEventNumber(eventNr);
+ //fLoader = AliRunLoader::Open("galice.root");
+ fLoader->GetEvent(eventNr);
+ //
+ AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
+ tpcl->LoadTracks();
+ fTreeRecTracks = tpcl->TreeT();
+
+ //
+ return 0;
+}
+
+
+
+////////////////////////////////////////////////////////////////////////
+void TPCCmpTr::Reset()
+{
+ fEventNr = 0;
+ fNEvents = 0;
+ fTreeCmp = 0;
+ //
+ fFileGenTracks = 0;
+ fDebug = 0;
+ //
+ fParamTPC = 0;
+ fTreeRecTracks = 0;
+ fTreePoints =0;
+
+ fTPCTrack = 0;
+ fTracks = 0;
+ fTrackPoints =0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
+{
+ fNEvents = nEvents;
+ fFirstEventNr = firstEventNr;
+ return Exec();
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t TPCCmpTr::Exec()
+{
+ TStopwatch timer;
+ timer.Start();
+
+ if (SetIO()==1)
+ return 1;
+
+ cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
+ for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
+ eventNr++) {
+ fEventNr=eventNr;
+ SetIO(fEventNr);
+ fNParticles = gAlice->GetEvent(fEventNr);
+
+ fIndexRecTracks = new Int_t[fNParticles*4]; //write at maximum 4 tracks corresponding to particle
+ fFakeRecTracks = new Int_t[fNParticles];
+ fMultiRecTracks = new Int_t[fNParticles];
+ for (Int_t i = 0; i<fNParticles; i++) {
+ fIndexRecTracks[4*i] = -1;
+ fIndexRecTracks[4*i+1] = -1;
+ fIndexRecTracks[4*i+2] = -1;
+ fIndexRecTracks[4*i+3] = -1;
+
+ fFakeRecTracks[i] = 0;
+ fMultiRecTracks[i] = 0;
+ }
+
+ cout<<"Start to process event "<<fEventNr<<endl;
+ cout<<"\tfNParticles = "<<fNParticles<<endl;
+ if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
+ if (TreeTLoop(eventNr)>0) return 1;
+
+ if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
+ if (TreeGenLoop(eventNr)>0) return 1;
+ if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
+
+ delete [] fIndexRecTracks;
+ delete [] fFakeRecTracks;
+ delete [] fMultiRecTracks;
+ }
+
+ CloseOutputFile();
+
+ cerr<<"Exec finished"<<endl;
+ timer.Stop();
+ timer.Print();
+ return 0;
+
+}
+////////////////////////////////////////////////////////////////////////
+Bool_t TPCCmpTr::ConnectGenTree()
+{
+//
+// connect all branches from the genTracksTree
+// use the same variables as for the new cmp tree, it may work
+//
+ fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
+ if (!fFileGenTracks) {
+ cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
+ return kFALSE;
+ }
+ fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
+ if (!fTreeGenTracks) {
+ cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
+ <<fFnGenTracks<<endl;
+ return kFALSE;
+ }
+ //
+ fMCInfo = new AliMCInfo;
+ fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
+
+
+ if (fDebug > 1) {
+ cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
+ }
+ return kTRUE;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+void TPCCmpTr::CreateTreeCmp()
+{
+ fFileCmp = TFile::Open(fFnCmp,"RECREATE");
+ if (!fFileCmp) {
+ cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
+ return;
+ }
+
+
+ fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
+ //
+ fMCInfo = new AliMCInfo;
+ fRecInfo = new AliTPCRecInfo;
+ //
+ fTPCTrack = new AliTPCtrack;
+ //
+ fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
+ fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
+ fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
+ //
+ fTreeCmp->AutoSave();
+
+}
+////////////////////////////////////////////////////////////////////////
+void TPCCmpTr::CloseOutputFile()
+{
+ if (!fFileCmp) {
+ cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
+ return;
+ }
+ fFileCmp->cd();
+ fTreeCmp->Write();
+ delete fTreeCmp;
+
+ fFileCmp->Close();
+ delete fFileCmp;
+ return;
+}
+////////////////////////////////////////////////////////////////////////
+
+TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
+ AliTPCParam *paramTPC) {
+
+ Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
+ Int_t index[4];
+ paramTPC->Transform0to1(x,index);
+ paramTPC->Transform1to2(x,index);
+ return TVector3(x);
+}
+////////////////////////////////////////////////////////////////////////
+
+Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
+{
+ //
+ // loop over all TPC reconstructed tracks and store info in memory
+ //
+ TStopwatch timer;
+ timer.Start();
+
+ if (!fTreeRecTracks) {
+ cerr<<"Can't get a tree with TPC rec. tracks "<<endl;
+ return 1;
+ }
+ //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
+
+ Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
+ if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
+ <<nEntries<<endl;
+ TBranch * br= fTreeRecTracks->GetBranch("tracks");
+ br->SetAddress(&fTPCTrack);
+ TBranch *brp = 0;
+ if (fTreePoints) brp = fTreePoints->GetBranch("debug");
+
+ if (fTracks){
+ fTracks->Delete();
+ delete fTracks;
+ }
+ if (fTrackPoints){
+ fTrackPoints->Delete();
+ delete fTrackPoints;
+ fTrackPoints =0;
+ }
+ fTracks = new TObjArray(nEntries);
+ if (brp){
+ fTrackPoints = new TObjArray(nEntries);
+ }
+ else fTrackPoints = 0;
+
+ //
+ //load tracks to the memory
+ for (Int_t i=0; i<nEntries;i++){
+ AliTPCtrack * track = new AliTPCtrack;
+ br->SetAddress(&track);
+ br->GetEntry(i);
+ fTracks->AddAt(track,i);
+ }
+ //
+ //load track points to the memory
+ if (brp) for (Int_t i=0; i<nEntries;i++){
+ TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
+ brp->SetAddress(&arr);
+ brp->GetEntry(i);
+ if (arr!=0)
+ for (Int_t j=0;j<arr->GetEntriesFast();j++){
+ AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
+ if (point && point->fID>=0){
+ fTrackPoints->AddAt(arr,point->fID);
+ break;
+ }
+ }
+ }
+ //
+
+ //
+ for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
+ //br->GetEntry(iEntry);
+ fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
+ //
+ Int_t label = fTPCTrack->GetLabel();
+ Int_t absLabel = abs(label);
+ if (absLabel < fNParticles) {
+ // fIndexRecTracks[absLabel] = iEntry;
+ if (label < 0) fFakeRecTracks[absLabel]++;
+
+ if (fMultiRecTracks[absLabel]>0){
+ if (fMultiRecTracks[absLabel]<4)
+ fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] = iEntry;
+ }
+ else
+ fIndexRecTracks[absLabel*4] = iEntry;
+ fMultiRecTracks[absLabel]++;
+ }
+ }
+ printf("Time spended in TreeTLoop\n");
+ timer.Print();
+
+ if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
+
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
+{
+//
+// loop over all entries for a given event, find corresponding
+// rec. track and store in the fTreeCmp
+//
+ TStopwatch timer;
+ timer.Start();
+ Int_t entry = 0;
+ Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
+ TBranch * branch = fTreeCmp->GetBranch("RC");
+ branch->SetAddress(&fRecInfo); // set all pointers
+
+ while (entry < nParticlesTR) {
+ fTreeGenTracks->GetEntry(entry);
+ entry++;
+ if (eventNr < fMCInfo->fEventNr) continue;
+ if (eventNr > fMCInfo->fEventNr) continue;
+ //
+ if (fDebug > 2 && fMCInfo->fLabel < 10) {
+ cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
+ }
+ fRecInfo->Reset();
+ //
+ if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
+ fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
+
+ // if (nBytes > 0) {
+ if (fTPCTrack) {
+ //
+ //
+ //if multifound find track with biggest pt - returning tracks - loos energy
+ //
+ if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
+ Double_t pt = TMath::Abs(1/fTPCTrack->Get1Pt());
+ for (Int_t i=1;i<4;i++){
+ if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
+ AliTPCtrack * track = (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
+ Double_t pt2 = TMath::Abs(1/track->Get1Pt());
+ if (pt2>pt)
+ fTPCTrack = track;
+ }
+ }
+ }
+ //
+ fRecInfo->fTP=0;
+ fRecInfo->fTrack =*fTPCTrack;
+ fRecInfo->fReconstructed = 1;
+ fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
+ fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
+ fRecInfo->Update(fMCInfo, fParamTPC, kTRUE, 1);
+ }
+ }
+ else{
+ fRecInfo->Update(fMCInfo, fParamTPC, kFALSE, 1);
+ }
+ fTreeCmp->Fill();
+ }
+ fTreeCmp->AutoSave();
+ fTracks->Delete();
+ fTPCTrack =0;
+ if (fTrackPoints){
+ fTrackPoints->Delete();
+ delete fTrackPoints;
+ fTrackPoints =0;
+ }
+ printf("Time spended in TreeGenLoop\n");
+ timer.Print();
+ if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
+ return 0;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
+
+void AliTPCComparisonDraw::SetIO(const char *fname)
+{
+ TFile *file = TFile::Open(fname);
+ if (!file){
+ printf("Comparison file couldn't be generated\n");
+ return;
+ }
+ //
+ fTree = (TTree*) file->Get("TPCcmpTracks");
+ if (!fTree) {
+ printf("no track comparison tree found\n");
+ file->Close();
+ delete file;
+ }
+
+}
+
+void AliTPCComparisonDraw::Eff()
+{
+
+}
+
+void AliTPCComparisonDraw::ResPt()
+{
+}
--- /dev/null
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class digitRow
+//
+////////////////////////////////////////////////////////////////////////
+const Int_t kgRowBytes = 32;
+
+class digitRow: public TObject {
+
+public:
+ digitRow();
+ virtual ~digitRow(){;}
+ void SetRow(Int_t row);
+ Bool_t TestRow(Int_t row);
+ digitRow & operator=(const digitRow &digOld);
+ Int_t RowsOn(Int_t upto=8*kgRowBytes);
+ Int_t Last();
+ Int_t First();
+ void Reset();
+
+//private:
+ UChar_t fDig[kgRowBytes];
+
+ ClassDef(digitRow,1) // container for digit pattern
+};
+ClassImp(digitRow)
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class AliTPCGenInfo
+//
+////////////////////////////////////////////////////////////////////////
+
+class AliTPCGenInfo: public TObject {
+
+public:
+ AliTPCGenInfo();
+ ~AliTPCGenInfo();
+
+ AliTrackReference fTrackRef; // track reference saved in the output tree
+ AliTrackReference fTrackRefOut; // decay track reference saved in the output tree
+ TParticle fParticle; // generated particle
+ Int_t fLabel; // track label
+ Int_t fEventNr; // event number
+
+ Float_t fDecayCoord[3]; // position of particle decay
+ Double_t fVDist[4]; //distance of the particle vertex from primary vertex
+
+ 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
+ Int_t fDigitsInSeed; // digits in the default seed rows
+ Float_t fPrim; // theoretical dedx in tpc according particle momenta and mass
+ digitRow fRow; // information about digits row pattern
+ TClonesArray * fReferences; //containner with all track references
+ ClassDef(AliTPCGenInfo,1) // container for
+};
+ClassImp(AliTPCGenInfo)
+
+
+
+class AliTPCGenV0Info: public TObject {
+public:
+ AliTPCGenInfo fMCd; //info about daughter particle
+ AliTPCGenInfo fMCm; //info about mother particle
+ void Update(); // put some derived info to special field
+ Double_t fDist1; //info about closest distance according closest MC - linear DCA
+ Double_t fDist2; //info about closest distance parabolic DCA
+ //
+ Double_t fPdr[3]; //momentum at vertex daughter - according approx at DCA
+ Double_t fPd[4]; //exact momentum from MC info
+ Double_t fX[3]; //exact position of the vertex
+ Double_t fXr[3]; //rec. position according helix
+ //
+ Double_t fPm[3]; //momentum at the vertex mother
+ Double_t fAngle[3]; //three angels
+ Double_t fRr; // rec position of the vertex
+ Double_t fR; //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(AliTPCGenV0Info,1) // container for
+};
+ClassImp(AliTPCGenV0Info)
+
+
+
+void AliTPCGenV0Info::Update()
+{
+ fPd[0] = fMCd.fParticle.Px();
+ fPd[1] = fMCd.fParticle.Py();
+ fPd[2] = fMCd.fParticle.Pz();
+ fPd[3] = fMCd.fParticle.P();
+ fX[0] = fMCd.fParticle.Vx();
+ fX[1] = fMCd.fParticle.Vy();
+ fX[2] = fMCd.fParticle.Vz();
+ fR = TMath::Sqrt( fX[0]*fX[0]+
+ fX[1]*fX[1]);
+ fPdg[0] = fMCd.fParticle.GetPdgCode();
+ fPdg[1] = fMCm.fParticle.GetPdgCode();
+ //
+ fLab[0] = fMCd.fParticle.GetUniqueID();
+ fLab[1] = fMCm.fParticle.GetUniqueID();
+
+}
+
+
+
+
+ /////////////////////////////////////////////////////////////////////////
+class AliTPCRecInfo: public TObject {
+
+public:
+ AliTPCRecInfo(){fTP = new TClonesArray("AliTPCTrackPoint2",0);}
+ ~AliTPCRecInfo(){if (fTP) {fTP->Delete();delete fTP;}}
+ //
+ AliTPCtrack fTPCTrack; // tpc track
+ Float_t fTRLocalCoord[3]; //local coordinates of the track ref.
+ Int_t fReconstructed; //flag if track was reconstructed
+ Double_t fRecPhi; // reconstructed phi angle (0;2*kPI)
+ Double_t fLambda; // reconstructed
+ Double_t fRecPt_1; // reconstructed
+ Float_t fdEdx; // reconstructed dEdx
+ Int_t fFake; // fake track
+ Int_t fMultiple; // number of reconstructions
+ TClonesArray *fTP; //container with track points
+ void Reset();
+ //
+ ClassDef(AliTPCRecInfo,1) // container for
+};
+ClassImp(AliTPCRecInfo)
+
+void AliTPCRecInfo::Reset()
+{
+ fMultiple =0;
+ fFake =0;
+ fReconstructed=0;
+ fRecPhi =0;
+ fLambda =0;
+}
+
+
+/////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////
+
+
+class AliTPCRecV0Info: public TObject {
+public:
+ AliTPCRecInfo fT1; //track1
+ AliTPCRecInfo 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 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
+ Int_t fLab[2]; //MC label of the partecle
+ ClassDef(AliTPCRecV0Info,1) // container for
+};
+
+ClassImp(AliTPCRecV0Info)
+
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class TPCFindGenTracks
+//
+////////////////////////////////////////////////////////////////////////
+
+class TPCFindGenTracks {
+
+public:
+ TPCFindGenTracks();
+ TPCFindGenTracks(char* fnHits,
+ char* fnDigits ="tpc.digits.root",
+ char* fnRes ="genTracks.root",
+ Int_t nEvents=1, Int_t firstEvent=0);
+ virtual ~TPCFindGenTracks();
+ 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();
+ //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 SetIO();
+ Float_t TR2LocalX(AliTrackReference *trackRef,
+ AliTPCParam *paramTPC);
+
+public:
+ AliTPCGenInfo* fMCInfo; //! information writen per particle
+ Int_t fDebug; //! debug flag
+ Int_t fEventNr; //! current event number
+ Int_t fLabel; //! track label
+ Int_t fNEvents; //! number of events to process
+ Int_t fFirstEventNr; //! first event to process
+ Int_t fNParticles; //! number of particles in TreeK
+ TTree *fTreeGenTracks; //! output tree with generated tracks
+ char *fFnRes; //! output file name with stored tracks
+ char *fFnHits; //! input file name with hits
+ char *fFnDigits; //! input file name with digits
+ TFile *fFileGenTracks; //! output file with stored fTreeGenTracks
+ TFile *fFileHits; //! input file with hits
+ TFile *fFileTreeD; //! input file with digits
+ //
+ TTree * fTreeD; //! current tree with digits
+ TTree * fTreeTR; //! current tree with TR
+ AliStack *fStack; //! current stack
+ //
+ digitRow *fContainerDigitRow; //! big container for partial information
+ //
+ AliTrackReference *fReferences; //! container with track references
+ Int_t *fReferenceIndex0; //! first index for given track
+ Int_t *fReferenceIndex1; //! last index for given track
+ //
+ AliTPCParam* fParamTPC; //! AliTPCParam
+ Double_t fVPrim[3]; //! primary vertex position
+ // the fVDist[3] contains size of the 3-vector
+
+private:
+
+// some constants for the original non-pareller tracking (by Y.Belikov)
+ static const Int_t seedRow11 = 158; // nRowUp - 1
+ static const Int_t seedRow12 = 139; // nRowUp - 1 - (Int_t) 0.125*nRowUp
+ static const Int_t seedRow21 = 149; // seedRow11 - shift
+ static const Int_t seedRow22 = 130; // seedRow12 - shift
+ static const Double_t kRaddeg = 180./kPI;
+
+ static const Int_t fgMaxIndexTR = 50000; // maximum number of tracks with a track ref
+ static const Int_t fgMaxTR = 1000000; // maximum number of track refs
+
+ static const Int_t fgMaxParticles = 2000000; // maximum number of generated particles
+ static const Double_t fgPtCut = .1; // do not store particles with generated pT less than this
+ static const Float_t fgTrackRefLocalXMax = 82.95;
+ static const Float_t fgTrackRefLocalXMaxDelta = 5.;
+
+ ClassDef(TPCFindGenTracks,1) // class which creates and fills tree with TPCGenTrack objects
+};
+ClassImp(TPCFindGenTracks)
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Start of implementation of the class TPCCmpTr
+//
+////////////////////////////////////////////////////////////////////////
+
+class TPCCmpTr {
+
+public:
+ TPCCmpTr();
+ TPCCmpTr(char* fnRecTracks,
+ char* fnGenTracks ="genTracks.root",
+ char* fnCmpRes ="cmpTracks.root",
+ char* fnGalice ="galice.root",
+ Int_t nEvents=1, Int_t firstEvent=0);
+ virtual ~TPCCmpTr();
+ void Reset();
+ Int_t Exec();
+ Int_t Exec(Int_t nEvents, Int_t firstEventNr);
+ void CreateTreeCmp();
+ void CloseOutputFile();
+ Bool_t ConnectGenTree();
+ Int_t TreeGenLoop(Int_t eventNr);
+ Int_t TreeTLoop(Int_t eventNr);
+ 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 TrackReferenceTPC
+ 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; //! output file name with cmp tracks
+ TFile *fFileCmp; //! output file with cmp tracks
+ TTree *fTreeCmp; //! output tree with cmp tracks
+ //
+ char *fFnGenTracks; //! input file name with gen tracks
+ TFile *fFileGenTracks;
+ TTree *fTreeGenTracks;
+ //
+ char *fFnHits; //! input file name with gAlice object (needed for B)
+ TFile *fFileHits; //! input file with gAlice
+ //
+ char *fFnRecTracks; //! input file name with tpc rec. tracks
+ TFile *fFileRecTracks; //! input file with reconstructed tracks
+ TTree *fTreeRecTracks; //! tree with reconstructed tracks
+ TTree *fTreePoints; //! tree with reconstructed points
+ //
+ Int_t *fIndexRecTracks; //! index of particle label in the TreeT_TPC
+ Int_t *fFakeRecTracks; //! number of fake tracks
+ Int_t *fMultiRecTracks; //! number of multiple reconstructions
+ //
+ TObjArray * fTracks; //!container with tracks
+ TObjArray * fTrackPoints; //! container with track points
+ //
+ 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
+ //
+ AliTPCGenInfo* fMCInfo; //! MC information writen per particle
+ AliTPCRecInfo* fRecInfo; //! Rec. information writen per particle
+ //
+ AliTPCtrack *fTPCTrack; //! pointer to TPC track to connect branch
+
+ ClassDef(TPCCmpTr,1) // class which creates and fills tree with TPCGenTrack objects
+};
+ClassImp(TPCCmpTr)
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff