123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692 |
- // C++ headers
- #include <iostream>
- #include <string>
- #include <utility>
- #include <set>
- #include <map>
- #include <functional>
- // ROOT headers
- #include <TROOT.h>
- #include <TFile.h>
- #include <TChain.h>
- #include <TTree.h>
- #include <TSystem.h>
- #include <TH1.h>
- #include <TH2.h>
- #include <TMath.h>
- #include <TString.h>
- #include <TStopwatch.h>
- // FemtoDst headers
- #include "StFemtoDstReader.h"
- #include "StFemtoDst.h"
- #include "StFemtoEvent.h"
- #include "StFemtoTrack.h"
- // AnalysisTree headers
- #include "AnalysisTree/Configuration.hpp"
- #include "AnalysisTree/DataHeader.hpp"
- #include "AnalysisTree/EventHeader.hpp"
- #include "AnalysisTree/Detector.hpp"
- #include "AnalysisTree/Matching.hpp"
- float get_beamP(float sqrtSnn, float m_target = 0.938, float m_beam = 0.938);
- float GetCentFromCent9(int cent9);
- float GetCentFromCent16(int cent16);
- TVector3 GetBbcPositionEast(int imod);
- TVector3 GetBbcPositionWest(int imod);
- TVector3 GetZdcPositionXEast(int strip);
- TVector3 GetZdcPositionYEast(int strip);
- TVector3 GetZdcPositionXWest(int strip);
- TVector3 GetZdcPositionYWest(int strip);
- int main(int argc, char **argv)
- {
- TString iFileName, oFileName;
- std::string system="";
- float sqrtSnn = -1.;
- bool isDataHeaderDefined = false;
- if (argc < 5)
- {
- std::cerr << "./MpdDst2AnalysisTree -i inputfile/inputlist -o outputfile [OPTIONAL: --sqrtSnn sqrtSnn --system colliding_system_name]" << std::endl;
- return 1;
- }
- for (int i = 1; i < argc; i++)
- {
- if (std::string(argv[i]) != "-i" &&
- std::string(argv[i]) != "-o" &&
- std::string(argv[i]) != "--sqrtSnn" &&
- std::string(argv[i]) != "--system")
- {
- std::cerr << "\n[ERROR]: Unknown parameter " << i << ": " << argv[i] << std::endl;
- return 1;
- }
- else
- {
- if (std::string(argv[i]) == "-i" && i != argc - 1)
- {
- iFileName = argv[++i];
- continue;
- }
- if (std::string(argv[i]) == "-i" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
- return 1;
- }
- if (std::string(argv[i]) == "-o" && i != argc - 1)
- {
- oFileName = argv[++i];
- continue;
- }
- if (std::string(argv[i]) == "-o" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
- return 1;
- }
- if (std::string(argv[i]) == "--system" && i != argc - 1)
- {
- system = std::string(argv[++i]);
- isDataHeaderDefined = true;
- continue;
- }
- if (std::string(argv[i]) == "--system" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
- return 1;
- }
- if (std::string(argv[i]) == "--sqrtSnn" && i != argc - 1)
- {
- sqrtSnn = std::stof(std::string(argv[++i]));
- isDataHeaderDefined = true;
- continue;
- }
- if (std::string(argv[i]) == "--sqrtSnn" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
- return 1;
- }
- }
- }
- const int NbbcModules = 32; // 0-15 - east, 16-31 - west
- const int NbbcModulesHalf = 16;
- const int NzdcStipsX = 7;
- const int NzdcStipsY = 8;
- TStopwatch timer;
- timer.Start();
- StFemtoDstReader* femtoReader = new StFemtoDstReader(iFileName.Data());
- femtoReader->Init();
- // This is a way if you want to spead up IO
- std::cout << "Explicit read status for some branches" << std::endl;
- femtoReader->SetStatus("*",0);
- femtoReader->SetStatus("Event",1);
- femtoReader->SetStatus("Track",1);
- std::cout << "Status has been set" << std::endl;
- if( !femtoReader->chain() ) {
- std::cout << "No chain has been found." << std::endl;
- }
- Long64_t eventsInTree = femtoReader->tree()->GetEntries();
- std::cout << "eventsInTree: " << eventsInTree << std::endl;
- Long64_t events2read = femtoReader->chain()->GetEntries();
- std::cout << "Number of events to read: " << events2read << std::endl;
- // Set up output dst
- TFile *outFile = new TFile(oFileName.Data(), "RECREATE");
- TTree *outTree = new TTree("aTree","AnalysisTree Dst at STAR");
- // Set up AnalysisTree data header
- AnalysisTree::DataHeader *dataHeader = new AnalysisTree::DataHeader;
- if (system != "") dataHeader->SetSystem(system);
- if (sqrtSnn > 0) dataHeader->SetBeamMomentum(get_beamP(sqrtSnn));
- auto &bbc_mod_pos = dataHeader->AddDetector();
- TVector3 modulePos;
- for (int imodule=0; imodule<NbbcModules; imodule++)
- {
- auto *module = bbc_mod_pos.AddChannel();
- if (imodule < NbbcModulesHalf) modulePos = GetBbcPositionEast(imodule);
- if (imodule >= NbbcModulesHalf) modulePos = GetBbcPositionWest(imodule-NbbcModulesHalf);
- module->SetPosition(modulePos);
- }
- // Set up AnalysisTree configureation
- AnalysisTree::Configuration *out_config = new AnalysisTree::Configuration;
- // Set up AnalysisTree configuration
- std::string str_reco_event_branch = "RecoEvent";
- std::string str_tpc_tracks_branch = "TpcTracks";
- std::string str_bbc_branch = "BbcModules";
- std::string str_zdc_horizontal_branch = "ZdcHitsHorizontal";
- std::string str_zdc_vertical_branch = "ZdcHitsVertical";
- AnalysisTree::BranchConfig reco_event_branch(str_reco_event_branch.c_str(), AnalysisTree::DetType::kEventHeader);
- AnalysisTree::BranchConfig tpc_tracks_branch(str_tpc_tracks_branch.c_str(), AnalysisTree::DetType::kTrack);
- AnalysisTree::BranchConfig bbc_branch(str_bbc_branch.c_str(), AnalysisTree::DetType::kModule);
- AnalysisTree::BranchConfig zdc_horizontal_branch(str_zdc_horizontal_branch.c_str(), AnalysisTree::DetType::kHit);
- AnalysisTree::BranchConfig zdc_vertical_branch(str_zdc_vertical_branch.c_str(), AnalysisTree::DetType::kHit);
- reco_event_branch.AddField<int>("runId");
- reco_event_branch.AddField<int>("refMult");
- reco_event_branch.AddField<int>("refMult2");
- reco_event_branch.AddField<float>("cent9");
- reco_event_branch.AddField<float>("cent16");
- reco_event_branch.AddField<int>("numberOfBTofHit");
- reco_event_branch.AddField<int>("numberOfTofMatched");
- reco_event_branch.AddField<float>("vpdVz");
- tpc_tracks_branch.AddField<int>("nHits");
- tpc_tracks_branch.AddField<int>("nHitsFit");
- tpc_tracks_branch.AddField<int>("nHitsPoss");
- tpc_tracks_branch.AddField<float>("nHits_Fit2Poss");
- tpc_tracks_branch.AddField<float>("nSigmaPion");
- tpc_tracks_branch.AddField<float>("nSigmaKaon");
- tpc_tracks_branch.AddField<float>("nSigmaProton");
- tpc_tracks_branch.AddField<float>("dEdx");
- tpc_tracks_branch.AddField<float>("tofMass2");
- tpc_tracks_branch.AddFields<float>({"dca_x","dca_y","dca_z"});
- tpc_tracks_branch.AddField<float>("dca");
- tpc_tracks_branch.AddField<bool>("isPrimary");
- tpc_tracks_branch.AddField<bool>("isTofTrack");
- auto hasher = std::hash<std::string>();
- // Initialize AnalysisTree Dst components
- out_config->AddBranchConfig(reco_event_branch);
- AnalysisTree::EventHeader *reco_event = new AnalysisTree::EventHeader( Short_t(hasher(reco_event_branch.GetName())) );
- out_config->AddBranchConfig(tpc_tracks_branch);
- AnalysisTree::TrackDetector *tpc_tracks = new AnalysisTree::TrackDetector( Short_t(hasher(tpc_tracks_branch.GetName())) );
- out_config->AddBranchConfig(bbc_branch);
- AnalysisTree::ModuleDetector *bbc_modules = new AnalysisTree::ModuleDetector( Short_t(hasher(bbc_branch.GetName())) );
- out_config->AddBranchConfig(zdc_horizontal_branch);
- AnalysisTree::HitDetector *zdc_horizontal_hits = new AnalysisTree::HitDetector( Short_t(hasher(zdc_horizontal_branch.GetName())) );
- out_config->AddBranchConfig(zdc_vertical_branch);
- AnalysisTree::HitDetector *zdc_vertical_hits = new AnalysisTree::HitDetector( Short_t(hasher(zdc_vertical_branch.GetName())) );
- reco_event->Init(reco_event_branch);
- // reco_event's additional field ids
- const int irunId = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("runId");
- const int irefMult = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("refMult");
- const int irefMult2 = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("refMult2");
- const int icent9 = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("cent9");
- const int icent16 = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("cent16");
- const int inumberOfBTofHit = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("numberOfBTofHit");
- const int inumberOfTofMatched = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("numberOfTofMatched");
- const int ivpdVz = out_config->GetBranchConfig(str_reco_event_branch).GetFieldId("vpdVz");
- // tpc_tracks' additional field ids
- const int inHits = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nHits");
- const int inHitsFit = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nHitsFit");
- const int inHitsPoss = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nHitsPoss");
- const int inHits_Fit2Poss = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nHits_Fit2Poss");
- const int inSigmaPion = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nSigmaPion");
- const int inSigmaKaon = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nSigmaKaon");
- const int inSigmaProton = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nSigmaProton");
- const int idEdx = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dEdx");
- const int itofMass2 = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("tofMass2");
- const int idca_x = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca_x");
- const int idca_y = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca_y");
- const int idca_z = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca_z");
- const int idca = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca");
- const int iisPrimary = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("isPrimary");
- const int iisTofTrack = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("isTofTrack");
- // Create branches in the output tree
- outTree->Branch(str_reco_event_branch.c_str(), "AnalysisTree::EventHeader", &reco_event, 32000, 99);
- outTree->Branch(str_tpc_tracks_branch.c_str(), "AnalysisTree::TrackDetector", &tpc_tracks, 256000, 99);
- outTree->Branch(str_bbc_branch.c_str(), "AnalysisTree::ModuleDetector", &bbc_modules, 128000, 99);
- outTree->Branch(str_zdc_horizontal_branch.c_str(), "AnalysisTree::HitDetector", &zdc_horizontal_hits, 128000, 99);
- outTree->Branch(str_zdc_vertical_branch.c_str(), "AnalysisTree::HitDetector", &zdc_vertical_hits, 128000, 99);
- // Printout basic configuration info
- out_config->Print();
- // Starting event loop
- TVector3 pVtx;
- // Loop over events
- for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
- {
- if (iEvent % 1000 == 0)
- {
- std::cout << "Working on event #[" << (iEvent+1)
- << "/" << events2read << "]" << std::endl;
- }
- Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
- if( !readEvent ) {
- std::cout << "Something went wrong! Nothing to analyze..." << std::endl;
- break;
- }
- // Retrieve femtoDst
- StFemtoDst *dst = femtoReader->femtoDst();
- // Retrieve event information
- StFemtoEvent *event = dst->event();
- if( !event ) {
- std::cout << "Something went wrong! Event is hiding from me..." << std::endl;
- break;
- }
- // Return primary vertex position
- pVtx = event->primaryVertex();
- // Reject vertices that are far from the central membrane along the beam
- if( TMath::Abs( pVtx.Z() ) > 200. ) continue;
- tpc_tracks->ClearChannels();
- bbc_modules->ClearChannels();
- zdc_horizontal_hits->ClearChannels();
- zdc_vertical_hits->ClearChannels();
- // Reading Reco Event
- reco_event->SetVertexPosition3(pVtx);
- reco_event->SetField(int(event->runId()), irunId);
- reco_event->SetField(int(event->refMult()), irefMult);
- reco_event->SetField(int(event->refMult2()), irefMult2);
- reco_event->SetField(float(GetCentFromCent9(event->cent9())), icent9);
- reco_event->SetField(float(GetCentFromCent16(event->cent16())), icent16);
- reco_event->SetField(int(event->numberOfBTofHit()), inumberOfBTofHit);
- reco_event->SetField(int(event->numberOfTofMatched()), inumberOfTofMatched);
- reco_event->SetField(float(event->vpdVz()), ivpdVz);
- // Track analysis
- Int_t nTracks = dst->numberOfTracks();
- tpc_tracks->Reserve(nTracks);
- // Track loop
- for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
- {
- // Retrieve i-th femto track
- StFemtoTrack *femtoTrack = dst->track(iTrk);
- if (!femtoTrack) continue;
- auto *track = tpc_tracks->AddChannel();
- track->Init(out_config->GetBranchConfig(str_tpc_tracks_branch));
- // Reading Reco Track
- track->SetMomentum(femtoTrack->pMom().X(),femtoTrack->pMom().Y(),femtoTrack->pMom().Z());
- track->SetField(int(femtoTrack->nHits()), inHits);
- track->SetField(int(femtoTrack->nHitsFit()), inHitsFit);
- track->SetField(int(femtoTrack->nHitsPoss()), inHitsPoss);
- track->SetField(float(femtoTrack->nHitsFit()/femtoTrack->nHitsPoss()), inHits_Fit2Poss);
- track->SetField(float(femtoTrack->nSigmaPion()), inSigmaPion);
- track->SetField(float(femtoTrack->nSigmaKaon()), inSigmaKaon);
- track->SetField(float(femtoTrack->nSigmaProton()), inSigmaProton);
- track->SetField(float(femtoTrack->dEdx()), idEdx);
- track->SetField(float(femtoTrack->massSqr()), itofMass2);
- track->SetField(float(femtoTrack->gDCA(pVtx).Mag()), idca);
- track->SetField(float(femtoTrack->gDCA(pVtx).X()), idca_x);
- track->SetField(float(femtoTrack->gDCA(pVtx).Y()), idca_y);
- track->SetField(float(femtoTrack->gDCA(pVtx).Z()), idca_z);
- track->SetField(bool(femtoTrack->isPrimary()), iisPrimary);
- track->SetField(bool(femtoTrack->isTofTrack()), iisTofTrack);
- } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
- // Fill BBC information
- bbc_modules->Reserve(NbbcModules);
- for (int imodule = 0; imodule<NbbcModules; imodule++)
- {
- auto *module = bbc_modules->AddChannel();
- module->Init(out_config->GetBranchConfig(str_bbc_branch));
- module->SetNumber(imodule);
- if (imodule < NbbcModulesHalf) module->SetSignal(event->bbcAdcEast(imodule));
- if (imodule >= NbbcModulesHalf) module->SetSignal(event->bbcAdcWest(imodule-NbbcModulesHalf));
- }
- // Fill ZDC-SMD information
- // There are 7 vertical and 8 horizontal strips of ZDC-SMD. To calculate Q-vectors from ZDC
- // one has to sum up signal from those strips separately:
- // Qx = sum zdcSmdVertical * position_of_the_vertical_strip, (sum over 7 strips)
- // Qy = sum zdcSmdHorizontal * position_of_the_Horizontal_strip, (sum over 8 strips)
- // To translate this information into AnalysisTree, HitDetector was chosen for horizontal
- // and vertical strips separately.
- zdc_horizontal_hits->Reserve(2*NzdcStipsY);
- zdc_vertical_hits->Reserve(2*NzdcStipsX);
- TVector3 geomZdcPos;
- // Fill ZDC-SMD East - horizontal
- for (int ihity=0; ihity<NzdcStipsY; ihity++)
- {
- auto *hit = zdc_horizontal_hits->AddChannel();
- hit->Init(out_config->GetBranchConfig(str_zdc_horizontal_branch));
- geomZdcPos = GetZdcPositionYEast(ihity);
- hit->SetPosition(geomZdcPos);
- hit->SetSignal(event->zdcSmdEastHorizontal(ihity));
- }
- // Fill ZDC-SMD East - vertical
- for (int ihitx=0; ihitx<NzdcStipsX; ihitx++)
- {
- auto *hit = zdc_vertical_hits->AddChannel();
- hit->Init(out_config->GetBranchConfig(str_zdc_vertical_branch));
- geomZdcPos = GetZdcPositionXEast(ihitx);
- hit->SetPosition(geomZdcPos);
- hit->SetSignal(event->zdcSmdEastVertical(ihitx));
- }
- // Fill ZDC-SMD West - horizontal
- for (int ihity=0; ihity<NzdcStipsY; ihity++)
- {
- auto *hit = zdc_horizontal_hits->AddChannel();
- hit->Init(out_config->GetBranchConfig(str_zdc_horizontal_branch));
- geomZdcPos = GetZdcPositionYWest(ihity);
- hit->SetPosition(geomZdcPos);
- hit->SetSignal(event->zdcSmdWestHorizontal(ihity));
- }
- // Fill ZDC-SMD West - vertical
- for (int ihitx=0; ihitx<NzdcStipsX; ihitx++)
- {
- auto *hit = zdc_vertical_hits->AddChannel();
- hit->Init(out_config->GetBranchConfig(str_zdc_vertical_branch));
- geomZdcPos = GetZdcPositionXWest(ihitx);
- hit->SetPosition(geomZdcPos);
- hit->SetSignal(event->zdcSmdWestVertical(ihitx));
- }
- outTree->Fill();
- } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
- outFile->cd();
- outTree->Print();
- outTree->Write();
- if (isDataHeaderDefined) dataHeader->Write("DataHeader");
- out_config->Write("Configuration");
- outFile->Close();
- femtoReader->Finish();
- timer.Stop();
- timer.Print();
- return 0;
- }
- float get_beamP(float sqrtSnn, float m_target, float m_beam)
- {
- return sqrt( pow((pow(sqrtSnn,2) - pow(m_target,2) - pow(m_beam,2))/(2*m_target), 2) - pow(m_beam,2) );
- }
- float GetCentFromCent9(int cent9)
- {
- if (cent9 == 0) return 75.;
- if (cent9 == 1) return 65.;
- if (cent9 == 2) return 55.;
- if (cent9 == 3) return 45.;
- if (cent9 == 4) return 35.;
- if (cent9 == 5) return 25.;
- if (cent9 == 6) return 15.;
- if (cent9 == 7) return 7.5;
- if (cent9 == 8) return 2.5;
- return -1.;
- }
- float GetCentFromCent16(int cent16)
- {
- if (cent16 == 0) return 77.5;
- if (cent16 == 1) return 72.5;
- if (cent16 == 2) return 67.5;
- if (cent16 == 3) return 62.5;
- if (cent16 == 4) return 57.5;
- if (cent16 == 5) return 52.5;
- if (cent16 == 6) return 47.5;
- if (cent16 == 7) return 42.5;
- if (cent16 == 8) return 37.5;
- if (cent16 == 9) return 32.5;
- if (cent16 == 10) return 27.5;
- if (cent16 == 11) return 22.5;
- if (cent16 == 12) return 17.5;
- if (cent16 == 13) return 12.5;
- if (cent16 == 14) return 7.5;
- if (cent16 == 15) return 2.5;
- return -1.;
- }
- //--------------------------------------------------------------------------------------------------------------------------//
- //this function simply connects the gain values read in to the BBC azimuthal distribution
- //since tiles 7 and 9 (+ 13 and 15) share a gain value it is ambiguous how to assign the geometry here
- //It is preferred sometimes assigning the angle between the tiles thus "greying out" the adcs.
- //Others have assigned all of the adc to one (exclusive) or the the other.
- TVector3 GetBbcPositionEast(int imod)
- {
- double z_pos = -375.; // in cm
- double radius=0.;;
- // Radius here is an artificial value just to distinguish between different rings of BBC modules.
- // It does not represent the true geometrical positions of the modules themself.
- if (imod >=0 && imod < 6) radius = 1.;
- if (imod >=6 && imod < 16) radius = 2.;
- //float GetPhiInBBC(int eastWest, int bbcN) { //imod=0 to 23
- const float Pi = TMath::Pi();
- const float Phi_div = Pi / 6;
- float bbc_phi = Phi_div;
- switch(imod) {
- case 0: bbc_phi = 3.*Phi_div;
- break;
- case 1: bbc_phi = Phi_div;
- break;
- case 2: bbc_phi = -1.*Phi_div;
- break;
- case 3: bbc_phi = -3.*Phi_div;
- break;
- case 4: bbc_phi = -5.*Phi_div;
- break;
- case 5: bbc_phi = 5.*Phi_div;
- break;
- //case 6: bbc_phi= (mRndm.Rndm() > 0.5) ? 2.*Phi_div:4.*Phi_div; //tiles 7 and 9 are gained together we randomly assign the gain to one XOR the other
- case 6: bbc_phi = 3.*Phi_div;
- break;
- case 7: bbc_phi = 3.*Phi_div;
- break;
- case 8: bbc_phi = Phi_div;
- break;
- case 9: bbc_phi = 0.;
- break;
- case 10: bbc_phi = -1.*Phi_div;
- break;
- //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div; //tiles 13 and 15 are gained together
- case 11: bbc_phi = -3.*Phi_div;
- break;
- case 12: bbc_phi = -3.*Phi_div;
- break;
- case 13: bbc_phi = -5.*Phi_div;
- break;
- case 14: bbc_phi = Pi;
- break;
- case 15: bbc_phi = 5.*Phi_div;
- break;
- }
- //if we're looking at the east BBC we need to flip around x in the STAR coordinates,
- //a line parallel to the beam would go through tile 1 on the W BBC and tile 3 on the
- // E BBC
- if (bbc_phi > -0.001){ //this is not a >= since we are talking about finite adcs -- not to important
- bbc_phi = Pi - bbc_phi;
- }
- else {
- bbc_phi= -Pi - bbc_phi;
- }
- if(bbc_phi < 0.0) bbc_phi += 2.*Pi;
- if(bbc_phi > 2.*Pi) bbc_phi -= 2.*Pi;
- TVector3 vresult;
- vresult.SetXYZ(radius*TMath::Cos(bbc_phi), radius*TMath::Sin(bbc_phi), z_pos);
- return vresult;
- }
- TVector3 GetBbcPositionWest(int imod)
- {
- double z_pos = 375.; // in cm
- double radius=0.;;
- // Radius here is an artificial value just to distinguish between different rings of BBC modules.
- // It does not represent the true geometrical positions of the modules themself.
- if (imod >=0 && imod < 6) radius = 1.;
- if (imod >=6 && imod < 16) radius = 2.;
- //float GetPhiInBBC(int eastWest, int bbcN) { //imod=0 to 23
- const float Pi = TMath::Pi();
- const float Phi_div = Pi / 6;
- float bbc_phi = Phi_div;
- switch(imod) {
- case 0: bbc_phi = 3.*Phi_div;
- break;
- case 1: bbc_phi = Phi_div;
- break;
- case 2: bbc_phi = -1.*Phi_div;
- break;
- case 3: bbc_phi = -3.*Phi_div;
- break;
- case 4: bbc_phi = -5.*Phi_div;
- break;
- case 5: bbc_phi = 5.*Phi_div;
- break;
- //case 6: bbc_phi= (mRndm.Rndm() > 0.5) ? 2.*Phi_div:4.*Phi_div; //tiles 7 and 9 are gained together we randomly assign the gain to one XOR the other
- case 6: bbc_phi = 3.*Phi_div;
- break;
- case 7: bbc_phi = 3.*Phi_div;
- break;
- case 8: bbc_phi = Phi_div;
- break;
- case 9: bbc_phi = 0.;
- break;
- case 10: bbc_phi = -1.*Phi_div;
- break;
- //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div; //tiles 13 and 15 are gained together
- case 11: bbc_phi = -3.*Phi_div;
- break;
- case 12: bbc_phi = -3.*Phi_div;
- break;
- case 13: bbc_phi = -5.*Phi_div;
- break;
- case 14: bbc_phi = Pi;
- break;
- case 15: bbc_phi = 5.*Phi_div;
- break;
- }
- if(bbc_phi < 0.0) bbc_phi += 2.*Pi;
- if(bbc_phi > 2.*Pi) bbc_phi -= 2.*Pi;
- TVector3 vresult;
- vresult.SetXYZ(radius*TMath::Cos(bbc_phi), radius*TMath::Sin(bbc_phi), z_pos);
- return vresult;
- }
- // Get position of each slat;strip starts from 0
- TVector3 GetZdcPositionXEast(int strip)
- {
- TVector3 vresult;
- double x_pos, y_pos;
- double z_pos = -1800.; // in cm
- Double_t zdcsmd_x[7] = {0.5,2,3.5,5,6.5,8,9.5};
- Double_t mZDCSMDCenterex = 4.72466;
- Double_t mZDCSMDCenterey = 5.53629;
- y_pos = mZDCSMDCenterey;
- if (strip>=0 && strip<7)
- {
- x_pos = zdcsmd_x[strip]-mZDCSMDCenterex;
- }
- else
- {
- x_pos = 0.;
- }
- vresult.SetXYZ(x_pos, y_pos, z_pos);
- return vresult;
- }
- TVector3 GetZdcPositionYEast(int strip)
- {
- TVector3 vresult;
- double x_pos, y_pos;
- double z_pos = -1800.; // in cm
- Double_t zdcsmd_y[8] = {1.25,3.25,5.25,7.25,9.25,11.25,13.25,15.25};
- Double_t mZDCSMDCenterex = 4.72466;
- Double_t mZDCSMDCenterey = 5.53629;
- x_pos = mZDCSMDCenterex;
- if (strip>=0 && strip<8)
- {
- y_pos = zdcsmd_y[strip]/TMath::Sqrt(2.)-mZDCSMDCenterey;
- }
- else
- {
- x_pos = 0.;
- }
- vresult.SetXYZ(x_pos, y_pos, z_pos);
- return vresult;
- }
- // Get position of each slat;strip starts from 0
- TVector3 GetZdcPositionXWest(int strip)
- {
- TVector3 vresult;
- double x_pos, y_pos;
- double z_pos = 1800.; // in cm
- Double_t zdcsmd_x[7] = {0.5,2,3.5,5,6.5,8,9.5};
- Double_t mZDCSMDCenterwx = 4.39604;
- Double_t mZDCSMDCenterwy = 5.19968;
- y_pos = mZDCSMDCenterwy;
- if (strip>=0 && strip<7)
- {
- x_pos = mZDCSMDCenterwx-zdcsmd_x[strip];
- }
- else
- {
- x_pos = 0.;
- }
- vresult.SetXYZ(x_pos, y_pos, z_pos);
- return vresult;
- }
- // Get position of each slat;strip starts from 0
- TVector3 GetZdcPositionYWest(int strip)
- {
- TVector3 vresult;
- double x_pos, y_pos;
- double z_pos = 1800.; // in cm
- Double_t zdcsmd_y[8] = {1.25,3.25,5.25,7.25,9.25,11.25,13.25,15.25};
- Double_t mZDCSMDCenterwx = 4.39604;
- Double_t mZDCSMDCenterwy = 5.19968;
- x_pos = mZDCSMDCenterwx;
- if (strip>=0 && strip<8)
- {
- y_pos = zdcsmd_y[strip]/sqrt(2.)-mZDCSMDCenterwy;
- }
- else
- {
- y_pos = 0.;
- }
- vresult.SetXYZ(x_pos, y_pos, z_pos);
- return vresult;
- }
|