From 54709d62a7d99bef6d883808e7a1544839187655 Mon Sep 17 00:00:00 2001 From: marian Date: Mon, 10 Aug 2009 08:05:32 +0000 Subject: [PATCH] Adding new AnalysisTask for evaluation of meterial budget For the moment only MC part implemented (test of principle) To be extended --- PWG1/AliMaterialBudget.cxx | 1052 ++++++++++++++++++++++++++++++++++++ PWG1/AliMaterialBudget.h | 93 ++++ PWG1/PWG1LinkDef.h | 2 + PWG1/libPWG1.pkg | 3 +- 4 files changed, 1149 insertions(+), 1 deletion(-) create mode 100644 PWG1/AliMaterialBudget.cxx create mode 100644 PWG1/AliMaterialBudget.h diff --git a/PWG1/AliMaterialBudget.cxx b/PWG1/AliMaterialBudget.cxx new file mode 100644 index 00000000000..6152c964e2c --- /dev/null +++ b/PWG1/AliMaterialBudget.cxx @@ -0,0 +1,1052 @@ + + + + +// +// + +// This class estiamtes the material budget of the inner detectors in ALICE based +// on the "upper/lower track matching"-method of the ALICE TPC. + +// +// +// + +// + + + + +// ROOT includes +#include +#include +#include +#include +#include +#include "TGeoGlobalMagField.h" + +// ALIROOT includes +#include +#include +#include +#include "AliStack.h" +#include "AliMCEvent.h" +#include "AliMCEventHandler.h" + +#include +#include "AliMaterialBudget.h" +#include "AliGenInfoMaker.h" +#include "AliHelix.h" + +// +#include "AliMCInfo.h" +#include "AliComparisonObject.h" +#include "AliESDRecInfo.h" +#include "AliTPCParamSR.h" +#include "AliTracker.h" +#include "AliTPCseed.h" + +// STL includes +#include + +using namespace std; + +ClassImp(AliMaterialBudget) + +//________________________________________________________________________ +AliMaterialBudget::AliMaterialBudget() : + AliAnalysisTask(), + fMCinfo(0), //! MC event handler + fESD(0), + fDebugStreamer(0), + fStreamLevel(0), + fDebugLevel(0), + fDebugOutputPath(), + fListHist(0), + fHistMult(0), + fCutMaxD(5), // maximal distance in rfi ditection + fCutMaxDz(40), // maximal distance in z ditection + fCutTheta(0.03), // maximal distan theta + fCutMinDir(-0.99) // direction vector products +{ + // + // Default constructor (should not be used) + // +} + +AliMaterialBudget::AliMaterialBudget(const AliMaterialBudget& /*info*/) : + AliAnalysisTask(), + fMCinfo(0), //! MC event handler + fESD(0), + // + fDebugStreamer(0), + fStreamLevel(0), + fDebugLevel(), + fDebugOutputPath(), + fListHist(0), + fHistMult(0), + fCutMaxD(5), // maximal distance in rfi ditection + fCutMaxDz(40), // maximal distance in z ditection + fCutTheta(0.03), // maximal distan theta + fCutMinDir(-0.99) // direction vector products +{ + // + // Default constructor + // +} + + + +//________________________________________________________________________ +AliMaterialBudget::AliMaterialBudget(const char *name) : + AliAnalysisTask(name, "AliMaterialBudget"), + fMCinfo(0), //! MC event handler + fESD(0), + fDebugStreamer(0), + fStreamLevel(0), + fDebugLevel(0), + fDebugOutputPath(), + fListHist(0), + fHistMult(0), + fCutMaxD(5), // maximal distance in rfi ditection + fCutMaxDz(40), // maximal distance in z ditection + fCutTheta(0.03), // maximal distan theta + fCutMinDir(-0.99) // direction vector products +{ + // + // Normal constructor + // + // Input slot #0 works with a TChain + DefineInput(0, TChain::Class()); + // Output slot #0 writes into a TList + DefineOutput(0, TList::Class()); + // + // +} + +AliMaterialBudget::~AliMaterialBudget(){ + // + // + // + if (fDebugLevel>0) printf("AliMaterialBudget::~AliMaterialBudget\n"); + if (fDebugStreamer) delete fDebugStreamer; + fDebugStreamer=0; +} + + +//________________________________________________________________________ +void AliMaterialBudget::ConnectInputData(Option_t *) +{ + // + // Connect the input data + // + if(fDebugLevel>3) + cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl; + + TTree* tree=dynamic_cast(GetInputData(0)); + if (!tree) { + //Printf("ERROR: Could not read chain from input slot 0"); + } + else { + AliESDInputHandler *esdH = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); + if (!esdH) { + //Printf("ERROR: Could not get ESDInputHandler"); + } + else { + esdH->SetActiveBranches("ESDfriend"); + fESD = esdH->GetEvent(); + //Printf("*** CONNECTED NEW EVENT ****"); + } + } + AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); + mcinfo->SetReadTR(kTRUE); + + fMCinfo = mcinfo->MCEvent(); + + +} + + + + + + +//________________________________________________________________________ +void AliMaterialBudget::CreateOutputObjects() +{ + // + // Connect the output objects + // + if(fDebugLevel>3) + cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl; + // + fListHist = new TList(); + fListHist->SetOwner(); + // + fHistMult = new TH1F("HistMult", "Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5); + fListHist->Add(fHistMult); + + +} + + +//________________________________________________________________________ +void AliMaterialBudget::Exec(Option_t *) { + // + // Execute analysis for current event + // + + if(fDebugLevel>3) + cout << "AliMaterialBudget::Exec()" << endl; + + fHistMult->Fill(fESD->GetNumberOfTracks()); + //FindPairs(fESD); // nearly everything takes place in find pairs... + + // If MC has been connected + + if (!fMCinfo){ + cout << "Not MC info\n" << endl; + }else{ + ProcessMCInfo(); + //mcinfo->Print(); + //DumpInfo(); + } + // + PostData(0, fListHist); +} + + + + +//________________________________________________________________________ +void AliMaterialBudget::Terminate(Option_t *) { + // + // Terminate loop + // + if(fDebugLevel>3) + printf("AliMaterialBudget: Terminate() \n"); + // + if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n"); + if (fDebugStreamer) delete fDebugStreamer; + fDebugStreamer = 0; + return; +} + + + +TTreeSRedirector *AliMaterialBudget::GetDebugStreamer(){ + // + // Get Debug streamer + // In case debug streamer not yet initialized and StreamLevel>0 create new one + // + if (fStreamLevel==0) return 0; + if (fDebugStreamer) return fDebugStreamer; + TString dsName; + dsName=GetName(); + dsName+="Debug.root"; + dsName.ReplaceAll(" ",""); + fDebugStreamer = new TTreeSRedirector(dsName.Data()); + return fDebugStreamer; +} + + + + + + +AliExternalTrackParam * AliMaterialBudget::MakeTrack(const AliTrackReference* ref, TParticle*part) +{ + // + // Make track out of the track ref + // part - TParticle used to determine chargr + // the covariance matrix - equal 0 - starting from ideal MC position + if (!part->GetPDG()) return 0; + Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()}; + Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()}; + Double_t charge = TMath::Nint(part->GetPDG()->Charge()/3.); + if (ref->X()*ref->Px()+ref->Y()*ref->Py() <0){ + pxyz[0]*=-1; + pxyz[1]*=-1; + pxyz[2]*=-1; + charge*=-1.; + } + Double_t cv[21]; + for (Int_t i=0; i<21;i++) cv[i]=0; + AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,charge); + return param; +} + +Bool_t AliMaterialBudget::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){ + // + // Propagate track to point xyz using + // AliTracker::PropagateTo functionality + // + // param - track parameters + // xyz - position to propagate + // mass - particle mass + // step - step to be used + Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); + AliTracker::PropagateTrackToBxByBz(param, radius+step, mass, step, kTRUE,0.99); + AliTracker::PropagateTrackToBxByBz(param, radius+0.5, mass, step*0.1, kTRUE,0.99); + Double_t sxyz[3]={0,0,0}; + AliESDVertex vertex(xyz,sxyz); + Bool_t isOK = param->PropagateToDCA(&vertex,AliTracker::GetBz(),10); + return isOK; +} + + +void AliMaterialBudget::ProcessMCInfo(){ + // + // + // + // + TParticle * particle= new TParticle; + TClonesArray * trefs = new TClonesArray("AliTrackReference"); + const Double_t kPcut=0.05; + const Double_t kMinDrITS = 2.; // minimal distance between references + const Double_t kMinDrTRD = 8.; // minimal distance between references + const Double_t kMinDrTOF = 10.; // minimal distance between references + // + // + // Process tracks + // + Int_t npart = fMCinfo->GetNumberOfTracks(); + if (npart==0) return; + Double_t vertex[4]={0,0,0,0}; + fMCinfo->GetParticleAndTR(0, particle, trefs); + if (particle){ + vertex[0]=particle->Vx(); + vertex[1]=particle->Vy(); + vertex[2]=particle->Vz(); + vertex[3]=particle->R(); + } + // + // + AliTrackReference dummy,*pdummy= &dummy; + AliExternalTrackParam epdummy,*pepdummy= &epdummy; + Int_t nRefITS =0; + Int_t nRefTPC =0; + Int_t nRefTRD =0; + Int_t nRefTOF =0; + AliTrackReference * refITS0, *refITS1; + AliTrackReference * refTPC0, *refTPC1; + AliTrackReference * refTPCIn0, *refTPCIn1; + AliTrackReference * refTRD0, *refTRD1; + AliTrackReference * refTOF0, *refTOF1; + AliTrackReference *refMinR; + // + for (Int_t ipart=0;ipartGetParticleAndTR(ipart, particle, trefs); + AliMCParticle * pp = fMCinfo->GetTrack(ipart); + if (!pp) continue; + if (particle->P()GetMass(); + + // RESET + nRefITS =0; + nRefTPC =0; + nRefTRD =0; + nRefTOF =0; + refITS0=pdummy; refITS1=pdummy; + refTPC0=pdummy; refTPC1=pdummy; + refTPCIn0=pdummy; refTPCIn1=pdummy; + refTRD0=pdummy; refTRD1=pdummy; + refTOF0=pdummy; refTOF1=pdummy; + refMinR = pdummy; + // + Int_t nref = pp->GetNumberOfTrackReferences(); + if (nref==0) continue; + for (Int_t iref = 0; iref < nref; iref++) { + AliTrackReference *ref = pp->GetTrackReference(iref); + if (ref->DetectorId()==AliTrackReference::kDisappeared) continue; + // if (ref.Px()*particle.Px()+ref.Py()*particle.Py()<0) break; // returning track + // + if (ref->DetectorId()==AliTrackReference::kITS){ + if (TMath::Abs(ref->R()-refITS1->R())>kMinDrITS) { + refITS1 = ref; + nRefITS++; + } + if (refITS0==pdummy) refITS0=ref; + } + if (ref->DetectorId()==AliTrackReference::kTPC){ + nRefTPC++; + refTPC1 = ref; + if (refTPC0==pdummy) refTPC0=ref; + } + if (ref->DetectorId()==AliTrackReference::kTRD){ + if (TMath::Abs(ref->R()-refTRD1->R())>kMinDrTRD) { + refTRD1 = ref; + nRefTRD++; + } + if (refTRD0==pdummy) refTRD0=ref; + } + if (ref->DetectorId()==AliTrackReference::kTOF){ + if (TMath::Abs(ref->X()-refTOF1->X()) + TMath::Abs(ref->Y()-refTOF1->Y()) + TMath::Abs(ref->Z()-refTOF1->Z())>kMinDrTOF) { + refTOF1 = ref; + nRefTOF++; + } + if (refTOF0==pdummy) refTOF0=ref; + } + // + // "find inner track ref" + if (ref->DetectorId()==AliTrackReference::kTPC){ + if (ref->Px()*ref->X()+ref->Py()*ref->Y()<0){ + // track in + if (refTPCIn0 == pdummy) refTPCIn0=ref; + if (refTPCIn0 != pdummy && refTPCIn0->R()>ref->R()) + refTPCIn0=ref; + } + if (ref->Px()*ref->X()+ref->Py()*ref->Y()>0){ + // track in + if (refTPCIn1 == pdummy) refTPCIn1=ref; + if (refTPCIn1 != pdummy && refTPCIn1->R()>ref->R()) + refTPCIn1=ref; + } + } + + + if (refMinR==pdummy && ref->P()>0 ){ + refMinR=ref; + } + if (refMinR->R()>ref->R() && ref->P()>0 ) refMinR=ref; + + } + // + AliExternalTrackParam * trackMC = pepdummy; + //track0->GetDZ(0,0,0,bz,dvertex0) + Float_t dist[2]={0,0}; + AliMagF* field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField(); + Double_t esdfield= fESD->GetMagneticField(); + Double_t xyz[3]={0,0,0}; + Double_t bxyz[3]={0,0,0}; + field->Field(xyz,bxyz); + if (refMinR->P()>0) { + trackMC = MakeTrack(refMinR,particle); + trackMC->GetDZ(0,0,0,bxyz[2],dist); + } + Double_t alphaTOF0 = TMath::ATan2(refTOF0->Y(),refTOF0->X()); + Double_t alphaTOF1 = TMath::ATan2(refTOF1->Y(),refTOF1->X()); + Int_t dsecTOF = TMath::Nint(180*(alphaTOF0-alphaTOF1)/(TMath::Pi()*20.)-0.5); + // + // make the two different TPC tracks and propagate them to their DCA + // + Double_t dP = 0; + Bool_t statusProp = kFALSE; + Double_t dY = 0; + Double_t dZ = 0; + AliExternalTrackParam * track0 = pepdummy; + AliExternalTrackParam * track1 = pepdummy; + AliExternalTrackParam * otrack0 = pepdummy; + AliExternalTrackParam * otrack1 = pepdummy; + if (refTPCIn0!=pdummy && refTPCIn1!=pdummy) { + track0 = MakeTrack(refTPCIn0,particle); + track1 = MakeTrack(refTPCIn1,particle); + otrack0 = MakeTrack(refTPCIn0,particle); + otrack1 = MakeTrack(refTPCIn1,particle); + dP = track0->P() - track1->P(); // momentum loss + statusProp = AliMaterialBudget::PropagateCosmicToDCA(track0,track1,mass); + if (statusProp) { + dY = track0->GetY() - track1->GetY(); + dZ = track0->GetZ() - track1->GetZ(); + } + } + // + TTreeSRedirector *pcstream = GetDebugStreamer(); + if (pcstream){ + char name[100]; + for (Int_t id=0; id<3;id++){ + + // if (id==0) sprintf(name,"mcAll"); // all tracks: inconvenient to cut if on is only interest in tracks which reach the TPC + if (id==0) continue; // require TPC + if (id==1) sprintf(name,"mcITS"); + if (id==2) sprintf(name,"mcTPC"); + if (id==1&& nRefITS==0) continue; + if (id==2&& nRefTPC==0) continue; + + (*pcstream)<X(),refIn->Y(), refIn->Z()}; + Double_t mass = part->GetMass(); + Double_t step=1; + // + param=MakeTrack(refOut,part); + paramMC=MakeTrack(refOut,part); + if (!param) return; + if (type<3) PropagateToPoint(param,xyzIn, mass, step); + if (type==3) { + AliTPCseed seed; + seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance()); + Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X()); + seed.Rotate(alpha-seed.GetAlpha()); + seed.SetMass(mass); + for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){ + seed.PropagateTo(xlayer); + } + seed.PropagateTo(refIn->R()); + param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance()); + } + TTreeSRedirector *pcstream = GetDebugStreamer(); + TVectorD gpos(3); + TVectorD gmom(3); + param->GetXYZ(gpos.GetMatrixArray()); + param->GetPxPyPz(gmom.GetMatrixArray()); + if (pcstream){ + (*pcstream)<<"MC"<< + "type="<GetEntries(); + if (nrefsAt(iref); + if (!ref) continue; + Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py(); + if (dir<0) break; + if (ref->DetectorId()!=AliTrackReference::kTPC) continue; + if (iref0<0) iref0 = iref; + iref1 = iref; + } + if (iref1-iref0At(iref0); + AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1); + AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part); + AliExternalTrackParam *paramUpdate = MakeTrack(refIn,part); + paramUpdate->AddCovariance(covar); + Double_t mass = part->GetMass(); + Double_t charge = part->GetPDG()->Charge()/3.; +/* + Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X()); + Float_t radiusIn= refIn->R(); + Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X()); + Float_t radiusOut= refOut->R(); +*/ + Bool_t isOKP=kTRUE; + Bool_t isOKU=kTRUE; + AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField(); + for (Int_t iref = iref0; iref<=iref1; iref++){ + AliTrackReference * ref = (AliTrackReference*)trefs->At(iref); + Float_t alphaC= TMath::ATan2(ref->Y(),ref->X()); + Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()}; + Double_t mag[3]; + field->Field(pos,mag); + isOKP&=paramPropagate->Rotate(alphaC); + isOKU&=paramUpdate->Rotate(alphaC); + for (Float_t xref= paramPropagate->GetX(); xrefR(); xref++){ + isOKP&=paramPropagate->PropagateTo(xref, mag[2]); + isOKU&=paramUpdate->PropagateTo(xref, mag[2]); + } + isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]); + isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]); + Double_t clpos[2] = {0, ref->Z()}; + Double_t clcov[3] = { 0.005,0,0.005}; + isOKU&= paramUpdate->Update(clpos, clcov); + } + TTreeSRedirector *pcstream = GetDebugStreamer(); + if (pcstream){ + TVectorD gposU(3); + TVectorD gmomU(3); + TVectorD gposP(3); + TVectorD gmomP(3); + paramUpdate->GetXYZ(gposU.GetMatrixArray()); + paramUpdate->GetPxPyPz(gmomU.GetMatrixArray()); + paramPropagate->GetXYZ(gposP.GetMatrixArray()); + paramPropagate->GetPxPyPz(gmomP.GetMatrixArray()); + + (*pcstream)<<"MCupdate"<< + "isOKU="<3) + cout << "AliMaterialBudget::FindPairs()" << endl; + + + AliESDfriend *ESDfriend=static_cast(event->FindListObject("AliESDfriend")); + Int_t ntracks=event->GetNumberOfTracks(); + TObjArray tpcSeeds(ntracks); + if (ntracks==0) return; + Double_t vtxx[3]={0,0,0}; + Double_t svtxx[3]={0.000001,0.000001,100.}; + AliESDVertex vtx(vtxx,svtxx); + // + //track loop + // + for (Int_t i=0;iGetTrack(i); + const AliExternalTrackParam * trackIn = track->GetInnerParam(); + const AliExternalTrackParam * trackOut = track->GetOuterParam(); + if (!trackIn) continue; + if (!trackOut) continue; + AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i); + TObject *calibObject; + AliTPCseed *seed = 0; + for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { + if ((seed=dynamic_cast(calibObject))) break; + } + if (seed) tpcSeeds.AddAt(seed,i); + } + + if (ntracks<2) return; + // + // Find pairs + // + for (Int_t i=0;iGetTrack(i); + // track0 - choosen upper part + if (!track0) continue; + if (!track0->GetOuterParam()) continue; + if (track0->GetOuterParam()->GetAlpha()<0) continue; + Double_t dir0[3]; + track0->GetDirection(dir0); + for (Int_t j=0;jGetTrack(j); + //track 1 lower part + if (!track1) continue; + if (!track1->GetOuterParam()) continue; + if (track1->GetOuterParam()->GetAlpha()>0) continue; + // + Double_t dir1[3]; + track1->GetDirection(dir1); + + AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i); + AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j); + if (! seed0) continue; + if (! seed1) continue; + // + Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]); + Float_t d0 = track0->GetLinearD(0,0); + Float_t d1 = track1->GetLinearD(0,0); + // + // conservative cuts - convergence to be guarantied + // applying before track propagation + if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0 + if (dir>fCutMinDir) continue; // direction vector product + Float_t bz = AliTracker::GetBz(); + Float_t dvertex0[2]; //distance to 0,0 + Float_t dvertex1[2]; //distance to 0,0 + track0->GetDZ(0,0,0,bz,dvertex0); + track1->GetDZ(0,0,0,bz,dvertex1); + if (TMath::Abs(dvertex0[1])>250) continue; + if (TMath::Abs(dvertex1[1])>250) continue; + // + // + // + Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1)); + AliExternalTrackParam param0(*track0); + AliExternalTrackParam param1(*track1); + // + // Propagate using Magnetic field and correct fo material budget + // + AliTracker::PropagateTrackToBxByBz(¶m0,dmax+1,0.0005,3,kTRUE); + AliTracker::PropagateTrackToBxByBz(¶m1,dmax+1,0.0005,3,kTRUE); + // + // Propagate rest to the 0,0 DCA - z should be ignored + // + Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000); + Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000); + // + param0.GetDZ(0,0,0,bz,dvertex0); + param1.GetDZ(0,0,0,bz,dvertex1); + if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue; + // + Double_t xyz0[3];//,pxyz0[3]; + Double_t xyz1[3];//,pxyz1[3]; + param0.GetXYZ(xyz0); + param1.GetXYZ(xyz1); + Bool_t isPair = IsPair(¶m0,¶m1); + // + // HERE WE WILL PUT THE ACCESS TO THE MC TRACKS AND MATCH THESE !!!! + // + Int_t label0 = TMath::Abs(track0->GetLabel()); + AliMCParticle *mcParticle0 = fMCinfo->GetTrack(label0); + TParticle *particle0 = mcParticle0->Particle(); + AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle0); // get the first TPC track reference + if (!ref0) continue; + AliExternalTrackParam *paramMC0 = 0; + paramMC0 = MakeTrack(ref0, particle0); + // + Int_t label1 = TMath::Abs(track1->GetLabel()); + AliMCParticle *mcParticle1 = fMCinfo->GetTrack(label1); + TParticle *particle1 = mcParticle1->Particle(); + AliTrackReference *ref1 = GetFirstTPCTrackRef(mcParticle1); // get the first TPC track reference + if (!ref1) continue; + AliExternalTrackParam *paramMC1 = 0; + paramMC1 = MakeTrack(ref1, particle1); + // + // ACCESS TOF INFORMATION + Int_t nTrackRefTOF0 = 0; + Int_t nTrackRefITS0 = 0; + AliTrackReference * refLastTOF0 = 0; + AliTrackReference * refFirstTOF0 = GetAllTOFinfo(mcParticle0, nTrackRefTOF0, nTrackRefITS0); + Float_t alphaTOF0 = 0; + if (refFirstTOF0) alphaTOF0 = refFirstTOF0->Alpha(); + // + Int_t nTrackRefTOF1 = 0; + Int_t nTrackRefITS1 = 0; + AliTrackReference * refLastTOF1 = 0; + AliTrackReference * refFirstTOF1 =GetAllTOFinfo(mcParticle1, nTrackRefTOF1, nTrackRefITS1); + Float_t alphaTOF1 = 0; + if (refFirstTOF1) alphaTOF1 = refFirstTOF1->Alpha(); + //cout << " STATUS "<0){ + TTreeSRedirector * cstream = GetDebugStreamer(); + AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam(); + AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam(); + AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam(); + AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam(); + // + // + // + if (cstream) { + (*cstream) << "Track0" << + "dir="<GetParameter(); + const Double_t *p1 = tr1->GetParameter(); + if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE; + if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE; + if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE; + + Double_t d0[3], d1[3]; + tr0->GetDirection(d0); + tr1->GetDirection(d1); + if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE; + // + return kTRUE; +} + + +AliTrackReference * AliMaterialBudget::GetFirstTPCTrackRef(AliMCParticle *mcParticle) +{ + // return first TPC track reference + if(!mcParticle) return 0; + + // find first track reference + // check direction to select proper reference point for looping tracks + Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences(); + AliTrackReference *ref = 0; + AliTrackReference *refIn = 0; + for (Int_t iref = 0; iref < nTrackRef; iref++) { + ref = mcParticle->GetTrackReference(iref); + if(ref && (ref->DetectorId()==AliTrackReference::kTPC)) + { + //Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py(); + //if(dir < 0.) break; + + refIn = ref; + break; + } + } + + return refIn; +} + + +AliTrackReference * AliMaterialBudget::GetAllTOFinfo(AliMCParticle *mcParticle, Int_t &nTrackRef, Int_t &nTrackRefITS, Int_t retValue) { + // + // + // + + if(!mcParticle) return 0; + Int_t counter = 0; + nTrackRef = 0; + nTrackRefITS = 0; + AliTrackReference *ref = 0; + for (Int_t iref = 0; iref < mcParticle->GetNumberOfTrackReferences(); iref++) { + ref = mcParticle->GetTrackReference(iref); + if(ref && (ref->DetectorId()==AliTrackReference::kTOF)) { + counter = iref; + nTrackRef++; + } + if(ref && (ref->DetectorId()==AliTrackReference::kITS)) nTrackRefITS++; + } + if (nTrackRef ==0) return 0; + if (retValue == 0) return mcParticle->GetTrackReference(counter - nTrackRef +1); + return mcParticle->GetTrackReference(counter); + +} + + +void AliMaterialBudget::FinishTaskOutput() +{ + // + // According description in AliAnalisysTask this method is call + // on the slaves before sending data + // + Terminate("slave"); + gSystem->Exec("pwd"); + RegisterDebugOutput(); + +} + + +void AliMaterialBudget::RegisterDebugOutput(){ + // + // + // + // + // store - copy debug output to the destination position + // currently ONLY for local copy + TString dsName; + dsName=GetName(); + dsName+="Debug.root"; + dsName.ReplaceAll(" ",""); + TString dsName2=fDebugOutputPath.Data(); + gSystem->MakeDirectory(dsName2.Data()); + dsName2+=gSystem->HostName(); + gSystem->MakeDirectory(dsName2.Data()); + dsName2+="/"; + dsName2+=gSystem->BaseName(gSystem->pwd()); + dsName2+="/"; + gSystem->MakeDirectory(dsName2.Data()); + dsName2+=dsName; + AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data())); + printf("copy %s\t%s\n",dsName.Data(),dsName2.Data()); + TFile::Cp(dsName.Data(),dsName2.Data()); +} + + +Bool_t AliMaterialBudget::PropagateCosmicToDCA(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass){ + // + // param0 - upper part of cosmic track + // param1 - lower part of cosmic track + // + // 0. Propagate both tracks to DCA to (0,0,0) + // 1. After propagation to DCA rotate track param1 to corrdinate system of track1 <-> rotate param0 to coordinate system of param 1 ???? + // 2. Propagate track 1 to refernce x from track0 + // + + // step 0. + + Float_t d0 = param0->GetLinearD(0,0); + Float_t d1 = param1->GetLinearD(0,0); + Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1)); + // + // propagate in the beginning taking all material into account + // + AliTracker::PropagateTrackToBxByBz(param0,dmax+1.,mass,0.5,kTRUE,0.99,-1.); + AliTracker::PropagateTrackToBxByBz(param1,dmax+1.,mass,0.5,kTRUE,0.99,1.); + // + Double_t vtxx[3]={0,0,0}; + Double_t svtxx[3]={0.000001,0.000001,100.}; + AliESDVertex vtx(vtxx,svtxx); + // + Bool_t b0 = param0->PropagateToDCA(&vtx,AliTracker::GetBz(),1000); + Bool_t b1 = param1->PropagateToDCA(&vtx,AliTracker::GetBz(),1000); + + if (!(b0 && b1)) return 0; + + // step 1. + + Float_t dAlpha = param0->GetAlpha(); + param1->Rotate(dAlpha); + + // step 2. + + Float_t refX = param0->GetX(); + param1->PropagateTo(refX,AliTracker::GetBz()); + + return kTRUE; + + +} + + diff --git a/PWG1/AliMaterialBudget.h b/PWG1/AliMaterialBudget.h new file mode 100644 index 00000000000..0fb1272e4c8 --- /dev/null +++ b/PWG1/AliMaterialBudget.h @@ -0,0 +1,93 @@ +#ifndef ALIMATERIALBUDGET_H +#define ALIMATERIALBUDGET_H + +// ROOT includes +#include +#include +#include + +// AliRoot includes +#include +#include +#include +#include +#include +#include +#include +class AliGenInfoMaker; +class TTreeSRedirector; +class AliMCEventHadnler; +class TParticle; +class AliMCInfo; +class AliESDRecInfo; +class AliESDEvent; +class AliMCEvent; +class AliComparisonObject; + +class AliMaterialBudget : public AliAnalysisTask { + public: + AliMaterialBudget(); + AliMaterialBudget(const char *name); + virtual ~AliMaterialBudget(); + + virtual void ConnectInputData(Option_t *); + virtual void CreateOutputObjects(); + virtual void Exec(Option_t *option); + virtual void Terminate(Option_t *); + virtual void FinishTaskOutput(); + void SetDebugOuputhPath(const char * name){fDebugOutputPath=name;} + + // + void FindPairs(AliESDEvent * event); + Bool_t IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1); + // + void ProcessMCInfo(); + void ProcessRefTracker(AliTrackReference* refIn, AliTrackReference* refOut, TParticle*part, Int_t type); + + void FitTrackRefs(TParticle * part, TClonesArray * trefs); + + // + // debug streamer part + // + TTreeSRedirector *GetDebugStreamer(); + void SetStreamLevel(Int_t streamLevel){fStreamLevel=streamLevel;} + void SetDebugLevel(Int_t level) {fDebugLevel = level;} + Int_t GetStreamLevel() const {return fStreamLevel;} + Int_t GetDebugLevel() const {return fDebugLevel;} + // + static Bool_t PropagateCosmicToDCA(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass); + static AliExternalTrackParam * MakeTrack(const AliTrackReference* ref, TParticle*part); + static Bool_t PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step); + // + AliTrackReference * GetFirstTPCTrackRef(AliMCParticle *mcParticle); + AliTrackReference * GetAllTOFinfo(AliMCParticle *mcParticle, Int_t & nTrackRef, Int_t &nTrackRefITS, Int_t retValue =0); + protected: + void RegisterDebugOutput(); + AliMaterialBudget(const AliMaterialBudget& /*info*/); + AliMaterialBudget& operator=(const AliMaterialBudget& /*info*/) { return *this;} + AliMCEvent * fMCinfo; //! MC event handler + AliESDEvent * fESD; //! current esd event + // + // + // + TTreeSRedirector *fDebugStreamer; //! debug streamer + Int_t fStreamLevel; // debug stream level + Int_t fDebugLevel; // debug level + TString fDebugOutputPath; // debug output path + // + // histogran + // + TList * fListHist; // list for histograms + TH1F * fHistMult; // track multiplicity histograms + // + // cuts + // + Float_t fCutMaxD; // maximal distance in rfi ditection + Float_t fCutMaxDz; // maximal distance in z ditection + Float_t fCutTheta; // maximal distance in theta ditection + Float_t fCutMinDir; // direction vector products + // + ClassDef(AliMaterialBudget, 1); // Analysis task base class for tracks +}; + +#endif diff --git a/PWG1/PWG1LinkDef.h b/PWG1/PWG1LinkDef.h index 8babf317024..4ec06fb18ad 100644 --- a/PWG1/PWG1LinkDef.h +++ b/PWG1/PWG1LinkDef.h @@ -44,6 +44,8 @@ #pragma link C++ class AliAlignmentDataFilterITS+; #pragma link C++ class AliAnalysisTaskV0QA+; +#pragma link C++ class AliMaterialBudget+; + #endif diff --git a/PWG1/libPWG1.pkg b/PWG1/libPWG1.pkg index 167c7b27d9c..f88b0ea5804 100644 --- a/PWG1/libPWG1.pkg +++ b/PWG1/libPWG1.pkg @@ -31,7 +31,8 @@ SRCS:= AliTreeDraw.cxx \ AliPerformanceTPC.cxx \ AliPerformanceMC.cxx \ AliAlignmentDataFilterITS.cxx \ - AliAnalysisTaskV0QA.cxx + AliAnalysisTaskV0QA.cxx \ + AliMaterialBudget.cxx HDRS:= $(SRCS:.cxx=.h) -- 2.43.0