#include "glpk.h" #include "PndSttMvdTracking.h" #include "PndSttHit.h" #include "PndSciTHit.h" #include "PndSttPoint.h" #include "PndSttTrack.h" #include "PndSttPoint.h" #include "PndSttHelixHit.h" #include "PndSttSingleStraw.h" #include "PndSttTube.h" #include "PndSttMapCreator.h" #include "PndSdsHit.h" #include "PndSdsMCPoint.h" #include "PndTrackCand.h" #include "PndTrackCandHit.h" #include "PndTrack.h" #include "FairRootManager.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "FairTrackParP.h" #include "TGeoManager.h" #include "TClonesArray.h" #include "TGeoVolume.h" #include "TVector3.h" #include "TRandom.h" #include "TH1F.h" #include "TMath.h" #include "TCanvas.h" #include "TGeoTube.h" #include #include using namespace std; // ----- Default constructor ------------------------------------------- PndSttMvdTracking::PndSttMvdTracking() : FairTask("STT Stt-Mvd Tracking") { fPersistence = kTRUE; fVerbose = 0; istampa = 0; iplotta = false; doMcComparison = false; YesClean = false; YesCleanMvd = true; YesSciTil = false ; MvdAloneTracking = true; Initialization_ClassVariables(); sprintf(fSttBranch,"STTHit"); sprintf(fMvdPixelBranch,"MVDHitsPixel"); sprintf(fMvdStripBranch,"MVDHitsStrip"); } // ------------------------------------------------------------------------- PndSttMvdTracking::PndSttMvdTracking(Int_t verbose) : FairTask("STT Stt-Mvd Tracking") { fPersistence = kTRUE; fVerbose = verbose; istampa = verbose; iplotta = false; doMcComparison = false; YesClean = false; YesCleanMvd = true; YesSciTil = false ; MvdAloneTracking = true; Initialization_ClassVariables(); sprintf(fSttBranch,"STTHit"); sprintf(fMvdPixelBranch,"MVDHitsPixel"); sprintf(fMvdStripBranch,"MVDHitsStrip"); } // ------------------------------------------------------------------------- PndSttMvdTracking::PndSttMvdTracking(int istamp, bool iplot, bool imc) : FairTask("STT Stt-Mvd Tracking") { fPersistence = kTRUE; istampa = istamp; iplotta = iplot; doMcComparison = imc; YesClean = false; YesCleanMvd = true; YesSciTil = false ; MvdAloneTracking = true; Initialization_ClassVariables(); sprintf(fSttBranch,"STTHit"); sprintf(fMvdPixelBranch,"MVDHitsPixel"); sprintf(fMvdStripBranch,"MVDHitsStrip"); } // ------------------------------------------------------------------------- PndSttMvdTracking::PndSttMvdTracking(int istamp, bool iplot, bool imc, bool doSciTil) : FairTask("STT Stt-Mvd Tracking") { fPersistence = kTRUE; istampa = istamp; iplotta = iplot; doMcComparison = imc; YesClean = false; YesCleanMvd = true; YesSciTil = doSciTil ; MvdAloneTracking = true; Initialization_ClassVariables(); sprintf(fSttBranch,"STTHit"); sprintf(fMvdPixelBranch,"MVDHitsPixel"); sprintf(fMvdStripBranch,"MVDHitsStrip"); } // ----- Destructor ---------------------------------------------------- PndSttMvdTracking::~PndSttMvdTracking() { } // ------------------------------------------------------------------------- //--------------- begin PndSttMvdTracking::Initialization_ClassVariables void PndSttMvdTracking::Initialization_ClassVariables() { // this is only for initializing the Class Variables. size_t len; // booleans : len = sizeof(InclusionListSciTil); memset (InclusionListSciTil,0,len); len = sizeof(InclusionListStt); memset (InclusionListStt,0,len); // char : len = sizeof(fSttBranch); memset (fSttBranch,0,len); len = sizeof(fMvdPixelBranch); memset (fMvdPixelBranch,0,len); len = sizeof(fMvdStripBranch); memset (fMvdStripBranch,0,len); // Short_t : nSciTilHits=0; len = sizeof(nTrackCandHit); memset (nTrackCandHit,0,len); len = sizeof(nSciTilHitsinTrack); memset (nSciTilHitsinTrack,0,len); len = sizeof(nSttParHitsinTrack); memset (nSttParHitsinTrack,0,len); len = sizeof(nSttSkewHitsinTrack); memset (nSttSkewHitsinTrack,0,len); len = sizeof(nMvdPixelHitsinTrack); memset (nMvdPixelHitsinTrack,0,len); len = sizeof(nMvdStripHitsinTrack); memset (nMvdStripHitsinTrack,0,len); len = sizeof(ListMvdPixelHitsinTrack); memset (ListMvdPixelHitsinTrack,0,len); len = sizeof(ListMvdStripHitsinTrack); memset (ListMvdStripHitsinTrack,0,len); len = sizeof(ListSciTilHitsinTrack); memset (ListSciTilHitsinTrack,0,len); len = sizeof(ListTrackCandHit); memset (ListTrackCandHit,0,len); len = sizeof(ListSttParHitsinTrack); memset (ListSttParHitsinTrack,0,len); len = sizeof(ListSttSkewHitsinTrack); memset (ListSttSkewHitsinTrack,0,len); len = sizeof(ListSttSkewHitsinTrackSolution); memset (ListSttSkewHitsinTrackSolution,0,len); // Short_t : nMvdPixelHit=0; nMvdStripHit=0; nMvdTrackCand=0; nMvdDSPixelHitNotTrackCand=0; nMvdUSPixelHitNotTrackCand=0; nMvdDSStripHitNotTrackCand=0; nMvdUSStripHitNotTrackCand=0; len = sizeof(nHitMvdTrackCand); memset (nHitMvdTrackCand,0,len); len = sizeof(ListHitMvdTrackCand); memset (ListHitMvdTrackCand,0,len); len = sizeof(ListHitTypeMvdTrackCand); memset (ListHitTypeMvdTrackCand,0,len); len = sizeof(ListMvdDSPixelHitNotTrackCand); memset (ListMvdDSPixelHitNotTrackCand,0,len); len = sizeof(ListMvdUSPixelHitNotTrackCand); memset (ListMvdUSPixelHitNotTrackCand,0,len); len = sizeof(ListMvdDSStripHitNotTrackCand); memset (ListMvdDSStripHitNotTrackCand,0,len); len = sizeof(ListMvdUSStripHitNotTrackCand); memset (ListMvdUSStripHitNotTrackCand,0,len); len = sizeof(ListTrackCandHitType); memset (ListTrackCandHitType,0,len); // Double_t : SEMILENGTH_STRAIGHT=0.; ZCENTER_STRAIGHT=0.; len = sizeof(CxMC); memset (CxMC,0,len); len = sizeof(CyMC); memset (CyMC,0,len); len = sizeof(R_MC); memset (R_MC,0,len); len = sizeof(MCtruthTrkInfo); memset (MCtruthTrkInfo,0,len); len = sizeof(MCSkewAloneX); memset (MCSkewAloneX,0,len); len = sizeof(MCSkewAloneY); memset (MCSkewAloneY,0,len); len = sizeof(XMvdPixel); memset (XMvdPixel,0,len); len = sizeof(YMvdPixel); memset (YMvdPixel,0,len); len = sizeof(ZMvdPixel); memset (ZMvdPixel,0,len); len = sizeof(sigmaXMvdPixel); memset (sigmaXMvdPixel,0,len); len = sizeof(sigmaYMvdPixel); memset (sigmaYMvdPixel,0,len); len = sizeof(sigmaZMvdPixel); memset (sigmaZMvdPixel,0,len); len = sizeof(refindexMvdPixel); memset (refindexMvdPixel,0,len); len = sizeof(XMvdStrip); memset (XMvdStrip,0,len); len = sizeof(YMvdStrip); memset (YMvdStrip,0,len); len = sizeof(ZMvdStrip); memset (ZMvdStrip,0,len); len = sizeof(sigmaXMvdStrip); memset (sigmaXMvdStrip,0,len); len = sizeof(sigmaYMvdStrip); memset (sigmaYMvdStrip,0,len); len = sizeof(sigmaZMvdStrip); memset (sigmaZMvdStrip,0,len); len = sizeof(refindexMvdStrip); memset (refindexMvdStrip,0,len); len = sizeof(ALFA); memset (ALFA,0,len); len = sizeof(BETA); memset (BETA,0,len); len = sizeof(GAMMA); memset (GAMMA,0,len); len = sizeof(posizSciTil); memset (posizSciTil,0,len); len = sizeof(SciTilHitsXwithTrack); memset (SciTilHitsXwithTrack,0,len); len = sizeof(SciTilHitsYwithTrack); memset (SciTilHitsYwithTrack,0,len); // puntatori : HANDLE=NULL; HANDLE2=NULL; HANDLEXYZ=NULL; PHANDLEX=NULL; PHANDLEY=NULL; PHANDLEZ=NULL; SHANDLEX=NULL; SHANDLEY=NULL; SHANDLEZ=NULL; hdeltaRPixel=NULL; hdeltaRStrip=NULL; hdeltaRPixel2=NULL; hdeltaRStrip2=NULL; fMCTrackArray=NULL; fSttTubeArray=NULL; fSttPointArray=NULL; fSttHitArray=NULL; fSttTrackArray=NULL; fSttTrackCandArray=NULL; fMvdPixelHitArray=NULL; fMvdStripHitArray=NULL; fMvdTrackCandArray=NULL; fSciTHitArray=NULL; fMvdMCPointArray=NULL; fSttMvdPndTrackCandArray=NULL; fSttMvdPndTrackArray=NULL; fSttParameters=NULL; } //--------------- end PndSttMvdTracking::Initialization_ClassVariables // ----- Public method Init -------------------------------------------- InitStatus PndSttMvdTracking::Init() { if(iplotta){ hdeltaRPixel = new TH1F("hdeltaRPixel", "distance MC Pixel point from trajectory in XY plane", 100, -1, 1); hdeltaRStrip = new TH1F("hdeltaRStrip", "distance MC Strip point from trajectory in XY plane", 100, -1, 1); hdeltaRPixel2 = new TH1F( "hdeltaRPixel2", "distance MC point from trajectory in XY plane (Pixels)", 100, -10, 10); hdeltaRStrip2 = new TH1F( "hdeltaRStrip2", "distance MC point from trajectory in XY plane (Strips)", 100, -10,10); } IVOLTE=-1; // --------------------------- opening files for special purposes if(istampa >=1 ){ //---- apertura file con info su Found tracce su cui si fa Helix fit dopo HANDLE2 = fopen("info_da_PndTrackFinderReal.txt","w"); // ---- open filehandle per statistica sugli hits etc. HANDLE = fopen("statistiche.txt","w"); // --------------- // HANDLEXYZ = fopen("infoPndTrackFinderRealXYZ.txt","w"); // --------------- open file delle info su deltaX, Y, Z degli hits in comune tra tracce trovate e MC PHANDLEX = fopen("deltaParXmio.txt","w"); PHANDLEY = fopen("deltaParYmio.txt","w"); PHANDLEZ = fopen("deltaParZmio.txt","w"); SHANDLEX = fopen("deltaSkewXmio.txt","w"); SHANDLEY = fopen("deltaSkewYmio.txt","w"); SHANDLEZ = fopen("deltaSkewZmio.txt","w"); } // end of if(istampa >=1) // ------------------------- // Get RootManager FairRootManager* ioman = FairRootManager::Instance(); if ( ! ioman ) { cout << "-E- PndSttMvdTracking::Init: " << "RootManager not instantiated, return!" << endl; return kFATAL; } // ----- maps of STT tubes // CHECK added PndSttMapCreator *mapper = new PndSttMapCreator(fSttParameters); fSttTubeArray = mapper->FillTubeArray(); //---------------------------------------------------- end map // get the MCTrack array fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack"); if ( ! fMCTrackArray) { cout << "-E- PndSttMvdTracking::Init: No MCTrack array, return!" << endl; return kERROR; } // open STTFoundTrack array fSttTrackArray = (TClonesArray*) ioman->GetObject("STTFoundTrack"); if(!fSttTrackArray) { cout << "-E- PndSttMvdTracking::Init: No STTFoundTrack array, return!" << endl; return kERROR; } // Get SttTrackCand array dal pattern recognition di STT fSttTrackCandArray = (TClonesArray*) ioman->GetObject("STTTrackCand"); if ( ! fSttTrackCandArray) { cout << "-E- PndSttMvdTracking::Init: No SttTrack Cand array, return!" << endl; return kERROR; } // Get input array questi sono i MC point di STT fSttPointArray = (TClonesArray*) ioman->GetObject("STTPoint"); if ( ! fSttPointArray ) { cout << "-W- PndSttHelixHitProducer::Init: " << "No STTPoint array, return!" << endl; return kERROR; } // Get input array hit di STT after digi fSttHitArray = (TClonesArray*) ioman->GetObject(fSttBranch); // fSttHitArray = (TClonesArray*) ioman->GetObject("STTHit"); if ( ! fSttHitArray ) { cout << "-W- PndSttMvdTracking::Init: " << "No STTHit array, return!" << endl; return kERROR; } // Create and register output array for PndTrackCand of Stt+Mvd combined fSttMvdPndTrackCandArray = new TClonesArray("PndTrackCand"); ioman->Register("SttMvdTrackCand","SttMvd",fSttMvdPndTrackCandArray, kTRUE); // Create and register output array for PndTrack of Stt+Mvd combined fSttMvdPndTrackArray = new TClonesArray("PndTrack"); ioman->Register("SttMvdTrack","SttMvd",fSttMvdPndTrackArray, kTRUE); // ------------------------- get the Mvd hits fMvdPixelHitArray = (TClonesArray*) ioman->GetObject(fMvdPixelBranch); // fMvdPixelHitArray = (TClonesArray*) ioman->GetObject("MVDHitsPixel"); if ( !fMvdPixelHitArray){ std::cout << "-W- PndSttMvdTracking::Init: " << "No MVD Pixel hitArray, return!" << std::endl; return kERROR; } fMvdStripHitArray = (TClonesArray*) ioman->GetObject(fMvdStripBranch); // fMvdStripHitArray = (TClonesArray*) ioman->GetObject("MVDHitsStrip"); if ( !fMvdStripHitArray){ std::cout << "-W- PndSttMvdTracking::Init: " << "No MVD Strip hitArray, return!" << std::endl; return kERROR; } // ------------------------- get the Mvd track candidates fMvdTrackCandArray = (TClonesArray*) ioman->GetObject("MVDRiemannTrackCand"); if ( !fMvdTrackCandArray){ std::cout << "-W- PndSttMvdTracking::Init: " << "No MVD TrackCand Array, return!" << std::endl; return kERROR; } cout << "-I- PndSttMvdTracking: Initialization successfull" << endl; // ------------------------- get the Mvd MC points fMvdMCPointArray = (TClonesArray*) ioman->GetObject("MVDPoint"); if ( !fMvdMCPointArray){ std::cout << "-W- PndSttMvdTracking::Init: " << "No MVD MC Point Array, return!" << std::endl; return kERROR; } // ------------------------- get the SciTil hits if(YesSciTil) { fSciTHitArray = (TClonesArray*) ioman->GetObject("SciTHit"); } else { fSciTHitArray = NULL; } //--------------------------- cout << "-I- PndSttMvdTracking: Initialization successfull" << endl; return kSUCCESS; } // ------------------------------------------------------------------------- void PndSttMvdTracking::SetParContainers() { FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb(); fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar"); } void PndSttMvdTracking::WriteHistograms(){ /* TFile* file = FairRootManager::Instance()->GetOutFile(); file->cd(); file->mkdir("PndSttHelixHit"); file->cd("PndSttHelixHit"); hx->Write(); delete hx; */ if(iplotta){ TFile* file = FairRootManager::Instance()->GetOutFile(); file->cd(); file->mkdir("PndSttMvdTracking"); file->cd("PndSttMvdTracking"); hdeltaRPixel->Write(); hdeltaRStrip->Write(); hdeltaRPixel2->Write(); hdeltaRStrip2->Write(); delete hdeltaRPixel; delete hdeltaRStrip; delete hdeltaRPixel2; delete hdeltaRStrip2; } } // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- void PndSttMvdTracking::Exec(Option_t* opt) { bool flag, intersect, Status, keepit[MAXTRACKSPEREVENT], status[MAXTRACKSPEREVENT], SttSZfit[MAXTRACKSPEREVENT], Mvdhits[MAXTRACKSPEREVENT]; Short_t Candidato, Charge, statusflag[MAXTRACKSPEREVENT], iHit, i, j, k, m, resultFitSZagain[MAXTRACKSPEREVENT], tubeID ; Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT], CHARGE[MAXTRACKSPEREVENT]; Short_t ipinco, kall, l, n, nalone, ncand, nhitsinfit, Nint, nMvdMCPoint, nRemainingCandidates, nSttHit, nSttParHit, nSttSkewHit, nSttHelixTrack, nSttMCPoint, nSttTrackCand, nTotalCandidates, nXYZhits, // TemporarynSkewHitsinTrack, FromHitToMCTrack[nmaxSttHits], ListAllParHits[nmaxSttHits], ListAllSkewHits[nmaxSttHits], nHitsInMCTrack[MAXTRACKSPEREVENT], nMCParalAlone[MAXTRACKSPEREVENT], nMCSkewAlone[MAXTRACKSPEREVENT], nParalCommon[MAXTRACKSPEREVENT], nSkewCommon[MAXTRACKSPEREVENT], nSkewHitsInMCTrack[MAXTRACKSPEREVENT], nSpuriParinTrack[MAXTRACKSPEREVENT], nSpuriSkewinTrack[MAXTRACKSPEREVENT], nSttHitsinTrack[MAXTRACKSPEREVENT], TemporarySkewList[2*nmaxSttHits][2], ListSttHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], ListSttHitsinTrackType[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], MCParalAloneList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], MCSkewAloneList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], ParalCommonList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], ParSpuriList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], SkewCommonList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], SkewSpuriList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack]; Int_t iaccept, nrounds0, nrounds1, ipunto, ListHits[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack]; Double_t Dist, Distance, Fi_final_helix_referenceframe, Ntras, posx, posy, Phi, Ptras, Rad, RR, SIGN, sqrtRR, TanL, Z, ddd, dis, dista, dista0, dista0_1, dista1, dista1_1, delta, emme, gap, highqualitycut, Pxini, px, Pyini, py, Pzini, qop, x, y, AloneX[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], AloneY[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], HoughFi[MAXTRACKSPEREVENT], KAPPA[MAXTRACKSPEREVENT], primoangolo[MAXTRACKSPEREVENT], s[2], Start[3], Trajectory_Start[MAXTRACKSPEREVENT][2], ultimoangolo[MAXTRACKSPEREVENT], versor[2], XintersectionList[2], YintersectionList[2], z[2], zeta0, zeta1, zdrift[2], zerror[2], WDX[nmaxSttHits], WDY[nmaxSttHits], WDZ[nmaxSttHits], Ox[MAXTRACKSPEREVENT], Oy[MAXTRACKSPEREVENT], R[MAXTRACKSPEREVENT], FI0[MAXTRACKSPEREVENT], Fifirst[MAXTRACKSPEREVENT], Fi_low_limit[MAXTRACKSPEREVENT], Fi_up_limit[MAXTRACKSPEREVENT], Posiz1[3], Posiz2[3], Px[MAXTRACKSPEREVENT], Py[MAXTRACKSPEREVENT], Pz[MAXTRACKSPEREVENT], X[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack+nmaxSciTilHitsinTrack], Y[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack+nmaxSciTilHitsinTrack], Sbis[2*nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack +nmaxMvdStripHitsInTrack+nmaxSciTilHitsinTrack], // multiplication by 2 in the ZEDbis[2*nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack +nmaxMvdStripHitsInTrack+nmaxSciTilHitsinTrack], // rather improbable chance that DriftRadiusbis[2*nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+ nmaxMvdStripHitsInTrack+nmaxSciTilHitsinTrack],// all skew hits have double ErrorDriftRadiusbis[2*nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack+2],// solutions SchosenPixel[MAXTRACKSPEREVENT][nmaxMvdPixelHits], SchosenStrip[MAXTRACKSPEREVENT][nmaxMvdStripHits], SchosenSkew[MAXTRACKSPEREVENT][nmaxSttHits], // NO multiplication by 2 here because for the // skew hits only one // solution is selected. XY[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack][2], ZchosenPixel[MAXTRACKSPEREVENT][nmaxMvdPixelHits], ZchosenStrip[MAXTRACKSPEREVENT][nmaxMvdStripHits], ZchosenSkew[MAXTRACKSPEREVENT][nmaxSttHits], ErrorchosenPixel[nmaxMvdPixelHits], ErrorchosenStrip[nmaxMvdStripHits], ErrorchosenSkew[nmaxSttHits], info[nmaxSttHits][7], TemporaryS[2*nmaxSttHits], TemporaryZ[2*nmaxSttHits], TemporaryZDrift[2*nmaxSttHits], TemporaryZErrorafterTilt[2*nmaxSttHits]; TVector3 Momentum,ErrMomentum,Position,ErrPosition; // FairMCPoint * pSttMCPoint[nmaxSttHits]; FairMCPoint * puntator; PndTrackCand * pMvdTrackCand, * pSttTrackCand; PndTrack * pSttTrack; PndSttTrack * pSttHelixTrack; PndSttTube * pSttTube; PndSttHit * pSttHit; PndSdsHit * pMvdPixelHit, * pMvdStripHit; // PndSdsMCPoint * pMvdMCPoint; PndMCTrack* pMCtr; PndTrackCandHit pndtrackcandhit; fSttMvdPndTrackCandArray->Delete(); fSttMvdPndTrackArray->Delete(); IVOLTE++; if(istampa>0) {cout<<"\n\nEntering in PndSttMvdTracking, IVOLTE = "<GetEntriesFast(); // num. tracce/evento if (istampa > 1 && IVOLTE<20) { cout<<"da PndSttMvdTracking : evento (partendo da 0) N. "<< IVOLTE<< "\n N. di MC truth tracks : "< MAXMCTRACKS){ cout<<"da PndSttMvdTracking : N. di MC truth tracks = "< MAXMCTRACKS = "<1) {cout<<"from PndSttMvdTracking, IVOLTE = "<At(ic); if ( !( fabs(pMC->GetStartVertex().X())<0.5 && fabs(pMC->GetStartVertex().Y())<0.5 && fabs(pMC->GetStartVertex().Z())<0.5 )) continue; double carica; int icode = pMC->GetPdgCode() ; // PDG code of track double Pxx = pMC->GetMomentum().X(); double Pyy = pMC->GetMomentum().Y(); double aaa = sqrt( Pxx*Pxx + Pyy*Pyy); double Rr = aaa*1000./(BFIELD*CVEL);// R (cm) of Helix of track projected // in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track cout<<"\tTraccia n. "<GetMomentum().X() <<", Py "<GetMomentum().Y() <<", Pz "<GetMomentum().Z() <<", carica = "<GetStartVertex().X() <<", Yvert "<GetStartVertex().Y() <<", Zvert "<GetStartVertex().Z()<GetEntriesFast(); if (nSttMCPoint ==0){ cout<<"da PndSttMvdTracking : N. di Stt MC points = 0\n"<nmaxSttHits){ cout<<"da PndSttMvdTracking : N. di Stt MC points = "< nmaxSttHits ("<1){ cout<<"from PndSttMvdTracking, printout dei "<At(i); cout<<"\tpoint n. "<GetX()<< ", Ymc "<GetY()<<", Zmc "<GetZ() <<", appartiene a MC Track "<< pSttMC->GetTrackID()<GetEntriesFast(); if (nSttHit ==0){ cout<<"da PndSttMvdTracking : N. di Stt Hits = 0, return!\n"< nmaxSttHits) { cout<<"da PndSttMvdTracking : N. di Stt Hits = "< nmaxSttHits (="<= 1 && IVOLTE<20) { cout<<"da PndSttMvdTracking : evento (partendo da 0) N. "<< IVOLTE<< "\n\tN. totale Hits in STT : "<GetEntriesFast(); nMvdStripHit = fMvdStripHitArray->GetEntriesFast(); if(nMvdPixelHit>nmaxMvdPixelHits){ cout<<"from PndSttMvdTracking, nMvdPixelHit is > maximum allowed (" <nmaxMvdStripHits){ cout<<"from PndSttMvdTracking, nMvdStripHit is > maximum allowed (" <GetEntriesFast(); if(nMvdMCPoint>nmaxMvdMCPoints) { cout<<"from PndSttMvdTracking, nMvdMCPoint = "< the maximum number allowed ("<2&& IVOLTE<20) cout<<"N. MC Points delle Mvd = "<GetEntriesFast(); if (nMvdTrackCand> MAXMVDTRACKSPEREVENT) { cout<<"da PndSttMvdTracking : N. of MvdTrackCand = "<< nMvdTrackCand<<" and it is > MAXMVDTRACKSPEREVENT (="<=3 && IVOLTE<20)cout<<"da PndSttMvdTracking : n. totale di Mvd track cand = " <GetEntriesFast(); if (nSttTrackCand ==0){ cout<<"da PndSttMvdTracking : N. of SttTrackCand = 0, return!\n"< MAXTRACKSPEREVENT) { cout<<"da PndSttMvdTracking : N. of nSttTrackCand = "<< nSttTrackCand <<" and it is > MAXTRACKSPEREVENT (="<GetEntriesFast(); if(istampa>0) cout<<"da PndSttMvdTracking, event "<0){ PndSciTHit *pPndSciTHit; TVector3 posiz; // the first SciTil hit; this cannot be duplicate hit by definition. pPndSciTHit = (PndSciTHit*) fSciTHitArray->At(0); posiz = pPndSciTHit->GetPosition(); if(istampa>0) cout<<"da PndSttMvdTracking SciTil non purgati, Xpos "<At(j); posiz = pPndSciTHit->GetPosition(); if(istampa>0) cout<<"da PndSttMvdTracking SciTil non purgati, Xpos " <0){ //-----------stampe. if(istampa>0){ cout<<"da PndSttMvdTracking, dopo purga di SciTil; n. hits = "<At(i); ipunto= pSttHit->GetRefIndex();// right way to extract the corrisponding MC point. tubeID = pSttHit->GetTubeID(); pSttTube = (PndSttTube *) fSttTubeArray->At(tubeID); TVector3 center = pSttTube->GetPosition(); // drift radius Double_t dradius = pSttHit->GetIsochrone(); // wire direction TVector3 wiredirection = pSttTube->GetWireDirection(); if(wiredirection.Z() >=0.) { WDX[i] = wiredirection.X(); WDY[i] = wiredirection.Y(); WDZ[i] = wiredirection.Z(); } else { WDX[i] = -wiredirection.X(); WDY[i] = -wiredirection.Y(); WDZ[i] = -wiredirection.Z(); } info[i][0]= pSttTube->GetPosition().X(); info[i][1]= pSttTube->GetPosition().Y(); info[i][2]= pSttTube->GetPosition().Z(); info[i][3]= dradius; info[i][4]= pSttTube->GetHalfLength(); if(ipunto>=0) { puntator = (FairMCPoint*) fSttPointArray->At(ipunto); info[i][6]= puntator->GetTrackID(); } else { info[i][6]= -10.; } if( fabs( WDX[i] )< 0.00001 && fabs( WDY[i] )< 0.00001 ){ info[i][5]= 1.; ListAllParHits[nSttParHit]=i; nSttParHit++; ZCENTER_STRAIGHT = info[i][2]; // this works because just few lines below there is the SEMILENGTH_STRAIGHT = info[i][4]; // requirement that Minclinations[0] > 2 (= at least 3 parallel straws) } else { info[i][5]= 99.; // to signal that it is a skew straw. ListAllSkewHits[nSttSkewHit]=i; nSttSkewHit++; } //--------------- inizio stampaggi, stampe di controllo if (istampa >= 1) { cout <<"da PndSttMvdTracking, Stt iHit "<< i << " e n. punto MC ottenuto con RefIndex = " <GetPosition().X() << " " << pSttTube->GetPosition().Y() << " " << pSttTube->GetPosition().Z() << "; R = "<GetPosition().X()*pSttTube->GetPosition().X()+ pSttTube->GetPosition().Y()*pSttTube->GetPosition().Y())<< ", suo drift radius = "<GetX() << " " << puntator->GetY() << " " << puntator->GetZ()<GetPosition().X() << " " << pSttTube->GetPosition().Y() << " " << pSttTube->GetPosition().Z() << "; R = "<GetPosition().X()*pSttTube->GetPosition().X()+ pSttTube->GetPosition().Y()*pSttTube->GetPosition().Y())<< ", suo drift radius = "<GetTrackID()<= //-------- fine stampaggi } // end of for( i= 0; i< nSttHit; i++) // fill the exclusion list for those straws with multiple hits // first the parallel straws for(i=0; i< nSttHit-1; i++){ if( !InclusionListStt[ i ] ) continue; for(j=i+1; j< nSttHit; j++){ if(InclusionListStt[ j ] && fabs(info[i][0] - info[j][0])<1.e-20 && fabs(info[i][1] - info[j][1])<1.e-20 ) { InclusionListStt[j]= false ; } } // end of for(j=i+1; j< Nhits;; j++) } // end of for(i=0; i< Nhits-1; i++) //----------------------------------- end of exclusion of straws with multiple hits //------------------------------------------------------------- fine hits STT // --------------------------------------------- estraggo gli HITS Pixel MVD for( i= 0; i< nMvdPixelHit; i++){ pMvdPixelHit = (PndSdsHit *) fMvdPixelHitArray->At(i); // ipunto = pMvdPixelHit->GetRefIndex(); TVector3 temp = pMvdPixelHit->GetPosition(); XMvdPixel[i] = temp.X(); YMvdPixel[i] = temp.Y(); ZMvdPixel[i] = temp.Z(); sigmaXMvdPixel[i] = pMvdPixelHit->GetDx(); sigmaYMvdPixel[i] = pMvdPixelHit->GetDy(); sigmaZMvdPixel[i] = pMvdPixelHit->GetDz(); refindexMvdPixel[i] = pMvdPixelHit->GetRefIndex(); // cout<<"\tPixel n. "<GetBranchId(fMvdPixelBranch) = "<< FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)<GetBranchId(fMvdStripBranch) = "<< FairRootManager::Instance()->GetBranchId(fMvdStripBranch)<=2 && IVOLTE<20) //------------------------------------------ stampaggi // ------------------------------------------ estraggo le altre info della trackcand del MVD bool inMvdTrackCandPixel[nmaxMvdPixelHits], inMvdTrackCandStrip[nmaxMvdStripHits]; for(j=0; jAt(i); TVector3 dirSeed=pMvdTrackCand->getDirSeed(); TVector3 posSeed=pMvdTrackCand->getPosSeed(); qop = pMvdTrackCand->getQoverPseed(); nHitMvdTrackCand[i] = pMvdTrackCand->GetNHits(); // n. hits in questa track cand if( nHitMvdTrackCand[i]>nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack){ cout<<"from PndSttMvdTracking, nHitMvdTrackCand[i] = "< nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack (=" <GetSortedHit(j); // the following case should never happen (in principle), but, just to be on // the safe side .... if(pndtrackcandhit.GetHitId()<0 || pndtrackcandhit.GetDetId()<0)continue; // this is a Pixel. if( pndtrackcandhit.GetDetId()== FairRootManager::Instance()->GetBranchId(fMvdPixelBranch) && kPixel < nmaxMvdPixelHitsInTrack){ // the following is a protection because the maximum Mvd hit n. cannot exceed // nmaxMvdPixelHits . if( pndtrackcandhit.GetHitId()>nmaxMvdPixelHits ){ cout<<"from PndSttMvdTracking, this Pixel Mvd hit has a number = " < nmaxMvdPixelHits (="<GetBranchId(fMvdStripBranch) && kStrip < nmaxMvdStripHitsInTrack){ // the following is a protection because the maximum Mvd hit n. cannot exceed // nmaxMvdStripHits . if( pndtrackcandhit.GetHitId()>nmaxMvdStripHits ){ cout<<"from PndSttMvdTracking, this Strip Mvd hit has a number = " < nmaxMvdStripHits (="<=0.){ ListMvdDSPixelHitNotTrackCand[nMvdDSPixelHitNotTrackCand] = i; nMvdDSPixelHitNotTrackCand++; } else { ListMvdUSPixelHitNotTrackCand[nMvdUSPixelHitNotTrackCand] = i; nMvdUSPixelHitNotTrackCand++; } } } for( i= 0; i< nMvdStripHit ; i++){ if( ! inMvdTrackCandStrip[i] ){ if( ZMvdStrip[i]>=0.){ ListMvdDSStripHitNotTrackCand[nMvdDSStripHitNotTrackCand] = i; nMvdDSStripHitNotTrackCand++; } else { ListMvdUSStripHitNotTrackCand[nMvdUSStripHitNotTrackCand] = i; nMvdUSStripHitNotTrackCand++; } } } if(istampa>=2){ for( i= 0; i< nMvdTrackCand ; i++){ cout<<" PndSttMvdTracking, MvdTrackCand n. "<0) { cout<<"PndSttMvdTracking : evt. "< nmaxSttHits (="<GetSortedHit(j); // necessary check to exclude hits whose number is > nmaxSttHits. if(pndtrackcandhit.GetDetId()==1001){// my SciTil detector ID. ListSciTilHitsinTrack[i][nSciTilHitsinTrack[i]]=pndtrackcandhit.GetHitId(); InclusionListSciTil[ ListSciTilHitsinTrack[i][nSciTilHitsinTrack[i]] ] = false; nSciTilHitsinTrack[i] ++; } else { // this is a Stt Det Id if( pndtrackcandhit.GetHitId() >= nmaxSttHits ) continue; igoodStt++; if(igoodStt==0) { // for the calculation of the approximate Fi of the first hit in those track x = info[ pndtrackcandhit.GetHitId() ][0];//in the middle of the tube y = info[ pndtrackcandhit.GetHitId() ][1];//in the middle of the tube } ListSttHitsinTrack[i][nSttParHitsinTrack[i]+nSttSkewHitsinTrack[i]] = pndtrackcandhit.GetHitId(); // # hit of Stt if ( fabs(info[ pndtrackcandhit.GetHitId() ][5]- 1.) < 0.0001) { // axial Stt hits ListSttHitsinTrackType[i][nSttParHitsinTrack[i]+nSttSkewHitsinTrack[i]] = 2; ListSttParHitsinTrack[i][nSttParHitsinTrack[i]] = pndtrackcandhit.GetHitId(); // # hit of Stt nSttParHitsinTrack[i]++; } else { // skew Stt hits ListSttHitsinTrackType[i][nSttParHitsinTrack[i]+nSttSkewHitsinTrack[i]] = 3; ListSttSkewHitsinTrack[i][nSttSkewHitsinTrack[i]] = pndtrackcandhit.GetHitId(); // # hit of Stt nSttSkewHitsinTrack[i]++; } } // end of if(pndtrackcandhit.GetDetId()==1001) } // end of for(j=0; j=2){ cout<<"da PndSttMvdTracking, evt. "<getDirSeed(); posSeed=pSttTrackCand->getPosSeed(); qop = pSttTrackCand->getQoverPseed(); daTrackFoundaTrackMC[i]=pSttTrackCand->getMcTrackId(); Px[i] = dirSeed.X(); Py[i] = dirSeed.Y(); Pz[i] = dirSeed.Z(); if( dirSeed.Z() < 9990.) { // normal case with skew hits information existing SttSZfit[i]=true; Ntras = sqrt(dirSeed.X()*dirSeed.X()+dirSeed.Y()*dirSeed.Y() ) ; R[i] = Ntras /(0.003*BFIELD*fabs(qop)); if(qop <0.){ Ox[i] = -R[i]*dirSeed.Y()/Ntras ; Oy[i] = R[i]*dirSeed.X()/Ntras ; KAPPA[i]=0.003*BFIELD*fabs(qop)/dirSeed.Z(); CHARGE[i]=-1; } else { Ox[i] = R[i]*dirSeed.Y()/Ntras ; Oy[i] = -R[i]*dirSeed.X()/Ntras ; KAPPA[i]=-0.003*BFIELD*fabs(qop)/dirSeed.Z(); CHARGE[i]=1; } } else if (dirSeed.Z()-10000.<0.) { // case without skew straws info present SttSZfit[i]=false; Ptras = sqrt(dirSeed.X()*dirSeed.X()+dirSeed.Y()*dirSeed.Y() ) ; R[i] = Ptras /(0.003*BFIELD); if(qop <0.){ Ox[i] = -R[i]*dirSeed.Y()/Ptras; Oy[i] = R[i]*dirSeed.X()/Ptras; CHARGE[i]=-1; } else { Ox[i] = R[i]*dirSeed.Y()/Ptras; Oy[i] = -R[i]*dirSeed.X()/Ptras; CHARGE[i]=1; } } else { // case in which KAPPA = 0 SttSZfit[i]=true; if(qop <0.){ Ox[i] = -R[i]*dirSeed.Y()/Ntras; Oy[i] = R[i]*dirSeed.X()/Ntras; CHARGE[i]=-1; } else { Ox[i] = R[i]*dirSeed.Y()/Ptras; Oy[i] = -R[i]*dirSeed.X()/Ptras; CHARGE[i]=1; } KAPPA[i]=0.; } FI0[i] = atan2(-Oy[i], -Ox[i]); Ox[i] += Trajectory_Start[i][0]; Oy[i] += Trajectory_Start[i][1]; Fifirst[i] = atan2( y-Oy[i], x-Ox[i]); HoughFi[i] = FI0[i] + PI; if(HoughFi[i]<0.) HoughFi[i]=0.; if( FI0[i] < 0. ) FI0[i]+= 2.*PI; if( FI0[i] < 0. ) FI0[i]=0.; if( Fifirst[i] < 0. ) Fifirst[i]+= 2.*PI;if( Fifirst[i] < 0. ) Fifirst[i]=0.; } // end of for( i= 0; i< nSttTrackCand; i++) //--------------- end fetching of the Stt PndTrackCand from PR of the STT //----------- treat here the special pathological case when there are no Mvd hits at all. //----------- Copy simply the quantities from the STT alone PndTrackCand. if(nMvdPixelHit+nMvdStripHit==0){ for( ncand= 0; ncand< nSttTrackCand; ncand++){ int nbuoni, nparbuoni, nskewbuoni; for(j=0, nbuoni=0,nparbuoni=0, nskewbuoni=0;j fabs(ddd) ){ dista0=ddd; zeta0 = z[0]-zdrift[0]; } else { zeta0 = z[0]+zdrift[0]; } dista1 = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]+zdrift[1],s[1],&nrounds1); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]-zdrift[1],s[1],&nrounds1); if( fabs(dista1)> fabs(ddd) ){ dista1=ddd; zeta1 = z[1]-zdrift[1]; } else { zeta1 = z[1]+zdrift[1]; } if( fabs(dista1) fabs(ddd) ){ ZchosenSkew[ncand][k]=z[0]-zdrift[0]; } else { ZchosenSkew[ncand][k]=z[0]+zdrift[0]; } ListTrackCandHit[ncand][nbuoni]= ListSttSkewHitsinTrack[ncand][nskewbuoni]= ListSttHitsinTrack[ncand][j]; ListTrackCandHitType[ncand][nbuoni]= ListSttHitsinTrackType[ncand][j]; nskewbuoni++; nbuoni++; } else if( z[1]<999998.){ dista = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand] ,z[1]+zdrift[1],s[1],&nrounds0); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]-zdrift[1],s[1],&nrounds1); SchosenSkew[ncand][k]=s[1]; if( fabs(dista)> fabs(ddd) ){ ZchosenSkew[ncand][k]=z[1]-zdrift[1]; } else { ZchosenSkew[ncand][k]=z[1]+zdrift[1]; } ListTrackCandHit[ncand][nbuoni]= ListSttSkewHitsinTrack[ncand][nskewbuoni]= ListSttHitsinTrack[ncand][j]; ListTrackCandHitType[ncand][nbuoni]= ListSttHitsinTrackType[ncand][j]; nskewbuoni++; nbuoni++; } // end of if( z[0]<999998. && z[1]<999998.) } else { // continuation of if(ListSttHitsinTrackType[ncand][j]==3) // so, here it must be a parallel hit. ListTrackCandHit[ncand][nbuoni]= ListSttParHitsinTrack[ncand][nparbuoni]= ListSttHitsinTrack[ncand][j]; ListTrackCandHitType[ncand][nbuoni]= ListSttHitsinTrackType[ncand][j]; nparbuoni++; nbuoni++; } // end of if(ListSttHitsinTrackType[ncand][j]==3) } // end of for(j=0; j Fi_low_limit and // possibly > 2PI. &statusflag[ncand],//it is a vector; =0, all well; =1, track contained completely between RMin // and RMax; = -1 track contained within RMin; =-2 track outside RMax. RStrawDetectorMin, RStrawDetectorMax ); if( statusflag[ncand] == -1) { Fi_low_limit[ncand] = -99999.; } else if (statusflag[ncand] == -2){ keepit[ncand] = false; }; } // end of for( ncand= 0; ncand< nSttTrackCand; ncand++) //--------------------- here call to the function that matches Mvd hits with Stt tracks delta=0.5; // parameter of proximity for associating Mvd hits to Stt tracks highqualitycut=0.2; // parameter of proximity for associating Mvd hits to Stt tracks // This method matches the Mvd hits to the found tracks. MatchMvdHitsToSttTracks2( keepit,// input and output. delta, highqualitycut, nSttTrackCand, Ox, Oy, R, FI0, Fi_low_limit, CHARGE, nMvdPixelHitsinTrack, // output ListMvdPixelHitsinTrack, // output nMvdStripHitsinTrack, // output ListMvdStripHitsinTrack // output ); //---------inizio stampe. if(istampa>=2&& IVOLTE<20){ cout<<"da PndSttMvdTracking, dopo MatchMvd..2 :\n"; stampetta(nSttTrackCand,keepit); } //end of if(istampa>=0) //---------------------------------- fine stampe. //----------------- end of section with match Mvd hits with Stt hits for( ncand= 0; ncand< nSttTrackCand; ncand++){ if(nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]==0){ Mvdhits[ncand]=false; ALFA[ncand]=-2.*Ox[ncand]; BETA[ncand]=-2.*Oy[ncand]; GAMMA[ncand]= Ox[ncand]*Ox[ncand]+Oy[ncand]*Oy[ncand]-R[ncand]*R[ncand]; // just copy and go to the next track cand. int nbuoni, nparbuoni, nskewbuoni; for(j=0, nbuoni=0,nparbuoni=0, nskewbuoni=0;j=2) cout<<"da PndSttMvdTracking, IVOLTE = "<< IVOLTE<<", ncand = " < fabs(ddd) ){ dista0=ddd; zeta0 = z[0]-zdrift[0]; } else { zeta0 = z[0]+zdrift[0]; } dista1 = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]+zdrift[1],s[1],&nrounds1); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]-zdrift[1],s[1],&nrounds1); if( fabs(dista1)> fabs(ddd) ){ dista1=ddd; zeta1 = z[1]-zdrift[1]; } else { zeta1 = z[1]+zdrift[1]; } if( fabs(dista1) fabs(ddd) ){ ZchosenSkew[ncand][k]=z[0]-zdrift[0]; } else { ZchosenSkew[ncand][k]=z[0]+zdrift[0]; } ListTrackCandHit[ncand][nbuoni]= ListSttSkewHitsinTrack[ncand][nskewbuoni]= ListSttHitsinTrack[ncand][j]; ListTrackCandHitType[ncand][nbuoni]= ListSttHitsinTrackType[ncand][j]; nskewbuoni++; nbuoni++; } else if( z[1]<999998.){ dista = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand] ,z[1]+zdrift[1],s[1],&nrounds0); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]-zdrift[1],s[1],&nrounds1); SchosenSkew[ncand][k]=s[1]; if( fabs(dista)> fabs(ddd) ){ ZchosenSkew[ncand][k]=z[1]-zdrift[1]; } else { ZchosenSkew[ncand][k]=z[1]+zdrift[1]; } ListTrackCandHit[ncand][nbuoni]= ListSttSkewHitsinTrack[ncand][nskewbuoni]= ListSttHitsinTrack[ncand][j]; ListTrackCandHitType[ncand][nbuoni]= ListSttHitsinTrackType[ncand][j]; nskewbuoni++; nbuoni++; } // end of if( z[0]<999998. && z[1]<999998.) } else { // continuation of if(ListSttHitsinTrackType[ncand][j]==3) // so, here it must be a parallel hit. ListTrackCandHit[ncand][nbuoni]= ListSttParHitsinTrack[ncand][nparbuoni]= ListSttHitsinTrack[ncand][j]; ListTrackCandHitType[ncand][nbuoni]= ListSttHitsinTrackType[ncand][j]; nparbuoni++; nbuoni++; } // end of if(ListSttHitsinTrackType[ncand][j]==3) } // end of for(j=0; j nmaxMvdPixelHitsInTrack){ cout<<"from PndSttMvdTracking, nMvdPixelHitsinTrack["< nmaxMvdPixelHitsInTrack (" < nmaxMvdStripHitsInTrack){ cout<<"from PndSttMvdTracking, nMvdStripHitsinTrack["< nmaxMvdStripHitsInTrack (" <=2&& IVOLTE<20){ cout<<"da PndSttMvdTracking, dopo riordinamento :\n"; stampetta(nSttTrackCand,keepit); } //end of if(istampa>=0) //---------------------------------- fine stampe. //--------------------- refit the Helix in XY plane using Stt + Mvd associated hits Short_t iexcl; Double_t d, diff, rotationangle, trajectory_vertex[2]; trajectory_vertex[0]=0.; trajectory_vertex[1]=0.; for(ncand=0; ncand< nTotalCandidates; ncand++){ if(!keepit[ncand]) continue; if( !Mvdhits[ncand]) continue; rotationangle= atan2( Py[ncand], Px[ncand]); // Py/Px = m in // v = m*u + q formula if( rotationangle<0. ) rotationangle += 2.*PI; //---------- translation of the reference system in the best Mvd hit position for(i=0, diff = 1.e20;i Pixel; = 1 --> Strip; = 2 Stt parallel; 3 = Stt Skew. // There are no -1 type hits at this point. if(!(ListTrackCandHitType[ncand][i]==1 || ListTrackCandHitType[ncand][i]==0) ) continue; if( ListTrackCandHitType[ncand][i] == 0 ){ d = fabs( sqrt( (Ox[ncand]-XMvdPixel[ListTrackCandHit[ncand][i]]) *(Ox[ncand]-XMvdPixel[ListTrackCandHit[ncand][i]]) +(Oy[ncand]-YMvdPixel[ListTrackCandHit[ncand][i]]) *(Oy[ncand]-YMvdPixel[ListTrackCandHit[ncand][i]]) ) - R[ncand]); if( d < diff ){ diff = d; iexcl=i; trajectory_vertex[0]= XMvdPixel[ListTrackCandHit[ncand][i]]; trajectory_vertex[1]= YMvdPixel[ListTrackCandHit[ncand][i]]; } } else { d = fabs( sqrt( (Ox[ncand]-XMvdStrip[ListTrackCandHit[ncand][i]]) *(Ox[ncand]-XMvdStrip[ListTrackCandHit[ncand][i]]) +(Oy[ncand]-YMvdStrip[ListTrackCandHit[ncand][i]]) *(Oy[ncand]-YMvdStrip[ListTrackCandHit[ncand][i]]) ) - R[ncand]); if( d < diff ){ diff = d; iexcl=i; trajectory_vertex[0]= XMvdStrip[ListTrackCandHit[ncand][i]]; trajectory_vertex[1]= YMvdStrip[ListTrackCandHit[ncand][i]]; } } } // end of for(i=0, diff = 1.e20;i<..... if( diff > 0.5 ) { trajectory_vertex[0]=trajectory_vertex[1]=0.; iexcl=-1; } RefitMvdStt( nTrackCandHit[ncand], //this is input &ListTrackCandHit[ncand][0],//this is both input and output &ListTrackCandHitType[ncand][0],//this is both input and output info, rotationangle, trajectory_vertex, iexcl, &ALFA[ncand], // output of the fit &BETA[ncand], // output of the fit &GAMMA[ncand],// output of the fit &status[ncand] // fit status; true = successful ); if(status[ncand]){ Ox[ncand] = -ALFA[ncand]/2.; Oy[ncand] = -BETA[ncand]/2.; R[ncand] = Ox[ncand]*Ox[ncand]+Oy[ncand]*Oy[ncand]-GAMMA[ncand]; if( R[ncand] > 0. ) R[ncand]=sqrt(R[ncand]) ; else R[ncand]=0.; FI0[ncand] = atan2(Trajectory_Start[ncand][1]-Oy[ncand], Trajectory_Start[ncand][0]-Ox[ncand]); if( FI0[ncand] < 0. ) FI0[ncand]+= 2.*PI; } if(R[ncand]=2&& IVOLTE<20){ cout<<"da PndSttMvdTracking, dopo XY Refit :\n"; stampetta(nSttTrackCand,keepit); } //end of if(istampa>=0) //---------------------------------- fine stampe. // find out if the SciTil hit associated with this track is still acceptable after the last // XY refit. for(i=0, iaccept=0;i Fi_low_limit and // possibly > 2PI. // Fi_low_limit and Fi_up_limit are undefined // when statusflag is negative. &statusflag[ncand],//it is a vector; =0, all well; =1, track contained completely between RMin // and RMax; = -1 track contained within RMin; =-2 track outside RMax. RStrawDetectorMin, RStrawDetectorMax ); // Fi_low_limit set to -99999. when it is contained in Mvd // region completely. if( statusflag[ncand] == -1 ) Fi_low_limit[ncand] = -99999.; if( statusflag[ncand] == -2) keepit[ncand] = false; } // end of for(ncand=0; ncand< nTotalCandidates; ncand++) //--------------------- end of refit the Helix in XY plane using Stt + Mvd associated hits if(nMvdPixelHit+nMvdStripHit>0){ //--------------------- here call to the function that matches Mvd hits with the new // circular trajectory in XY found for the second time, after first refit delta=0.5; // parameter of proximity for associating Mvd hits to Stt tracks // highqualitycut=0.3; // parameter of proximity for associating Mvd hits to Stt tracks highqualitycut=0.5; // parameter of proximity for associating Mvd hits to Stt tracks MatchMvdHitsToSttTracksagain( keepit, Mvdhits, delta, highqualitycut, nSttTrackCand, Ox, Oy, R, FI0, Fi_low_limit, // because here we deal with hits in Mvd region CHARGE, nMvdPixelHitsinTrack, // input and output ListMvdPixelHitsinTrack, // input and output nMvdStripHitsinTrack, // input and output ListMvdStripHitsinTrack // input and output ); if(istampa>=2){ cout<<"\n\nda PndSttMvdTracking, dopo MatchMvdHitsToSttTracksagain :"<0) //--------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // use the risult just obtained from the fit in XY to redo the association of the Skew Straw hits for(ncand=0; ncand< nTotalCandidates; ncand++) { if(!keepit[ncand]) continue; if( ! Mvdhits[ncand]) { if(YesCleanMvd){ // reject the candidate if it is NOT contained in the pipe and // therefore it should have at least 1 Mvd hit but it has none. if( (!IsInTargetPipe( Ox[ncand], Oy[ncand], R[ncand], FI0[ncand], KAPPA[ncand], CHARGE[ncand], VERTICALGAP/2.) ) ) keepit[ncand]=false; } // end of (YesCleanMvd) continue; } if( statusflag[ncand] == -1 ) continue; // this is when the XY circle is contained in the // the Mvd region completely; skip the association of // the Skews. nSttSkewHitsinTrack[ncand]= AssociateSkewHitsToXYTrack( InclusionListStt, // hit is excluded only if it multiple hit. nSttSkewHit, ListAllSkewHits, Ox[ncand], // input : X of center of XY plane circle Oy[ncand], // input : Y of center of XY plane circle R[ncand], // input : Radius of XY plane circle info, WDX, WDY, WDZ, Fi_low_limit[ncand], // Fi (in Helix XY frame) lower limit using the Stt detector minimum/maximum radius Fi_up_limit[ncand], // Fi (in Helix XY frame) upper limit using the Stt detector maximum/minimum radius CHARGE[ncand], TemporarySkewList, // output, list of selected skew hits (in skew numbering) TemporaryS, // output, S coordinate of selected Skew hit TemporaryZ, // output, Z coordinate of selected Skew hit (center wire) TemporaryZDrift, // output, drift distance IN Z DIRECTION only, of selected Skew hit TemporaryZErrorafterTilt // output, Radius taking into account the tilt, IN Z DIRECTION only, of selected Skew hit ); for(j=0;j=2){ cout<<"\n\nda PndSttMvdTracking, dopo AssociateSkewHitsToXYTrack :"<1.e-10) { ErrorDriftRadiusbis[kall]=TemporaryZErrorafterTilt[j]; } else { ErrorDriftRadiusbis[kall]=0.5; } } // end of for(j=0;j=2){ cout<<"\n\nda PndSttMvdTracking, prima di FitSZspace :"<=2){ cout<<"\n\nda PndSttMvdTracking, dopo FitSZspace :"<=2){ cout<<"\n\nda PndSttMvdTracking, dopo isintargetpipe (Ox="<0){ //---------inizio stampe. if(istampa>=2&& IVOLTE<20){ cout<<"da PndSttMvdTracking, prima di CollectParSttHitsagain :\n"; stampetta(nSttTrackCand,keepit); } //end of if(istampa>=0) //---------------------------------- fine stampe. CollectParSttHitsagain( keepit, Mvdhits, info, nSttParHit, ListAllParHits, 0, // starting candidate n. (included). nTotalCandidates, // ending candidate n. (excluded). Ox, Oy, R, KAPPA, FI0, Fi_low_limit, Fi_up_limit, nSttParHitsinTrack, // input and output ListSttParHitsinTrack // input and output ); //---------inizio stampe. if(istampa>=2&& IVOLTE<20){ cout<<"da PndSttMvdTracking, dopo CollectParSttHitsagain :\n"; stampetta(nSttTrackCand,keepit); } //end of if(istampa>=0) //---------------------------------- fine stampe. // ordering all the hits belonging to the candidate track, by increasing R (large // trajectories) or Conformal variables (better for small trajectories); // from candidate n. 0 to candidate n. nTotalCandidates-1; loading ListTrackCandHit. // the array ordered are : // ListTrackCandHit, ListTrackCandHitType, ListSttParHitsinTrack, ListSttSkewHitsinTrack // and also at the end the SciTil hit (if present) is added. Ordering_Loading_ListTrackCandHit( keepit, 0, nTotalCandidates, info, Ox, Oy, R, Trajectory_Start, CHARGE, SchosenSkew ); // adding at the end the SciTil hits (if present). for(ncand=0; ncand< nTotalCandidates; ncand++){ i=nMvdPixelHitsinTrack[ncand]+nMvdStripHitsinTrack[ncand]+ nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand]; for(j=0;j=2&&IVOLTE<20){ cout<<"\n------------------------------------------------\n" <<"\tIVOLTE = "<1) cout<<"PndSttMvdTracking, entra in TrackCleanup tracce normali, IVOLTE "<0) { // load the UsedPixel and UsedStrip vectors. bool UsedPixel[nmaxMvdPixelHits], UsedStrip[nmaxMvdStripHits]; Short_t List[nmaxMvdPixelHits+nmaxMvdStripHits], ListType[nmaxMvdPixelHits+nmaxMvdStripHits]; Double_t centrex, centrey, radius; for(i=0;iGetBranchId(fMvdPixelBranch)){ if(UsedPixel[ ListHitMvdTrackCand[i][j] ]) continue; AloneX[nalone] = XMvdPixel[ ListHitMvdTrackCand[i][j] ]; AloneY[nalone] = YMvdPixel[ ListHitMvdTrackCand[i][j] ]; List[nalone] = ListHitMvdTrackCand[i][j]; ListType[nalone] = 0; nalone++; }else {// at this point this is a Strip hit; already made sure // earlier in the code that there is no third possibility. // if (ListHitTypeMvdTrackCand[i][j] == // this condition in principle at this // // point is always satisfied. // FairRootManager::Instance()->GetBranchId(fMvdStripBranch)){ if(UsedStrip[ ListHitMvdTrackCand[i][j] ]) continue; AloneX[nalone] = XMvdStrip[ ListHitMvdTrackCand[i][j] ]; AloneY[nalone] = YMvdStrip[ ListHitMvdTrackCand[i][j] ]; List[nalone] = ListHitMvdTrackCand[i][j]; ListType[nalone] = 1; nalone++; } } // end of for(j=0;j RStrawDetectorMin/2.){ Ox[nTotalCandidates] = centrex; Oy[nTotalCandidates] = centrey; R[nTotalCandidates] = radius; FI0[nTotalCandidates] = atan2(-Oy[nTotalCandidates], -Ox[nTotalCandidates]); if( FI0[nTotalCandidates] < 0. ) FI0[nTotalCandidates]+= 2.*PI; nSttParHitsinTrack[nTotalCandidates]=0; nSttSkewHitsinTrack[nTotalCandidates]=0; nMvdPixelHitsinTrack[nTotalCandidates]=0; nMvdStripHitsinTrack[nTotalCandidates]=0; for(j=0; j2){ cout<<"from PndSttMvdTracking, evt. "<1){ cout<<"PndSttMvdTracking, entra in TrackCleanup tracce Mvd, IVOLTE "<1){ cout<<"\tPndSttMvdTracking, uscito da TrackCleanup tracce Mvd, keepit = "; if(keepit[ncand]) cout<<" true!\n"; else cout<<" false!\n"; } } // end of if(YesClean) //----------- end of cleanup of new tracks found. } // end of if(MvdAloneTracking) // ----------------------------- end of finding Tracks starting with the Mvd tracks // ------------------------------------------------------------------------------------- // ------------------------------------------------------------------------------------- // ------------------------------------------------------------------------------------- // skipping: ; //----------stampaggi if(istampa>=2){ cout<<"Evt. n. "<0){ cout<<"\tscitil hit n. "<0 && nTotalCandidates > 0 && doMcComparison ){ int nmid,nn; Double_t X1[MAXTRACKSPEREVENT], Y1[MAXTRACKSPEREVENT], X2[MAXTRACKSPEREVENT], Y2[MAXTRACKSPEREVENT], X3[MAXTRACKSPEREVENT], Y3[MAXTRACKSPEREVENT]; for(i=0; i0 ){ // clockwise. if( FI0[i] < angle ) angle -= 2.*PI; if( FI0[i] < angle ) angle =0.; } else{ // counterclockwise. if( FI0[i] > angle ) angle += 2.*PI; if( FI0[i] > angle ) angle = FI0[i]; } middle = (FI0[i]+angle)/2.; X2[i] = Ox[i] + R[i]*cos( middle ); Y2[i] = Oy[i] + R[i]*sin( middle ); } // end of for(i=0; i=2){ // for(i=0;i=3){ // for(ncand=0; ncand< nSttTrackCand; ncand++){ int nbuone=-1; for(ncand=0; ncand< nTotalCandidates; ncand++){ if(!keepit[ncand]) continue; nbuone++; //----------------- ora la traccia MC corrispondente a questa traccia Stt // cout<<"\n\nda PndSttMvdTracking : la PndTrackCand n. "<-1 ) { Int_t icode, im; Double_t aaa,Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Pxx, Pyy, carica ; PndMCTrack* pMC; im=daTrackFoundaTrackMC[ncand]; pMC = (PndMCTrack*) fMCTrackArray->At(im); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Pxx = pMC->GetMomentum().X(); Pyy = pMC->GetMomentum().Y(); aaa = sqrt( Pxx*Pxx + Pyy*Pyy); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if(fabs(carica)<1.e-5) { cout<<"...associata ad una traccia neutra, assurdo!\n"; continue ;} Cx = Oxx + Pyy*1000./(BFIELD*CVEL*carica); Cy = Oyy - Pxx*1000./(BFIELD*CVEL*carica); cout<<" e' associata alla traccia MC n. "<=2) //-------------fine stampaggi // ora il confronto per il meeting di Groningen //---------- conteggio delle tracce MC accettabili!! int citata; int nMCTracksaccettabili=0; int ListaMCTracksaccettabili[nMCTracks]; for (i=0;iAt(i); if ( ! pMCtr ) continue; Double_t aaa, carica, Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Pxx, Pyy ; Int_t icode; icode = pMCtr->GetPdgCode() ; // PDG code of track Oxx = pMCtr->GetStartVertex().X(); // X of starting point track Oyy = pMCtr->GetStartVertex().Y(); // Y of starting point track Pxx = pMCtr->GetMomentum().X(); Pyy = pMCtr->GetMomentum().Y(); aaa = sqrt( Pxx*Pxx + Pyy*Pyy); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla if(istampa>2){ TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if(fabs(carica)<1.e-5) continue; Cx = Oxx + Pyy*1000./(BFIELD*CVEL*carica); Cy = Oyy - Pxx*1000./(BFIELD*CVEL*carica); cout<<"da PndSttMvdTracking, evento (cominciando da 0) n. "<=1 ) fprintf(HANDLE, "\n Evento %d NTotaleTracceMC %d ------\n",IVOLTE, nMCTracks); if(istampa>=1 ) { fprintf(HANDLE, "\n Evento %d NTotaleTracceMC %d ------\n",IVOLTE, nMCTracksaccettabili); int ibene=0; if(nMCTracksaccettabili>0){ for(int ii=0; ii0) fprintf(HANDLE,"\tn. volte almeno 1 traccia MC accettabile e' ricostruita %d\n" ,ibene); bool flaggo; int ii, ibuone=-1; Double_t HoughFiii; for (ii=0; ii1) {cout<<"\tevt. n "<At(i); if ( ! pMCtr ){ fprintf(HANDLE, " MC track n. %d doesn't have pointer to MC Track TClones Array\n", i); continue; } // controllo che la traccia associata MC sia una delle tracce MC 'ragionevoli'. flaggo = true; for(int g=0; gGetPdgCode() ; // PDG code of track Oxx = pMCtr->GetStartVertex().X(); // X of starting point track Oyy = pMCtr->GetStartVertex().Y(); // Y of starting point track Pxx = pMCtr->GetMomentum().X(); Pyy = pMCtr->GetMomentum().Y(); aaa = sqrt( Pxx*Pxx + Pyy*Pyy); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if(fabs(carica)<1.e-5) fprintf(HANDLE," MC track n. %d e' neutra, assurdo!\n",i); Cx = Oxx + Pyy*1000./(BFIELD*CVEL*carica); Cy = Oyy - Pxx*1000./(BFIELD*CVEL*carica); Fifi = atan2(Cy, Cx); // MC truth Fifi angle of circle of Helix trajectory if(Fifi<0.) Fifi += 2.*PI; Double_t Kakka ; if( fabs( pMCtr->GetMomentum().Z() )< 1.e-20) Kakka = 99999999.; else Kakka = -carica*0.001*BFIELD*CVEL/pMCtr->GetMomentum().Z(); fprintf(HANDLE, " R_MC %g R %g Fi_MC %g Fi %g KAPPA_MC %g KAPPA %g FI0_MC %g FI0 %g\n", Rr, R[ ii ], Fifi, HoughFiii, Kakka, KAPPA[ ii ], fmod(Fifi+ PI, 2.*PI), // FI0 da MC truth FI0[ ii ] ); //------------------------ // HoughFi = atan2(Oy[i],Ox[i]); // if(HoughFi<0.) HoughFi += 2.*PI; fprintf(HANDLE2,"Evento n. %d Found track %d messa in PndTrackCand",IVOLTE, ii); fprintf(HANDLE2, " R_MC %g R %g Fi_MC %g Fi %g KAPPA_MC %g KAPPA %g FI0_MC %g FI0 %g\n", Rr, R[ ii ], Fifi, HoughFiii, Kakka, KAPPA[ ii ], fmod(Fifi+ PI, 2.*PI), FI0[ ii ] ); //---------------------------------- } // end of for (ii=0; ii=1) //--------------ghosts // fa il conto delle ghost solo sugli eventi che hanno almeno 1 traccia MC accettabile. int NParghost=0, NParhitsghost=0,icc,ibuone=-1; if( istampa>=1 && nMCTracksaccettabili>0){ for(icc=0; icc=1) } // end of if( nMCTracks >0 && nTotalCandidates > 0 && doMcComparison ) //------------------ end section with comparison MC Mvd hits - associated hits to a certain track //------- load the new PndTrackCand ; each track has the STT and the Mvd hits associated //------- also load the new PndTrack ; each track has the STT and the Mvd hits associated Double_t Oxx, Oyy; // Int_t iflaggo; for(ncand=0, ipinco = 0; ncand< nTotalCandidates; ncand++){ if(!keepit[ncand]) continue; // case in which there was no Skew hits and no KAPPA info and that // candidate could not be associated to any Mvd hits --> no KAPPA information! if(ncand1.e-20 ){ Pzini = -CHARGE[ncand]*0.003*BFIELD/KAPPA[ncand]; if(fabs(Pzini) > PMAX) continue; } else { continue; } // PndTrackCand Array loading new((*fSttMvdPndTrackCandArray)[ipinco]) PndTrackCand; PndTrackCand *pTrckCand = (PndTrackCand*) fSttMvdPndTrackCandArray->At(ipinco); TVector3 dirSeed(Pxini,Pyini,Pzini); // momentum direction in starting point qop = CHARGE[ncand]/dirSeed.Mag(); dirSeed.SetMag(1.); pTrckCand->setTrackSeed(posSeed, dirSeed, qop); if(doMcComparison ){ pTrckCand->setMcTrackId( daTrackFoundaTrackMC[ncand] ); }else{ pTrckCand->setMcTrackId(-1); } // pTrckCand->sorted=true; for(j=0; j< nTrackCandHit[ncand]; j++){ switch (ListTrackCandHitType[ncand][j]){ case 0: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId(fMvdPixelBranch), (Int_t)ListTrackCandHit[ncand][j],j); break; case 1: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId(fMvdStripBranch), (Int_t)ListTrackCandHit[ncand][j],j); break; case 2: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId(fSttBranch), (Int_t)ListTrackCandHit[ncand][j],j); break; case 3: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId(fSttBranch), (Int_t)ListTrackCandHit[ncand][j],j); break; } } // PndTrack Array loading // last = last hit in found track; first = first hit in found track; // pTrckCand = pointer to the corresponding PndTrackCand. // the first hit if (ListTrackCandHitType[ncand][0] == 0){ // Mvd Pixel Posiz1[0] = XMvdPixel[ ListTrackCandHit[ncand][0] ]; Posiz1[1] = YMvdPixel[ ListTrackCandHit[ncand][0] ]; Posiz1[2] = ZMvdPixel[ ListTrackCandHit[ncand][0] ]; ErrPosition.SetX(sigmaXMvdPixel[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); ErrPosition.SetY(sigmaXMvdPixel[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); ErrPosition.SetZ(sigmaXMvdPixel[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); } else if (ListTrackCandHitType[ncand][0] == 1){ // Mvd Strip Posiz1[0] = XMvdStrip[ ListTrackCandHit[ncand][0] ]; Posiz1[1] = YMvdStrip[ ListTrackCandHit[ncand][0] ]; Posiz1[2] = ZMvdStrip[ ListTrackCandHit[ncand][0] ]; ErrPosition.SetX(sigmaXMvdStrip[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); ErrPosition.SetY(sigmaXMvdStrip[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); ErrPosition.SetZ(sigmaXMvdStrip[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); } else if( ListTrackCandHitType[ncand][0] == 2 ){ // it is a parallel straw hit PndSttInfoXYZParal( info, ListTrackCandHit[ncand][0], Ox[ncand], Oy[ncand], R[ncand], KAPPA[ncand], FI0[ncand], CHARGE[ncand], Posiz1 ); // cases in which the calculation of the position failed, see // PndSttInfoXYZParal method. if(Posiz1[2]<-888888887. || Posiz1[0] < -999999998.) continue; ErrPosition.SetX(0.02); // 200 microns ErrPosition.SetY(0.02); // 200 microns ErrPosition.SetZ(1.); // 1 cm } else if ( ListTrackCandHitType[ncand][0] == 3 ){ // it is a skew straw hit Posiz1[0] = Ox[ncand]+R[ncand]*cos(SchosenSkew[ncand][ ListTrackCandHit[ncand][0] ]); Posiz1[1] = Oy[ncand]+R[ncand]*sin(SchosenSkew[ncand][ ListTrackCandHit[ncand][0] ]); Posiz1[2] = ZchosenSkew[ncand][ ListTrackCandHit[ncand][0] ]; ErrPosition.SetX(0.02); // 200 microns ErrPosition.SetY(0.02); // 200 microns ErrPosition.SetZ(1.); // 1 cm } Position.SetX( Posiz1[0] ); Position.SetY( Posiz1[1] ); Position.SetZ( Posiz1[2] ); versor[0] = Ox[ncand]-Posiz1[0]; versor[1] = Oy[ncand]-Posiz1[1]; Distance = sqrt(versor[0]*versor[0]+versor[1]*versor[1]); versor[0] /= Distance; versor[1] /= Distance; px = -CHARGE[ncand]*Ptras*versor[1]; py = CHARGE[ncand]*Ptras*versor[0]; Momentum.SetX(px); Momentum.SetY(py); Momentum.SetZ(Pzini); ErrMomentum.SetX(0.05*Ptras); // set at 5% all the times. ErrMomentum.SetY(0.05*Ptras); // set at 5% all the times. ErrMomentum.SetZ(0.05*Pzini); // set at 5% all the times. // the plane of this FairTrackParP better is perpendicular to // the momentum direction ddd = Ptras*sqrt(Ptras*Ptras+Pzini*Pzini); if(istampa>=3){ cout<<" evento = "<=2) //------------------------------ plottamenti -------------------------------------------- if(iplotta && IVOLTE=2&&IVOLTE<20){ cout<<"\n\n---------------------------------------------\n"; cout<<"PndSttMvdTracking, plottamenti, evt. n. "<0.){ if( ultimoangolo[i]> primoangolo[i]) ultimoangolo[i]-=2.*PI; primo=ultimoangolo[i]*180./PI; ultimoangolo[i]=primoangolo[i]*180./PI; primoangolo[i]=primo; }else{ if( ultimoangolo[i]0 && doMcComparison) { for( j=0;jAt(MCSkewAloneList[i][j]); MCSkewAloneX[ MCSkewAloneList[i][j] ]=puntator->GetX(); MCSkewAloneY[ MCSkewAloneList[i][j] ]=puntator->GetY(); } if(istampa>1) cout<<"PndSttMvdTracking, evt. "< xmax) xmax = posizSciTil[i][0]; if (posizSciTil[i][1] < ymin) ymin = posizSciTil[i][1]; if (posizSciTil[i][1] > ymax) ymax = posizSciTil[i][1]; } //------------- for( ii=0; ii< Nhits; ii++) { i = ListHitsinTrack[iTrack][ii] ; if (info[i][0]-info[i][3] < xmin) xmin = info[i][0]-info[i][3]; if (info[i][0]+info[i][3] > xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } for( i=0; i< nSkewHitsinTrack; i++) { ii = ListSkewHitsinTrack[iTrack][i] ; aaa=Ox+R*cos(SchosenSkew[ii]); bbb=Oy+R*sin(SchosenSkew[ii]); if (aaa < xmin) xmin = aaa; if (aaa > xmax) xmax = aaa; if (bbb < ymin) ymin = bbb; if (bbb > ymax) ymax = bbb; } // for( ii=0; ii< nMvdStripHit; ii++) { for( i=0; i< nMvdStripHitsAssociatedToSttTra; i++) { ii = ListStripHitsinTrack[iTrack][i]; if (XMvdStrip[ii] < xmin) xmin = XMvdStrip[ii]; if (XMvdStrip[ii] > xmax) xmax = XMvdStrip[ii] ; if (YMvdStrip[ii] < ymin) ymin = YMvdStrip[ii]; if (YMvdStrip[ii] > ymax) ymax = YMvdStrip[ii]; } // for( ii=0; ii< nMvdPixelHit; ii++) { for( i=0; i< nMvdPixelHitsAssociatedToSttTra; i++) { ii = ListPixelHitsinTrack[iTrack][i]; if (XMvdPixel[ii] < xmin) xmin = XMvdPixel[ii]; if (XMvdPixel[ii] > xmax) xmax = XMvdPixel[ii] ; if (YMvdPixel[ii] < ymin) ymin = YMvdPixel[ii]; if (YMvdPixel[ii] > ymax) ymax = YMvdPixel[ii]; } for( i=0; i< nMCMvdPixelAlone; i++) { ii = MCMvdPixelAloneList[i]; if (XMvdPixel[ii] < xmin) xmin = XMvdPixel[ii]; if (XMvdPixel[ii] > xmax) xmax = XMvdPixel[ii] ; if (YMvdPixel[ii] < ymin) ymin = YMvdPixel[ii]; if (YMvdPixel[ii] > ymax) ymax = YMvdPixel[ii]; } for( i=0; i< nMCMvdStripAlone; i++) { ii = MCMvdStripAloneList[i]; if (XMvdStrip[ii] < xmin) xmin = XMvdStrip[ii]; if (XMvdStrip[ii] > xmax) xmax = XMvdStrip[ii] ; if (YMvdStrip[ii] < ymin) ymin = YMvdStrip[ii]; if (YMvdStrip[ii] > ymax) ymax = YMvdStrip[ii]; } for( ii=0; ii< nMCParalAlone[iTrack]; ii++) { i = MCParalAloneList[iTrack][ii] ; if( info[i][0] < xmin) xmin = info[i][0]; if( info[i][0] > xmax) xmax = info[i][0]; if( info[i][1] < ymin) ymin = info[i][1]; if( info[i][1] > ymax) ymax = info[i][1]; } for( ii=0; ii< nMCSkewAlone[iTrack]; ii++) { i = MCSkewAloneList[iTrack][ii] ; if( info[i][0] < xmin) xmin = info[i][0]; if( info[i][0] > xmax) xmax = info[i][0]; if( info[i][1] < ymin) ymin = info[i][1]; if( info[i][1] > ymax) ymax = info[i][1]; } if( xmin > 0. ) xmin = 0.; if( xmax < 0.) xmax = 0.; if( ymin > 0. ) ymin = 0.; if( ymax < 0.) ymax = 0.; deltax = xmax-xmin; deltay = ymax - ymin; if( deltax > deltay) { ymin -= 0.5*(deltax-deltay); ymax = ymin+ deltax; delta = deltax; } else { xmin -= 0.5*(deltay-deltax); xmax = xmin+ deltay; delta= deltay; } xmax = xmax + delta*0.15; xmin = xmin - delta*0.15; ymax = ymax + delta*0.15; ymin = ymin - delta*0.15; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); // disegna il BiHexagon destro e sinistro delle inner parallel straws. char myname[100]; sprintf(myname, "InnerPar"); DrawBiHexagonInMacro( VERTICALGAP, MACRO, RStrawDetectorMin, ApotemaMaxInnerParStraw, 4, // color code, 4= blue. myname ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. sprintf(myname, "Skew"); DrawBiHexagonInMacro( VERTICALGAP, MACRO, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, 2, // color code. myname ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. sprintf(myname, "OuterPar"); DrawHexagonCircleInMacro( VERTICALGAP, MACRO, ApotemaMinOuterParStraw, RStrawDetectorMax, 4, // color code. myname ); //-------------- //---- disegna gli Scitil. for( ii=0; ii< nSciTilHitsinTrack[iTrack]; ii++) { i = ListSciTilHitsinTrack[iTrack][ii] ; disegnaSciTilHit( MACRO, i, posizSciTil[i][0], posizSciTil[i][1], 0 ); } //------------------------ fprintf(MACRO,"TEllipse* FoundTrack = new TEllipse(%f,%f,%f,%f,%f,%f);\n" ,Ox,Oy,R,R,primoangolo,ultimoangolo); fprintf(MACRO, "FoundTrack->SetLineColor(2);\nFoundTrack->SetFillStyle(0);\nFoundTrack->Draw(\"only\");\n"); disegnaAssiXY(MACRO,xmin,xmax,ymin,ymax); //------------- hits paralleli in comune con traccia MC for( ii=0; ii< nParalCommon[iTrack]; ii++) { i = ParalCommonList[iTrack][ii] ; fprintf(MACRO, "TEllipse* CommonParalHit%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nCommonParalHit%d->SetFillStyle(0);\nCommonParalHit%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } //------------- hits paralleli spuri for( ii=0; ii< nSpuriParinTrack[iTrack]; ii++) { i = ParSpuriList[iTrack][ii] ; fprintf(MACRO, "TEllipse* SpurParalHit%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nSpurParalHit%d->SetFillStyle(0);\nSpurParalHit%d->SetLineColor(2);\nSpurParalHit%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i,i); } //------------- hits paralleli MC 'alone' for( ii=0; ii< nMCParalAlone[iTrack]; ii++) { i = MCParalAloneList[iTrack][ii] ; fprintf(MACRO, "TEllipse* AloneParalHit%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nAloneParalHit%d->SetFillStyle(0);\nAloneParalHit%d->SetLineColor(4);\nAloneParalHit%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i,i); } //------------- //------------- hits skew in comune e spuri di traccia MC for( i=0; i< nSkewHitsinTrack; i++) { ii = ListSkewHitsinTrack[iTrack][i]; aaa=Ox+R*cos(SchosenSkew[ii]); bbb=Oy+R*sin(SchosenSkew[ii]); bool flaggo = true; for( int k=0; kSetMarkerColor(1);\nCommonSkewHit%d->Draw();\n" ,ii,ii); flaggo = false; break; } } if(flaggo){ fprintf(MACRO,"TMarker* SpurSkewHit%d = new TMarker(%f,%f,%d);\n", ii,aaa,bbb,28); fprintf(MACRO,"SpurSkewHit%d->SetMarkerColor(2);\nSpurSkewHit%d->Draw();\n",ii,ii); } } //------------- hits paralleli MC 'alone' for( ii=0; ii< nMCSkewAlone[iTrack]; ii++) { i = MCSkewAloneList[iTrack][ii] ; fprintf(MACRO, "TMarker* AloneSkewHit%d = new TMarker(%f,%f,%d);\nAloneSkewHit%d->SetMarkerColor(4);\nAloneSkewHit%d->Draw();\n", i,MCSkewAloneX[i],MCSkewAloneY[i],28,i,i); } //------------- now the Strips for( i=0; i< nMvdStripHitsAssociatedToSttTra; i++) { ii = ListStripHitsinTrack[iTrack][i]; x1= XMvdStrip[ii]-sigmaXMvdStrip[ii]; x2= XMvdStrip[ii]+sigmaXMvdStrip[ii]; y1= YMvdStrip[ii]-sigmaYMvdStrip[ii]; y2= YMvdStrip[ii]+sigmaYMvdStrip[ii]; bool flaggo=true; for( int k=0; kSetMarkerColor(1);\nCommonMvdStrip%d->Draw();\n", ii,ii); flaggo=false; break; } } if(flaggo){ fprintf(MACRO,"TMarker* SpurMvdStrip%d = new TMarker(%f,%f,%d);\n", ii,XMvdStrip[ii],YMvdStrip[ii],25); fprintf(MACRO,"SpurMvdStrip%d->SetMarkerColor(2);\nSpurMvdStrip%d->Draw();\n",ii,ii); } } //-------------- hit Mvd Strips 'Alone' for( ii=0; ii< nMCMvdStripAlone; ii++) { i = MCMvdStripAloneList[ii] ; fprintf(MACRO, "TMarker* AloneMvdStrip%d = new TMarker(%f,%f,%d);\nAloneMvdStrip%d->SetMarkerColor(4);\nAloneMvdStrip%d->Draw();\n", i,XMvdStrip[i],YMvdStrip[i],25,i,i); } //------------- // for( ii=0; ii< nMvdPixelHit; ii++) { for( i=0; i< nMvdPixelHitsAssociatedToSttTra; i++) { ii = ListPixelHitsinTrack[iTrack][i]; x1= XMvdPixel[ii]-sigmaXMvdPixel[ii]; x2= XMvdPixel[ii]+sigmaXMvdPixel[ii]; y1= YMvdPixel[ii]-sigmaYMvdPixel[ii]; y2= YMvdPixel[ii]+sigmaYMvdPixel[ii]; // fprintf(MACRO,"TBox* BP%d = new TBox(%f,%f,%f,%f);\nBP%d->SetFillColor(4);\nBP%d->Draw();\n", // ii,x1,y1,x2,y2,ii,ii); bool flaggo=true; for( int k=0; kSetMarkerColor(1);\nCommonMvdPixel%d->Draw();\n", ii,ii); flaggo=false; break; } } if(flaggo){ fprintf(MACRO,"TMarker* SpurMvdPixel%d = new TMarker(%f,%f,%d);\n", ii,XMvdPixel[ii],YMvdPixel[ii],26); fprintf(MACRO,"SpurMvdPixel%d->SetMarkerColor(2);\nSpurMvdPixel%d->Draw();\n",ii,ii); } } //-------------- hit Mvd Pixel 'Alone' for( ii=0; ii< nMCMvdPixelAlone; ii++) { i = MCMvdPixelAloneList[ii] ; fprintf(MACRO, "TMarker* AloneMvdPixel%d = new TMarker(%f,%f,%d);\nAloneMvdPixel%d->SetMarkerColor(4);\nAloneMvdPixel%d->Draw();\n", i,XMvdPixel[i],YMvdPixel[i],26,i,i); } //------------- //----------------- ora la traccia MC corrispondente a questa traccia Stt if( daSttTrackaMCTrack>-1 ) { Int_t icode, im; Double_t Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Px, Py, carica ; PndMCTrack* pMC; im=daSttTrackaMCTrack; pMC = (PndMCTrack*) fMCTrackArray->At(im); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if (fabs(carica)>=0.1 ){ Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); fprintf(MACRO, "TEllipse* MC%d = new TEllipse(%f,%f,%f,%f,%f,%f);\nMC%d->SetFillStyle(0);\nMC%d->SetLineColor(3);\nMC%d->Draw(\"only\");\n", im,Cx,Cy,Rr,Rr,0.,360.,im,im,im); } // end of if (fabs(carica)>=0.1 ) } // end if ( pMC ) };// end of if( daSttTrackaMCTrack>-1 //----------- fine parte del MC fprintf(MACRO,"}\n"); fclose(MACRO); return ; } //--end of function PndSttMvdTracking::WriteMacroSttParallelAssociatedHitsandMvdwithMC //----------start of function PndSttMvdTracking::WriteMacroParallelHitsGeneral void PndSttMvdTracking::WriteMacroParallelHitsGeneral( Int_t Nhits, Double_t info[][7], Short_t nTracksFoundSoFar, bool *keepit, Double_t *Ox, Double_t *Oy, Double_t *Radius, Double_t *FI0, Double_t *ultimoangolo, Double_t *primoangolo ) { Int_t i, j, i1, ii, index, Kincl, nlow, nup, STATUS; Double_t xmin , xmax, ymin, ymax, xl, xu, yl, yu, gamma, dx, dy, diff, d1, d2, delta, deltax, deltay, deltaz, deltaS, factor,ff, zmin, zmax, Smin, Smax, S1, S2, z1, z2, y1, y2,x1,x2, vx1, vy1, vz1, C0x1, C0y1, C0z1, aaa, bbb, ccc, rrr, angle, minor, major, distance, Rx, Ry, LL, Aellipsis1, Bellipsis1,fi1, fmin, fmax, offset, step, SkewInclWithRespectToS, zpos, zpos1, zpos2, Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; char nome[300], nome2[300]; //---------- parallel straws Macro now sprintf(nome,"MacroSttMvdAllHitsEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); xmin=1.e20; xmax=-1.e20; ymin=1.e20; ymax=-1.e20; for( i=0; i< Nhits; i++) { // all straws, anche le skew if (info[i][0]-info[i][3] < xmin) xmin = info[i][0]-info[i][3]; if (info[i][0]+info[i][3] > xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } for( ii=0; ii< nMvdPixelHit; ii++) { if (XMvdPixel[ii] < xmin) xmin = XMvdPixel[ii]; if (XMvdPixel[ii] > xmax) xmax = XMvdPixel[ii] ; if (YMvdPixel[ii] < ymin) ymin = YMvdPixel[ii]; if (YMvdPixel[ii] > ymax) ymax = YMvdPixel[ii]; } for( ii=0; ii< nMvdStripHit; ii++) { if (XMvdStrip[ii] < xmin) xmin = XMvdStrip[ii]; if (XMvdStrip[ii] > xmax) xmax = XMvdStrip[ii] ; if (YMvdStrip[ii] < ymin) ymin = YMvdStrip[ii]; if (YMvdStrip[ii] > ymax) ymax = YMvdStrip[ii]; } // SciTil hits. for( ii=0; ii< nSciTilHits; ii++) { if (posizSciTil[ii][0] < xmin) xmin = posizSciTil[ii][0]; if (posizSciTil[ii][0] > xmax) xmax = posizSciTil[ii][0] ; if (posizSciTil[ii][1] < ymin) ymin = posizSciTil[ii][1]; if (posizSciTil[ii][1] > ymax) ymax = posizSciTil[ii][1]; } // if( xmin > 0. ) xmin = 0.; // if( xmax < 0.) xmax = 0.; // if( ymin > 0. ) ymin = 0.; // if( ymax < 0.) ymax = 0.; deltax = xmax-xmin; deltay = ymax - ymin; if( deltax > deltay) { ymin -= 0.5*(deltax-deltay); ymax = ymin+ deltax; delta = deltax; } else { xmin -= 0.5*(deltay-deltax); xmax = xmin+ deltay; delta= deltay; } xmax = xmax + delta*0.15; xmin = xmin - delta*0.15; ymax = ymax + delta*0.15; ymin = ymin - delta*0.15; if( xmin>-1.05*RStrawDetectorMax) xmin=-1.05*RStrawDetectorMax; if( ymin>-1.05*RStrawDetectorMax) ymin=-1.05*RStrawDetectorMax; if( xmax<1.05*RStrawDetectorMax) xmax=1.05*RStrawDetectorMax; if( ymax<1.05*RStrawDetectorMax) ymax=1.05*RStrawDetectorMax; // ymin=xmin=-1.05*RStrawDetectorMax; // ymax=xmax= 1.05*RStrawDetectorMax; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); disegnaAssiXY(MACRO,xmin,xmax,ymin,ymax); // disegna il BiHexagon destro e sinistro delle inner parallel straws. char myname[100]; sprintf(myname,"InnerPar"); DrawBiHexagonInMacro( VERTICALGAP, MACRO, RStrawDetectorMin, ApotemaMaxInnerParStraw, 4, // color code, 4= blue. myname ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. sprintf(myname,"Skew"); DrawBiHexagonInMacro( VERTICALGAP, MACRO, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, 2, // color code. myname ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. sprintf(myname,"OuterPar"); DrawHexagonCircleInMacro( VERTICALGAP, MACRO, ApotemaMinOuterParStraw, RStrawDetectorMax, 4, // color code. myname ); //-------------- for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO,"TEllipse* Paral%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nParal%d->SetFillStyle(0);\nParal%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } else { // skew straws. fprintf(MACRO,"TMarker* Skew%d = new TMarker(%f,%f,%d);\nSkew%d->SetMarkerColor(1);\nSkew%d->Draw();\n", i,info[i][0],info[i][1],28,i,i,i); } } for( ii=0; ii< nMvdStripHit; ii++) { x1= XMvdStrip[ii]-sigmaXMvdStrip[ii]; x2= XMvdStrip[ii]+sigmaXMvdStrip[ii]; y1= YMvdStrip[ii]-sigmaYMvdStrip[ii]; y2= YMvdStrip[ii]+sigmaYMvdStrip[ii]; fprintf(MACRO, "TMarker* Strip%d = new TMarker(%f,%f,%d);\nStrip%d->SetMarkerColor(1);\nStrip%d->Draw();\n", ii,XMvdStrip[ii],YMvdStrip[ii],25,ii,ii,ii); } for( ii=0; ii< nMvdPixelHit; ii++) { x1= XMvdPixel[ii]-sigmaXMvdPixel[ii]; x2= XMvdPixel[ii]+sigmaXMvdPixel[ii]; y1= YMvdPixel[ii]-sigmaYMvdPixel[ii]; y2= YMvdPixel[ii]+sigmaYMvdPixel[ii]; fprintf(MACRO, "TMarker* Pixel%d = new TMarker(%f,%f,%d);\nPixel%d->SetMarkerColor(1);\nPixel%d->Draw();\n", ii,XMvdPixel[ii],YMvdPixel[ii],26,ii,ii); } //---- disegna gli Scitil. for( i=0; i< nSciTilHits; i++) { disegnaSciTilHit( MACRO, i, posizSciTil[i][0], posizSciTil[i][1], 0 // 0 --> disegna in XY. ); } //------------------------ //------------------------------- plotting all the tracks found for(i=0, ii=0; iSetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", // i,aaa,bbb,rrr,rrr,i,i,i); "TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,%f,%f);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw(\"only\");\n", ii,aaa,bbb,rrr,rrr,primoangolo[i],ultimoangolo[i],ii,ii,ii); } // ----------- fprintf(MACRO,"}\n"); fclose(MACRO); //------------------------------------------------------------------------------------------------------------ //---------- parallel straws Macro now con anche le tracce MC if(doMcComparison){ sprintf(nome,"MacroSttMvdAllHitswithMCEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); disegnaAssiXY(MACRO,xmin,xmax,ymin,ymax); // disegna il BiHexagon destro e sinistro delle inner parallel straws. sprintf(myname,"InnerPar"); DrawBiHexagonInMacro( VERTICALGAP, MACRO, RStrawDetectorMin, ApotemaMaxInnerParStraw, 4, // color code, 4= blue. myname ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. sprintf(myname,"Skew"); DrawBiHexagonInMacro( VERTICALGAP, MACRO, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, 2, // color code. myname ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. sprintf(myname,"OuterPar"); DrawHexagonCircleInMacro( VERTICALGAP, MACRO, ApotemaMinOuterParStraw, RStrawDetectorMax, 4, // color code. myname ); //-------------- for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO,"TEllipse* Paral%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nParal%d->SetFillStyle(0);\nParal%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } else { // skew straws. fprintf(MACRO,"TMarker* Skew%d = new TMarker(%f,%f,%d);\nSkew%d->SetMarkerColor(1);\nSkew%d->Draw();\n", i,info[i][0],info[i][1],28,i,i,i); } } for( ii=0; ii< nMvdStripHit; ii++) { x1= XMvdStrip[ii]-sigmaXMvdStrip[ii]; x2= XMvdStrip[ii]+sigmaXMvdStrip[ii]; y1= YMvdStrip[ii]-sigmaYMvdStrip[ii]; y2= YMvdStrip[ii]+sigmaYMvdStrip[ii]; fprintf(MACRO, "TMarker* Strip%d = new TMarker(%f,%f,%d);\nStrip%d->SetMarkerColor(1);\nStrip%d->Draw();\n", ii,XMvdStrip[ii],YMvdStrip[ii],25,ii,ii,ii); } for( ii=0; ii< nMvdPixelHit; ii++) { x1= XMvdPixel[ii]-sigmaXMvdPixel[ii]; x2= XMvdPixel[ii]+sigmaXMvdPixel[ii]; y1= YMvdPixel[ii]-sigmaYMvdPixel[ii]; y2= YMvdPixel[ii]+sigmaYMvdPixel[ii]; // fprintf(MACRO,"TBox* BP%d = new TBox(%f,%f,%f,%f);\nBP%d->SetFillColor(4);\nBP%d->Draw();\n", // ii,x1,y1,x2,y2,ii,ii); fprintf(MACRO, "TMarker* Pixel%d = new TMarker(%f,%f,%d);\nPixel%d->SetMarkerColor(1);\nPixel%d->Draw();\n", ii,XMvdPixel[ii],YMvdPixel[ii],26,ii,ii); } //---- disegna gli Scitil. for( i=0; i< nSciTilHits; i++) { disegnaSciTilHit( MACRO, i, posizSciTil[i][0], posizSciTil[i][1], 0 ); } //------------------------ //------------------------------- plotting all the tracks found for(i=0, ii=0; iSetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", // i,aaa,bbb,rrr,rrr,i,i,i); "TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,%f,%f);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw(\"only\");\n", ii,aaa,bbb,rrr,rrr,primoangolo[i],ultimoangolo[i],ii,ii,ii); } // ----------- //----------------- ora le traccia MC for(i=0; iAt(i); TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); if( fabs(carica) > 0.1 ){ Rr = aaa*1000./(BFIELD*CVEL); //R(cm) of Helix of track projected in XY plane; // B = 2 Tesla Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); // calcolo per plottare solo la parte rilevante della traccia MC. primo=alfa0 = atan2(Oyy-Cy, Oxx-Cx); bool flaggo=true; for(j=0;j<90;j++){ newalfa = alfa0 - carica*j*PI/45; newx = Cx + Rr*cos(newalfa); newy = Cy + Rr*sin(newalfa); if(newx > xmax || newx < xmin || newy>ymax||newy ultimo ) { primo = ultimo; ultimo = alfa0;}; flaggo=false; break; } } if(flaggo){ primo = 0.; ultimo = 2.*PI; } fprintf(MACRO, "TEllipse* MC%d = new TEllipse(%f,%f,%f,%f,%f,%f);\nMC%d->SetFillStyle(0);\nMC%d->SetLineColor(3);\nMC%d->Draw(\"only\");\n", i,Cx,Cy,Rr,Rr,primo*180./PI,ultimo*180./PI,i,i,i); } else { // continuation of if( fabs(carica) > 0.1 ) continue; // per ora non plotto le neutre. // this is a neutral particle (green dashed straight line). if( Oxx < xmin || Oxx > xmax || Oyy < ymin || Oyy > ymax) continue; double time, time1, time2; if(fabs(Px) >1.e-10 ){ time2 = (xmax-Oxx)/Px; time1 = (xmin-Oxx)/Px; if(time1<0. && time2<0.) { continue; }else if(time2>0.&&time1<0.){ time = time2; x2 = xmax; }else if (time1>0. && time2<0.) { time = time1; x2 = xmin; } else{ if( time2>time1 ) { time=time2; x2 = xmax; } else { time=time1; x2 = xmin; } } y2 = time*Py + Oyy; fprintf(MACRO,"TLine* MCneut%d = new TLine(%f,%f,%f,%f);\n", i,Oxx,Oyy,x2,y2); fprintf(MACRO,"MCneut%d->SetLineStyle(2);\n",i); fprintf(MACRO,"MCneut%d->SetLineColor(3);\nMCneut%d->Draw(\"only\");\n" ,i,i); } else if( fabs(Px) >1.e-10 ) { if( fabs(Py) >1.e-10 ) continue; time2 = (ymax-Oyy)/Py; time1 = (ymin-Oyy)/Py; if(time1<0. && time2<0.) { continue; }else if(time2>0.&&time1<0.){ y2 = ymax; }else if (time1>0. && time2<0.) { y2 = ymin; } else{ y2= time2>time1 ? ymax : ymin ; } fprintf(MACRO,"TLine* MCneut%d = new TLine(%f,%f,%f,%f);\n", i,Oxx,Oyy,Oxx,y2); fprintf(MACRO,"MCneut%d->SetLineStyle(2);\n",i); fprintf(MACRO,"MCneut%d->SetLineColor(3);\nMCneut%d->Draw(\"only\");\n" ,i,i); } } // end of if( fabs(carica) > 0.1 ) } // end of if ( pMC ) } // end of for(i=0; i=0.){ i = (Int_t) time; sprintf(nome,"MacroSttMvdAllHitsTime%dEvent%d",i, IVOLTE); } else { i = (Int_t) -time; sprintf(nome,"MacroSttMvdAllHitsTime-%dEvent%d",i, IVOLTE); } sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); // fprintf(MACRO,"void %s()\n{\n",nome); fprintf(MACRO,"{\n"); xmin=1.e20; xmax=-1.e20; ymin=1.e20; ymax=-1.e20; for( i=0; i< Nhits; i++) { // all straws, anche le skew if (info[i][0]-info[i][3] < xmin) xmin = info[i][0]-info[i][3]; if (info[i][0]+info[i][3] > xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } for( ii=0; ii< nMvdPixelHit; ii++) { if (XMvdPixel[ii] < xmin) xmin = XMvdPixel[ii]; if (XMvdPixel[ii] > xmax) xmax = XMvdPixel[ii] ; if (YMvdPixel[ii] < ymin) ymin = YMvdPixel[ii]; if (YMvdPixel[ii] > ymax) ymax = YMvdPixel[ii]; } for( ii=0; ii< nMvdStripHit; ii++) { if (XMvdStrip[ii] < xmin) xmin = XMvdStrip[ii]; if (XMvdStrip[ii] > xmax) xmax = XMvdStrip[ii] ; if (YMvdStrip[ii] < ymin) ymin = YMvdStrip[ii]; if (YMvdStrip[ii] > ymax) ymax = YMvdStrip[ii]; } if( xmin > 0. ) xmin = 0.; if( xmax < 0.) xmax = 0.; if( ymin > 0. ) ymin = 0.; if( ymax < 0.) ymax = 0.; deltax = xmax-xmin; deltay = ymax - ymin; if( deltax > deltay) { ymin -= 0.5*(deltax-deltay); ymax = ymin+ deltax; delta = deltax; } else { xmin -= 0.5*(deltay-deltax); xmax = xmin+ deltay; delta= deltay; } xmax = xmax + delta*0.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); if( time>=0.){ i = (Int_t) time; fprintf(MACRO,"TPaveLabel* pt= new TPaveLabel(%f,%f,%f,%f,\"Bkg; Time shift = %d nsec\");\n", xmin+delta*0.05,ymax-delta*0.15,xmin+delta*0.35,ymax- delta*0.05,i); } else { i = (Int_t) -time; fprintf(MACRO,"TPaveLabel* pt= new TPaveLabel(%f,%f,%f,%f,\"Bkg; Time shift = -%d nsec\");\n", xmin+delta*0.05,ymax-delta*0.15,xmin+delta*0.35,ymax- delta*0.05,i); } fprintf(MACRO,"pt->Draw();\n"); disegnaAssiXY(MACRO, xmin, xmax, ymin, ymax); for( i=0; i< Nhits; i++) { if( info[i][3]+dR >STRAWRADIUS || info[i][3]+dR <0.) continue; if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO, "TEllipse* Paral%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nParal%d->SetFillStyle(0);\nParal%d->Draw();\n", i,info[i][0],info[i][1],info[i][3]+dR,info[i][3]+dR,i,i); } } if(fabs(time)<5.){ for( ii=0; ii< nMvdStripHit; ii++) { x1= XMvdStrip[ii]-sigmaXMvdStrip[ii]; x2= XMvdStrip[ii]+sigmaXMvdStrip[ii]; y1= YMvdStrip[ii]-sigmaYMvdStrip[ii]; y2= YMvdStrip[ii]+sigmaYMvdStrip[ii]; fprintf(MACRO, "TMarker* Strip%d = new TMarker(%f,%f,%d);\nStrip%d->SetMarkerColor(1);\nStrip%d->Draw();\n", ii,XMvdStrip[ii],YMvdStrip[ii],25,ii,ii,ii); } for( ii=0; ii< nMvdPixelHit; ii++) { x1= XMvdPixel[ii]-sigmaXMvdPixel[ii]; x2= XMvdPixel[ii]+sigmaXMvdPixel[ii]; y1= YMvdPixel[ii]-sigmaYMvdPixel[ii]; y2= YMvdPixel[ii]+sigmaYMvdPixel[ii]; fprintf(MACRO, "TMarker* Pixel%d = new TMarker(%f,%f,%d);\nPixel%d->SetMarkerColor(1);\nPixel%d->Draw();\n", ii,XMvdPixel[ii],YMvdPixel[ii],26,ii,ii); } } // end of if(fabs(time)>5.) fprintf(MACRO,"}\n"); fclose(MACRO); //------------------------------------------------------------------------------------------------------------ return ; } //----------end of function PndSttMvdTracking::WriteMacroParallelHitsGeneralspecial //----------start of function PndSttMvdTracking::WriteMacroSkewAssociatedHitswithMC void PndSttMvdTracking::WriteMacroSkewAssociatedHitswithMC( Double_t KAPPA, Double_t FI0, Double_t Ox, Double_t Oy, Double_t R, Short_t charge, Double_t info[][7], Double_t WDX[nmaxSttHits], Double_t WDY[nmaxSttHits], Double_t WDZ[nmaxSttHits], Short_t iTrack, Int_t iNome, // questo e' per il nome delle Macro solamente. Short_t nSkewHitsinTrack, Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Short_t nSkewCommon, Short_t SkewCommonList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Short_t daTrackFoundaTrackMC, Short_t nMCSkewAlone, Short_t MCSkewAloneList[MAXMCTRACKS][nmaxSttHitsInTrack], Short_t nPixelHitsinTrack[MAXTRACKSPEREVENT], // output Short_t ListPixelHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdPixelHitsInTrack], // output Short_t nStripHitsinTrack[MAXTRACKSPEREVENT], // output Short_t ListStripHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdStripHitsInTrack], // output Short_t nMvdPixelCommon, Short_t *MvdPixelCommonList, Short_t nMvdPixelSpuriinTrack, Short_t *MvdPixelSpuriList, Short_t nMCMvdPixelAlone, Short_t *MCMvdPixelAloneList, Short_t nMvdStripCommon, Short_t *MvdStripCommonList, Short_t nMvdStripSpuriinTrack, Short_t *MvdStripSpuriList, Short_t nMCMvdStripAlone, Short_t *MCMvdStripAloneList, Double_t *ESSE ) { TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); Int_t i, j, i1, ii, iii, index, Kincl, nlow, nup, STATUS, imc, Nmin, Nmax; Double_t xmin , xmax, ymin, ymax, dx, dy, diff, d1, d2, delta, deltax, deltay, deltaz, deltaS, esse,factor, zmin, zmax, zmin2, zmax2, Smin, Smax, S1, S2, z1, z2, y1, y2, vx1, vy1, vz1, C0x1, C0y1, C0z1, aaa, bbb, ccc, angle, minor, major, distance, Rx, Ry, LL, Aellipsis1, Bellipsis1,fi1, fmin, fmax, offset, step, SkewInclWithRespectToS, zpos, zpos1, zpos2, dist[2], Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; //------------------- skew straws hits Macro now char nome2[300],nome[300]; FILE *MACRO; sprintf(nome, "MacroSttMvdSZwithMCEvent%dT%d", IVOLTE,iNome); sprintf(nome2, "%s.C",nome); MACRO = fopen(nome2,"w"); fprintf(MACRO,"{\n"); //KAPPA = 1./166.67 ; FI0 = 1.5*PI; index=0; //--------------- ricerca del minimo e massimo. if( nSciTilHitsinTrack[iTrack]+nSkewHitsinTrack+ nPixelHitsinTrack[iTrack]+ nStripHitsinTrack[iTrack] == 1) // solo 1 punto da disegnare, in questo caso aggiunge { // un punto finto in (0, FI0 ). Smax=Smin= FI0; zmax =zmin =0.; } else { Smin=zmin = 1.e10; Smax=zmax = -zmin; } //--------------------- // prima lo (gli) hits SciTil; assumo al massimo 1 SciTil per traccia; // quindi nSciTilHits[iTrack] puo' essere 0 o 1. for(i=0; iSmax ) Smax=ESSE[i]; if( ESSE[i]zmax ) zmax=posizSciTil[j][2]; if( posizSciTil[j][2]= info[i][4] + Aellipsis1) continue; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; } // end of for( ii=0; ii<2; ii++) } // end of for( iii=0; iii< nSkewHitsinTrack; iii++) //------ aggiungo in blu eventuali punti della traccia MC che sono non mecciati for( iii=0; iii< nMCSkewAlone; iii++) { i = MCSkewAloneList[iTrack][iii]; aaa = sqrt(WDX[i]*WDX[i]+WDY[i]*WDY[i]+ WDZ[i]*WDZ[i]); vx1 = WDX[i]/aaa; vy1 = WDY[i]/aaa; vz1 = WDZ[i]/aaa; C0x1 = info[i][0]; C0y1 = info[i][1]; C0z1 = info[i][2]; calculateintersections(Ox,Oy,R,C0x1,C0y1,C0z1,info[i][3], vx1,vy1,vz1, &STATUS,POINTS1); if(STATUS < 0 ) continue ; for( ii=0; ii<2; ii++){ j=3*ii; dist[ii] = sqrt( (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) ); } if(dist[0]>dist[1]) ii=1; else ii=0; distance = dist[ii]; j=3*ii; distance = sqrt( (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) ); Rx = POINTS1[j]-Ox ; // x component Radial vector of cylinder of trajectory Ry = POINTS1[1+j]-Oy ; // y direction Radial vector of cylinder of trajectory aaa = sqrt(Rx*Rx+Ry*Ry); SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ; SkewInclWithRespectToS /= R; bbb = sqrt( SkewInclWithRespectToS*SkewInclWithRespectToS + vz1*vz1); // the tilt direction of this ellipse is (1,0) when major axis along Z direction LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; // } // end of for( ii=0; ii<2; ii++) } // end of for( iii=0; iii< nMCSkewAlone[imaxima]; iii++) //------------------------------- //------ fine aggiunta in blu eventuali punti della traccia MC che sono non mecciati // ora la parte delle Mvd // prima i pixel for(i=0; i ZMvdPixel[ ii ] ) zmin = ZMvdPixel[ ii ]; if( zmax < ZMvdPixel[ii ] ) zmax = ZMvdPixel[ii ]; esse = atan2( YMvdPixel[ ii ]-Oy, XMvdPixel[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; } // poi le strip for(i=0; i ZMvdStrip[ ii ] ) zmin = ZMvdStrip[ ii ]; if( zmax < ZMvdStrip[ii ] ) zmax = ZMvdStrip[ii ]; esse = atan2( YMvdStrip[ ii ]-Oy, XMvdStrip[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; } // ora i pixel 'Alone' for(i=0; i ZMvdPixel[ ii ] ) zmin = ZMvdPixel[ ii ]; if( zmax < ZMvdPixel[ii ] ) zmax = ZMvdPixel[ii ]; esse = atan2( YMvdPixel[ ii ]-Oy, XMvdPixel[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; } // ora le strip 'Alone' for(i=0; i ZMvdStrip[ ii ] ) zmin = ZMvdStrip[ ii ]; if( zmax < ZMvdStrip[ii ] ) zmax = ZMvdStrip[ii ]; esse = atan2( YMvdStrip[ ii ]-Oy, XMvdStrip[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; } if( zmax >= zmin && Smax >= Smin ) { aaa = Smax-Smin; Smin -= aaa*0.2; Smax += aaa*0.2; aaa = zmax-zmin; zmin -= aaa*0.05; zmax += aaa*0.05; if(Smax > 2.*PI) Smax = 2.*PI; if( Smin < 0.) Smin = 0.; deltaz = zmax-zmin; deltaS = Smax-Smin; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n", zmin-0.1*deltaz,R*(Smin-.1*deltaS),zmax+0.1*deltaz,R*(Smax+0.1*deltaS)); // fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",zmin,Smin,zmax,Smax); for( ii=0; ii< index; ii++) { fprintf(MACRO,"E%d->Draw();\n",ii); } fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", // zmin+0.05*deltaz,Smin+0.05*deltaS,zmax-0.05*deltaz,Smin+0.05*deltaS,zmin+0.05*deltaz,zmax-0.05*deltaz); zmin-0.01*deltaz,R*(Smin+0.05*deltaS),zmax+0.01*deltaz,R*(Smin+0.05*deltaS), zmin-0.01*deltaz,zmax+0.01*deltaz); fprintf(MACRO,"Assex->SetTitle(\"Z (cm)\");\n"); fprintf(MACRO,"Assex->Draw();\n"); fprintf(MACRO,"TGaxis *Assey = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", zmin+0.05*deltaz,R*(Smin-0.01*deltaS),zmin+0.05*deltaz,R*(Smax+0.01*deltaS), R*(Smin-0.01*deltaS),R*(Smax+0.01*deltaS)); fprintf(MACRO,"Assey->SetTitle(\"Helix crf (cm)\");\n"); fprintf(MACRO,"Assey->Draw();\n"); //------------ // plot di eventuali hits SciTil; for(i=0;i 1.e-10){ Tiltdirection1[0] = vz1/bbb; Tiltdirection1[1] = SkewInclWithRespectToS/bbb; } else { Tiltdirection1[0] = 1.; Tiltdirection1[1] = 0.; } LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; if( distance >= info[i][4] + Aellipsis1) continue; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; Double_t rotation1 = 180.*atan2(Tiltdirection1[1],Tiltdirection1[0])/PI; fprintf(MACRO,"TEllipse* Skew%d_%d = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\nSkew%d_%d->SetFillStyle(0);\n", // i,ii,POINTS1[j+2],fi1,Aellipsis1,Bellipsis1,rotation1,i,ii); i,ii,POINTS1[j+2],R*fi1,Aellipsis1,R*Bellipsis1,rotation1,i,ii); // ------ se lo hit e' spurio marcalo in rosso bool flaggo=true; for( i1=0; i1SetLineColor(2);\n",i,ii); fprintf(MACRO,"Skew%d_%d->Draw();\n",i,ii); index++; } // end of for( ii=0; ii<2; ii++) } // end of for( iii=0; iii< nSkewHitsinTrack; iii++) //------ aggiungo in blu eventuali punti della traccia MC che sono non mecciati for( iii=0; iii< nMCSkewAlone; iii++) { i = MCSkewAloneList[iTrack][iii]; aaa = sqrt(WDX[i]*WDX[i]+WDY[i]*WDY[i]+ WDZ[i]*WDZ[i]); vx1 = WDX[i]/aaa; vy1 = WDY[i]/aaa; vz1 = WDZ[i]/aaa; C0x1 = info[i][0]; C0y1 = info[i][1]; C0z1 = info[i][2]; calculateintersections(Ox,Oy,R,C0x1,C0y1,C0z1,info[i][3], vx1,vy1,vz1, &STATUS,POINTS1); if(STATUS < 0 ) continue ; for( ii=0; ii<2; ii++){ j=3*ii; dist[ii] = sqrt( (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) ); } if(dist[0]>dist[1]) ii=1; else ii=0; distance = dist[ii]; j=3*ii; Rx = POINTS1[j]-Ox ; // x component Radial vector of cylinder of trajectory Ry = POINTS1[1+j]-Oy ; // y direction Radial vector of cylinder of trajectory aaa = sqrt(Rx*Rx+Ry*Ry); SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ; SkewInclWithRespectToS /= R; bbb = sqrt( SkewInclWithRespectToS*SkewInclWithRespectToS + vz1*vz1); // the tilt direction of this ellipse is (1,0) when major axis along Z direction if( bbb > 1.e-10){ Tiltdirection1[0] = vz1/bbb; Tiltdirection1[1] = SkewInclWithRespectToS/bbb; } else { Tiltdirection1[0] = 1.; Tiltdirection1[1] = 0.; } LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; Double_t rotation1 = 180.*atan2(Tiltdirection1[1],Tiltdirection1[0])/PI; fprintf(MACRO,"TEllipse* AloneSkew%d_%d = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\nAloneSkew%d_%d->SetFillStyle(0);\n", // i,ii,POINTS1[j+2],fi1,Aellipsis1,Bellipsis1,rotation1,i,ii); i,ii,POINTS1[j+2],R*fi1,Aellipsis1,R*Bellipsis1,rotation1,i,ii); // ------ marca lo hit in blu fprintf(MACRO,"AloneSkew%d_%d->SetLineColor(4);\n",i,ii); fprintf(MACRO,"AloneSkew%d_%d->Draw();\n",i,ii); index++; } // end of for( iii=0; iii< nMCSkewAlone; iii++) //------------------------------- //------ fine aggiunta in blu eventuali punti della traccia MC che sono non mecciati // ora la parte delle Mvd // prima i pixel for(i=0; i ZMvdPixel[ ii ] ) zmin = ZMvdPixel[ ii ]; if( zmax < ZMvdPixel[ii ] ) zmax = ZMvdPixel[ii ]; esse = atan2( YMvdPixel[ ii ]-Oy, XMvdPixel[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; bool flaggo=true; for( int k=0; kSetMarkerColor(1);\n", ii,ZMvdPixel[ii],R*esse,26,ii); fprintf(MACRO,"CommonPixel%d->Draw();\n",ii); flaggo=false; break; } } if(flaggo){ fprintf(MACRO,"TMarker* SpuriousPixel%d = new TMarker(%f,%f,%d);\nSpuriousPixel%d->SetMarkerColor(2);\n", ii,ZMvdPixel[ii],R*esse,26,ii); fprintf(MACRO,"SpuriousPixel%d->Draw();\n",ii); } } // poi le strip for(i=0; i ZMvdStrip[ ii ] ) zmin = ZMvdStrip[ ii ]; if( zmax < ZMvdStrip[ii ] ) zmax = ZMvdStrip[ii ]; esse = atan2( YMvdStrip[ ii ]-Oy, XMvdStrip[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; bool flaggo=true; for( int k=0; kSetMarkerColor(1);\n", ii,ZMvdStrip[ii],R*esse,25,ii); fprintf(MACRO,"CommonStrip%d->Draw();\n",ii); flaggo=false; break; } } if(flaggo){ fprintf(MACRO,"TMarker* SpuriousStrip%d = new TMarker(%f,%f,%d);\nSpuriousStrip%d->SetMarkerColor(2);\n", ii,ZMvdStrip[ii],R*esse,25,ii); fprintf(MACRO,"SpuriousStrip%d->Draw();\n",ii); } } // ora i pixel 'Alone' for(i=0; i ZMvdPixel[ ii ] ) zmin = ZMvdPixel[ ii ]; if( zmax < ZMvdPixel[ii ] ) zmax = ZMvdPixel[ii ]; esse = atan2( YMvdPixel[ ii ]-Oy, XMvdPixel[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; fprintf(MACRO,"TMarker* AlonePixel%d = new TMarker(%f,%f,%d);\nAlonePixel%d->SetMarkerColor(4);\n", ii,ZMvdPixel[ii],R*esse,26,ii); fprintf(MACRO,"AlonePixel%d->Draw();\n",ii); } // ora le strip 'Alone' for(i=0; i ZMvdStrip[ ii ] ) zmin = ZMvdStrip[ ii ]; if( zmax < ZMvdStrip[ii ] ) zmax = ZMvdStrip[ii ]; esse = atan2( YMvdStrip[ ii ]-Oy, XMvdStrip[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; fprintf(MACRO,"TMarker* AloneStrip%d = new TMarker(%f,%f,%d);\nAloneStrip%d->SetMarkerColor(4);\n", ii,ZMvdStrip[ii],R*esse,25,ii); fprintf(MACRO,"AloneStrip%d->Draw();\n",ii); } // plot della traccia trovata dal finder zmin2=zmin; zmax2=zmax; bool flaggo=true; if( -KAPPA*charge>0.) { // Pz>0. if( zmax <0.) { cout<<"da WriteMacroSkewAssociatedHitswithMC, questa traccia" <<" e' inconsistente col proprio Pz, non plottata!\n"; flaggo=false; } else { zmin = 0.; } } else { // Pz<0. if( zmin >0.) { cout<<"da WriteMacroSkewAssociatedHitswithMC, questa traccia" <<" e' inconsistente col proprio Pz, non plottata!\n"; flaggo=false; } else{ zmax = 0.; } } if(flaggo){ if ( KAPPA >= 0.) { fmin = KAPPA*zmin + FI0; fmax = KAPPA*zmax + FI0; } else { fmax = KAPPA*zmin + FI0; fmin = KAPPA*zmax + FI0; } if( fmax>=0.) { Nmax = (int) (0.5*fmax/ PI); } else { Nmax = ( (int) (0.5*fmax/ PI) ) -1; } if( fmin>=0.) { Nmin = (int) (0.5*fmin/ PI); } else { Nmin = ((int) (0.5*fmin/ PI) )-1; } if(fabs(KAPPA)<1.e-10) { cout<<"da WriteMacroSkewAssociatedHitswithMC, questa traccia Found da PR non plottata" <<" perche' ha fabs(KAPPA)<1.e-10.\n"; } else { for(i=Nmin; i<= Nmax;i++){ offset = 2.*PI*i; z1 = (i*2.*PI-FI0)/KAPPA; z2 = ((i+1)*2.*PI-FI0)/KAPPA; fprintf(MACRO, "TLine* FOUND%d = new TLine(%f,%f,%f,%f);\nFOUND%d->SetLineColor(2);\nFOUND%d->Draw();\n", i-Nmin,z1,0.,z2, R*2.*PI,i-Nmin,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) } // end of if(fabs(KAPPA)<1.e-10) } // end of if(flaggo) zmin=zmin2; zmax=zmax2; //----------------- ora la traccia MC corrispondente a questa traccia Stt // for(imc=0; imc-1 ) { Int_t icode ; Double_t Rr, Dd, Fifi, Kakka, Oxx, Oyy, Cx, Cy, Px, Py, carica ; PndMCTrack* pMC; pMC = (PndMCTrack*) fMCTrackArray->At(imc); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if (fabs(carica)>=0.1 ){ Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); Fifi = atan2(Cy, Cx); // MC truth Fifi angle of circle of Helix trajectory if(Fifi<0.) Fifi += 2.*PI; if( fabs( pMC->GetMomentum().Z() )< 1.e-20) Kakka = 99999999.; else Kakka = -carica*0.001*BFIELD*CVEL/pMC->GetMomentum().Z(); KAPPA=Kakka; FI0 = fmod(Fifi+ PI, 2.*PI); if ( KAPPA >= 0.) { fmin = KAPPA*zmin + FI0; fmax = KAPPA*zmax + FI0; } else { fmax = KAPPA*zmin + FI0; fmin = KAPPA*zmax + FI0; } if( fmax>=0.) { Nmax = (int) (0.5*fmax/ PI); } else { Nmax = ( (int) (0.5*fmax/ PI) ) -1; } if( fmin>=0.) { Nmin = (int) (0.5*fmin/ PI); } else { Nmin = ((int) (0.5*fmin/ PI) )-1; } for(i=Nmin; i<= Nmax;i++){ offset = 2.*PI*i; z1 = (i*2.*PI-FI0)/KAPPA; z2 = ((i+1)*2.*PI-FI0)/KAPPA; fprintf(MACRO,"TLine* MC%d_%d = new TLine(%f,%f,%f,%f);\nMC%d_%d->SetLineColor(3);\nMC%d_%d->Draw();\n", imc,i-Nmin,z1,0.,z2,R* 2.*PI,imc,i-Nmin,imc,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) } // end of if (fabs(carica)>=0.1 ) } // end of if ( pMC ) } // end of if( imc>-1 ) } // end of if( zmax >= zmin && Smax >= Smin ) fprintf(MACRO,"}\n"); fclose(MACRO); } //----------end of function PndSttMvdTracking::WriteMacroSkewAssociatedHitswithMC //----------begin of function PndSttMvdTracking::WriteMacroAllHitsRestanti void PndSttMvdTracking::WriteMacroAllHitsRestanti( Short_t nSttHit, Short_t nSttParHit, Short_t nSttSkewHit, Double_t info[][7], bool *keepit, Short_t nSttTrackCand, Short_t nCandHit[MAXTRACKSPEREVENT], Short_t ListCandHit[MAXTRACKSPEREVENT][nmaxSttHitsInTrack+ nmaxMvdPixelHitsInTrack+ nmaxMvdStripHitsInTrack+nmaxSciTilHitsinTrack], Short_t ListCandHitType[MAXTRACKSPEREVENT][nmaxSttHitsInTrack+ nmaxMvdPixelHitsInTrack+ nmaxMvdStripHitsInTrack+nmaxSciTilHitsinTrack] ) { // nSttHit = parallel+skew. bool exclusionStt[nmaxSttHits], exclusionPixel[nmaxMvdPixelHits], exclusionStrip[nmaxMvdStripHits]; char nome[300], nome2[300]; int i,j,k; Double_t delta, deltax, deltay, xmin, xmax, ymin, ymax; for(i=0;i xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } for( i=0; i< nMvdPixelHit; i++) { if (XMvdPixel[i] < xmin) xmin = XMvdPixel[i]; if (XMvdPixel[i] > xmax) xmax = XMvdPixel[i]; if (YMvdPixel[i] < ymin) ymin = YMvdPixel[i]; if (YMvdPixel[i] > ymax) ymax = YMvdPixel[i]; } for( i=0; i< nMvdStripHit; i++) { if (XMvdStrip[i] < xmin) xmin = XMvdStrip[i]; if (XMvdStrip[i] > xmax) xmax = XMvdStrip[i]; if (YMvdStrip[i] < ymin) ymin = YMvdStrip[i]; if (YMvdStrip[i] > ymax) ymax = YMvdStrip[i]; } for( i=0; i< nSciTilHits; i++) { if (posizSciTil[i][0] < xmin) xmin = posizSciTil[i][0]; if (posizSciTil[i][0] > xmax) xmax = posizSciTil[i][0]; if (posizSciTil[i][1] < ymin) ymin = posizSciTil[i][1]; if (posizSciTil[i][1] > ymax) ymax = posizSciTil[i][1]; } if( xmin > 0. ) xmin = 0.; if( xmax < 0.) xmax = 0.; if( ymin > 0. ) ymin = 0.; if( ymax < 0.) ymax = 0.; deltax = xmax-xmin; deltay = ymax - ymin; if( deltax > deltay) { ymin -= 0.5*(deltax-deltay); ymax = ymin+ deltax; delta = deltax; } else { xmin -= 0.5*(deltay-deltax); xmax = xmin+ deltay; delta= deltay; } xmax = xmax + delta*0.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; sprintf(nome,"MacroSttMvdHitsRestantiEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); disegnaAssiXY(MACRO,xmin,xmax,ymin,ymax); // disegna il BiHexagon destro e sinistro delle inner parallel straws. char myname[100]; sprintf(myname,"InnerPar"); DrawBiHexagonInMacro( VERTICALGAP, MACRO, RStrawDetectorMin, ApotemaMaxInnerParStraw, 4, // color code, 4= blue. myname ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. sprintf(myname,"Skew"); DrawBiHexagonInMacro( VERTICALGAP, MACRO, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, 2, // color code. myname ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. sprintf(myname,"OuterPar"); DrawHexagonCircleInMacro( VERTICALGAP, MACRO, ApotemaMinOuterParStraw, RStrawDetectorMax, 4, // color code. myname ); //-------------- for( i=0; i< nSttHit; i++) { if( !exclusionStt[i]) { // all straws if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nE%d->SetFillStyle(0);\nE%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } else { // skew straws. fprintf(MACRO,"TMarker* SS%d = new TMarker(%f,%f,%d);\nSS%d->SetMarkerColor(1);\nSS%d->Draw();\n", i,info[i][0],info[i][1],28,i,i,i); } } } for( i=0; i< nMvdPixelHit; i++) { if( !exclusionPixel[i]) { // all Pixels fprintf(MACRO, "TMarker* Pixel%d = new TMarker(%f,%f,%d);\nPixel%d->SetMarkerColor(1);\nPixel%d->Draw();\n", i,XMvdPixel[i],YMvdPixel[i],26,i,i); } } for( i=0; i< nMvdStripHit; i++) { if( !exclusionStrip[i]) { // all Pixels fprintf(MACRO, "TMarker* Strip%d = new TMarker(%f,%f,%d);\nStrip%d->SetMarkerColor(1);\nStrip%d->Draw();\n", i,XMvdStrip[i],YMvdStrip[i],25,i,i,i); } } for( i=0; i< nSciTilHits; i++) { if( InclusionListSciTil[i]) { // all SciTil hit never used. fprintf(MACRO, "TMarker* Strip%d = new TMarker(%f,%f,%d);\nStrip%d->SetMarkerColor(1);\nStrip%d->Draw();\n", i,posizSciTil[i][0],posizSciTil[i][1],30,i,i,i); } } fprintf(MACRO,"}\n"); fclose(MACRO); return; } //----------end of function PndSttMvdTracking::WriteMacroAllHitsRestanti //----------begin of function PndSttMvdTracking::DrawBiHexagonInMacro void PndSttMvdTracking::DrawBiHexagonInMacro( Double_t vgap, FILE * MACRO, Double_t Ami, Double_t Ama, Short_t color, char *name ) { Short_t iside; // these are the points defining the sides of the biexhagon on the left(inner). Double_t side_x[] = { -vgap/2., -Ama , -Ama, -vgap/2., -vgap/2., -Ami, -Ami, -vgap/2., -vgap/2.}, side_y[] = {(-0.5*vgap+2.*Ama)/sqrt(3.), Ama/sqrt(3.), -Ama/sqrt(3.), -(-0.5*vgap+2.*Ama)/sqrt(3.), -(-0.5*vgap+2.*Ami)/sqrt(3.), -Ami/sqrt(3.), Ami/sqrt(3.), (-0.5*vgap+2.*Ami)/sqrt(3.), (-0.5*vgap+2.*Ama)/sqrt(3.) }; // Inner straws left for(iside=0;iside<8;iside++){ fprintf(MACRO, "TLine* %sL%d = new TLine(%f,%f,%f,%f);\n",name, iside,side_x[iside],side_y[iside],side_x[iside+1],side_y[iside+1]); fprintf(MACRO,"%sL%d->SetLineColor(%d);\n",name,iside,color); fprintf(MACRO,"%sL%d->SetLineStyle(2);\n",name,iside); fprintf(MACRO,"%sL%d->Draw();\n",name,iside); } // Inner straws right for(iside=0;iside<8;iside++){ fprintf(MACRO, "TLine* %sR%d = new TLine(%f,%f,%f,%f);\n",name, iside,-side_x[iside],side_y[iside],-side_x[iside+1],side_y[iside+1]); fprintf(MACRO,"%sR%d->SetLineColor(%d);\n",name,iside,color); fprintf(MACRO,"%sR%d->SetLineStyle(2);\n",name,iside); fprintf(MACRO,"%sR%d->Draw();\n",name,iside); } } //----------end of function PndSttMvdTracking::DrawBiHexagonInMacro //----------begin of function PndSttMvdTracking::DrawHexagonCircleInMacro void PndSttMvdTracking::DrawHexagonCircleInMacro( Double_t GAP, FILE * MACRO, Double_t ApotemaMin, Double_t Rma, Short_t color, char *name ) { Short_t iside; Double_t angle1, angle2; // these are the points defining the sides of the exhagon on the left(outer). Double_t side_x[] = { -GAP/2., -GAP/2. , -ApotemaMin, -ApotemaMin, -GAP/2., -GAP/2. }, side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-0.5*GAP)/sqrt(3.), ApotemaMin/sqrt(3.), -ApotemaMin/sqrt(3.), -(2.*ApotemaMin-0.5*GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)}; // Outer straws left for(iside=0;iside<5;iside++){ fprintf(MACRO, "TLine* %sL%d = new TLine(%f,%f,%f,%f);\n",name, iside,side_x[iside],side_y[iside],side_x[iside+1],side_y[iside+1]); fprintf(MACRO,"%sL%d->SetLineColor(%d);\n",name,iside,color); fprintf(MACRO,"%sL%d->SetLineStyle(2);\n",name,iside); fprintf(MACRO,"%sL%d->Draw();\n",name,iside); } // Outer straws right for(iside=0;iside<5;iside++){ fprintf(MACRO, "TLine* %sR%d = new TLine(%f,%f,%f,%f);\n",name, iside,-side_x[iside],side_y[iside],-side_x[iside+1],side_y[iside+1]); fprintf(MACRO,"%sR%d->SetLineColor(%d);\n",name,iside,color); fprintf(MACRO,"%sR%d->SetLineStyle(2);\n",name,iside); fprintf(MACRO,"%sR%d->Draw();\n",name,iside); } // drawing the left circle. angle1 = atan2 ( side_y[0], side_x[0])*180./PI; angle2 = 360. + atan2 ( side_y[4], side_x[4])*180./PI; fprintf(MACRO,"TEllipse* %sCircleL = new TEllipse(0.,0.,%f,%f,%f,%f);\n", name,Rma,Rma,angle1,angle2); fprintf(MACRO,"%sCircleL->SetFillStyle(0);\n",name); fprintf(MACRO,"%sCircleL->SetLineColor(%d);\n",name,color); fprintf(MACRO,"%sCircleL->Draw(\"only\");\n",name,color); // drawing the right circle. angle2 = atan2 ( side_y[0], -side_x[0])*180./PI; angle1 = atan2 ( side_y[4], -side_x[4])*180./PI; fprintf(MACRO,"TEllipse* %sCircleR = new TEllipse(0.,0.,%f,%f,%f,%f);\n", name,Rma,Rma,angle1,angle2); fprintf(MACRO,"%sCircleR->SetFillStyle(0);\n",name); fprintf(MACRO,"%sCircleR->SetLineColor(%d);\n",name,color); fprintf(MACRO,"%sCircleR->Draw(\"only\");\n",name,color); } //----------end of function PndSttMvdTracking::DrawHexagonCircleInMacro //----------begin of function PndSttMvdTracking::DrawSttDetectorInMacro void PndSttMvdTracking::AssociateFoundTrackstoMC( Double_t info[][7], Short_t nTracksFoundSoFar, Short_t nHitsinTrack[MAXTRACKSPEREVENT], Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], Short_t nSkewHitsinTrack[MAXTRACKSPEREVENT], Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool flaggo, inclusionMC[MAXTRACKSPEREVENT][nmaxSttHits], inclusionExp[MAXTRACKSPEREVENT]; Short_t ntoMCtrack[MAXTRACKSPEREVENT], toMCtracklist[MAXTRACKSPEREVENT][nmaxSttHits], toMCtrackfrequency[MAXTRACKSPEREVENT][nmaxSttHits]; Short_t i, j, jtemp,jexp; Short_t itemp, massimo; Int_t enne; for(i=0; i -1){ itemp=-1; massimo = -1; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if( !inclusionExp[jexp]) continue; for(i=0; i< ntoMCtrack[jexp]; i++){ if( !inclusionMC[jexp][i]) continue; if( toMCtrackfrequency[jexp][i]>massimo){ massimo=toMCtrackfrequency[jexp][i]; itemp = toMCtracklist[jexp][i]; jtemp = jexp; } } } if( itemp>-1 ){ daTrackFoundaTrackMC[jtemp]=itemp; inclusionExp[jtemp]=false; for(jexp=0; jexp -1) return; } //----------end of function PndSttMvdTracking::AssociateFoundTrackstoMC //----------begin of function PndSttMvdTracking::AssociateFoundTrackstoMCbis void PndSttMvdTracking::AssociateFoundTrackstoMCbis( Double_t info[][7], Short_t nTracksFoundSoFar, Short_t nHitsinTrack[MAXTRACKSPEREVENT], Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], Short_t nSkewHitsinTrack[MAXTRACKSPEREVENT], Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool flaggo, inclusionMC[MAXTRACKSPEREVENT][nmaxSttHits], inclusionExp[MAXTRACKSPEREVENT]; Short_t ntoMCtrack[MAXTRACKSPEREVENT], toMCtracklist[MAXTRACKSPEREVENT][nmaxSttHits], toMCtrackfrequency[MAXTRACKSPEREVENT][nmaxSttHits]; Short_t i, j, enne, jtemp,jexp; Short_t itemp, massimo; for(i=0; i -1){ itemp=-1; massimo = -1; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if( !inclusionExp[jexp]) continue; for(i=0; i< ntoMCtrack[jexp]; i++){ if( !inclusionMC[jexp][i]) continue; if( toMCtrackfrequency[jexp][i]>massimo){ massimo=toMCtrackfrequency[jexp][i]; itemp = toMCtracklist[jexp][i]; jtemp = jexp; } } } if( itemp>-1 ){ daTrackFoundaTrackMC[jtemp]=itemp; inclusionExp[jtemp]=false; for(jexp=0; jexp -1) return; } //----------end of function PndSttMvdTracking::AssociateFoundTrackstoMCbis //----------begin of function PndSttMvdTracking::AssociateFoundTrackstoMCtris void PndSttMvdTracking::AssociateFoundTrackstoMCtris( Double_t info[][7], Double_t Ox[MAXTRACKSPEREVENT], Double_t Oy[MAXTRACKSPEREVENT], Double_t R[MAXTRACKSPEREVENT], Short_t nTracksFoundSoFar, Short_t nHitsinTrack[MAXTRACKSPEREVENT], Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], Short_t nSkewHitsinTrack[MAXTRACKSPEREVENT], Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool firstime, flaggo, inclusionMC[MAXTRACKSPEREVENT][nmaxSttHits], inclusionExp[MAXTRACKSPEREVENT]; Short_t ntoMCtrack[MAXTRACKSPEREVENT], toMCtrackfrequency[MAXTRACKSPEREVENT][nmaxSttHits]; Short_t i, j, jtemp,jexp , nmid; Short_t itemp, massimo, toMCtracklist[MAXTRACKSPEREVENT][nmaxSttHits]; Int_t enne; Double_t dx, Cx, Cy, Rr, alfa, beta, gamma, minimo, tanlow[MAXTRACKSPEREVENT], tanmid[MAXTRACKSPEREVENT], tanup[MAXTRACKSPEREVENT], toMCtrackdistance[MAXTRACKSPEREVENT][nmaxSttHits]; int nevento=4; for(i=0; i 2) { dx = info[ ListHitsinTrack[i][0] ][0] - Ox[i]; if(fabs(dx)> 1.e-10){ tanlow[i] = ( info[ ListHitsinTrack[i][0] ][1] - Oy[i] )/dx; } else { tanlow[i] = 999999.; } nmid = (Short_t) (nHitsinTrack[i]+1)/2; dx = info[ ListHitsinTrack[i][nmid] ][0] - Ox[i]; if(fabs(dx)> 1.e-10){ tanmid[i] = ( info[ ListHitsinTrack[i][nmid] ][1] - Oy[i] )/dx; } else { tanmid[i] = 999999.; } dx = info[ ListHitsinTrack[i][nHitsinTrack[i]-1] ][0] - Ox[i]; if(fabs(dx)> 1.e-10){ tanup[i] = ( info[ ListHitsinTrack[i][nHitsinTrack[i]-1] ][1] - Oy[i] )/dx; } else { tanup[i] = 999999.; } } // end of if( nHitsinTrack[i] > 2) //----------------------------- } // end of for(i=0; i -1){ itemp=-1; massimo = -1; minimo = 999999999999.; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if( !inclusionExp[jexp]) continue; for(i=0; i< ntoMCtrack[jexp]; i++){ if( !inclusionMC[jexp][i]) continue; // if( toMCtrackfrequency[jexp][i]>massimo){ if( toMCtrackdistance[jexp][i]<-0.5) continue; if( toMCtrackdistance[jexp][i]-1 ){ daTrackFoundaTrackMC[jtemp]=itemp; inclusionExp[jtemp]=false; for(jexp=0; jexp -1) return; } //----------end of function PndSttMvdTracking::AssociateFoundTrackstoMCtris //----------begin of function PndSttMvdTracking::AssociateFoundTrackstoMCquater void PndSttMvdTracking::AssociateFoundTrackstoMCquater( bool *keepit, Double_t info[][7], Double_t Ox[MAXTRACKSPEREVENT], Double_t Oy[MAXTRACKSPEREVENT], Double_t R[MAXTRACKSPEREVENT], Double_t X1[MAXTRACKSPEREVENT], Double_t Y1[MAXTRACKSPEREVENT], Double_t X2[MAXTRACKSPEREVENT], Double_t Y2[MAXTRACKSPEREVENT], Double_t X3[MAXTRACKSPEREVENT], Double_t Y3[MAXTRACKSPEREVENT], Short_t nTracksFoundSoFar, Short_t nHitsinTrack[MAXTRACKSPEREVENT], Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Short_t nSkewHitsinTrack[MAXTRACKSPEREVENT], Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Short_t nPixelHitsinTrack[MAXTRACKSPEREVENT], Short_t ListPixelHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdPixelHitsInTrack], Int_t *FromPixeltoMCTrack, Short_t nStripHitsinTrack[MAXTRACKSPEREVENT], Short_t ListStripHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdStripHitsInTrack], Int_t *FromStriptoMCTrack, Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool firstime, flaggo, inclusionMC[nTracksFoundSoFar][nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], inclusionExp[nTracksFoundSoFar]; Short_t ntoMCtrack[nTracksFoundSoFar], toMCtrackfrequency[nTracksFoundSoFar][nmaxSttHitsInTrack]; Short_t i, j, jtemp,jexp , nmid; Short_t itemp, massimo, toMCtracklist[nTracksFoundSoFar][nmaxSttHitsInTrack]; Int_t enne; Double_t dx, Cx, Cy, Rr, alfa, beta, gamma, minimo, tanlow[nTracksFoundSoFar], tanmid[nTracksFoundSoFar], tanup[nTracksFoundSoFar], toMCtrackdistance[nTracksFoundSoFar][nmaxSttHitsInTrack]; // tanlow[MAXTRACKSPEREVENT], // tanmid[MAXTRACKSPEREVENT], // tanup[MAXTRACKSPEREVENT], // toMCtrackdistance[MAXTRACKSPEREVENT][nmaxSttHits]; for(i=0; i 2 || nn == 1) { dx = X1[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanlow[i] = (Y1[i] - Oy[i] )/dx; } else { tanlow[i] = 999999.; } dx = X2[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanmid[i] = ( Y2[i] - Oy[i] )/dx; } else { tanmid[i] = 999999.; } dx = X3[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanup[i] = (Y3[i] - Oy[i] )/dx; } else { tanup[i] = 999999.; } } else if (nn==2) { // continuation of if( nn > 2 || nn == 1) dx = X1[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanlow[i] = (Y1[i] - Oy[i] )/dx; } else { tanlow[i] = 999999.; } dx = X2[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanmid[i] = ( Y2[i] - Oy[i] )/dx; } else { tanmid[i] = 999999.; } tanup[i] = tanmid[i]; } /* else{ // continuation of if( nn > 2), case in which there is only 1 hit in the Track cand. dx = X1[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanlow[i] = (Y1[i] - Oy[i] )/dx; } else { tanlow[i] = 999999.; } tanmid[i] = tanlow[i]; tanup[i] = tanlow[i]; } // end of if( nn > 2) */ //----------------------------- } // end of for(i=0; i -1){ itemp=-1; massimo = -1; minimo = 999999999999.; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if(!keepit[jexp]) continue; if( !inclusionExp[jexp]) continue; for(i=0; i< ntoMCtrack[jexp]; i++){ if( !inclusionMC[jexp][i]) continue; // if( toMCtrackfrequency[jexp][i]>massimo){ if( toMCtrackdistance[jexp][i]<-0.5) continue; if( toMCtrackdistance[jexp][i]-1 ){ daTrackFoundaTrackMC[jtemp]=itemp; inclusionExp[jtemp]=false; for(jexp=0; jexp -1) return; } //----------end of function PndSttMvdTracking::AssociateFoundTrackstoMCquater //------------------------- begin of function PndSttMvdTracking::SttMatchedSpurious void PndSttMvdTracking::SttMatchedSpurious( bool *keepit, Short_t ntotalHits, Double_t info[][7], Short_t nTracksFoundSoFar, // quelle trovate dal PR Short_t nHitsinTrack[], // n. hits PARALLELI+SKEW, dal PR Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], // dal PR Short_t nSkewHitsinTrack[MAXTRACKSPEREVENT], // n. hits skew, dal PR Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], // dal PR Short_t nParalCommon[MAXTRACKSPEREVENT], Short_t ParalCommonList[MAXMCTRACKS][nmaxSttHitsInTrack], Short_t nSpuriParinTrack[MAXTRACKSPEREVENT], Short_t ParSpuriList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Short_t nSkewCommon[MAXTRACKSPEREVENT], Short_t SkewCommonList[MAXMCTRACKS][nmaxSttHitsInTrack], Short_t nSpuriSkewinTrack[MAXTRACKSPEREVENT], Short_t SkewSpuriList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Short_t nHitsInMCTrack[MAXTRACKSPEREVENT], Short_t nSkewHitsInMCTrack[MAXTRACKSPEREVENT], Short_t nMCParalAlone[MAXTRACKSPEREVENT], Short_t MCParalAloneList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Short_t nMCSkewAlone[MAXTRACKSPEREVENT], Short_t MCSkewAloneList[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool flaggo; Short_t i, jexp, exphit, iHit, enne[MAXTRACKSPEREVENT][nmaxSttHits]; Short_t emme; for(jexp=0; jexp 2. || (emme != daTrackFoundaTrackMC[jexp]) ) continue; if( !InclusionListStt[i]) continue; // escludo gli hits con multiple hits flaggo=true; for(exphit=0; exphit=3 && IVOLTE<20) cout<<"MatchMvdHitsToSttTracks, i= "<=4){ cout<<"n. punti nel fit "<=3 && IVOLTE <= 20){ sprintf(stringa, "/home/boca/panda/glpk/glpk-4.39/examples/glpsol --min -o soluztrack%dEvent%dstep%d GeneralParallelHitsConformeTraccia%dEvent%d.mcs", 0,IVOLTE,1,0,IVOLTE); } else { sprintf(stringa, "/home/boca/panda/glpk/glpk-4.39/examples/glpsol --min -o soluztrack%dEvent%dstep%d GeneralParallelHitsConformeTraccia%dEvent%d.mcs >& /dev/null", 0, IVOLTE,1,0,IVOLTE,1); } */ m1_result = final_values[0]; m2_result = final_values[1]; q1_result = final_values[2]; q2_result = final_values[3]; if(istampa>2) cout<<"Results : m1 = "< 1.e-10&& DriftRadius[ i ]>0.){ n++; ave += (S[ i ] - FInot)/Z[ i ]; avex += Z[ i ]; avey += (S[ i ] - FInot); } } if( n>0) { ave /=n; avex /=n; avey /=n; rotationangle = atan2(avey,avex); } else { rotationangle=PI/2.; } */ // use the trick of increasing the rotation angle by 10 degrees in order to obtain always a positive m // rotationangle -= PI/18.; rotationangle = PI/2.; cose = cos(rotationangle); sine = sin(rotationangle); nSttHits = nMvdHits = 0; for(i=0;i 1.e-10) { *emme=((*emme)*cose+sine)/(cose-(*emme)*sine); return 1; } else { // in this case the equation is 0 = U in the Conformal plane --> x=0 in the XY plane. return -99; } } //----------end of function PndSttMvdTracking::FitSZspace //------------------ begin function PndSttMvdTracking::RefitMvdStt void PndSttMvdTracking::RefitMvdStt( Short_t nCandHit, Short_t *ListCandHit, Short_t *ListCandHitType, Double_t info[][7], Double_t rotationangle, // this is between 0. and 2*PI Double_t tv[2], Short_t iexcl, Double_t *pAlfa, // output of the fit Double_t *pBeta, // output of the fit Double_t *pGamma,// set at zero always for now bool *status // fit status; true = successful ) { bool Type; Short_t i, iparallel; Short_t MAXIMUMHITSINFIT = 20; if(MAXIMUMHITSINFIT>nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack) MAXIMUMHITSINFIT = nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack; Short_t exitstatus; Double_t dist2, mindis, emme, factor, gamma, qu, ErrorStraw = 0.03, ErrorMvd = 0.01, Xconformal[nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], Yconformal[nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], DriftRadiusconformal[nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], ErrorDriftRadiusconformal[nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack]; *status= false; factor=3.; mindis=0.5; // trajectory_vertex[0]=trajectory_vertex[1]=0.; for(i=0, iparallel=0; i2){ cout<<"PndSttMvdTracking::RefitMvdStt, n. Hits (Mvd+Stt || ) = "<At(j); TVector3 position; pMvdMCPoint->Position(position); XMvdMCPoint=position.X(); YMvdMCPoint=position.Y(); ZMvdMCPoint=position.Z(); MCPointtoMCTrackID= pMvdMCPoint->GetTrackID(); if( !inclusionMCPoint[j]) continue; distance = (XMvdMCPoint-XMvdStrip[i])*(XMvdMCPoint-XMvdStrip[i])+ (YMvdMCPoint-YMvdStrip[i])*(YMvdMCPoint-YMvdStrip[i])+ (ZMvdMCPoint-ZMvdStrip[i])*(ZMvdMCPoint-ZMvdStrip[i]); if(istampa>2) cout<<"distanza**2 di Strip hit n. "< -1){ if( daTrackFoundaTrackMC[i] == FromPixeltoMCTrack[ ListPixelHitsinTrack[i][j] ] ){ index = i*nMvdPixelHit+nMvdPixelCommon[i]; *(MvdPixelCommonList+index) = ListPixelHitsinTrack[i][j]; nMvdPixelCommon[i]++; } else { index = i*nMvdPixelHit+nMvdPixelSpuriinTrack[i]; *(MvdPixelSpuriList+index) = ListPixelHitsinTrack[i][j]; nMvdPixelSpuriinTrack[i]++; } } } for(j=0;j -1){ if( daTrackFoundaTrackMC[i] == FromStriptoMCTrack[ ListStripHitsinTrack[i][j] ] ){ index = i*nMvdStripHit+nMvdStripCommon[i]; *(MvdStripCommonList+index) = ListStripHitsinTrack[i][j]; nMvdStripCommon[i]++; } else { index = i*nMvdStripHit+nMvdStripSpuriinTrack[i]; *(MvdStripSpuriList+index) = ListStripHitsinTrack[i][j]; nMvdStripSpuriinTrack[i]++; } } } } // end of for(i=0; i -1){ for(i=0;i -1){ for(i=0;i0){ // track must rotate clockwise looking into the beam. anglemax = FI0[i]; anglemin = Fifirst[i]; } else { anglemin = FI0[i]; anglemax = Fifirst[i]; } if(anglemax < anglemin) anglemax += 2.*PI; if(anglemax < anglemin) anglemax=anglemin; // this is just to be super-sure. } // end of if( Fifirst[i] < -99998. ) //-------------------- ngoodmix=0; nn[0]=0; for( imvdcand=0; imvdcandGetBranchId(fMvdPixelBranch)){ ncont++; angle = atan2( YMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]]-Oy[i], XMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]]-Ox[i] ); if(angle<0.) angle += 2.*PI; if( angle>anglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax) { dist=fabs( sqrt( (Ox[i]-XMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]])* (Ox[i]-XMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]]) +(Oy[i]-YMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]])* (Oy[i]-YMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]]))-R[i]); if(distGetBranchId(fMvdPixelBranch); Dist += dist; if( dist anglemin) } else {// at this point this is a Strip hit; already made sure // earlier in the code that there is no third possibility. ncont++; angle = atan2( YMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]]-Oy[i], XMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]]-Ox[i] ); if(angle<0.) angle += 2.*PI; if( angle>anglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]])* (Ox[i]-XMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]]) +(Oy[i]-YMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]])* (Oy[i]-YMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]])) -R[i]); if(distGetBranchId(fMvdStripBranch); Dist += dist; if( dist anglemin) } // end of if(ListHitTypeMvdTrackCand[imvdcand][jmvdhit] } // end of for( jmvdhit=0; jmvdhit0) { DIST[ngoodmix]=Dist/nn[ngoodmix]; ngoodmix++; //--------- stampaggi if(istampa>=3 ){cout<<"\tquesto Mvd candidato (n. ngoodmix = "<GetBranchId(fMvdPixelBranch)) { cout<<"\tPixel hit n. "<GetBranchId(fMvdStripBranch)){ cout<<"\tStrip hit n. "<anglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdPixel[ListMvdDSPixelHitNotTrackCand[jmvdhit]])* (Ox[i]-XMvdPixel[ListMvdDSPixelHitNotTrackCand[jmvdhit]]) +(Oy[i]-YMvdPixel[ListMvdDSPixelHitNotTrackCand[jmvdhit]])* (Oy[i]-YMvdPixel[ListMvdDSPixelHitNotTrackCand[jmvdhit]]) ) -R[i]); if(distGetBranchId(fMvdPixelBranch); DIST[ngoodmix] += dist; if( dist anglemin ) } // end of for( jmvdhit=0; jmvdhitanglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdStrip[ListMvdDSStripHitNotTrackCand[jmvdhit]])* (Ox[i]-XMvdStrip[ListMvdDSStripHitNotTrackCand[jmvdhit]]) +(Oy[i]-YMvdStrip[ListMvdDSStripHitNotTrackCand[jmvdhit]])* (Oy[i]-YMvdStrip[ListMvdDSStripHitNotTrackCand[jmvdhit]]) ) -R[i]); if(distGetBranchId(fMvdStripBranch); DIST[ngoodmix] += dist; if( dist anglemin ) } // end of for( jmvdhit=0; jmvdhit0) { DIST[ngoodmix] /= nn[ngoodmix]; ngoodmix++; //--------- stampaggi if(istampa>=3){cout<<"\tevento n. "<GetBranchId(fMvdPixelBranch)) { cout<<"\tDS Pixel hit n. "<GetBranchId(fMvdStripBranch)){ cout<<"\tDS Strip hit n. "<0) //-------- now the US (upstream) Mvd hits nn[ngoodmix]=0; DIST[ngoodmix] = 0.; nHighQuality[ngoodmix]=0; for( jmvdhit=0; jmvdhitanglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdPixel[ListMvdUSPixelHitNotTrackCand[jmvdhit]])* (Ox[i]-XMvdPixel[ListMvdUSPixelHitNotTrackCand[jmvdhit]]) +(Oy[i]-YMvdPixel[ListMvdUSPixelHitNotTrackCand[jmvdhit]])* (Oy[i]-YMvdPixel[ListMvdUSPixelHitNotTrackCand[jmvdhit]]) ) -R[i]); if(distGetBranchId(fMvdPixelBranch); DIST[ngoodmix] += dist; if( dist anglemin ) } // end of for( jmvdhit=0; jmvdhitanglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdStrip[ListMvdUSStripHitNotTrackCand[jmvdhit]])* (Ox[i]-XMvdStrip[ListMvdUSStripHitNotTrackCand[jmvdhit]]) +(Oy[i]-YMvdStrip[ListMvdUSStripHitNotTrackCand[jmvdhit]])* (Oy[i]-YMvdStrip[ListMvdUSStripHitNotTrackCand[jmvdhit]]) ) -R[i]); if(distGetBranchId(fMvdStripBranch); DIST[ngoodmix] += dist; if( dist anglemin ) } // end of for( jmvdhit=0; jmvdhit0) { DIST[ngoodmix] /= nn[ngoodmix]; ngoodmix++; //--------- stampaggi if(istampa>=3){cout<<"\tevento n. "<0) { for(int icc=0; iccGetBranchId(fMvdPixelBranch)) { cout<<"\tUS Pixel hit n. "<GetBranchId(fMvdStripBranch)){ cout<<"\tUS Strip hit n. "<0) //------- end of using the Mvd which are in no Mvd Track Candidate if(istampa>=3 ){cout<<"da PndSttMvdTracking : appena prima arbitration, IVOLTE = "<< IVOLTE<<", Stt track cand = "< MaxTurnofTracks && abs(nr2)>MaxTurnofTracks) { continue; } else if ( abs(nr2)>MaxTurnofTracks) { zeta0 = ZED[i]+DriftRadius[i]; nchosen=nrounds0; } else if ( abs(nrounds0) > MaxTurnofTracks ) { dista = ddd; zeta0 = ZED[i]-DriftRadius[i]; nchosen=nr2; } else { if( fabs(dista)>fabs(ddd) ) { dista = ddd; zeta0 = ZED[i]-DriftRadius[i]; nchosen=nr2; } else { zeta0 = ZED[i]+DriftRadius[i]; nchosen=nrounds0; } } error = DriftRadius[i]; if(nchosen>=0)SchosenSkew[ListSkewHitsinTrack[j]]=S[i]+nchosen*2.*PI; else SchosenSkew[ListSkewHitsinTrack[j]]=S[i]+(nchosen-1)*2.*PI; ZchosenSkew[ListSkewHitsinTrack[j]]=zeta0; //-----------------stampe. if(istampa>=3){ cout<<"in EliminateSpuriousSZ : insomma, dista prima della selezione = "<< dista<< ", ed e' da comparare\n\tcon 4*error = "<< 4.*error< SEMILENGTH_STRAIGHT && fabs(Zlast2 - ZCENTER_STRAIGHT) > SEMILENGTH_STRAIGHT ) continue; */ NAssociated++; } // end of for( ii=0; ii<2; ii++) } // for( iii=0; iii< NSkewhits; iii++) return NAssociated; } //----------end of function PndSttMvdTracking::AssociateSkewHitsToXYTrack //----------begin of function PndSttMvdTracking::CalculateCircleThru3Points bool PndSttMvdTracking::CalculateCircleThru3Points( Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3, Double_t *Ox, Double_t *Oy, Double_t *R) { // inspiration from http://mathworld.wolfram.com/Circle.html Double_t a, d, e, q1, q2, q3; a = x1*y2 + y1*x3 + x2*y3 - x3*y2 - x2*y1 - x1*y3; if( fabs(a)<1.e-10 ) return false; q1 = x1*x1 + y1*y1; q2 = x2*x2 + y2*y2; q3 = x3*x3 + y3*y3; d = -(q1*y2 + y1*q3 + q2*y3 - q3*y2 - q2*y1 - q1*y3); e = -(x1*q2 + q1*x3 + x2*q3 - x3*q2 - x2*q1 - x1*q3); *Ox = -0.5*d/a; *Oy = -0.5*e/a; *R = sqrt( (x1-*Ox)*(x1-*Ox) + (y1-*Oy)*(y1-*Oy) ); return true; } //----------end of function PndSttMvdTracking::CalculateCircleThru3Points //----------begin of function PndSttMvdTracking::Ordering_Loading_ListTrackCandHit void PndSttMvdTracking::Ordering_Loading_ListTrackCandHit( bool *keepit, Short_t FirstCandidate, Short_t LastCandidate, Double_t info[][7], Double_t * Ox, Double_t * Oy, Double_t * Rr, Double_t Trajectory_Start[MAXTRACKSPEREVENT][2], Short_t *CHARGE, Double_t SchosenSkew[][nmaxSttHits] ) { Short_t ncand; for(ncand=FirstCandidate; ncand< LastCandidate; ncand++){ // for small radius trajectory better the ordering with conformal. if( Rr[ncand]< RStrawDetectorMax/2.){ OrderingConformal_Loading_ListTrackCandHit( keepit, ncand, info, Ox, Oy, Rr, Trajectory_Start, CHARGE, SchosenSkew ); } else { // otherwise it is better distance from (0,0) method. OrderingR_Loading_ListTrackCandHit( keepit, ncand, info ); } } // end of for(ncand=FirstCandidate; ncand< LastCandidate; ncand++) return; } //----------end of function PndSttMvdTracking::Ordering_Loading_ListTrackCandHit //----------begin of function PndSttMvdTracking::OrderingR_Loading_ListTrackCandHit void PndSttMvdTracking::OrderingR_Loading_ListTrackCandHit( bool *keepit, Short_t ncand, Double_t info[][7] ) { Short_t i, j, ipar, iskew; // ordering all the hits belonging to the candidate track, by increasing R; // forming the new track with Mvd+Stt hits if(!keepit[ncand]) return; nTrackCandHit[ncand] =nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand]+ nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; Short_t tempmvdindex[nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] ], tempmvdtype[nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] ]; Int_t auxIndex[nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] ]; Double_t auxR[nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] ]; // adding the Mvd hits (Pixel and Strips) for(i=0; i< nMvdPixelHitsinTrack[ncand]; i++){ auxR[i] = XMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]* XMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]+ YMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]* YMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]; tempmvdindex[i]=ListMvdPixelHitsinTrack[ncand][i]; tempmvdtype[i]=0; auxIndex[i] = i; } for(i=0; i< nMvdStripHitsinTrack[ncand]; i++){ auxR[i+nMvdPixelHitsinTrack[ncand]] = XMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]* XMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]+ YMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]* YMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]; tempmvdindex[i+nMvdPixelHitsinTrack[ncand]]= ListMvdStripHitsinTrack[ncand][i]; tempmvdtype[i+nMvdPixelHitsinTrack[ncand]]=1; auxIndex[i+nMvdPixelHitsinTrack[ncand]]= i+nMvdPixelHitsinTrack[ncand]; } // ordering the Mvd Hits if( nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] >0){ Merge_Sort( nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand], auxR, auxIndex); // constructing the first part of the ordered new Track Candidate for(i=0; i< nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; i++){ ListTrackCandHit[ncand][i] = tempmvdindex[ auxIndex[i] ]; ListTrackCandHitType[ncand][i] = tempmvdtype[ auxIndex[i] ]; } } // end of if( nMvdPixelHitsinTrack[ncand]+ // construction of the second part of the ordered new Track Candidate Short_t tempmvdindex2[nSttParHitsinTrack[ncand]+ nSttSkewHitsinTrack[ncand] ], tempmvdtype2[nSttParHitsinTrack[ncand]+ nSttSkewHitsinTrack[ncand] ]; Int_t auxIndex2[nSttParHitsinTrack[ncand]+ nSttSkewHitsinTrack[ncand] ]; Double_t auxR2[nSttParHitsinTrack[ncand]+ nSttSkewHitsinTrack[ncand] ]; for(i=0; i0){ Merge_Sort( nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand], auxR2, auxIndex2); for(j=0,ipar=0,iskew=0; j< nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand];j++){ i = j+nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; ListTrackCandHit[ncand][i] = tempmvdindex2[ auxIndex2[j] ]; ListTrackCandHitType[ncand][i] = tempmvdtype2[ auxIndex2[j] ]; if( ListTrackCandHitType[ncand][i]==2) { ListSttParHitsinTrack[ncand][ipar]=tempmvdindex2[auxIndex2[j]]; ipar++; } else { ListSttSkewHitsinTrack[ncand][iskew]=tempmvdindex2[auxIndex2[j]]; iskew++; } } } // end of if( nSttParHitsinTrack[ncand]+ return; } //----------end of function PndSttMvdTracking::OrderingR_Loading_ListTrackCandHit //----------begin of function PndSttMvdTracking::OrderingConformal_Loading_ListTrackCandHit void PndSttMvdTracking::OrderingConformal_Loading_ListTrackCandHit( bool *keepit, Short_t ncand, Double_t info[][7], Double_t * Ox, Double_t * Oy, Double_t * Rr, Double_t Trajectory_Start[MAXTRACKSPEREVENT][2], Short_t *CHARGE, Double_t SchosenSkew[][nmaxSttHits] ) { Short_t i, j, ipar, iskew; // ordering all the hits belonging to the candidate track, by increasing R; // forming the new track with Mvd+Stt hits // arrays used to store temporarily the info of Mvd hits to be ordered. Int_t ListHits[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], ListHits2[nmaxSttHitsInTrack]; Double_t XY[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack][2], XY2[nmaxSttHitsInTrack][2]; // for(ncand=FirstCandidate; ncand< LastCandidate; ncand++){ if(!keepit[ncand]) return; nTrackCandHit[ncand] =nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand]+ nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; // adding the Mvd hits (Pixel and Strips) if( nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] >0){ for(i=0; i< nMvdPixelHitsinTrack[ncand]; i++){ XY[i][0] = XMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]; XY[i][1] = YMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]; ListHits[i] = ListMvdPixelHitsinTrack[ncand][i]; } for(i=0; i< nMvdStripHitsinTrack[ncand]; i++){ XY[i+nMvdPixelHitsinTrack[ncand]][0] = XMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]; XY[i+nMvdPixelHitsinTrack[ncand]][1] = YMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]; // to distinguish between Pixels and Strips, add a number // to the original Strip hit number. ListHits[i+nMvdPixelHitsinTrack[ncand]] = ListMvdStripHitsinTrack[ncand][i]+ (nmaxMvdPixelHits+nmaxMvdStripHits)*10 ; } // ordering the Mvd Hits OrderingUsingConformal( Ox[ncand], Oy[ncand], &Trajectory_Start[ncand][0], nMvdPixelHitsinTrack[ncand]+nMvdStripHitsinTrack[ncand], XY, // XY[*][0] = X position, XY[*][0] = Y position. CHARGE[ncand], // input ListHits // output ); // constructing the ordered new Track Candidate now for(i=0; i< nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; i++){ if(ListHits[i]<(nmaxMvdPixelHits+nmaxMvdStripHits)*10){//Pixel. ListTrackCandHit[ncand][i] = ListHits[i]; ListTrackCandHitType[ncand][i] = 0; } else { // Strip hits. ListTrackCandHit[ncand][i] = ListHits[i]- (nmaxMvdPixelHits+nmaxMvdStripHits)*10; ListTrackCandHitType[ncand][i] = 1; } } // end of for(i=0; i< nMvdPixelHitsinTrack[ncand]+ } // end of if( nMvdPixelHitsinTrack[ncand]+ // construction of the second part of the ordered new Track Candidate if( nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand] >0){ for(i=0; i=3) for(int ica=0; ica 0.9 * RminStrawDetector/2 and consequently Ox and Oy are not both 0. // The scheme for the ordering of the hit is as follows : // 1) order hits by increasing U or V of the conformal mapping; see Gianluigi's Logbook page 283; // 2) find the charge of the track by checking if it is closest to the center in XY // the first or the last of the ordered hits. // 3) in case, invert the ordering of U, V and ListHits such that the first hits in the // list are always those closer to the Trajectory_Start. // ordering of the hits aaa = atan2( oY-Traj_Sta[1], oX-Traj_Sta[0]); // atan2 defined between -PI and PI. // the following statement is necessary since for unknown reason the root interpreter // gives a weird error when using PI directly in the if statement below!!!!!!! I lost // 2 hours trying to figure this out! b1 = PI/4.; if((aaa>b1&&aaa<3.*b1) || (aaa>-3.*b1&&aaa<-b1)){//use U as ordering variable; //[case 1 or 3 Gianluigi's Logbook page 285]. for (j = 0; j< nHits; j++){ bbb = XY[j][0]-Traj_Sta[0]; ccc = XY[j][1]-Traj_Sta[1]; U[j]= bbb/(bbb*bbb+ccc*ccc); } Merge_Sort( nHits, U, ListHits); if((aaa>b1&&aaa<3.*b1)){ // case #1; if( Charge == -1){ // inverting the order of the hits. for(i=0;ib1&&aaa<3.*b1)) } else { // use V as ordering variable [case 2 or 4 Gianluigi's Logbook page 285]. for (j = 0; j< nHits; j++){ bbb = XY[j][0]-Traj_Sta[0]; ccc = XY[j][1]-Traj_Sta[1]; V[j]= ccc/(bbb*bbb+ccc*ccc); } Merge_Sort( nHits, V, ListHits); if((aaa<=-3.*b1 || aaa>=3.*b1)){ // case #2; if( Charge == -1){ // inverting the order of the hits. for(i=0;ib1&& .... return; } //----------end of function PndSttMvdTracking::OrderingUsingConformal //----------begin of function PndSttMvdTracking::TrackCleanup bool PndSttMvdTracking::TrackCleanup( Double_t GAP, Double_t Oxx, Double_t Oyy, Double_t Rr, Double_t KAPPA, Double_t FI0, Short_t Charge, Double_t Start[3], Short_t &nHitsPar, // n. hits parall Stt Short_t *ListHitsPar, Short_t &nHitsSkew, // n. hits parall Stt Short_t *ListHitsSkew, Double_t *auxS, Double_t info[][7], Double_t RStrawDetMin, Double_t ApotemaMaxInnerPar, Double_t ApotemaMinSkew, Double_t ApotemaMaxSkew, Double_t ApotemaMinOuterPar, Double_t RStrawDetMax ) { // this method does 3 things : // // 1) finds the entrance and exit points in the STT parallel and skew volumes of the current track; // 2) eliminates from the track hit list possible spurious hits that are not encompassed // by the entrance and exit point; // 3) eliminates the tracks if the hit sequence is not continuous enough. bool ConsiderLastHit; Short_t flagInnerStt, flagOuterStt; Short_t ihit, nInnerHits, nOuterHits, nIntersections[2], ListInnerHits[nmaxSttHits], ListOuterHits[nmaxSttHits]; Double_t FiLimitAdmissible, r, Xcross[2], Ycross[2], XintersectionList[12][2], // second index =0 --> inner Hexagon, =1 --> outer. YintersectionList[12][2]; // first index : all the possible intersections // (up to 12 intersections). //------------------------ // calculation of the Maximum FI angle possible (if it is a +ve charge) of the Minimum // for this track, taking into account // that the maximum possible Z of a hit is ZCENTER_STRAIGHT + SEMILENGTH_STRAIGHT; the minimum // Z of a hit is ZCENTER_STRAIGHT - SEMILENGTH_STRAIGHT. if(Charge<0){ if( KAPPA>0.){ FiLimitAdmissible = FI0 + KAPPA*(ZCENTER_STRAIGHT + SEMILENGTH_STRAIGHT) ; } else { FiLimitAdmissible = FI0 + KAPPA*(ZCENTER_STRAIGHT - SEMILENGTH_STRAIGHT) ; } } else { if( KAPPA>0.){ FiLimitAdmissible = FI0 + KAPPA*(ZCENTER_STRAIGHT - SEMILENGTH_STRAIGHT) ; } else { FiLimitAdmissible = FI0 + KAPPA*(ZCENTER_STRAIGHT + SEMILENGTH_STRAIGHT) ; } } // end of if(Charge<0) //----------------------------------------------------------------------------------------------------- // parallel cleanup. //----------------stampe if(istampa>=2){ cout<<" IVOLTE = "<1) cout<<"\tentra in SttParalCleanup\n"; // if(nHitsPar>0 && !SttParalCleanup( if(!SttParalCleanup( GAP, Oxx, Oyy, Rr, Charge, Start, FI0, FiLimitAdmissible, nHitsPar, // it doesn't get modify for now. ListHitsPar, // input only for now. info, RStrawDetMin, ApotemaMaxInnerPar, ApotemaMinOuterPar, RStrawDetMax ) ){ if(istampa>1) cout<<"uscito da SttParalCleanup : false\n"; return false; } if(istampa>1) cout<<"uscito da : SttParalCleanup true\n"; //---------------------------------------------------------------------------- // skew cleanup. if(istampa>1) cout<<"\tentra in SttSkewCleanup\n"; if ( ! (SttSkewCleanup( GAP, Oxx, Oyy, Rr, Charge, Start, // strarting point of trajectory. FI0, FiLimitAdmissible, nHitsSkew, // it doesn't get modify for now. ListHitsSkew, // it doesn't get modify for now. auxS, info, ApotemaMinSkew, //distance hexagon side from (0,0) ApotemaMaxSkew, //distance hexagon side from (0,0) 3., // cut distance 1 // max number of failures allowed. ) ) ) { if(istampa>1) cout<<"uscito da SttSkewCleanup false\n"; return false; } if(istampa>1) cout<<"uscito da SttSkewCleanup true\n"; return true; }; //----------end of function PndSttMvdTracking::TrackCleanup //----------begin of function PndSttMvdTracking::SttParalCleanup bool PndSttMvdTracking::SttParalCleanup( Double_t GAP, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t FI0, Double_t FiLimitAdmissible, Short_t nHits, Short_t *Listofhits, Double_t info[][7], Double_t RStrawDetMin, Double_t ApotemaInnerParMax, Double_t ApotemaOuterParMin, Double_t RStrawDetMax ) { // this method does 3 things : // // 1) finds the entrance and exit points in the STT parallel volumes of the current track; // 2) eliminates from the track hit list possible spurious hits that are not encompassed // by the entrance and exit point; // 3) eliminates the tracks if the hit sequence is not continuous enough. bool flaggo; Short_t flagInnerSttR, flagOuterSttR, flagInnerSttL, flagOuterSttL, flagOutStt; Short_t enne, i, ihit, ipurged, islack, nintersections, nnn, nInnerHits, nInnerHitsLeft, nInnerHitsRight, nOuter, nOuterHits, nOuterHitsLeft, nOuterHitsRight, nIntersections[2], ListHits[nHits], ListInnerHits[nHits], ListInnerHitsLeft[nHits], ListInnerHitsRight[nHits], ListOuterHits[nHits], ListOuterHitsLeft[nHits], ListOuterHitsRight[nHits]; Double_t epsilonTheta, fi, r, aux[2], Xcross[2], Ycross[2], XcrossL[2], YcrossL[2], XcrossR[2], YcrossR[2], XcrossOut[2], YcrossOut[2], XintersectionList[7], // there is also the last boundary FiLimitAdmissible YintersectionList[7]; // take into account and the two possible // intersections with the external circle. islack=1; // uncertainty allowed in the # of straws that should be hit in a given part // of the Stt detector. //------------------------ // elimination of hits outside the physical FI range (FiLimitAdmissible) due to finite length of // Straws. if(istampa>1) { cout<<"SttParalCleanup, evento n. "< FI0){ if( fi>FiLimitAdmissible+epsilonTheta) continue; } else { fi += 2.*PI; if( fi >FiLimitAdmissible+epsilonTheta ) continue; } // end of if( fi > FI0) } else { // continuation of if(Charge <0) if( fi > FI0){ fi -= 2.*PI; } // end of if( fi > FI0) if (fi < FiLimitAdmissible-epsilonTheta) continue; } // end of if(Charge <0) ListHits[ipurged]=Listofhits[i]; ipurged++; } // end of for(i=0, ipurged=0; i< nHits; i++) nHits = ipurged; if(istampa>1) { cout<<"SttParalCleanup, evento n. "<1) { cout<<"SttParalCleanup, evento n. "< track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; } // end of if(nHits==0) // first of all, find possible intersection points with outer circle encompassing // the Stt system. flagOutStt = FindIntersectionsOuterCircle( Oxx, Oyy, Rr, RStrawDetectorMax, XcrossOut, YcrossOut ); // intersection with Inner Section. flagInnerSttL=FindTrackEntranceExitbiHexagonLeft( GAP, Oxx, Oyy, Rr, Charge, Start, RStrawDetMin, ApotemaInnerParMax, XcrossL, YcrossL ); // find the entrance and exit of the track in the Inner Right Parallel Straw region. // This region is bounded by two Hexagons, and it has the target gap in the middle. flagInnerSttR=FindTrackEntranceExitbiHexagonRight( GAP, Oxx, Oyy, Rr, Charge, Start, RStrawDetMin, ApotemaInnerParMax, XcrossR, YcrossR ); if(istampa>1) { cout<<"SttParalCleanup, evento n. "<1) { cout<<"SttParalCleanup, evento n. "< cut. islack // uncertainty allowed as far as the n. of hits that should be present. ) ){ return false; } if(istampa>1) cout<<"uscito da BadTrack_ParStt.\n"; //----------------------------------------------------- } // end of if(flaggo) } // end of if( flagInnerSttL == -1 && flagInnerSttR == -1 ) //outer: ; islack=1; // reset the extra uncertainty in the # Stt. //------------ Outer Parallel Stt hits section. // find the entrance and exit of the track in the Outer Parallel Straw region, Left side. // This region is bounded by a Hexagon (inner), a Circle (outer) and it has // the target gap in the middle. // Returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections. flagOuterSttL=FindTrackEntranceExitHexagonCircleLeft( Oxx, Oyy, Rr, Charge, Start, ApotemaMinOuterParStraw, RStrawDetMax, GAP, XcrossL, YcrossL ); //------------ // find the entrance and exit of the track in the Outer Parallel Straw region, Right side. // This region is bounded by a Hexagon (inner), a Circle (outer) and it has // the target gap in the middle. flagOuterSttR=FindTrackEntranceExitHexagonCircleRight( Oxx, Oyy, Rr, Charge, Start, ApotemaMinOuterParStraw, RStrawDetMax, GAP, XcrossR, YcrossR ); if(istampa>1) { cout<<"SttParalCleanup, evento n. "<1) { cout<<"in SttParalCleanup Outer, caso traccia entra in L and R outer. Dopo scelta, flagOuterSttR " <1) { cout<<"in SttParalCleanup Outer, nhit considerati "<< nnn <1) { cout<<"in SttParalCleanup Outer, caso in cui FiLimitAdmissible = "<< FiLimitAdmissible <<" conta!" <1) { cout<<"in SttParalCleanup Outer, prima di ChooseEntranceExitbis, nintersections " << nintersections<<" e loro lista :"<=2&&IVOLTE<20){ cout<<"SttParalCleanup, OUTER, caso R || L true, IVOLTE = "< cut. islack // uncertainty allowed as far as the n. of hits that should be present. ) ) return false; // } // end of if(nOuterHits==0) //---------------------------------------------------------------------------- // finito: ; // if the code comes here it means that the track is acceptable. // nHits = nOuterHits+nInnerHits; return true; }; //----------end of function PndSttMvdTracking::SttParalCleanup //----------begin of function PndSttMvdTracking::SttSkewCleanup bool PndSttMvdTracking::SttSkewCleanup( Double_t GAP, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t FI0, Double_t FiLimitAdmissible, Short_t nHits, Short_t *Listofhits, Double_t *S, Double_t info[][7], Double_t RminStrawSkew, Double_t RmaxStrawSkew, Double_t cut, // cut distance (in cm). Short_t maxnum // max number allowed of failures to pass the cut. ) { bool ConsiderLastHit; Short_t flagSttL, flagSttR, flagOutStt; Short_t i, ipurged, ibad, islack, nHitsLeft, nHitsRight, nintersections, ninside, nnn, nIntersections[2], ListHits[nHits], ListHitsRight[nHits], ListHitsLeft[nHits]; Double_t cut2, epsilonTheta, fi, FiStart, length, r, Sprevious, aux[2], Distance[nmaxSttHits+1], Xcross[2], Ycross[2], XcrossL[2], YcrossL[2], XcrossR[2], YcrossR[2], XcrossOut[2], YcrossOut[2], XintersectionList[5], // second index =0 --> inner Hexagon, =1 --> outer. YintersectionList[5]; // first index : all the possible intersections // (up to 12 intersections). cut2=cut*cut; islack=1;// uncertainty allowed as far as // the n. of hits that should be present in a given section of the Stt track. //------------------------ // elimination of hits outside the physical FI range (FiLimitAdmissible) due to finite length of // Straws. epsilonTheta = STRAWRADIUS/Rr; // some extra slac for being conservative. if(istampa>1&&IVOLTE<20) cout<<"\n\nevt "<1&&IVOLTE<20)cout<<"\thit // n. "< FI0){ if( fi>FiLimitAdmissible+epsilonTheta) continue; } else { fi += 2.*PI; if( fi > FiLimitAdmissible+epsilonTheta ) continue; } // end of if( fi > FI0) } else { // continuation of if(Charge <0) if( fi > FI0){ fi -= 2.*PI; } // end of if( fi > FI0) if (fi < FiLimitAdmissible-epsilonTheta) continue; } // end of if(Charge <0) if(istampa>1&&IVOLTE<20)cout<<"in SttSkewCleanup : hit preso!"<1&&IVOLTE<20)cout<<"in SttSkewCleanup : hit skew prima di purga = " <1&&IVOLTE<20)cout<<"in SttSkewCleanup : n. hit skew Left = " <1)cout<<"in SttSkewCleanup : flagLeft (-1,0,1) = "<1)cout<<"in SttSkewCleanup : distanza entrata-uscita<4*STRAWRADIUS,flagSttR set at -1!\n"; } if( flagSttL == 0 && (XcrossL[0]-XcrossL[1])*(XcrossL[0]-XcrossL[1])+ (YcrossL[0]-YcrossL[1])*(YcrossL[0]-YcrossL[1]) < 16.*STRAWRADIUS*STRAWRADIUS ){ flagSttR=-1; if(istampa>1)cout<<"in SttSkewCleanup : distanza entrata-uscita<4*STRAWRADIUS,flagSttL set at -1!\n"; } if (flagSttR != 0 && flagSttL != 0 ) { //nHits=0; if(istampa>1)cout<<"in SttSkewCleanup : flagSttR = "<=2&&IVOLTE<20){ cout<<"in SttSkewCleanup, IVOLTE = "<1)cout<<"in SttSkewCleanup :hit n. "<< ListHits[i] <<" is NOT inside the arc between entrance and exit; hit excluded!\n"; } ninside++; Distance[i] = 2.*Rr*Rr*(1.-cos(S[i]-Sprevious)); // this is the usual // distance**2 formula: (x1-x2)**2+(y1-y2)**2; // it is already 'protected' against S[i] jumps // around 2PI/0. Sprevious = S[i]; if(Distance[i]<0.) Distance[i]=0.; // rounding errors protection. if(istampa>=2)cout<<"in SttSkewCleanup, Hit n. "<< ListHits[i]<<" has Distance " <cut2){ if(Distance[i]>16.*cut2){ if(istampa>=2)cout<<"in SttSkewCleanup, Hit n. "<< ListHits[i]<<" has Distance " <4.*cut [="<=2)cout<<"in SttSkewCleanup, Hit n. "<< ListHits[i]<<" has Distance " < cut [="<1){ cout<<"in SttSkewCleanup, n. Hits inside = "<=2)cout<<"in SttSkewCleanup, last Hit n. "<< ListHits[nHits-1]<<" has Distance from boundary " <cut2 ){ if( Distance[nHits]>16.*cut2){ if(istampa>=2)cout<<"in SttSkewCleanup, last Hit n. "<< ListHits[nHits-1]<<" has Distance from boundary " <4 .*cut [="<=2)cout<<"in SttSkewCleanup, last Hit n. "<< ListHits[nHits-1]<<" has Distance from boundary " < cut [="< maxnum){ if(istampa>=2&&IVOLTE<20)cout<<"in SttSkewCleanup, reject this track because ibad = "<< ibad <<" and it is > maxnum [="<1) {cout<<"in BadTrack_ParStt : Xingresso "<1) {cout<<"in BadTrack_ParStt :hit || n. "<1)cout<<"in BadTrack_ParStt, Stt || hit n. (original notation) "<< ListHits[nHits-1]<<", Distance to boundary = "<cut2 ){ if( Distance[nHits]>16.*cut2){ return true; } ibad++; } } // end of if( IsInsideArc if(istampa>1)cout<<"in BadTrack_ParStt, ibad "< length_segmentq || (x-P2x)*(x-P2x)+(y-P2y)*(y-P2y) > length_segmentq ) continue; status = true; XintersectionList[*Nintersections] = x; YintersectionList[*Nintersections] = y; (*Nintersections)++; } // end of for(ipossibility=-1; ... return status; } //----------end of function PndSttMvdTracking::IntersectionCircle_Segment //----------begin of function PndSttMvdTracking::IntersectionsWithGapSemicircle Short_t PndSttMvdTracking::IntersectionsWithGapSemicircle( Double_t Oxx, Double_t Oyy, Double_t Rr, Double_t GAP, bool left, // if true --> Left semicircle; false --> Right semicircle. Double_t Rma, Double_t *XintersectionList, Double_t *YintersectionList ) { Short_t nIntersectionsCircle; Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; nIntersectionsCircle=0; aaa = sqrt(Oxx*Oxx+Oyy*Oyy); // preliminary condition for having intersections between trajectory and Circle. if( !( aaa >= Rr + Rma || Rr >= aaa + Rma) && !( aaa + Rr <= Rma || aaa - Rr >= Rma ) ) { // now the calculation FI0 = atan2(-Oyy,-Oxx); cosFi = (aaa*aaa + Rr*Rr - Rma*Rma)/(2.*Rr*aaa); if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.; Fi = acos(cosFi); // (x1, y1) and (x2,y2) are the intersections between the trajectory // and the circle, in the laboratory reference frame. x1 = Oxx+Rr*cos(FI0 - Fi); y1 = Oyy+Rr*sin(FI0 - Fi); theta1 = atan2(y1,x1); // in this way theta1 is between -PI and PI. x2 = Oxx+Rr*cos(FI0 + Fi); y2 = Oyy+Rr*sin(FI0 + Fi); theta2 = atan2(y2,x2); // in this way theta2 is between -PI and PI. // Theta1, Theta2 = angle of the edges of the outer circle + Gap in the laboratory frame. // Theta1, Theta2 are angles between -PI/2 and +PI/2. if(!left){ // Right (looking into the beam) Semicircle. Theta2 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.); Theta1 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.); if( Theta1<= theta1 && theta1 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x1; YintersectionList[nIntersectionsCircle]=y1; nIntersectionsCircle++; } if( Theta1<= theta2 && theta2 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x2; YintersectionList[nIntersectionsCircle]=y2; nIntersectionsCircle++; } } else { // Left (looking into the beam) Semicircle. Theta2 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.); Theta1 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.); if( Theta1<= theta1 || theta1 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x1; YintersectionList[nIntersectionsCircle]=y1; nIntersectionsCircle++; } if( Theta1<= theta2 || theta2 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x2; YintersectionList[nIntersectionsCircle]=y2; nIntersectionsCircle++; } } } // end of if (!( a >= Rr + Rma || ..... //---------------------------- end of calculation intersection with outer circle. return nIntersectionsCircle; } //----------end of function PndSttMvdTracking::IntersectionsWithGapSemicircle //----------star of function PndSttMvdTracking::IsInternal bool PndSttMvdTracking::IsInternal( Double_t Px, // point Double_t Py, Double_t Xtraslation, Double_t Ytraslation, Double_t Theta ) { // for explanations see Gianluigi's logbook on page 278-280. if( (Xtraslation-Px)*sin(Theta) + (Py-Ytraslation)*cos(Theta) >= 0. ) return true; else return false; } //----------end of function PndSttMvdTracking::IsInternal //----------star of function PndSttMvdTracking::ChooseEntranceExitbis void PndSttMvdTracking::ChooseEntranceExitbis( Double_t Oxx, Double_t Oyy, Short_t Charge, Double_t FiStart, Short_t nIntersections, Double_t *XintersectionList, // second index =1 -->inner polygon; Double_t *YintersectionList, // second index =2 -->outer polygon. Double_t Xcross[2], // output Double_t Ycross[2] // output ) { Short_t i, j; // this method works under the hypothesis that there are at least 2 intersections. if (nIntersections<2) return; if(Charge > 0) { Int_t auxIndex[100]; Double_t fi[100]; for( i=0;i FiStart) fi[i] -= 2.*PI; if( fi[i] > FiStart) fi[i] = FiStart; auxIndex[i]=i; } // end of for( i=0;i2.*ApotemaMaxInnerParStraw/sqrt(3.) ){ // outer Parallel hit. ListOuterHits[ *nOuterHits ] = ListHits[ihit]; (*nOuterHits)++; if(info[ListHits[ihit]][0]<0.){ ListOuterHitsLeft[ *nOuterHitsLeft ] = ListHits[ihit]; (*nOuterHitsLeft)++; } else { ListOuterHitsRight[ *nOuterHitsRight ] = ListHits[ihit]; (*nOuterHitsRight)++; } }else{ ListInnerHits[ *nInnerHits ] = ListHits[ihit]; (*nInnerHits)++; if(info[ListHits[ihit]][0]<0.){ ListInnerHitsLeft[ *nInnerHitsLeft ] = ListHits[ihit]; (*nInnerHitsLeft)++; } else { ListInnerHitsRight[ *nInnerHitsRight ] = ListHits[ihit]; (*nInnerHitsRight)++; } } } } //----------end of function PndSttMvdTracking::SeparateInnerOuterParallel //----------begin of function PndSttMvdTracking::FindTrackEntranceExitbiHexagonLeft Short_t PndSttMvdTracking::FindTrackEntranceExitbiHexagonLeft( Double_t vgap, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t ApotemaMax, Double_t Xcross[2], Double_t Ycross[2] ) { Short_t flag; Short_t nIntersections; Double_t FiStart, XintersectionList[16], // all the possible intersections YintersectionList[16]; // (up to 16 intersections). // The following is the form of the Left (looking from downstream into the beam) biHexagon // geometrical shape considered in this method : // /* /| / | / | / / / / / / / / | | | | | | | | \ \ \ \ \ \ \ \ \ | \ | \| */ // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; flag=IntersectionsWithClosedbiHexagonLeft( vgap, Oxx, Oyy, Rr, ApotemaMin, // Apotema of the inner Hexagon, ApotemaMax,// Apotema of the outer Hexagon. &nIntersections, XintersectionList, // XintersectionList[..][0] --> inner polygon, // XintersectionList[..][1] --> outer polygon. YintersectionList ); // IMPORTANT : // this is true because here it is assumed that the track comes from (0,0,0) // otherwise the code must be changed! if (!(flag == 0)) return flag; if( nIntersections<2) return -1; //------- among all possible intersection find the entrance point (Xcross[0], Ycross[0]) // and the exit point (Xcross[1], Ycross[1]) of this track. //-------- the starting point of the track. FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); if(FiStart<0.) FiStart+= 2.*PI; if(FiStart<0.) FiStart =0.; // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. // so at this point, the intersections are at least 2. ChooseEntranceExitbis( Oxx, Oyy, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); //---------------------- return flag; } //----------end of function PndSttMvdTracking::FindTrackEntranceExitbiHexagonLeft //----------begin of function PndSttMvdTracking::FindTrackEntranceExitbiHexagonRight Short_t PndSttMvdTracking::FindTrackEntranceExitbiHexagonRight( Double_t vgap, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t ApotemaMax, Double_t Xcross[2], Double_t Ycross[2] ) { Short_t flag; Short_t nIntersections; Double_t FiStart, XintersectionList[16], // all the possible intersections YintersectionList[16]; // (up to 16 intersections). // The following is the form of the Left (looking from downstream into the beam) biHexagon // geometrical shape considered in this method : // /* |\ | \ | \ \ \ \ \ | | | | / / / / / / | / | / |/ */ // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; flag=IntersectionsWithClosedbiHexagonRight( vgap, Oxx, Oyy, Rr, ApotemaMin, // Apotema of the inner Hexagon, ApotemaMax,// Apotema of the outer Hexagon. &nIntersections, XintersectionList, YintersectionList ); // IMPORTANT : // this is true because here it is assumed that the track comes from (0,0,0) // otherwise the code must be changed! if (!(flag == 0)) return flag; if( nIntersections<2) return -1; //------- among all possible intersection find the entrance point (Xcross[0], Ycross[0]) // and the exit point (Xcross[1], Ycross[1]) of this track. //-------- the starting point of the track. FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); if(FiStart<0.) FiStart+= 2.*PI; if(FiStart<0.) FiStart =0.; // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. // so at this point, the intersections are at least 2. ChooseEntranceExitbis( Oxx, Oyy, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); //---------------------- return flag; } //----------end of function PndSttMvdTracking::FindTrackEntranceExitbiHexagonRight //----------begin of function PndSttMvdTracking::FindTrackEntranceExitHexagonCircleLeft Short_t PndSttMvdTracking::FindTrackEntranceExitHexagonCircleLeft( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t Rma, // outer radius of the Circle. Double_t GAP, Double_t Xcross[2], Double_t Ycross[2] ) { // This methods finds the intersections between a trajectory coming from (0,0) and // parameters Oxx, Oyy, Rr with the closed geometrical figure (in XY) formed by the STT // external Left semicircle and the Outer STT parallel Left straw semi-Hexagon + Gap for // the pellet target target. // It returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections. Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; //------------------ Short_t nIntersectionsCircle, nIntersections; Double_t FiStart, XintersectionList[12], YintersectionList[12]; // all the possible intersections (up to 12 intersections). // finding all possible intersections with inner parallel straw region. Double_t Side_x[] = { -GAP/2., -GAP/2. , -ApotemaMin, -ApotemaMin, -GAP/2., -GAP/2. }, Side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-0.5*GAP)/sqrt(3.), ApotemaMin/sqrt(3.), -ApotemaMin/sqrt(3.), -(2.*ApotemaMin-0.5*GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)}, a[] = {1., -1./sqrt(3.), 1., 1./sqrt(3.), 1.}, b[] = {0., 1., 0., 1., 0.}, c[] = {GAP/2., -2.*ApotemaMin/sqrt(3.),ApotemaMin, 2.*ApotemaMin/sqrt(3.), GAP/2.}; nIntersections=IntersectionsWithOpenPolygon( Oxx, Oyy, Rr, 5, // n. Sides of open Polygon. a, // coefficient of formula : aX + bY + c = 0 defining the Polygon sides. b, c, Side_x, // X,Y coordinate of the Sides vertices (in sequence, following Side_y, // the Polygon along. XintersectionList, // XintersectionList YintersectionList // YintersectionList. ); //------------------------------------------------------------------------- // finding intersections of trajectory [assumed to originate from (0,0) ] // with outer semicircle, the Left part. nIntersectionsCircle=IntersectionsWithGapSemicircle( Oxx, // input from trajectory Oyy, // input from trajectory Rr, // input from trajectory GAP, // input, vertical gap in XY plane of STT detector. true, // true --> Left semicircle, false --> Right semicircle. Rma, // radius of external Stt containment. &XintersectionList[nIntersections], // output, X list of intersections (maximal 2). &YintersectionList[nIntersections] ); nIntersections += nIntersectionsCircle; //-------- the starting point of the track. if(nIntersections<2) return -1; FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); if(FiStart<0.) FiStart+= 2.*PI; if(FiStart<0.) FiStart =0.; // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. ChooseEntranceExitbis( Oxx, Oyy, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); return 0; } //----------end of function PndSttMvdTracking::FindTrackEntranceExitHexagonCircleLeft //----------begin of function PndSttMvdTracking::FindTrackEntranceExitHexagonCircleRight Short_t PndSttMvdTracking::FindTrackEntranceExitHexagonCircleRight( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t Rma, // outer radius of the Circle. Double_t GAP, Double_t Xcross[2], Double_t Ycross[2] ) { // This methods finds the intersections between a trajectory coming from (0,0) and // parameters Oxx, Oyy, Rr with the closed geometrical figure (in XY) formed by the STT // external Left semicircle and the Outer STT parallel Left straw semi-Hexagon + Gap for // the pellet target target. // It returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections. Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; //------------------ Short_t nIntersectionsCircle, nIntersections; Double_t FiStart, XintersectionList[12], YintersectionList[12]; // all the possible intersections (up to 10 intersections). // finding all possible intersections with inner parallel straw region. Double_t Side_x[] = { GAP/2., GAP/2. , ApotemaMin, ApotemaMin, GAP/2., GAP/2. }, Side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-GAP)/sqrt(3.), ApotemaMin/sqrt(3.), -ApotemaMin/sqrt(3.), -(2.*ApotemaMin-GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)}, a[] = {1., 1./sqrt(3.), 1., -1./sqrt(3.), 1.}, b[] = {0., 1., 0., 1., 0.}, c[] = {-GAP/2.,-2.*ApotemaMin/sqrt(3.),-ApotemaMin, 2.*ApotemaMin/sqrt(3.), -GAP/2.}; nIntersections=IntersectionsWithOpenPolygon( Oxx, Oyy, Rr, 5, // n. Sides of open Polygon. a, // coefficient of formula : aX + bY + c = 0 defining the Polygon sides. b, c, Side_x, // X,Y coordinate of the Sides vertices (in sequence, following Side_y, // the Polygon along. XintersectionList, // XintersectionList YintersectionList // YintersectionList. ); //------------------------------------------------------------------------- // finding intersections of trajectory [assumed to originate from (0,0) ] // with outer semicircle, the Left part. nIntersectionsCircle=IntersectionsWithGapSemicircle( Oxx, // input from trajectory Oyy, // input from trajectory Rr, // input from trajectory GAP, // input, vertical gap in XY plane of STT detector. false, // true --> Left semicircle, false --> Right semicircle. Rma, // radius of external Stt containment. &XintersectionList[nIntersections], // output, X list of intersections (maximal 2). &YintersectionList[nIntersections] ); nIntersections += nIntersectionsCircle; //-------- the starting point of the track. if(nIntersections<2) return -1; FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); if(FiStart<0.) FiStart+= 2.*PI; if(FiStart<0.) FiStart =0.; // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. ChooseEntranceExitbis( Oxx, Oyy, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); return 0; } //----------end of function PndSttMvdTracking::FindTrackEntranceExitHexagonCircleRight //----------begin of function PndSttMvdTracking::FindIntersectionsOuterCircle Short_t PndSttMvdTracking::FindIntersectionsOuterCircle( Double_t oX, Double_t oY, Double_t R, Double_t Rma, Double_t Xcross[2], Double_t Ycross[2] ) { // return -1 --> non intersection; // return 0 --> 2 intersections. Double_t a, cosFi, Fi, FI0; a = sqrt(oX*oX+oY*oY); // case with no intersections. if( a >= R + Rma || R >= a + Rma || a + R <= Rma) return -1; FI0 = atan2(-oY,-oX); cosFi = (a*a + R*R - Rma*Rma)/(2.*R*a); if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.; Fi = acos(cosFi); Xcross[0] = oX + R*cos(FI0+Fi); Ycross[0] = oY + R*sin(FI0+Fi); Xcross[1] = oX + R*cos(FI0-Fi); Ycross[1] = oY + R*sin(FI0-Fi); return 0; } //----------end of function PndSttMvdTracking::FindIntersectionsOuterCircle //----------begin of function PndSttMvdTracking::IsInsideCircle bool PndSttMvdTracking::IsInsideArc( Double_t Oxx, Double_t Oyy, Short_t Charge, Double_t Xcross[2], Double_t Ycross[2], Double_t f // f should be between 0 and 2PI. ) { Double_t f1, f2; // Xcross[0],Ycross[0] is the point of entrance. f1 = atan2(Ycross[0]-Oyy, Xcross[0]-Oxx); if(f1<0.) f1+= 2.*PI; if(f1<0.) f1= 0.; f2 = atan2(Ycross[1]-Oyy, Xcross[1]-Oxx); if(f2<0.) f2+= 2.*PI; if(f2<0.) f2= 0.; if(Charge<0){ if(f1 > f2 ) f2 +=2.*PI; if(f1 > f2 ) f2 = f1; if( f>f1){ if(f 0. if(f1 < f2 ) f1 +=2.*PI; if(f1 < f2 ) f1 = f2; if( f>f2){ if(f0.) { XintersectionList[0] = gap; YintersectionList[0] = sqrt(delta) + Oyy; XintersectionList[1] = gap; YintersectionList[1] = -sqrt(delta) + Oyy; nintersections+=2; } // intersections of trajectory with plane X = -gap; delta = Rr*Rr - (-gap-Oxx)*(-gap-Oxx); if(delta>0.) { XintersectionList[nintersections] = -gap; YintersectionList[nintersections] = sqrt(delta) + Oyy; XintersectionList[nintersections+1] = -gap; YintersectionList[nintersections+1] = -sqrt(delta) + Oyy; nintersections+=2; } // intersections of trajectory with plane Z = gap; XintersectionList[nintersections] = Oxx + Rr*cos(fi0+kappa*gap); YintersectionList[nintersections] = Oyy + Rr*sin(fi0+kappa*gap); nintersections++; // intersections of trajectory with plane Z = -gap; XintersectionList[nintersections] = Oxx + Rr*cos(fi0-kappa*gap); YintersectionList[nintersections] = Oyy + Rr*sin(fi0-kappa*gap); nintersections++; ChooseEntranceExitbis( Oxx, Oyy, Charge, fi0, nintersections,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); // here I suppose that the Mvd system CONSERVATIVELY ends at Rmax=5 cm. if( fabs(Ycross[0]) > 4. ) return true; else return false; } //----------end of function PndSttMvdTracking::IsInTargetPipe //----------begin of function PndSttMvdTracking::CalculateArcLength Double_t PndSttMvdTracking::CalculateArcLength( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t charge, Double_t Xcross[2], // entrance-exit point Double_t Ycross[2] // entrance-exit point ) { Double_t dis, theta1, theta2; theta1 = atan2(Ycross[0]-Oyy, Xcross[0]- Oxx); theta2 = atan2(Ycross[1]-Oyy, Xcross[1]- Oxx); if(charge>0){ // the rotation was clockwise. dis = theta1-theta2; } else { // the rotation was counterclockwise. dis = theta2-theta1; } if(dis<0.) dis += 2.*PI; if(dis<0.) dis =0.; dis *= Rr; return dis; } //----------end of function PndSttMvdTracking::CalculateArcLength //----------begin of function PndSttMvdTracking::FindCharge void PndSttMvdTracking::FindCharge( Double_t oX, Double_t oY, Short_t nParallelHits, Double_t *X, Double_t *Y, Short_t * Charge ) { Short_t ihit, nleft, nright; Double_t cross, disq, minl, minr; // this methods works with the hypothesis that this track comes // from (0,0) for(ihit=0, nleft=0, nright=0, minr = 9999999., minl = 9999999.; ihit0 hits stays 'on the left' (which means clockwise to go from the origin // to the hit following the smaller path) otherwise it stays 'on the right'. if (cross>0.) { disq = X[ihit]*X[ihit]+Y[ihit]*Y[ihit]; nleft++; } else { nright++; } } // end of for(ihit=0, nleft=0, nright=0;.... if( nright> nleft) { *Charge = -1; } else if ( nleft > nright) { *Charge = 1; } else { // then choose according the closest hit ti the center if( minr < minl ) *Charge = -1; else *Charge = 1; } } //----------end of function PndSttMvdTracking::FindCharge //----------begin of function PndSttMvdTracking::AssociateSciTilHit Short_t PndSttMvdTracking::AssociateSciTilHit( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t *List, // output Double_t *esse // output ) { bool intersect; Short_t igoodScit, iScitHit, Nint; Double_t distance, QQ, sqrtRR, SIGN, XintersectionList[2], YintersectionList[2]; igoodScit=0; for(iScitHit=0; iScitHitSetLineColor(1);\n",ScitilHit); if(tipo==0){ // disegna in XY. fprintf(MACRO,"Tile%d->SetLineWidth(2);\n",ScitilHit); } else if (tipo==1) { // disegna in SZ. fprintf(MACRO,"Tile%d->SetLineWidth(3);\n",ScitilHit); } else { // disegna in UV. fprintf(MACRO,"Tile%d->SetLineWidth(3);\n",ScitilHit); } fprintf(MACRO,"Tile%d->Draw();\n",ScitilHit); /* fprintf(MACRO,"TMarker* SciT%d = new TMarker(%f,%f,%d);\n", ScitilHit,posx,posy,30); fprintf(MACRO,"SciT%d->SetMarkerSize(1.5);\n",ScitilHit); fprintf(MACRO,"SciT%d->SetMarkerColor(1);\nSciT%d->Draw();\n" ,ScitilHit,ScitilHit); */ } //----------end of function PndSttMvdTracking::disegnaSciTilHit //----------end of function PndSttMvdTracking::disegnaassiXY void PndSttMvdTracking::disegnaAssiXY( FILE * MACRO, double xmin, double xmax, double ymin, double ymax ) { fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n",xmin,0.,xmax,0.,xmin,xmax); fprintf(MACRO,"Assex->SetTitle(\"X\");\n"); fprintf(MACRO,"Assex->SetTitleOffset(1.5);\n"); fprintf(MACRO,"Assex->Draw();\n"); fprintf(MACRO,"TGaxis *Assey = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", 0.,ymin,0.,ymax,ymin,ymax); fprintf(MACRO,"Assey->SetTitle(\"Y\");\n"); fprintf(MACRO,"Assey->SetTitleOffset(1.5);\n"); fprintf(MACRO,"Assey->Draw();\n"); } //----------end of function PndSttMvdTracking::disegnaassiXY ClassImp(PndSttMvdTracking)