+++ /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
-/**************************************************************************
- * 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:
-
-
-
-gSystem->Load("libPWG1.so");
-//
-AliESDComparisonMI *t2 = new AliESDComparisonMI("genTracks.root","cmpESDTracks.root","galice.root",-1,0,0);
-t2->Exec();
-
-
-TFile f("cmpESDTracks.root");
-AliTreeDraw comp;
-TTree * tree = (TTree*)f.Get("ESDcmpTracks");
-comp.SetTree(tree)
-
-
-
-//
-//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.")
-
-
-//
-//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);
-
-
-
-
-*/
-
-
-#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 "AliTPCseed.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"
-//
-#include "AliTreeDraw.h"
-#include "AliGenInfo.h"
-#include "AliESDComparisonMI.h"
-
-
-
-ClassImp(AliESDRecInfo)
-ClassImp(AliESDRecV0Info)
-ClassImp(AliESDRecKinkInfo)
-ClassImp(AliESDComparisonMI)
-
-
-
-
-
-void MakeAliases(AliTreeDraw&comp)
-{
- //
- // aliases definition
- //
- comp.T()->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
- comp.T()->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
- comp.T()->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
- comp.T()->SetAlias("theta","MC.fTrackRef.Theta()");
- comp.T()->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
- comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
- comp.T()->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
-
- comp.T()->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.T()->SetAlias("dtofmc2","fESDtrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
- comp.T()->SetAlias("dtofrc2","(fESDtrack.fTrackTime[2]-fESDtrack.fTOFsignal)");
-
- comp.T()->SetAlias("psum","fESDtrack.fTOFr[4]+fESDtrack.fTOFr[3]+fESDtrack.fTOFr[2]+fESDtrack.fTOFr[1]+fESDtrack.fTOFr[0]");
- comp.T()->SetAlias("P0","fESDtrack.fTOFr[0]/psum");
- comp.T()->SetAlias("P1","fESDtrack.fTOFr[1]/psum");
- comp.T()->SetAlias("P2","fESDtrack.fTOFr[2]/psum");
- comp.T()->SetAlias("P3","fESDtrack.fTOFr[3]/psum");
- comp.T()->SetAlias("P4","fESDtrack.fTOFr[4]/psum");
- comp.T()->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
-}
-
-
-AliESDRecInfo::AliESDRecInfo():
- fITSOn(0), // ITS refitted inward
- fTRDOn(0), // ITS refitted inward
- fDeltaP(0), //delta of momenta
- fSign(0), // sign
- fReconstructed(0), //flag if track was reconstructed
- fFake(0), // fake track
- fMultiple(0), // number of reconstructions
- fTPCOn(0), // TPC refitted inward
- fBestTOFmatch(0), //best matching between times
- fESDtrack(0), // esd track
- fTrackF(0), // friend track
- fTPCtrack(0), // tpc track
- fITStrack(0), // its track
- fTRDtrack(0) // trd track
-{
- //
- // default constructor
- //
-}
-
-
-AliESDRecInfo::AliESDRecInfo(const AliESDRecInfo& recinfo):
- TObject()
-{
- //
- //
- //
- memcpy(this,&recinfo, sizeof(recinfo));
- fESDtrack=0; fTrackF=0; fTPCtrack=0;fITStrack=0;fTRDtrack=0;
- SetESDtrack(recinfo.GetESDtrack());
-}
-
-
-AliESDRecInfo::~AliESDRecInfo()
-
-{
- //
- // destructor
- //
- if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
- if (fTrackF) { delete fTrackF; fTrackF=0;}
- if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
- if (fITStrack) { delete fITStrack; fITStrack=0;}
- if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
-
-}
-
-
-
-void AliESDRecInfo::Reset()
-{
- //
- // reset info
- //
- fMultiple =0;
- fFake =0;
- fReconstructed=0;
- if (fESDtrack) { delete fESDtrack; fESDtrack=0;}
- if (fTrackF) { delete fTrackF; fTrackF=0;}
- if (fTPCtrack) { delete fTPCtrack; fTPCtrack=0;}
- if (fITStrack) { delete fITStrack; fITStrack=0;}
- if (fTRDtrack) { delete fTRDtrack; fTRDtrack=0;}
-}
-
-void AliESDRecInfo::SetESDtrack(const AliESDtrack *track){
- //
- //
- //
- if (fESDtrack) delete fESDtrack;
- fESDtrack = (AliESDtrack*)track->Clone();
- if (track->GetFriendTrack()){
- if (fTrackF) delete fTrackF;
- fTrackF = (AliESDfriendTrack*)track->GetFriendTrack()->Clone();
- }
-
-}
-
-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;
-
-
-}
-
-//
-//
-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.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
- fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
- AliExternalTrackParam paramm(x,alpha,param,cov);
- //
- fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
- fT2.GetESDtrack()->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.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
- fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
- AliExternalTrackParam paramm(x,alpha,param,cov);
- //
- fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
- fT2.GetESDtrack()->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);
- */
- }
-
- //
- //
- }
-
-}
-
-
-////////////////////////////////////////////////////////////////////////
-AliESDComparisonMI::AliESDComparisonMI()
-{
- Reset();
-}
-
-////////////////////////////////////////////////////////////////////////
-AliESDComparisonMI::AliESDComparisonMI(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);
-
-}
-
-
-////////////////////////////////////////////////////////////////////////
-AliESDComparisonMI::~AliESDComparisonMI()
-{
- if (fLoader) {
- delete fLoader;
- }
-}
-
-//////////////////////////////////////////////////////////////
-Int_t AliESDComparisonMI::SetIO()
-{
- //
- //
- CreateTreeCmp();
- if (!fTreeCmp) return 1;
- fParamTPC = GetTPCParam();
- //
- if (!ConnectGenTree()) {
- cerr<<"Cannot connect tree with generated tracks"<<endl;
- return 1;
- }
- return 0;
-}
-
-//////////////////////////////////////////////////////////////
-
-Int_t AliESDComparisonMI::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 AliESDComparisonMI::Reset()
-{
- fEventNr = 0;
- fNEvents = 0;
- fTreeCmp = 0;
- fTreeCmpKinks =0;
- fTreeCmpV0 =0;
- // fFnCmp = "cmpTracks.root";
- fFileGenTracks = 0;
- fDebug = 0;
- //
- fParamTPC = 0;
- fEvent =0;
-}
-
-////////////////////////////////////////////////////////////////////////
-Int_t AliESDComparisonMI::Exec(Int_t nEvents, Int_t firstEventNr)
-{
- fNEvents = nEvents;
- fFirstEventNr = firstEventNr;
- return Exec();
-}
-
-////////////////////////////////////////////////////////////////////////
-Int_t AliESDComparisonMI::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); // no V0 info for a moment
- 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 AliESDComparisonMI::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 AliESDComparisonMI::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 AliESDComparisonMI::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 AliESDComparisonMI::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 AliESDComparisonMI::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;
- }
- //
- for (Int_t i=0; i<nV0MIs;i++){
- AliV0 * v0MI = (AliV0*)fEvent->GetV0(i);
- fSignedV0[i]=0;
- }
-
- //
- //
- 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->GetLabel(0));
-// Int_t label1 = TMath::Abs(v0MI->GetLabel(1));
-// //
-// for (Int_t i=0;i<2;i++){
-// Int_t absLabel = TMath::Abs(v0MI->GetLabel(i));
-// 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 AliESDComparisonMI::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) {
- fRecInfo->SetESDtrack(track);
- }else{
- fRecInfo->SetESDtrack(&dummytrack);
- }
- //
-
- fRecInfo->fReconstructed = 1;
- fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
- fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
- //
- fRecInfo->Update(fMCInfo,fParamTPC,kTRUE, fEvent);
- }
- else{
- fRecInfo->SetESDtrack(&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 AliESDComparisonMI::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->GetLabel(0))==label &&TMath::Abs(kink2->GetLabel(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->GetLabel(1))==label &&TMath::Abs(kink2->GetLabel(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->GetLabel(0)));
- AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->GetLabel(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 AliESDComparisonMI::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.GetESDtrack()->GetInnerParam() && fRecV0Info->fT2.GetESDtrack()->GetInnerParam()){
- fRecV0Info->fT1.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
- fRecV0Info->fT1.GetESDtrack()->GetInnerExternalCovariance(cov);
- AliExternalTrackParam paramP(x,alpha,param,cov);
- //
- fRecV0Info->fT2.GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
- fRecV0Info->fT2.GetESDtrack()->GetInnerExternalCovariance(cov);
- AliExternalTrackParam paramM(x,alpha,param,cov);
- //
- fRecV0Info->fV0tpc->SetParamN(paramM);
- fRecV0Info->fV0tpc->SetParamP(paramP);
- Double_t pid1[5],pid2[5];
- fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid1);
- fRecV0Info->fT1.GetESDtrack()->GetESDpid(pid2);
- //
- //fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
- fRecV0Info->fV0tpc->Update(vertex);
-
- //
- //
- fRecV0Info->fT1.GetESDtrack()->GetExternalParameters(x,param);
- fRecV0Info->fT1.GetESDtrack()->GetExternalCovariance(cov);
- alpha = fRecV0Info->fT1.GetESDtrack()->GetAlpha();
- new (¶mP) AliExternalTrackParam(x,alpha,param,cov);
- //
- fRecV0Info->fT2.GetESDtrack()->GetExternalParameters(x,param);
- fRecV0Info->fT2.GetESDtrack()->GetExternalCovariance(cov);
- alpha = fRecV0Info->fT2.GetESDtrack()->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->GetLabel(0))==label &&TMath::Abs(v0MI2->GetLabel(1))==label2) {
-// v0MI =v0MI2;
-// fRecV0Info->fV0Multiple++;
-// fSignedV0[index]=1;
-// }
-// if (TMath::Abs(v0MI2->GetLabel(1))==label &&TMath::Abs(v0MI2->GetLabel(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->GetLabel(0)));
-// AliESDRecInfo* fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->GetLabel(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;
-}
-////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////
-
-
+++ /dev/null
-#ifndef ALIESDCOMPARISONMI_H
-#define ALIESDCOMPARISONMI_H
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice */
-
-
-
-//////////////////////////////////////////////////////////////////////////////
-// Class AliESDComparisonMI //
-// collect together MC info and Rec info for comparison purposes
-// - effieciency studies and so on// //
-// marian.ivanov@cern.ch //
-//////////////////////////////////////////////////////////////////////////////
-
-
-#include "TObject.h"
-#include "AliESD.h"
-
-#include "AliESDtrack.h"
-#include "AliV0.h"
-#include "AliESDfriendTrack.h"
-#include "AliITStrackMI.h"
-#include "AliTRDtrack.h"
-class AliTPCseed;
-
-/////////////////////////////////////////////////////////////////////////
-class AliESDRecInfo: public TObject {
-
-public:
- AliESDRecInfo();
- AliESDRecInfo(const AliESDRecInfo& recinfo);
- ~AliESDRecInfo();
- void UpdatePoints(AliESDtrack* track);
- void Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed, AliESD *event);
- void Reset();
- //
- void SetESDtrack(const AliESDtrack *track);
- AliESDtrack *GetESDtrack() const { return fESDtrack;}
- AliESDfriendTrack *GetTrackF() const { return fTrackF;}
- AliTPCseed *GetTPCtrack() const { return fTPCtrack;}
- AliITStrackMI *GetITStrack() const { return fITStrack;}
- AliTRDtrack *GetTRDtrack() const { return fTRDtrack;}
- //
- 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
- Float_t fTRLocalCoord[3]; //local coordinates of the track ref.
- Int_t fStatus[4]; // status -0 not found - 1 -only in - 2 -in-out -3 -in -out-refit
- Int_t fLabels[2]; // labels
-
- Bool_t fITSOn; // ITS refitted inward
- Bool_t fTRDOn; // ITS refitted inward
- Float_t fDeltaP; //delta of momenta
- Double_t fSign; // sign
- 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
- Float_t fBestTOFmatch; //best matching between times
-
-private:
- AliESDtrack *fESDtrack; // esd track
- AliESDfriendTrack *fTrackF; // friend track
- AliTPCseed *fTPCtrack; // tpc track
- AliITStrackMI *fITStrack; // its track
- AliTRDtrack *fTRDtrack; // trd track
-
- ClassDef(AliESDRecInfo,2) // container for
-};
-
-
-
-/////////////////////////////////////////////////////////
-/////////////////////////////////////////////////////////
-/////////////////////////////////////////////////////////
-
-
-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
-};
-
-
-
-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
-};
-
-
-
-
-////////////////////////////////////////////////////////////////////////
-//
-// Start of implementation of the class AliESDComparisonMI
-//
-////////////////////////////////////////////////////////////////////////
-
-class AliESDComparisonMI {
-
-public:
- AliESDComparisonMI();
- AliESDComparisonMI(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 ~AliESDComparisonMI();
- 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(AliESDComparisonMI,1) // class which creates and fills tree with ESDGenTrack objects
-};
-
-
-#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff