From: mivanov Date: Thu, 6 Jun 2013 05:53:19 +0000 (+0000) Subject: for Jens and Martin X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;ds=sidebyside;h=526ddf0e8f51e0f015a38f96dca585ef1548e831;p=u%2Fmrichter%2FAliRoot.git for Jens and Martin Continuous readout toy simulation: loadlibs.C //modified version of the one used for GEMtest code. I haven't looked into this so it probably contains a lot of stuff that is not needed ToyMCTrack, ToyMCEvent //classes for storing event and track parameters ToyMCEventGenerator //base class for event generators, contains DistortTrack method to distort track using method borrowed from AliTPCCorrection ToyMCEventGeneratorSimple //simple first MC event generator, borrowing from Marians macro. Inherits from ToyMCEventGenerator to get DistortTrack method. Generate(Double_t time) generates and returns an event at time time makeTree.C //short macro to generate events from a poisson distr and storing in a ROOT tree. to run : root -l loadlibs.C .L makeTree.C+ makeTree(Double_t collFreq/*kHz*/, Double_t bunchFreq/*MHz*/, Int_t nEvents) --- diff --git a/TPC/Upgrade/README b/TPC/Upgrade/README new file mode 100644 index 00000000000..f04ab93609f --- /dev/null +++ b/TPC/Upgrade/README @@ -0,0 +1,16 @@ +Continuous readout toy simulation + +loadlibs.C //modified version of the one used for GEMtest code. I haven't looked into this so it probably contains a lot of stuff that is not needed + +ToyMCTrack, ToyMCEvent //classes for storing event and track parameters + +ToyMCEventGenerator //base class for event generators, contains DistortTrack method to distort track using method borrowed from AliTPCCorrection + +ToyMCEventGeneratorSimple //simple first MC event generator, borrowing from Marians macro. Inherits from ToyMCEventGenerator to get DistortTrack method. Generate(Double_t time) generates and returns an event at time time + +makeTree.C //short macro to generate events from a poisson distr and storing in a ROOT tree. + + +to run : root -l loadlibs.C + .L makeTree.C+ + makeTree(Double_t collFreq/*kHz*/, Double_t bunchFreq/*MHz*/, Int_t nEvents) \ No newline at end of file diff --git a/TPC/Upgrade/ToyMCEvent.cxx b/TPC/Upgrade/ToyMCEvent.cxx new file mode 100644 index 00000000000..df808c194d9 --- /dev/null +++ b/TPC/Upgrade/ToyMCEvent.cxx @@ -0,0 +1,45 @@ +#include "ToyMCEvent.h" + +ClassImp(ToyMCEvent); +Int_t ToyMCEvent::evCounter = 0; + +ToyMCEvent::ToyMCEvent() + :TObject() + ,fT0(-1.) + ,fX(-1000.) + ,fY(-1000.) + ,fZ(-1000.) + ,fTracks("ToyMCTrack") +{ + fEventNumber = evCounter; + evCounter++; +} + +//____________________________________________________ +ToyMCEvent::ToyMCEvent(const ToyMCEvent &event) + : TObject(event) + ,fEventNumber(event.fEventNumber) + ,fT0(event.fT0) + ,fX(event.fX) + ,fY(event.fY) + ,fZ(event.fZ) + ,fTracks(event.fTracks) +{ + // +} + +//_____________________________________________________ +ToyMCEvent& ToyMCEvent::operator = (const ToyMCEvent &event) +{ + //assignment operator + if (&event == this) return *this; + new (this) ToyMCEvent(event); + + return *this; +} +//_____________________________________________________ +ToyMCTrack* ToyMCEvent::AddTrack(const ToyMCTrack &track) +{ + return new(fTracks[fTracks.GetEntriesFast()]) ToyMCTrack(track); +} + diff --git a/TPC/Upgrade/ToyMCEvent.h b/TPC/Upgrade/ToyMCEvent.h new file mode 100644 index 00000000000..53e0a9de370 --- /dev/null +++ b/TPC/Upgrade/ToyMCEvent.h @@ -0,0 +1,48 @@ +#ifndef ToyMCEvent_H +#define ToyMCEvent_H + + +#include +#include "ToyMCTrack.h" + +class ToyMCEvent : public TObject { + public: + ToyMCEvent(); + ToyMCEvent(const ToyMCEvent &event); + virtual ~ToyMCEvent() {} + ToyMCEvent& operator = (const ToyMCEvent &event); + + + ToyMCTrack* AddTrack(const ToyMCTrack &track); + + Int_t GetNumberOfTracks() const { return fTracks.GetEntriesFast(); } + const ToyMCTrack* GetTrack(Int_t track) const { return static_cast(fTracks.At(track)); } + + void SetT0 (Float_t time) { fT0 = time; } + void SetX(Float_t var) { fX = var; } + void SetY(Float_t var) { fY = var; } + void SetZ(Float_t var) { fZ = var; } + + UInt_t GetEventNumber() const {return fEventNumber; } + Float_t GetT0() const {return fT0; } + Float_t GetX() const {return fX; } + Float_t GetY() const {return fY; } + Float_t GetZ() const {return fZ; } + + private: + static Int_t evCounter; + + UInt_t fEventNumber; + + Float_t fT0; + Float_t fX; + Float_t fY; + Float_t fZ; + + TClonesArray fTracks; + + ClassDef(ToyMCEvent, 1); + +}; + +#endif diff --git a/TPC/Upgrade/ToyMCEventGenerator.cxx b/TPC/Upgrade/ToyMCEventGenerator.cxx new file mode 100644 index 00000000000..309db5db62b --- /dev/null +++ b/TPC/Upgrade/ToyMCEventGenerator.cxx @@ -0,0 +1,191 @@ +#include +#include "ToyMCEventGenerator.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +ClassImp(ToyMCEventGenerator); + + +ToyMCEventGenerator::ToyMCEventGenerator() + :TObject() + ,fTPCParam(0x0) +{ + //TODO: Add method to set custom space charge file + fSpaceCharge = new AliTPCSpaceCharge3D(); + fSpaceCharge->SetSCDataFileName("SC_NeCO2_eps5_50kHz.root"); + fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2 + //fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2 + fSpaceCharge->InitSpaceCharge3DDistortion(); + fSpaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1"); + fSpaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz"); + const char* ocdb="local://$ALICE_ROOT/OCDB/"; + AliCDBManager::Instance()->SetDefaultStorage(ocdb); + AliCDBManager::Instance()->SetRun(0); + TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG)); + AliGeomManager::LoadGeometry(""); + fTPCParam = AliTPCcalibDB::Instance()->GetParameters(); + fTPCParam->ReadGeoMatrices(); +} +//________________________________________________________________ +ToyMCEventGenerator::ToyMCEventGenerator(const ToyMCEventGenerator &gen) + :TObject(gen) + ,fSpaceCharge(gen.fSpaceCharge) +{ + // +} +//________________________________________________________________ +ToyMCEventGenerator::~ToyMCEventGenerator() +{ + delete fSpaceCharge; + delete fTPCParam; +} +//________________________________________________________________ +Bool_t ToyMCEventGenerator::DistortTrack(ToyMCTrack &trackIn, Double_t t0) { + + if(!fTPCParam) { + fTPCParam = AliTPCcalibDB::Instance()->GetParameters(); + fTPCParam->ReadGeoMatrices(); + std::cout << "init tpc param" << std::endl; + } + const Int_t kNRows=AliTPCROC::Instance()->GetNRows(0)+AliTPCROC::Instance()->GetNRows(36); + // const Double_t kRTPC0 =AliTPCROC::Instance()->GetPadRowRadii(0,0); + // const Double_t kRTPC1 =AliTPCROC::Instance()->GetPadRowRadii(36,AliTPCROC::Instance()->GetNRows(36)-1); + const Double_t kMaxSnp = 0.85; + const Double_t kSigmaY=0.1; + const Double_t kSigmaZ=0.1; + const Double_t kDriftVel = fTPCParam->GetDriftV()/1000000; + const Double_t kMaxZ0=fTPCParam->GetZLength(); + const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); + + const Double_t iFCRadius= 83.5; //radius constants found in AliTPCCorrection.cxx + const Double_t oFCRadius= 254.5; + + ToyMCTrack track(trackIn); + const Int_t nMaxPoints = kNRows+50; + AliTrackPointArray pointArray0(nMaxPoints); + AliTrackPointArray pointArray1(nMaxPoints); + Double_t xyz[3]; + + if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp)) return 0; + + Int_t npoints=0; + Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame + + + //Simulate track from inner field cage radius to outer and save points along trajectory + Int_t nMinPoints = 40; + for (Double_t radius=iFCRadius; radiusGaus(0,0.000005); + xyz[1]+=gRandom->Gaus(0,0.000005); + xyz[2]+=gRandom->Gaus(0,0.000005); + + if (TMath::Abs(track.GetZ())>kMaxZ0) continue; + if (TMath::Abs(track.GetX())oFCRadius) continue; + + AliTrackPoint pIn0; // space point + AliTrackPoint pIn1; + Int_t sector= (xyz[2]>0)? 0:18; + pointArray0.GetPoint(pIn0,npoints); + pointArray1.GetPoint(pIn1,npoints); + Double_t alpha = TMath::ATan2(xyz[1],xyz[0]); + Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]}; + fSpaceCharge->DistortPoint(distPoint, sector); + pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]); + pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]); + track.Rotate(alpha); + AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point + AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point + prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint); + prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint); + pIn0=prot0.Rotate(-alpha); // rotate back to global frame + pIn1=prot1.Rotate(-alpha); // rotate back to global frame + pointArray0.AddPoint(npoints, &pIn0); + pointArray1.AddPoint(npoints, &pIn1); + npoints++; + if (npoints>=nMaxPoints) break; + } + + if (npointsSetX(xArr[iPoint]); + tempCl->SetY(yArr[iPoint]); + tempCl->SetZ(zArr[iPoint]); + + Float_t xyz[3] = {xArrDist[iPoint],yArrDist[iPoint],zArrDist[iPoint]}; + Int_t index[3] = {0}; + + Float_t padRow = fTPCParam->GetPadRow(xyz,index); + Int_t sector = index[1]; + + if(padRow < 0 || (sector < 36 && padRow>62) || (sector > 35 &&padRow > 95) ) continue; //outside sensitive area + + + if(lastRow!=padRow) { + + + //make and save distorted cluster + if(pntsInCurrentRow>0){ + xt/=pntsInCurrentRow; + yt/=pntsInCurrentRow; + zt/=pntsInCurrentRow; + AliTPCclusterMI *tempDistCl = trackIn.AddDistortedSpacePoint(AliTPCclusterMI()); + tempDistCl->SetX(xt); + tempDistCl->SetY(yt); + tempDistCl->SetZ(zt); + Float_t clxyz[3] = {xt,yt,zt}; + Int_t clindex[3] = {0}; + Int_t clRow = fTPCParam->GetPadRow(clxyz,clindex); + Int_t nPads = fTPCParam->GetNPads(clindex[1], clRow); + Int_t pad = TMath::Nint(clxyz[1] + nPads/2); //taken from AliTPC.cxx + tempDistCl->SetPad(pad); + tempDistCl->SetRow(clRow); + tempDistCl->SetDetector(clindex[1]); + tempDistCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(zt))/kDriftVel); // set time as t0 + drift time from dist z coordinate to pad plane + } + xt = 0; + yt = 0; + zt = 0; + pntsInCurrentRow = 0; + } + + xt+=xArrDist[iPoint]; + yt+=yArrDist[iPoint]; + zt+=zArrDist[iPoint]; + pntsInCurrentRow++; + lastRow = padRow; + } + + return 1; + + + +} diff --git a/TPC/Upgrade/ToyMCEventGenerator.h b/TPC/Upgrade/ToyMCEventGenerator.h new file mode 100644 index 00000000000..af5112ff9a5 --- /dev/null +++ b/TPC/Upgrade/ToyMCEventGenerator.h @@ -0,0 +1,30 @@ +#ifndef ToyMCEventGenerator_H +#define ToyMCEventGenerator_H + + +#include "ToyMCEvent.h" +#include "ToyMCTrack.h" +#include +#include +class ToyMCEventGenerator : public TObject { + public: + ToyMCEventGenerator(); + ToyMCEventGenerator(const ToyMCEventGenerator &gen); + virtual ~ToyMCEventGenerator(); + + virtual ToyMCEvent* Generate(Double_t time) = 0; + + Bool_t DistortTrack(ToyMCTrack &trackIn, Double_t t0); + + protected: + AliTPCParam *fTPCParam; + + private: + AliTPCSpaceCharge3D *fSpaceCharge; + + ClassDef(ToyMCEventGenerator, 1) + +}; + +#endif + diff --git a/TPC/Upgrade/ToyMCEventGeneratorSimple.cxx b/TPC/Upgrade/ToyMCEventGeneratorSimple.cxx new file mode 100644 index 00000000000..de073ae0187 --- /dev/null +++ b/TPC/Upgrade/ToyMCEventGeneratorSimple.cxx @@ -0,0 +1,88 @@ +#include +#include "ToyMCEventGeneratorSimple.h" +#include +#include +#include +#include "ToyMCEvent.h" + +ClassImp(ToyMCEventGeneratorSimple); + + +ToyMCEventGeneratorSimple::ToyMCEventGeneratorSimple() + :ToyMCEventGenerator() + ,fVertexMean(0.) + ,fVertexSigma(7.) +{ + // + +} +//________________________________________________________________ +ToyMCEventGeneratorSimple::ToyMCEventGeneratorSimple(const ToyMCEventGeneratorSimple &gen) + :ToyMCEventGenerator(gen) + ,fVertexMean(gen.fVertexMean) + ,fVertexSigma(gen.fVertexSigma) +{ +} +//________________________________________________________________ +ToyMCEventGeneratorSimple::~ToyMCEventGeneratorSimple() +{ +} + //________________________________________________________________ +ToyMCEventGeneratorSimple& ToyMCEventGeneratorSimple::operator = (const ToyMCEventGeneratorSimple &gen) +{ + //assignment operator + if (&gen == this) return *this; + new (this) ToyMCEventGeneratorSimple(gen); + + return *this; +} +//________________________________________________________________ +void ToyMCEventGeneratorSimple::SetParameters(Double_t vertexMean, Double_t vertexSigma) { + fVertexMean = vertexMean; + fVertexSigma = vertexSigma; +} +//________________________________________________________________ +ToyMCEvent* ToyMCEventGeneratorSimple::Generate(Double_t time) { + + ToyMCEvent *retEvent = new ToyMCEvent(); + retEvent->SetT0(time); + retEvent->SetX(0); + retEvent->SetX(0); + retEvent->SetZ(gRandom->Gaus(fVertexMean,fVertexSigma)); + + Double_t etaCuts=.9; + Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); + TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10); + fpt.SetParameters(7.24,0.120); + fpt.SetNpx(10000); + Int_t nTracks = 400; //TODO: draw from experim dist + + for(Int_t iTrack=0; iTrackUniform(0.0, 2*TMath::Pi()); + Double_t eta = gRandom->Uniform(-etaCuts, etaCuts); + Double_t pt = fpt.GetRandom(); // momentum for f1 + Short_t sign=1; + if(gRandom->Rndm() < 0.5){ + sign =1; + }else{ + sign=-1; + } + + Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.; + Double_t pxyz[3]; + pxyz[0]=pt*TMath::Cos(phi); + pxyz[1]=pt*TMath::Sin(phi); + pxyz[2]=pt*TMath::Tan(theta); + Double_t vertex[3]={0,0,retEvent->GetZ()}; + Double_t cv[21]={0}; + ToyMCTrack *tempTrack = new ToyMCTrack(vertex,pxyz,cv,sign); + + Bool_t trackDist = DistortTrack(*tempTrack, time); + if(trackDist) retEvent->AddTrack(*tempTrack); + delete tempTrack; + } + + return retEvent; +} + + diff --git a/TPC/Upgrade/ToyMCEventGeneratorSimple.h b/TPC/Upgrade/ToyMCEventGeneratorSimple.h new file mode 100644 index 00000000000..b924950f381 --- /dev/null +++ b/TPC/Upgrade/ToyMCEventGeneratorSimple.h @@ -0,0 +1,32 @@ +#ifndef ToyMCEventGeneratorSimple_H +#define ToyMCEventGeneratorSimple_H + + +#include "ToyMCEvent.h" +#include "ToyMCTrack.h" +#include "ToyMCEventGenerator.h" +class ToyMCEventGeneratorSimple : public ToyMCEventGenerator { + public: + ToyMCEventGeneratorSimple(); + ToyMCEventGeneratorSimple(const ToyMCEventGeneratorSimple &gen); + virtual ~ToyMCEventGeneratorSimple(); + ToyMCEventGeneratorSimple & operator = (const ToyMCEventGeneratorSimple &gen); + + ToyMCEvent* Generate(Double_t time); + void SetParameters(Double_t vertexMean, Double_t vertexSigma); + private: + + Double_t fVertexMean; + Double_t fVertexSigma; + + ClassDef(ToyMCEventGeneratorSimple, 1) + +}; + + + + + + +#endif + diff --git a/TPC/Upgrade/ToyMCTrack.cxx b/TPC/Upgrade/ToyMCTrack.cxx new file mode 100644 index 00000000000..3be217dbf17 --- /dev/null +++ b/TPC/Upgrade/ToyMCTrack.cxx @@ -0,0 +1,56 @@ +#include "ToyMCTrack.h" + +ClassImp(ToyMCTrack); + +ToyMCTrack::ToyMCTrack() + :AliExternalTrackParam() + ,fSpacePoints("AliTPCclusterMI") + ,fDistortedSpacePoints("AliTPCclusterMI") +{ + //default constructor +} +//____________________________________________________ +ToyMCTrack::ToyMCTrack(const ToyMCTrack &track) + : AliExternalTrackParam(track) + ,fSpacePoints(track.fSpacePoints) + ,fDistortedSpacePoints(track.fDistortedSpacePoints) +{ + //copy constructor +} +//_____________________________________________________ +ToyMCTrack& ToyMCTrack::operator = (const ToyMCTrack &track) +{ + //assignment operator + if (&track == this) return *this; + new (this) ToyMCTrack(track); + + return *this; +} +//________________________________________________________________ +ToyMCTrack::ToyMCTrack(Double_t x, Double_t alpha, + const Double_t param[5], + const Double_t covar[15]) + :AliExternalTrackParam(x,alpha,param,covar) + ,fSpacePoints("AliTPCclusterMI") + ,fDistortedSpacePoints("AliTPCclusterMI") +{ + //create external track parameters from given arguments +} +//________________________________________________________________ +ToyMCTrack::ToyMCTrack(Double_t xyz[3],Double_t pxpypz[3], + Double_t cv[21],Short_t sign) + :AliExternalTrackParam(xyz,pxpypz,cv,sign) + ,fSpacePoints("AliTPCclusterMI") + ,fDistortedSpacePoints("AliTPCclusterMI") +{ +} +//________________________________________________________________ +AliTPCclusterMI* ToyMCTrack::AddSpacePoint(const AliTPCclusterMI &spoint) +{ + return new(fSpacePoints[fSpacePoints.GetEntriesFast()]) AliTPCclusterMI(spoint); +} +//________________________________________________________________ +AliTPCclusterMI* ToyMCTrack::AddDistortedSpacePoint(const AliTPCclusterMI &spoint) +{ + return new(fDistortedSpacePoints[fDistortedSpacePoints.GetEntriesFast()]) AliTPCclusterMI(spoint); +} diff --git a/TPC/Upgrade/ToyMCTrack.h b/TPC/Upgrade/ToyMCTrack.h new file mode 100644 index 00000000000..80b5efbc89a --- /dev/null +++ b/TPC/Upgrade/ToyMCTrack.h @@ -0,0 +1,39 @@ +#ifndef TOYMCTRACK_H +#define TOYMCTRACK_H + +#include +#include +#include + +class ToyMCTrack : public AliExternalTrackParam { + + public: + ToyMCTrack(); + ToyMCTrack(Double_t x, Double_t alpha, + const Double_t param[5], + const Double_t covar[15]); + ToyMCTrack(Double_t xyz[3],Double_t pxpypz[3], + Double_t cv[21],Short_t sign); + ToyMCTrack(const ToyMCTrack &track); + ToyMCTrack& operator=(const ToyMCTrack &track); + virtual ~ToyMCTrack() {} + + Int_t GetNumberOfSpacePoints() const { return fSpacePoints.GetEntriesFast(); } + Int_t GetNumberOfDistSpacePoints() const { return fDistortedSpacePoints.GetEntriesFast(); } + + const AliTPCclusterMI* GetSpacePoint(Int_t spoint) const { return static_cast (fSpacePoints.At(spoint)); } + const AliTPCclusterMI* GetDistortedSpacePoint(Int_t dspoint) const { return static_cast (fDistortedSpacePoints.At(dspoint)); } + + AliTPCclusterMI* AddSpacePoint(const AliTPCclusterMI &spoint); + AliTPCclusterMI* AddDistortedSpacePoint(const AliTPCclusterMI &dspoint); + + private: + + TClonesArray fSpacePoints; + TClonesArray fDistortedSpacePoints; + + + ClassDef(ToyMCTrack,1) +}; + +#endif