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