--- /dev/null
+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
--- /dev/null
+#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);
+}
+
--- /dev/null
+#ifndef ToyMCEvent_H
+#define ToyMCEvent_H
+
+
+#include <TClonesArray.h>
+#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<const ToyMCTrack*>(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
--- /dev/null
+#include <iostream>
+#include "ToyMCEventGenerator.h"
+#include <AliTPCROC.h>
+#include <AliTrackPointArray.h>
+#include <TDatabasePDG.h>
+#include <TRandom.h>
+#include <TH2F.h>
+#include <AliTrackerBase.h>
+#include <AliCDBManager.h>
+#include <AliTPCParam.h>
+#include <AliGeomManager.h>
+#include <AliTPCcalibDB.h>
+#include <TGeoGlobalMagField.h>
+#include <AliTPCclusterMI.h>
+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; radius<oFCRadius; radius++){
+
+ if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) return 0;
+ track.GetXYZ(xyz);
+ xyz[0]+=gRandom->Gaus(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())<iFCRadius) 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 (npoints<nMinPoints) return 0;
+
+ //save space points and make clusters of distorted points
+ Int_t lastRow = 0;
+ Int_t pntsInCurrentRow = 0;
+ Double_t xt=0, yt=0, zt=0;
+ for(Int_t iPoint = 0; iPoint<npoints; iPoint++){
+
+ AliTPCclusterMI *tempCl = trackIn.AddSpacePoint(AliTPCclusterMI());
+
+ const Float_t *xArr = pointArray0.GetX();
+ const Float_t *yArr = pointArray0.GetY();
+ const Float_t *zArr = pointArray0.GetZ();
+
+ const Float_t *xArrDist = pointArray1.GetX();
+ const Float_t *yArrDist = pointArray1.GetY();
+ const Float_t *zArrDist = pointArray1.GetZ();
+
+ tempCl->SetX(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;
+
+
+
+}
--- /dev/null
+#ifndef ToyMCEventGenerator_H
+#define ToyMCEventGenerator_H
+
+
+#include "ToyMCEvent.h"
+#include "ToyMCTrack.h"
+#include <AliTPCSpaceCharge3D.h>
+#include <AliTPCParam.h>
+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
+
--- /dev/null
+#include <iostream>
+#include "ToyMCEventGeneratorSimple.h"
+#include <TDatabasePDG.h>
+#include <TRandom.h>
+#include <TF1.h>
+#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; iTrack<nTracks; iTrack++){
+ Double_t phi = gRandom->Uniform(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;
+}
+
+
--- /dev/null
+#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
+
--- /dev/null
+#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);
+}
--- /dev/null
+#ifndef TOYMCTRACK_H
+#define TOYMCTRACK_H
+
+#include <AliExternalTrackParam.h>
+#include <TClonesArray.h>
+#include <AliTPCclusterMI.h>
+
+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<AliTPCclusterMI*> (fSpacePoints.At(spoint)); }
+ const AliTPCclusterMI* GetDistortedSpacePoint(Int_t dspoint) const { return static_cast<AliTPCclusterMI*> (fDistortedSpacePoints.At(dspoint)); }
+
+ AliTPCclusterMI* AddSpacePoint(const AliTPCclusterMI &spoint);
+ AliTPCclusterMI* AddDistortedSpacePoint(const AliTPCclusterMI &dspoint);
+
+ private:
+
+ TClonesArray fSpacePoints;
+ TClonesArray fDistortedSpacePoints;
+
+
+ ClassDef(ToyMCTrack,1)
+};
+
+#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff