#include "glpk.h" #include "PndTracking.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 // definition of constants used in this class // integer constants #define MAXHITSINCELL 50 // the following MAXHITSINFIT cannot be too large because the fit // crashes 'silently' for too much memory consumption in the character arrays // or takes too long time; #define MAXHITSINFIT 12 #define MAXMVDMCPOINTS 2000 #define MINIMUMHITSPERTRACK 3 #define MINOUTERHITSPERTRACK 5 #define TIMEOUT 60 // floating point constants #define APOTEMAMAXINNERPARSTRAW 23.246827 #define APOTEMAMAXSKEWSTRAW 31.517569 // delimitation of the skew area #define APOTEMAMINOUTERPARSTRAW 31.863369 #define APOTEMAMINSKEWSTRAW 23.246827 // delimitation of the skew area #define BFIELD 2. // in Tesla #define CVEL 2.99792 // velocity of light #define DELTAnR 2. //range of nR in TrkAssociatedParallelHitsToHelixquater #define DIAMETERSTRAWTUBE 1. #define DIMENSIONSCITIL 2.85 // cm #define ERRORPIXEL 0.02611 #define ERRORSTRIP 0.02611 #define ERRORSQPIXEL 0.00068175 #define ERRORSQSTRIP 0.00068175 #define PI 3.141592654 #define PMAX 100. #define RSTRAWDETECTORMAX 40.73 // maximum radius of the Stt detector in cm #define RSTRAWDETECTORMIN 16.119 // minimum radius of the Stt detector in cm #define STRAWRADIUS 0.5 #define STRAWRESOLUTION 0.015 #define STRAW_SKEW_INCLINATION 3. #define STTDRIFTVEL 0.0025 // in cm/nsec #define VERTICALGAP 4. // (cm) gap between Left and Right sections of detector. using namespace std; // ----- Default constructor ------------------------------------------- PndTracking::PndTracking() : FairTask("Tracking") { fPersistence = kTRUE; 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"); } // ------------------------------------------------------------------------- PndTracking::PndTracking(Int_t verbose) : FairTask("Tracking") { fPersistence = kTRUE; 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"); } // ------------------------------------------------------------------------- PndTracking::PndTracking(int istamp, bool iplot, bool imc) : FairTask("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"); } // ------------------------------------------------------------------------- PndTracking::PndTracking(int istamp, bool iplot, bool imc, bool doSciTil) : FairTask("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 ---------------------------------------------------- PndTracking::~PndTracking() { } // ------------------------------------------------------------------------- //--------------- begin PndTracking::Initialization_ClassVariables void PndTracking::Initialization_ClassVariables() { // this is only for initializing the Class Variables. size_t len; // booleans : len = sizeof(SingleHitListStt); memset (SingleHitListStt,true,len); len = sizeof(InclusionListSciTil); memset (InclusionListSciTil,true,len); len = sizeof(InclusionListStt); memset (InclusionListStt,true,len); len = sizeof(inMvdTrackCandPixel); memset (inMvdTrackCandPixel,false,len); len = sizeof(inMvdTrackCandStrip); memset (inMvdTrackCandStrip,false,len); len = sizeof(TypeConf); memset (TypeConf,false,len); // char : len = sizeof(fSttBranch); memset (fSttBranch,0,len); len = sizeof(fMvdPixelBranch); memset (fMvdPixelBranch,0,len); len = sizeof(fMvdStripBranch); memset (fMvdStripBranch,0,len); // UShort_t : nMCTracks=0; nSciTilHits=0; len = sizeof(ListMvdPixelHitsinTrack); memset (ListMvdPixelHitsinTrack,0,len); len = sizeof(ListMvdStripHitsinTrack); memset (ListMvdStripHitsinTrack,0,len); len = sizeof(ListSciTilHitsinTrack); memset (ListSciTilHitsinTrack,0,len); len = sizeof(ListSttParHitsinTrack); memset (ListSttParHitsinTrack,0,len); len = sizeof(ListSttSkewHitsinTrack); memset (ListSttSkewHitsinTrack,0,len); len = sizeof(ListSttSkewHitsinTrackSolution); memset (ListSttSkewHitsinTrackSolution,0,len); len = sizeof(ListTrackCandHit); memset (ListTrackCandHit,0,len); len = sizeof(nMvdPixelHitsinTrack); memset (nMvdPixelHitsinTrack,0,len); len = sizeof(nMvdStripHitsinTrack); memset (nMvdStripHitsinTrack,0,len); 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); // Short_t : nMvdDSPixelHitNotTrackCand=0; nMvdDSStripHitNotTrackCand=0; nMvdPixelHit=0; nMvdStripHit=0; nMvdTrackCand=0; nMvdUSPixelHitNotTrackCand=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 : Fimin=0.; SEMILENGTH_STRAIGHT=75.; ZCENTER_STRAIGHT=35.; nRdivConformalEffective=0; len = sizeof(ALFA); memset (ALFA,0,len); len = sizeof(BETA); memset (BETA,0,len); len = sizeof(GAMMA); memset (GAMMA,0,len); 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(radiaConf); memset (radiaConf,0,len); len = sizeof(Ox); memset (Ox,0,len); len = sizeof(Oy); memset (Oy,0,len); len = sizeof(R); memset (R,0,len); len = sizeof(refindexMvdPixel); memset (refindexMvdPixel,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(XMvdPixel); memset (XMvdPixel,0,len); len = sizeof(YMvdPixel); memset (YMvdPixel,0,len); len = sizeof(ZMvdPixel); memset (ZMvdPixel,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(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 PndTracking::Initialization_ClassVariables // ----- Public method Init -------------------------------------------- InitStatus PndTracking::Init() { IVOLTE=-1; SEMILENGTH_STRAIGHT = 75.; ZCENTER_STRAIGHT = 35.; 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); } // --------------------------- 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- PndTracking::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- PndTracking::Init: No MCTrack array, return!" << endl; return kERROR; } // ------------------------- get the SciTil hits if(YesSciTil) { fSciTHitArray = (TClonesArray*) ioman->GetObject("SciTHit"); } else { fSciTHitArray = NULL; } //--------------------------- // 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- PndTracking::Init: " << "No STTHit array, return!" << endl; return kERROR; } // ------------------------- get the Mvd hits fMvdPixelHitArray = (TClonesArray*) ioman->GetObject(fMvdPixelBranch); // fMvdPixelHitArray = (TClonesArray*) ioman->GetObject("MVDHitsPixel"); if ( !fMvdPixelHitArray){ cout << "-W- PndTracking::Init: " << "No MVD Pixel hitArray, return!" <GetObject(fMvdStripBranch); // fMvdStripHitArray = (TClonesArray*) ioman->GetObject("MVDHitsStrip"); if ( !fMvdStripHitArray){ cout << "-W- PndTracking::Init: " << "No MVD Strip hitArray, return!" <GetObject("MVDRiemannTrackCand"); if ( !fMvdTrackCandArray){ cout << "-W- PndTracking::Init: " << "No MVD TrackCand Array, return!" <GetObject("MVDPoint"); if ( !fMvdMCPointArray){ cout << "-W- PndTracking::Init: " << "No MVD MC Point Array, return!" <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); //----------------------- // calculate the boundaries of the Box in Conformal Space, see Gianluigi logbook on pag. 210-211 UShort_t i; Double_t A, r1, r2; radiaConf[0] = 1./RSTRAWDETECTORMAX; r1 = RSTRAWDETECTORMIN; A = (RSTRAWDETECTORMAX - r1)/NRDIVCONFORMAL; if ( NRDIVCONFORMAL > 1 ) { for(i = 1; i< NRDIVCONFORMAL ; i++){ r2 = r1 + A; radiaConf[NRDIVCONFORMAL-i] = 1./r2; r1=r2; } } nRdivConformalEffective = NRDIVCONFORMAL; return kSUCCESS; } // ------------------------------------------------------------------------- void PndTracking::SetParContainers() { FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb(); fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar"); } void PndTracking::WriteHistograms(){ TFile* file = FairRootManager::Instance()->GetOutFile(); file->cd(); file->mkdir("PndTracking"); file->cd("PndTracking"); hdeltaRPixel->Write(); hdeltaRStrip->Write(); hdeltaRPixel2->Write(); hdeltaRStrip2->Write(); delete hdeltaRPixel; delete hdeltaRStrip; delete hdeltaRPixel2; delete hdeltaRStrip2; } // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- void PndTracking::Exec(Option_t* opt) { bool flag, intersect, keepit[MAXTRACKSPEREVENT], Mvdhits[MAXTRACKSPEREVENT], outcome, status[MAXTRACKSPEREVENT], SttSZfit[MAXTRACKSPEREVENT]; UShort_t nalone, ncand, nhitsinfit, Nint, NNN, nRemainingCandidates, nXYZhits, nHitsInMCTrack[MAXTRACKSPEREVENT], nMCParalAlone[MAXTRACKSPEREVENT], nMCSkewAlone[MAXTRACKSPEREVENT], nParalCommon[MAXTRACKSPEREVENT], nSkewCommon[MAXTRACKSPEREVENT], nSkewHitsInMCTrack[MAXTRACKSPEREVENT], nSpuriParinTrack[MAXTRACKSPEREVENT], nSpuriSkewinTrack[MAXTRACKSPEREVENT], // given a Hit number it gives its radial box number RConformalIndex[MAXSTTHITS], // given a Hit number it gives its azimuthal box number FiConformalIndex[MAXSTTHITS], nTotalCandidates, tempore[MAXSTTHITS], TemporarySkewList[2*MAXSTTHITS][2], BigList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], MCParalAloneList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], MCSkewAloneList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], // nBoxConformal, first index -> radial divisions, // 2nd index -> azimuthal divisions; n. of hits falling in this cell. nBoxConformal[NRDIVCONFORMAL][NFIDIVCONFORMAL], ParalCommonList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], ParSpuriList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], SkewCommonList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], SkewSpuriList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], HitsinBoxConformal[MAXHITSINCELL][NRDIVCONFORMAL][NFIDIVCONFORMAL]; Short_t i, iParHit, ipunto, j, k, kall, l, nFicell, nRcell, tubeID, Charge[MAXTRACKSPEREVENT], daTrackFoundaTrackMC[MAXTRACKSPEREVENT], resultFitSZagain[MAXTRACKSPEREVENT], statusflag[MAXTRACKSPEREVENT] ; UInt_t iaccept, nSttHit, nSttMCPoint, nSttParHit, nSttSkewHit, nSttTrackCand ; Int_t nrounds0, nrounds1, ListHits[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK]; Double_t Distance, Fi, Ptras, ddd, delta, dis, dista, dista0, dista1, emme, gap, highqualitycut, Pxini, px, Pyini, py, Pzini, qop, x, y, // AloneX[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK], AloneY[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK], DriftRadiusbis[2*MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+ MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK],// all skew hits have double ErrorchosenPixel[MAXMVDPIXELHITS], ErrorchosenStrip[MAXMVDSTRIPHITS], ErrorchosenSkew[MAXSTTHITS], ErrorDriftRadiusbis[2*MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+ MAXMVDSTRIPHITSINTRACK+2],// solutions FI0[MAXTRACKSPEREVENT], Fi_final_helix_referenceframe[MAXTRACKSPEREVENT], Fi_initial_helix_referenceframe[MAXTRACKSPEREVENT], Fi_low_limit[MAXTRACKSPEREVENT], Fi_up_limit[MAXTRACKSPEREVENT], KAPPA[MAXTRACKSPEREVENT], Posiz1[3], primoangolo[MAXTRACKSPEREVENT], s[2], Sbis[2*MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK +MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK], // multiplication by 2 in the SchosenPixel[MAXTRACKSPEREVENT][MAXMVDPIXELHITS], SchosenStrip[MAXTRACKSPEREVENT][MAXMVDSTRIPHITS], SchosenSkew[MAXTRACKSPEREVENT][MAXSTTHITS], // NO multiplication by 2 here because for the // skew hits only one // solution is selected. Start[3], TemporaryS[2*MAXSTTHITS], TemporaryZ[2*MAXSTTHITS], TemporaryZDrift[2*MAXSTTHITS], TemporaryZErrorafterTilt[2*MAXSTTHITS], temporeZErrorafterTilt[MAXSTTHITSINTRACK], temporeS[MAXSTTHITSINTRACK], temporeZ[MAXSTTHITSINTRACK], temporeZDrift[MAXSTTHITSINTRACK], tmpErrorZDrift[MAXSTTHITSINTRACK+MAXSCITILHITSINTRACK], tmpS[MAXSTTHITSINTRACK+MAXSCITILHITSINTRACK], tmpZ[MAXSTTHITSINTRACK+MAXSCITILHITSINTRACK], tmpZDrift[MAXSTTHITSINTRACK+MAXSCITILHITSINTRACK], trajectory_vertex[2], Trajectory_Start[MAXTRACKSPEREVENT][2], ultimoangolo[MAXTRACKSPEREVENT], versor[2], WDX[MAXSTTHITS], WDY[MAXSTTHITS], WDZ[MAXSTTHITS], X[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK], XintersectionList[2], Y[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK], YintersectionList[2], z[2], ZEDbis[2*MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK +MAXMVDSTRIPHITSINTRACK+MAXSCITILHITSINTRACK], // rather improbable chance that zdrift[2], zerror[2], ZDrift[2*MAXSTTHITS], ZErrorafterTilt[2*MAXSTTHITS], zeta0, zeta1, // info[MAXSTTHITS][7], infoparalConformal[MAXSTTHITS][5], Sfinal[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], XY[MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK][2], ZchosenPixel[MAXTRACKSPEREVENT][MAXMVDPIXELHITS], ZchosenStrip[MAXTRACKSPEREVENT][MAXMVDSTRIPHITS], ZchosenSkew[MAXTRACKSPEREVENT][MAXSTTHITS], Zfinal[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], ZDriftfinal[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], ZErrorafterTiltfinal[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK]; TVector3 ErrMomentum, ErrPosition, Momentum, Position; FairMCPoint *puntator; PndMCTrack* pMCtr; PndSttHit *pSttHit; PndSttTube *pSttTube; PndSdsHit *pMvdPixelHit, *pMvdStripHit; //------------------------------------ IVOLTE++; if(istampa>0) cout<<"\nEntering in PndTrack : evt (starting from 0) n. "<GetEntriesFast(); nMvdStripHit = fMvdStripHitArray->GetEntriesFast(); if(nMvdPixelHit>MAXMVDPIXELHITS){ cout<<"from PndTracking, nMvdPixelHit is > maximum allowed (" <MAXMVDSTRIPHITS){ cout<<"from PndTracking, nMvdStripHit is > maximum allowed (" <At(i); 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(); } // ------------------------------------------- extract info from HITS Strip MVD for( i= 0; i< nMvdStripHit; i++){ pMvdStripHit = (PndSdsHit *) fMvdStripHitArray->At(i); TVector3 temp = pMvdStripHit->GetPosition(); XMvdStrip[i] = temp.X(); YMvdStrip[i] = temp.Y(); ZMvdStrip[i] = temp.Z(); sigmaXMvdStrip[i] = pMvdStripHit->GetDx(); sigmaYMvdStrip[i] = pMvdStripHit->GetDy(); sigmaZMvdStrip[i] = pMvdStripHit->GetDz(); refindexMvdStrip[i] = pMvdStripHit->GetRefIndex(); } if(istampa>=2) stampaMvdHits(); // ------------------------------------------ get info from trackcand of MVD nMvdTrackCand = fMvdTrackCandArray->GetEntriesFast(); if (nMvdTrackCand> MAXMVDTRACKSPEREVENT) { cout<<"da PndTracking : N. of MvdTrackCand = "<< nMvdTrackCand<<" and it is > MAXMVDTRACKSPEREVENT (="<=3)cout<<"da PndTracking : n. totale di Mvd track cand = " <GetEntriesFast(); if (nSttMCPoint ==0){ cout<<"da PndTracking : N. di Stt MC points = 0\n"<MAXSTTHITS){ cout<<"da PndTracking : N. di Stt MC points = "< MAXSTTHITS ("<1) stampaMCSttPoints(nSttMCPoint); //---------- fetching the STT hits nSttHit = fSttHitArray->GetEntriesFast(); if (nSttHit ==0){ cout<<"da PndTracking : N. di Stt Hits = 0, return!\n"< MAXSTTHITS) { cout<<"da PndTracking : N. di Stt Hits = "< MAXSTTHITS (="<= 1) { cout<<"da PndTracking : evento (partendo da 0) N. "<< IVOLTE<< "\n\tN. totale Hits in STT : "<At(i); // right way to extract the corrisponding MC point. ipunto= pSttHit->GetRefIndex(); 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.; ListSttParHits[nSttParHit]=i; nSttParHit++; } else { info[i][5]= 99.;// to signal that it is a skew straw. ListSttSkewHits[nSttSkewHit]=i; nSttSkewHit++; } if (istampa >= 1) { stampaSttHits( i, ipunto, dradius, WDX, WDY, WDZ, puntator, pSttTube ); } } // end of for( i= 0; i< nSttHit; i++) // fill the inclusion list for Stt, include only first hit for those straws with // multiple hits. MakeInclusionListStt(nSttHit, info); //----------------------------------- end of exclusion of straws with multiple hits //-------------------------------------------- fetch the SciTil hits if( fSciTHitArray != NULL){ // number SciTil hits/event nSciTilHits = fSciTHitArray->GetEntriesFast(); if(istampa>0) cout<<"da PndTracking, 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 PndTracking SciTil non purgati, Xpos "<< posiz.X()<<", Ypos "<< posiz.Y()<<", Zpos "<At(j); posiz = pPndSciTHit->GetPosition(); if(istampa>0) cout<<"da PndTracking SciTil non purgati, Xpos " <0){ //-----------stampe. if(istampa>0){ cout<<"da PndTracking, dopo purga di SciTil; n. hits = "<0){ cout<<"da PndTracking, dopo SciTil section.\n"; stampetta(nSttTrackCand, keepit); } //-------fine stampe. //----- loop over the parallel hits // begins the first iteration with more severe cuts on the # hits in track candidate for(iParHit=0; iParHit MAXTRACKSPEREVENT) continue; if( ! InclusionListStt[ListSttParHits[iParHit]] ) continue; nRcell = RConformalIndex[ListSttParHits[iParHit]]; nFicell = FiConformalIndex[ListSttParHits[iParHit]]; outcome = FindTrackInXYProjection( iParHit,//seed hit; it is positive for STT Hits. nRcell, nFicell, nSttParHit, info, InclusionListStt, RConformalIndex, FiConformalIndex, nBoxConformal, HitsinBoxConformal, nSttTrackCand, nSttParHitsinTrack, ListSttParHitsinTrack, trajectory_vertex, infoparalConformal, 1., // dummy value, there is no SciTil info in this case; 1., // dummy value, there is no SciTil info in this case NULL, // dummy pointer, there is no // S_SciTilHitsinTrack[nSttTrackCand][0] info. Fi_low_limit, Fi_up_limit, Fi_initial_helix_referenceframe, Fi_final_helix_referenceframe, Charge, &U[nSttTrackCand][0], &V[nSttTrackCand][0] ); if(!outcome) continue; // -------- here the track and its hits were found, filling the Inclusion list for(j=0; j0){ cout<<"da PndTracking, dopo || hits loop.\n"; stampetta(nSttTrackCand, keepit); } //-------fine stampe. nSttTrackCand++; } // end of for(iParHit=0; iParHit0) { zmin = -FI0/KAPPA; zmax = (2.*PI-FI0)/KAPPA; } else { zmax = -FI0/KAPPA; zmin = (2.*PI-FI0)/KAPPA; } deltaz = zmax-zmin; for(i=0; i zmax ){ tempZ[sign]=fmod( tempZ[sign]-zmax, deltaz) + zmin; } else if (tempZ[sign] -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 PndTracking::AssociateFoundTrackstoMC //----------begin of function PndTracking::AssociateFoundTrackstoMCbis void PndTracking::AssociateFoundTrackstoMCbis( Double_t info[][7], UShort_t nTracksFoundSoFar, UShort_t nHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][MAXSTTHITS], UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][MAXSTTHITS], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool flaggo, inclusionMC[MAXTRACKSPEREVENT][MAXSTTHITS], inclusionExp[MAXTRACKSPEREVENT]; UShort_t ntoMCtrack[MAXTRACKSPEREVENT], toMCtracklist[MAXTRACKSPEREVENT][MAXSTTHITS], toMCtrackfrequency[MAXTRACKSPEREVENT][MAXSTTHITS]; UShort_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 PndTracking::AssociateFoundTrackstoMCbis //----------begin of function PndTracking::AssociateFoundTrackstoMCtris void PndTracking::AssociateFoundTrackstoMCtris( Double_t info[][7], UShort_t nTracksFoundSoFar, UShort_t nHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][MAXSTTHITS], UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][MAXSTTHITS], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool firstime, flaggo, inclusionMC[MAXTRACKSPEREVENT][MAXSTTHITS], inclusionExp[MAXTRACKSPEREVENT]; UShort_t ntoMCtrack[MAXTRACKSPEREVENT], toMCtrackfrequency[MAXTRACKSPEREVENT][MAXSTTHITS]; UShort_t i, j, jtemp,jexp , nmid; Short_t itemp, massimo, toMCtracklist[MAXTRACKSPEREVENT][MAXSTTHITS]; Int_t enne; Double_t dx, Cx, Cy, Rr, alfa, beta, gamma, minimo, tanlow[MAXTRACKSPEREVENT], tanmid[MAXTRACKSPEREVENT], tanup[MAXTRACKSPEREVENT], toMCtrackdistance[MAXTRACKSPEREVENT][MAXSTTHITS]; 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 = (UShort_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 PndTracking::AssociateFoundTrackstoMCtris //----------begin of function PndTracking::AssociateFoundTrackstoMCquater void PndTracking::AssociateFoundTrackstoMCquater( bool *keepit, Double_t info[][7], Double_t X1[MAXTRACKSPEREVENT], Double_t Y1[MAXTRACKSPEREVENT], Double_t X2[MAXTRACKSPEREVENT], Double_t Y2[MAXTRACKSPEREVENT], Double_t X3[MAXTRACKSPEREVENT], Double_t Y3[MAXTRACKSPEREVENT], UShort_t nTracksFoundSoFar, UShort_t nHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], UShort_t nPixelHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListPixelHitsinTrack[MAXTRACKSPEREVENT][MAXMVDPIXELHITSINTRACK], Int_t *FromPixeltoMCTrack, UShort_t nStripHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListStripHitsinTrack[MAXTRACKSPEREVENT][MAXMVDSTRIPHITSINTRACK], Int_t *FromStriptoMCTrack, Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool firstime, flaggo, inclusionMC[nTracksFoundSoFar][MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK], inclusionExp[nTracksFoundSoFar]; UShort_t ntoMCtrack[nTracksFoundSoFar], toMCtrackfrequency[nTracksFoundSoFar][MAXSTTHITSINTRACK]; UShort_t i, j, jtemp,jexp , nmid; Short_t itemp, massimo, toMCtracklist[nTracksFoundSoFar][MAXSTTHITSINTRACK]; Int_t enne; Double_t dx, Cx, Cy, Rr, alfa, beta, gamma, minimo, tanlow[nTracksFoundSoFar], tanmid[nTracksFoundSoFar], tanup[nTracksFoundSoFar], toMCtrackdistance[nTracksFoundSoFar][MAXSTTHITSINTRACK]; 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 PndTracking::AssociateFoundTrackstoMCquater //----------begin of function PndTracking::AssociateSciTilHit UShort_t PndTracking::AssociateSciTilHit( Double_t Oxx, Double_t Oyy, Double_t Rr, UShort_t *List, // output Double_t *esse // output ) { bool intersect; UShort_t igoodScit, iScitHit, Nint; Double_t distance, QQ, sqrtRR, SIGN, XintersectionList[2], YintersectionList[2]; igoodScit=0; for(iScitHit=0; iScitHit 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]/Rr; if( distance >= info[i][4]+Aellipsis1 ){ continue; } S[NAssociated] = atan2(POINTS1[j+1]-Oyy, POINTS1[j]-Oxx) ; // atan2 returns radians in (-pi and +pi] if( S[NAssociated] < 0.) S[NAssociated] += 2.*PI; if( S[NAssociated] < Fi_low_limit) { if( S[NAssociated]+2.*PI > Fi_up_limit) continue; } else if( S[NAssociated] > Fi_up_limit) { if( S[NAssociated]- 2.*PI < Fi_low_limit) continue; } //--------------------------- end check Z[NAssociated] = POINTS1[j+2]; ZDrift[NAssociated] = Aellipsis1*Tiltdirection1[0]; ZErrorafterTilt[NAssociated] = 0.02*aaa*Tiltdirection1[0]/LL; SkewList[NAssociated][0] = i; // n. skew hit in original hit numbering SkewList[NAssociated][1] = ii; // solution 0 or solution 1 were accepted // check if this skew hit doesn't "push out" the most external parallel hit // (see Gianluigi's logbook on page 251) /* Double_t Zh1 = Z[NAssociated] - ZDrift[NAssociated]; Double_t Zh2 = Z[NAssociated] + ZDrift[NAssociated]; Double_t Sh1 = S[NAssociated] - Aellipsis1*Tiltdirection1[1]; Double_t Sh2 = S[NAssociated] + Aellipsis1*Tiltdirection1[1]; Double_t Zlast1 = (Fi_final_helix_referenceframe-Fi_initial_helix_referenceframe)*Zh1 /(Sh1-Fi_initial_helix_referenceframe); Double_t Zlast2 = (Fi_final_helix_referenceframe-Fi_initial_helix_referenceframe)*Zh2 /(Sh2-Fi_initial_helix_referenceframe); if( fabs(Zlast1 - ZCENTER_STRAIGHT) > 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 PndTracking::AssociateSkewHitsToXYTrack //----------begin of function PndTracking::BadTrack_ParStt bool PndTracking::BadTrack_ParStt( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Xcross[2], // Xcross[0]=point of entrance; // Xcross[1]=point of exit. Double_t Ycross[2], UShort_t nHits, UShort_t* ListHits, Double_t info[][7], Double_t cut, UShort_t maxnum, UShort_t islack// uncertainty allowed as far as // the n. of hits that should be present in a given section of the Stt track. ) { UShort_t ibad, ihit, ninside; Double_t cut2, length, Xprevious, Yprevious, S, Distance[MAXSTTHITS+1]; cut2=cut*cut; ibad=0; Xprevious=Xcross[0]; Yprevious=Ycross[0]; length=CalculateArcLength(Oxx, Oyy, Rr, Charge, Xcross, Ycross ); if(istampa>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 "< MAXMCTRACKS){ cout<<"da PndTracking : N. di MC truth tracks = "< MAXMCTRACKS = "<=2) stampaMCTracks(); // return; // --------------------------------------------- get MC Points of MVD nMvdMCPoint = fMvdMCPointArray->GetEntriesFast(); if(nMvdMCPoint>MAXMVDMCPOINTS) { cout<<"from PndTracking, nMvdMCPoint = "< the maximum number allowed ("<=2) cout<<"N. MC Points delle Mvd = "<0 && nTotalCandidates > 0){ 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; i0) AssociateFoundTrackstoMCquater( keepit, info, X1, Y1, X2, Y2, X3, Y3, nTotalCandidates, nSttParHitsinTrack, ListSttParHitsinTrack, nSttSkewHitsinTrack, ListSttSkewHitsinTrack, nMvdPixelHitsinTrack, ListMvdPixelHitsinTrack, FromPixeltoMCTrack, nMvdStripHitsinTrack, ListMvdStripHitsinTrack, FromStriptoMCTrack, daTrackFoundaTrackMC ); //---------- stampe. if(istampa>=3){ 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 PndTracking : 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 PndTracking, 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 ] ); //------------------------ 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 if( nMCTracks >0 && nTotalCandidates > 0) } //------------------------- end function PndTracking::ComparisonwithMC //------------------------- begin of function PndTracking::Dist_SZ Double_t PndTracking::Dist_SZ( Double_t Rr, Double_t KAPPA, Double_t FI0, Double_t ZED, Double_t S, Int_t *nrounds ) { // Defining : ZZ = (S-FI0)/KAPPA // this method returns the distance (WITH ITS SIGN ) : ZZ - ZED. Therefore this number // can be negative. // Care is taken to calculate this distance properly taking into // account that we are dealing with the function FI = mod(KAPPA*Z + FI0, 2*3.14). // the limits of an Int_t are : -2147483648 <= n <= 2147483647 Double_t aaa, ABSdis1, dis1, dis2, dis_segments, gap; if(fabs(KAPPA) < 1.e-10){ return -999999999.; } else if (fabs(KAPPA)>1.e10) { return -ZED; } // gap = fabs(2.*PI/KAPPA); aaa = (KAPPA*ZED)/(2.*PI); if( aaa<-2147483648.) { *nrounds = -2147483647; } else if (aaa > 2147483647.) { *nrounds = 2147483647; } else { *nrounds = (Int_t) aaa; } dis_segments = 2.*PI*Rr/sqrt(1.+KAPPA*KAPPA*Rr*Rr); // distance between // two consecutive segments // of trajectory. dis1 = fabs( Rr*S -Rr*KAPPA*ZED - Rr*FI0)/sqrt( 1.+KAPPA*KAPPA*Rr*Rr); dis1 = fmod(dis1,dis_segments); dis2 = dis_segments-dis1; if(dis2<0.)dis2=0.; if( dis1 < dis2 ) { return dis1; } else { return dis2; } } //------------------------- end of function PndTracking::Dist_SZ //----------end of function PndTracking::disegnaassiXY void PndTracking::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 PndTracking::disegnaassiXY //----------begin of function PndTracking::disegnaSciTilHit void PndTracking::disegnaSciTilHit( FILE * MACRO, int ScitilHit, double posx, double posy, int tipo ) { double x1,x2,y1,y2,L,Rr, RR; L=DIMENSIONSCITIL/2.; if(tipo==0){ // SciTil disegnate in XY. Rr = sqrt(posx*posx+posy*posy); x1 = posx + posy*L/Rr; x2 = posx - posy*L/Rr; y1 = posy - posx*L/Rr; y2 = posy + posx*L/Rr; } else if (tipo==1) {// SciTil disegnate in SZ. x1 = posx - L; x2 = posx + L; y1 = posy; y2 = posy; } else { // SciTil disegnate in UV. RR = posx*posx+posy*posy; Rr = sqrt(RR); x1 = posx + posy*L/Rr; x1 /= RR; x2 = posx - posy*L/Rr; x2 /= RR; y1 = posy - posx*L/Rr; y1 /= RR; y2 = posy + posx*L/Rr; y2 /= RR; } fprintf(MACRO,"TLine *Tile%d = new TLine(%f,%f,%f,%f);\n",ScitilHit,x1,y1,x2,y2); fprintf(MACRO,"Tile%d->SetLineColor(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); } //----------end of function PndTracking::disegnaSciTilHit //----------begin of function PndTracking::DrawBiHexagonInMacro void PndTracking::DrawBiHexagonInMacro( Double_t vgap, FILE * MACRO, Double_t Ami, Double_t Ama, UShort_t color, char *name ) { UShort_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 PndTracking::DrawBiHexagonInMacro //----------begin of function PndTracking::DrawHexagonCircleInMacro void PndTracking::DrawHexagonCircleInMacro( Double_t GAP, FILE * MACRO, Double_t ApotemaMin, Double_t Rma, UShort_t color, char *name ) { UShort_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 PndTracking::DrawHexagonCircleInMacro //------------------------- begin of function PndTracking::EliminateSpuriousSZ void PndTracking::EliminateSpuriousSZ( UShort_t MaxTurnofTracks, UShort_t *nPixelHitsinTrack, UShort_t *ListPixelHitsinTrack, UShort_t *nStripHitsinTrack, UShort_t *ListStripHitsinTrack, UShort_t *nSkewHitsinTrack, UShort_t *ListSkewHitsinTrack, Double_t *S, Double_t *ZED, Double_t *DriftRadius, Double_t *ErrorDriftRadius, Double_t *SchosenPixel, Double_t *SchosenStrip, Double_t *SchosenSkew, Double_t *ZchosenPixel, Double_t *ZchosenStrip, Double_t *ZchosenSkew, Double_t *ErrorchosenPixel, Double_t *ErrorchosenStrip, Double_t *ErrorchosenSkew, Double_t KAPPA, Double_t FI0, Double_t Rr ) { UShort_t i, j, k, Kmin, auxnMvdPixel, auxListMvdPixel[MAXMVDPIXELHITS], auxnMvdStrip, auxListMvdStrip[MAXMVDSTRIPHITS], auxnSttSkew, auxListSttSkew[MAXSTTHITS]; Int_t nr2, nrounds0, nrounds1, nchosen, Nround[4]; const Double_t MvdCut=0.8, // const Double_t MvdCut=0.3, // const Double_t MvdCut=0.1, minimumSttDriftError = 1.; Double_t ddd, dista, dista1, dista0, d_min, error, zeta0, zeta1, Dista[4], Errore[4], Esse[4], Zeta[4]; auxnMvdPixel=0; auxnMvdStrip=0; auxnSttSkew=0; for(i=0;i<*nPixelHitsinTrack;i++){ if( fabs(Dist_SZ(Rr,KAPPA,FI0,ZED[i],S[i],&nrounds0)) < MvdCut && abs(nrounds0)<=MaxTurnofTracks){ auxListMvdPixel[auxnMvdPixel]=ListPixelHitsinTrack[i]; if(nrounds0>=0)SchosenPixel[ListPixelHitsinTrack[i]]=S[i]+nrounds0*2.*PI; else SchosenPixel[ListPixelHitsinTrack[i]]=S[i]+(nrounds0-1)*2.*PI; ZchosenPixel[ListPixelHitsinTrack[i]]=ZED[i]; ErrorchosenPixel[ListPixelHitsinTrack[i]]=ERRORPIXEL; auxnMvdPixel++; } } // end of for(i=0;i<*nPixelHitsinTrack;i++) for(j=0;j<*nStripHitsinTrack;j++){ i=j+(*nPixelHitsinTrack); if( fabs(Dist_SZ(Rr,KAPPA,FI0,ZED[i],S[i],&nrounds0)) < MvdCut && abs(nrounds0)<=MaxTurnofTracks){ auxListMvdStrip[auxnMvdStrip]=ListStripHitsinTrack[j]; if(nrounds0>=0)SchosenStrip[ListStripHitsinTrack[j]]=S[i]+nrounds0*2.*PI; else SchosenStrip[ListStripHitsinTrack[j]]=S[i]+(nrounds0-1)*2.*PI; ZchosenStrip[ListStripHitsinTrack[j]]=ZED[i]; ErrorchosenStrip[ListStripHitsinTrack[j]]=ERRORSTRIP; auxnMvdStrip++; } } // end of for(j=0;j<*nStripHitsinTrack;j++) //-----------------stampe. if(istampa>=3){ cout<<"in EliminateSpuriousSZ : nSkew hit = "<<*nSkewHitsinTrack <<", K = "<< KAPPA<<", FI0 = "<< FI0 <1){ cout<<"in EliminateSpuriousSZ : skew hit n. "< 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< MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK (=" <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 < MAXMVDPIXELHITSINTRACK){ // the following is a protection because the maximum Mvd hit n. // cannot exceed MAXMVDPIXELHITS . if( pndtrackcandhit.GetHitId()>MAXMVDPIXELHITS ){ cout<<"from PndTracking, this Pixel Mvd hit has a number = " < MAXMVDPIXELHITS (="<< MAXMVDPIXELHITS<<"), rejected!\n"; continue; } inMvdTrackCandPixel[ pndtrackcandhit.GetHitId() ]= true; ListHitTypeMvdTrackCand[i][k] = pndtrackcandhit.GetDetId(); // this in reality is the Branch name of // Hit; this Branch name is used to identify // a Pixel. // If it is -1 is supposed to be noise but // a in which Pixel or Strip?? Mistery. kPixel++; // this is a Strip. } else if( pndtrackcandhit.GetDetId()== FairRootManager::Instance()->GetBranchId(fMvdStripBranch) && kStrip < MAXMVDSTRIPHITSINTRACK){ // the following is a protection because the maximum Mvd hit n. cannot // exceed MAXMVDSTRIPHITS . if( pndtrackcandhit.GetHitId()>MAXMVDSTRIPHITS ){ cout<<"from PndTracking, this Strip Mvd hit has a number = " < MAXMVDSTRIPHITS (="<< MAXMVDSTRIPHITS<<"), rejected!\n"; continue; } inMvdTrackCandStrip[ pndtrackcandhit.GetHitId() ]= true; ListHitTypeMvdTrackCand[i][k] = pndtrackcandhit.GetDetId(); // this in reality is the Branch name of // Hit; this Branch name is used to identify // a Pixel. // If it is -1 is supposed to be noise but // a in which Pixel or Strip?? Mistery. kStrip++; } else { // this is the case should (in principle) never happen. continue; // ignore this hit. } // the following is the native Hit number that one can use // to extract all the info. If it is -1 I won't consider it (noise). ListHitMvdTrackCand[i][k] = pndtrackcandhit.GetHitId(); k++; } // end of for(j=0; j=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) stampaMvdTrackCandInfo(); } //----------------------- end of function PndTracking::ExtractInfoFromMvdTrackCand //----------begin of function PndTracking::FindCharge void PndTracking::FindCharge( Double_t oX, Double_t oY, UShort_t nParallelHits, Double_t *X, Double_t *Y, Short_t * Charge ) { UShort_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 PndTracking::FindCharge //------------------------- begin of function PndTracking::FindDistance Double_t PndTracking::FindDistance( Double_t Oxx, // center from wich distance is calculated Double_t Oyy, // center from wich distance is calculated Double_t Rr, Double_t tanlow, Double_t tanmid, Double_t tanup, Double_t alfa, // intersection circumference parameter Double_t beta, // intersection circumference parameter Double_t gamma // intersection circumference parameter ) { UShort_t i, n; Double_t Delta, m[3], q, dist, dist1, dist2, distlow, distmid, distup, totaldist, x1, x2, y1, y2; int nevento=1; m[0] = tanlow; m[1] = tanmid; m[2] = tanup; n=0; totaldist=0.; for(i=0;i<3;i++){ if( tanlow<999998.) { q = Oyy-m[i]*Oxx; Delta = alfa*alfa + beta*beta*m[i]*m[i] - 4.*gamma*m[i]*m[i] - 4.*q*q + 4.*m[i]*q*alfa + + 2.*alfa*beta*m[i] - 4.*beta*q - 4.*gamma; if( Delta < 0.){ dist = -1.; } else if (Delta==0.){ x1 = 0.5*(-alfa - 2.*m[i]*q - beta*m[i] )/(1.+m[i]*m[i]); y1 = m[i]*x1+q; dist = fabs( sqrt( (Oxx-x1)*(Oxx-x1) + (Oyy-y1)*(Oyy-y1) ) - Rr); } else { Delta = sqrt(Delta); x1 = 0.5*(-alfa - 2.*m[i]*q - beta*m[i] - Delta)/(1.+m[i]*m[i]); x2 = 0.5*(-alfa - 2.*m[i]*q - beta*m[i] + Delta)/(1.+m[i]*m[i]); y1 = m[i]*x1+q; y2 = m[i]*x2+q; dist1 = fabs( sqrt( (Oxx-x1)*(Oxx-x1) + (Oyy-y1)*(Oyy-y1) ) - Rr); dist2 = fabs( sqrt( (Oxx-x2)*(Oxx-x2) + (Oyy-y2)*(Oyy-y2) ) - Rr); if(dist1-0.5) { totaldist+=dist; n++; } } // end of for(i=0;i<3;i++) if(n!=3) totaldist = -1.; else totaldist = totaldist / n; return totaldist; } //------------------------- end of function PndTracking::FindDistance //----------start function PndTracking::FindingParallelTrackAngularRange void PndTracking::FindingParallelTrackAngularRange( Double_t oX, Double_t oY, Double_t Rr, Short_t Charge, Double_t *Fi_low_limit, Double_t *Fi_up_limit, Short_t * status, Double_t Rmi, // Rmin of cylindrical volume intersected by track; Double_t Rma // Rmax of cylindrical volume intersected by track; ) { // -------------- calculate the maximum fi and minimum fi spanned by this track, // see logbook pag.270; by using the Rmin and Rmax of the straw detector. // this function works also when circular trajectory in XY doesn't pass through (0,0). bool intersection_inner, intersection_outer; Double_t teta1, teta2, tetavertex, a, cosT, cost, cosFi, cosfi, Fi, fi, FI0, Px, Py, tmp; Rma += 1. ; // add a safety margin. Rmi -= 1. ; // add a safety margin. a = sqrt(oX*oX+oY*oY); // preliminary condition if(a + Rr <= Rmi ) // in this case there might be hits at radius < Rmi. { *status = -1 ;return;} if( a >= Rr + Rma || Rr >= a + Rma) // in this case there can be no hits with radius < Rma. { *status = -2;return;} if( a - Rr >= Rmi ) intersection_inner = false; else intersection_inner = true; if( a + Rr <= Rma || a - Rr >= Rma ) intersection_outer = false; else intersection_outer = true; if( (! intersection_inner) && (! intersection_outer) ){ *Fi_low_limit = 0.; *Fi_up_limit = 2.*PI; *status = 1; return; } // now the calculation FI0 = atan2(-oY,-oX); if( intersection_outer ){ cosFi = (a*a + Rr*Rr - Rma*Rma)/(2.*Rr*a); if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.; Fi = acos(cosFi); } if( intersection_inner ){ cosfi = (a*a + Rr*Rr - Rmi*Rmi)/(2.*Rr*a); if(cosfi<-1.) cosfi=-1.; else if(cosfi>1.) cosfi=1.; fi = acos(cosfi); } if( Charge < 0.){ // this particle rotates counterclockwise when looking into the beam if( intersection_outer && intersection_inner){ *Fi_low_limit=FI0 + fi; *Fi_up_limit= FI0 +Fi; } else if (intersection_inner) { *Fi_low_limit=FI0 + fi; *Fi_up_limit= FI0 - fi; } else { *Fi_low_limit=FI0 - Fi; *Fi_up_limit= FI0 + Fi; } // end of if( intersection_outer && intersection_inner } else { // continuation of if( Charge < 0.) if( intersection_outer && intersection_inner){ *Fi_low_limit=FI0 - Fi; *Fi_up_limit= FI0 - fi; } else if (intersection_inner) { *Fi_low_limit=FI0 + fi; // must invert because low limit must be < up limit *Fi_up_limit= FI0 - fi; } else { *Fi_low_limit=FI0 - Fi; *Fi_up_limit= FI0 + Fi; } // end of if( intersection_outer && intersection_inner } // end of if( Charge < 0.) if(*Fi_low_limit<0.) { *Fi_low_limit=fmod(*Fi_low_limit,2.*PI); *Fi_low_limit += 2.*PI; } else if (*Fi_low_limit>=2.*PI){ *Fi_low_limit=fmod(*Fi_low_limit,2.*PI); } if(*Fi_up_limit<0.) { *Fi_up_limit=fmod(*Fi_up_limit,2.*PI); *Fi_up_limit += 2.*PI; } else if (*Fi_up_limit>=2.*PI){ *Fi_up_limit=fmod(*Fi_up_limit,2.*PI); } // Modify *Fi_up_limit by adding // 2PI if it is the case, in order to make *Fi_up_limit > *Fi_low_limit. if( *Fi_up_limit < *Fi_low_limit ) *Fi_up_limit += 2.*PI; if( *Fi_up_limit < *Fi_low_limit ) *Fi_up_limit = *Fi_low_limit; *status = 0; return; } //---------- end of function PndTracking::FindingParallelTrackAngularRange //----------begin of function PndTracking::FindIntersectionsOuterCircle Short_t PndTracking::FindIntersectionsOuterCircle( Double_t oX, Double_t oY, Double_t Rr, 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 >= Rr + Rma || Rr >= a + Rma || a + Rr <= Rma) return -1; FI0 = atan2(-oY,-oX); cosFi = (a*a + Rr*Rr - Rma*Rma)/(2.*Rr*a); if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.; Fi = acos(cosFi); Xcross[0] = oX + Rr*cos(FI0+Fi); Ycross[0] = oY + Rr*sin(FI0+Fi); Xcross[1] = oX + Rr*cos(FI0-Fi); Ycross[1] = oY + Rr*sin(FI0-Fi); return 0; } //----------end of function PndTracking::FindIntersectionsOuterCircle //----------begin of function PndTracking::FindTrackEntranceExitbiHexagonLeft Short_t PndTracking::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; UShort_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 PndTracking::FindTrackEntranceExitbiHexagonLeft //----------begin of function PndTracking::FindTrackEntranceExitbiHexagonRight Short_t PndTracking::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; UShort_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 PndTracking::FindTrackEntranceExitbiHexagonRight //----------begin function PndTracking::FindTrackEntranceExitHexagonCircleLeft Short_t PndTracking::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; //------------------ UShort_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 PndTracking::FindTrackEntranceExitHexagonCircleLeft //----------begin of function PndTracking::FindTrackEntranceExitHexagonCircleRight Short_t PndTracking::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; //------------------ UShort_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 PndTracking::FindTrackEntranceExitHexagonCircleRight //----------begin of function PndTracking::FindTrackStrictCollection Short_t PndTracking::FindTrackStrictCollection( UShort_t NFiCELLDISTANCE, // seed track (original notation) as far as the Fi angle is concerned UShort_t iSeed, // n. of hits to search in ListHitsinTrackinWhichToSearch UShort_t NParallelToSearch, UShort_t *ListHitsinTrackinWhichToSearch, bool *InclusionList, UShort_t *FiConformalIndex, UShort_t *OutputListHitsinTrack ) { UShort_t i, iFi, iFiseed, iR, j, nHitsinTrack; Double_t auxRvalues[MAXSTTHITS]; // iSeed is the hit number in the ORIGINAL number scheme iFiseed = FiConformalIndex[iSeed]; nHitsinTrack=0; for(i=0; i0){ cout<<"PndSttTracking, evt. "<MAXSTTHITSINTRACK) { return false; } //----------------------- // find among the ListHitsinTrack if there are at least // a minimum # of hits belonging to the outer part of the STT system. // At this point of the code the Stt hits are already ordered from // the outermost to the innermost. for(j=0, Nouter =0; j= MINOUTERHITSPERTRACK) { for(i=0; i< Nouter;i++){ ListHitsinTrackinWhichToSearch[i]= ListHitsinTrack[nTracksFoundSoFar][i]; } for(i=0; i< Nouter;i++){ Naux = FindTrackPatterninBoxConformalSpecial( 3, // NRCELLDISTANCE 1, // NFiCELLDISTANCE nsttparhit, Nouter, ListHitsinTrackinWhichToSearch[i], // seed hit. ListHitsinTrackinWhichToSearch, info, InclusionList, RConformalIndex, FiConformalIndex, nBoxConformal, HitsinBoxConformal, OutputListHitsinTrack); if( Naux >= MINOUTERHITSPERTRACK && Naux > 0.7 * Nouter ) break; if( Naux >= MINOUTERHITSPERTRACK) { // further collection of hits in the inner region but this time strictly connected // to the outer ones // first the list of non outer hits for(j=Nouter;j= MINIMUMHITSPERTRACK && nHitsinTrack[nTracksFoundSoFar]<=MAXSTTHITSINTRACK) { for(j=0;j= .... } // end of if( Naux >= MINOUTERHITSPERTRACK) } // end of for(i=0; i< Nouter;i++) } // end of if( Nouter >= MINOUTERHITSPERTRACK) if( nHitsinTrack[nTracksFoundSoFar] < MINIMUMHITSPERTRACK || nHitsinTrack[nTracksFoundSoFar]>MAXSTTHITSINTRACK) { return false; } //--------------------------- // finding the rotation angle for best utilization of the MILP procedure for(j=0, rotationcos=0., rotationsin=0.; j0){ cout<<"PndSttTrackFinderReal, in FindTrackInXYProjection, evt. " <0){ nSciTilHitsinTrack[nTracksFoundSoFar]=nSciT; for(j=0;j=3) { cout<<"prima di TrkAssociatedParallelHitsToHelixQuater, evt. " <MAXSTTHITSINTRACK) return false; nHitsinTrack[nTracksFoundSoFar]=NN; for(i=0; i 0 always and it is between 0. and // 2 PI here (later, FixDiscontinuitiesFiangleinSZplane may change it adding +2PI or // -2PI if necessary). // Fi_final_helix_referenceframe is made such that it is < Fi_initial_helix_referenceframe // when the track is positive, and it is > Fi_initial_helix_referenceframe for negative // tracks. OrderingParallel( Ox[nTracksFoundSoFar], Oy[nTracksFoundSoFar], info, nHitsinTrack[nTracksFoundSoFar], &ListHitsinTrack[nTracksFoundSoFar][0], Charge[nTracksFoundSoFar], &Fi_initial_helix_referenceframe[nTracksFoundSoFar],//output &Fi_final_helix_referenceframe[nTracksFoundSoFar],// output U, V ); // finding the FI angular range (in the laboratory frame) spanned by this parallel track FindingParallelTrackAngularRange( Ox[nTracksFoundSoFar], Oy[nTracksFoundSoFar], R[nTracksFoundSoFar], Charge[nTracksFoundSoFar], &Fi_low_limit[nTracksFoundSoFar], &Fi_up_limit[nTracksFoundSoFar], &flagStt, RSTRAWDETECTORMIN, RSTRAWDETECTORMAX ); if( flagStt == -2 ) return false; // track outside cylinder RSTRAWDETECTORMAX. if( flagStt == -1 ) return false; // track inside cylinder RSTRAWDETECTORMIN. if( flagStt == 1 ) return false; // track comprised in cylinder : discard, // because at this stage only tracks from // vertex are searched. return true; }; //----------end of function PndTracking::FindTrackInXYProjection //----------begin of function PndTracking::FindTrackPatterninBoxConformal Short_t PndTracking::FindTrackPatterninBoxConformal( UShort_t NRCELLDISTANCE, UShort_t NFiCELLDISTANCE, UShort_t Nparal, Short_t ihit, // seed hit; Short_t nRcell, // R cell of the seed hit; // can be negative beacuse of SciTil hits; UShort_t nFicell, // Fi cell of the seed hit; Double_t info[][7], bool *InclusionList, UShort_t *RConformalIndex, UShort_t *FiConformalIndex, UShort_t nBoxConformal[][NFIDIVCONFORMAL], UShort_t HitsinBoxConformal[][NRDIVCONFORMAL][NFIDIVCONFORMAL], UShort_t *ListHitsinTrack ) { bool status, TemporaryInclusionList[MAXSTTHITS]; Short_t i; UShort_t j, iFi, iR, nRmin, nRmax, nRemainingHits, nHitsinTrack, auxIndex[MAXSTTHITS], Remaining[MAXSTTHITS]; Short_t iFi2; Double_t auxRvalues[MAXSTTHITS]; // ihit is the hit number in the PARALLEL number scheme for(i=0, nRemainingHits=0; i=0){ nHitsinTrack=1; ListHitsinTrack[0]= ihit ; i = 0; } else { nHitsinTrack=0; i = -1; } status=true; while( nRemainingHits > 0 && i < nHitsinTrack && status) { if (nRcell - NRCELLDISTANCE < 0 ) { nRmin = 0; } else { nRmin = nRcell - NRCELLDISTANCE; } if (nRcell + NRCELLDISTANCE >= nRdivConformalEffective ) { nRmax = nRdivConformalEffective-1; } else { nRmax = nRcell + NRCELLDISTANCE; } if(istampa>0) {cout<<"\tin FindTrackPatterninBoxConformal, nRmin = "<=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 PndTracking::FitSZspace //----------start of function PndTracking::FixDiscontinuitiesFiangleinSZplane void PndTracking::FixDiscontinuitiesFiangleinSZplane( UShort_t TemporarynSkewHitsinTrack, Double_t *S, Double_t *Fi_initial_helix_referenceframe, Short_t Charge ) { UShort_t i; Double_t max, min; if( Charge >0 ) { for(i=0 ; i *Fi_initial_helix_referenceframe ) S[i]-= 2.*PI; } } else { for(i=0 ; iAt(MCTrack); 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) { *Rr = -3.; return;} *Cx = Oxx + Pyy*1000./(BFIELD*CVEL*carica); *Cy = Oyy - Pxx*1000./(BFIELD*CVEL*carica); } else { *Rr = -1.; } return; } //------------------------- end of function PndTracking::getMCInfo //---------- begin of function PndTracking::InfoXYZParal void PndTracking::InfoXYZParal( Double_t info[][7], UShort_t infopar, Double_t Oxx, Double_t Oyy, Double_t Rr, Double_t KAPPA, Double_t FI0, Short_t Charge, Double_t *Posiz ) { Double_t fi, norm, vers[2]; vers[0] = Oxx - info[infopar][0]; vers[1] = Oyy - info[infopar][1]; norm = sqrt( vers[0]*vers[0] + vers[1]*vers[1] ); if(norm < 1.e-20) { Posiz[0] = -999999999.; return; } if( fabs( Rr - fabs( norm - info[infopar][3] ) ) // distance trajectory-drift radius < fabs( Rr - (norm + info[infopar][3]) ) ) { Posiz[0] = info[infopar][0] + info[infopar][3]*vers[0]/norm; Posiz[1] = info[infopar][1] + info[infopar][3]*vers[1]/norm; } else { Posiz[0] = info[infopar][0] - info[infopar][3]*vers[0]/norm; Posiz[1] = info[infopar][1] - info[infopar][3]*vers[1]/norm; } // end of if ( fabs( Rr - fabs( Distance - info[infopar][3] ) )..... // Posiz[0] = info[infopar][0] + info[infopar][3]*vers[0]/norm; // Posiz[1] = info[infopar][1] + info[infopar][3]*vers[1]/norm; if( fabs(KAPPA)<1.e-20 ){ Posiz[2] = -888888888.; return; } fi = atan2(-vers[1],-vers[0]); if(fi<0.) fi += 2.*PI; if ( Charge > 0){ if(fi > FI0 ) FI0 += 2.*PI; // Posiz[2] = (FI0-fi)/KAPPA; } else { if(fi < FI0 ) fi += 2.*PI; } Posiz[2] = (fi-FI0)/KAPPA; return; } //----------end of function PndTracking::InfoXYZParal //----------begin of function PndTracking::IntersectionCircle_Segment bool PndTracking::IntersectionCircle_Segment( Double_t a, // coefficients implicit equation. Double_t b, // of segment : a*x + b*y + c =0. Double_t c, Double_t P1x, // point delimiting the segment. Double_t P2x, // point delimiting the segment. Double_t P1y, // point delimiting the segment. Double_t P2y, // point delimiting the segment. Double_t Oxx, // center of circle. Double_t Oyy, Double_t Rr, // Radius of circle. UShort_t * Nintersections, Double_t XintersectionList[2], Double_t YintersectionList[2], Double_t *distance ) { // this method finds the intersection of a circle with a segment. If the circle // passes through an endpoint of the segment, that is considered an intersection also. bool status; Short_t ipossibility; Double_t aperp, bperp, cperp, det, distq, dist1, dist2, length, length_segmentq, Rq, x, y, Xintersection, Yintersection; *Nintersections=0; det = a*a+b*b ; if(det < 1.e-20){ cout<<"from PndTracking::IntersectionCircle_Segment :" <<" this is not the equation of a segment, return!\n"; return false; } // find if intersection circle - segment is possible. // check if there is an intersection (with the squares, it's the same!). distq = ( a*Oxx+ b*Oyy + c )*( a*Oxx+ b*Oyy + c )/det; *distance = sqrt(distq); Rq=Rr*Rr; length = Rq - distq; if(length <= 0. ) return false; // no intersection between trajectory and this // segment. length = sqrt(length); // coefficients of line perpendicular to input segment and passing for (Oxx, Oyy). aperp = -b; bperp = a; cperp = - aperp*Oxx - bperp*Oyy; // find intersection of segment with perpendicular : no need to check if // det is different from 0. Xintersection = (-bperp*c + b*cperp)/det; Yintersection = (-a*cperp + aperp*c)/det; det = sqrt(det); length_segmentq = (P1x-P2x)*(P1x-P2x) + (P1y-P2y)*(P1y-P2y); status = false; for(ipossibility=-1;ipossibility<2;ipossibility +=2){ x = Xintersection + ipossibility*length*b/det ; y = Yintersection - ipossibility*length*a/det ; if ( (x-P1x)*(x-P1x)+(y-P1y)*(y-P1y) > 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 PndTracking::IntersectionCircle_Segment //----------begin of function PndTracking::IntersectionSciTil_Circle bool PndTracking::IntersectionSciTil_Circle( Double_t posizSciTilx, Double_t posizSciTily, Double_t Oxx, // center of circle. Double_t Oyy, Double_t Rr, // Radius of circle. UShort_t * Nintersections, Double_t XintersectionList[2], Double_t YintersectionList[2] ) { bool intersect; Double_t distance, QQ, sqrtRR, SIGN; QQ = posizSciTilx*posizSciTilx+posizSciTily*posizSciTily; sqrtRR=sqrt(QQ); if( posizSciTily<0.) SIGN=-1.; else SIGN=1.; intersect = IntersectionCircle_Segment( posizSciTilx, posizSciTily, -QQ, posizSciTilx-SIGN*0.5*DIMENSIONSCITIL*posizSciTily/sqrtRR, posizSciTilx+SIGN*0.5*DIMENSIONSCITIL*posizSciTily/sqrtRR, posizSciTily-SIGN*posizSciTilx*0.5*DIMENSIONSCITIL/sqrtRR, posizSciTily+SIGN*posizSciTilx*0.5*DIMENSIONSCITIL/sqrtRR, Oxx, Oyy, Rr, Nintersections, // output XintersectionList, // output YintersectionList, // output &distance // output ); return intersect; } //----------end of function PndTracking::IntersectionSciTil_Circle //----------start function PndTracking::IntersectionsWithClosedbiHexagonLeft Short_t PndTracking::IntersectionsWithClosedbiHexagonLeft( //-------- inputs Double_t vgap, Double_t Oxx, Double_t Oyy, Double_t Rr, Double_t Ami, // min Apotema of hexagonal volume intersected by track; Double_t Ama, // max Apotema of hexagonal volume intersected by track; //-------- outputs UShort_t *nIntersections, Double_t *XintersectionList, Double_t *YintersectionList ) { // return integer convention : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon; // 1 --> track contained completely between the two polygons; // inner Hexagon --> 0 // outer Hexagon --> 1 bool internal, AtLeast1; UShort_t i, is, j, Nintersections; Double_t aaa, maxdistq, distance, //------------------- // a,b,c == coefficients of the implicit equations of the 3 sides of the half inner Hexagon // plus 3 sides of the half outer Hexagon plus 2 vertical sides corresponding to the Gap : // a*x + b*y +c =0; the numbering of these 8 sides foolows the convention of Gianluigi's // logbook on page 286. a[] = {-1./sqrt(3.) , 1., 1./sqrt(3.), 1., 1./sqrt(3.), 1., -1./sqrt(3.), 1.}, b[] = {1., 0., 1., 0., 1., 0., 1., 0.}, c[] = {-2.*Ama/sqrt(3.),Ama, 2.*Ama/sqrt(3.),vgap/2.,2.*Ami/sqrt(3.),Ami, -2.*Ami/sqrt(3.),vgap/2.}, //---------------------- tempX[2], tempY[2]; // side_x and side_y are ordered as side1, side2, ..... in Gianluigi's logbook on page 286. 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.) }; //----------------------- // find intersections (maximum 16) with the 8 sides. AtLeast1 = false; *nIntersections =0; internal = true; maxdistq=-9999.; for(is=0; is<8; is++){ aaa = (side_x[is]-Oxx)*(side_x[is]-Oxx)+(side_y[is]-Oyy)*(side_y[is]-Oyy); if(aaa>maxdistq) maxdistq=aaa; aaa = (side_x[is+1]-Oxx)*(side_x[is+1]-Oxx)+(side_y[is+1]-Oyy)*(side_y[is+1]-Oyy); if(aaa>maxdistq) maxdistq=aaa; if ( IntersectionCircle_Segment( a[is], b[is], c[is], side_x[is], side_x[is+1], side_y[is], side_y[is+1], Oxx, Oyy, Rr, &Nintersections, tempX, tempY, &distance // distance of (Oxx,Oyy) from line // defined by a*x+b*y+c=0. ) ){ AtLeast1=true; for(j=0;j track outside outer perimeter; // 0 --> at least 1 intersection with polygon; // 1 --> track contained completely between the two polygons; // inner Hexagon --> 0 // outer Hexagon --> 1 bool internal, AtLeast1; UShort_t i, is, j, Nintersections; Double_t aaa, maxdistq, distance, //------------------- // a,b,c == coefficients of the implicit equations of the 3 sides of the half inner Hexagon // plus 3 sides of the half outer Hexagon plus 2 vertical sides corresponding to the Gap : // a*x + b*y +c =0; the numbering of these 8 sides foolows the convention of Gianluigi's // logbook on page 286. a[] = {1./sqrt(3.) , 1., -1./sqrt(3.), 1., -1./sqrt(3.), 1., 1./sqrt(3.), 1.}, b[] = {1., 0., 1., 0., 1., 0., 1., 0.}, c[] = {-2.*Ama/sqrt(3.),-Ama, 2.*Ama/sqrt(3.),-vgap/2.,2.*Ami/sqrt(3.),-Ami, -2.*Ami/sqrt(3.),-vgap/2.}, //---------------------- tempX[2], tempY[2]; // side_x and side_y are ordered as side1, side2, ..... in Gianluigi's logbook on page 286. 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.) }; //----------------------- // find intersections (maximum 16) with the 8 sides. AtLeast1 = false; *nIntersections =0; internal = true; maxdistq=-9999.; for(is=0; is<8; is++){ aaa = (side_x[is]-Oxx)*(side_x[is]-Oxx)+(side_y[is]-Oyy)*(side_y[is]-Oyy); if(aaa>maxdistq) maxdistq=aaa; aaa = (side_x[is+1]-Oxx)*(side_x[is+1]-Oxx)+(side_y[is+1]-Oyy)*(side_y[is+1]-Oyy); if(aaa>maxdistq) maxdistq=aaa; if ( IntersectionCircle_Segment( a[is], b[is], c[is], side_x[is], side_x[is+1], side_y[is], side_y[is+1], Oxx, Oyy, Rr, &Nintersections, tempX, tempY, &distance // distance of (Oxx,Oyy) from line // defined by a*x+b*y+c=0. ) ){ AtLeast1=true; for(j=0;j track outside outer polygon; // -1 --> track contained completely within inner polygon; // 0 --> at least 1 intersection with inner polygon, at least 1 with outer polygon; // 1 --> track contained completely between the two polygons; // 2 --> track contained completely by larger polygons, with intersections in the smaller; // 3 --> track completely outsiede the small polygon with intersections in the bigger. // inner Hexagon --> 0 // outer Hexagon --> 1 bool internal[2], AtLeast1[2]; UShort_t i, is, j, Nintersections; Double_t mindist[2], distance, //------------------- // a,b,c == coefficients of the implicit equations of the six sides of the Hexagon // centered at 0 : a*x + b*y +c =0; see Gianluigi's logbook on page 277; // the coefficient c has to be multiplied by Erre. // The first side is a[] = { 1./sqrt(3.) , 1., -1./sqrt(3.), 1./sqrt(3.) , 1. , -1./sqrt(3.) } , b[] = { 1., 0. , 1., 1., 0. , 1. }, c[] = { -2./sqrt(3.) , -1., 2./sqrt(3.), 2./sqrt(3.), 1., -2./sqrt(3.)}, //---------------------- Erre[] = {Rmi, Rma}, // this is the distance from (0,0) // of the Verteces of the Hexagon delimiting the Skew area tempX[2], tempY[2]; // both hexagon_side_xlow and hexagon_side_xup must be multiplied by appropriate Erre; // sides of Hexagon ordered as a, b, c ..... in Gianluigi's logbook on page 280. Double_t hexagon_side_x[] = { 0., 1. , 1., 0., -1., -1., 0. }, hexagon_side_y[] = { 2./sqrt(3.),1./sqrt(3.),-1./sqrt(3.), -2./sqrt(3.),-1./sqrt(3.), 1./sqrt(3.), 2./sqrt(3.)}; //----------------------- // find intersection with the 6 sides of the small exhagon delimiting the skew straws zone // see on page 277 of Gianluigi's logbook. // status =0, at least 1 intersection with inner Hexagon, at least 1 intersection // with the outer Hexagon; =1, track contained completely between the // two Hexagons; = -1 track contained within inner Hexagon; // =-2 track outside outer Hexagon. for(i=0;i<2;i++){ // i=0 --> inner Hexagon, i= 1 --> outer Hexagon. AtLeast1[i] = false; nIntersections[i]=0; internal[i] = true; mindist[i]=999999.; for(is=0; is<6; is++){ if ( IntersectionCircle_Segment(a[is], b[is], c[is]*Erre[i], hexagon_side_x[is]*Erre[i], hexagon_side_x[is+1]*Erre[i], hexagon_side_y[is]*Erre[i], hexagon_side_y[is+1]*Erre[i], Oxx, Oyy, Rr, &Nintersections, tempX, tempY, &distance // distance of (Oxx,Oyy) from line // defined by a*x+b*y+c=0. ) ){ AtLeast1[i]=true; for(j=0;jdistance) mindist[i]=distance; // the definition of 'internal' here is when the given Point // stays at the same side of the origin (0,0) with respect to // the given line of equation a*x+b*y+c=0. internal[i] = internal[i] && IsInternal(Oxx, Oyy, a[is], b[is], c[is]*Erre[i] ); } // end of for(is=0; is<6; is++) } // end of for(i=0;i<2;i++) if( (!AtLeast1[0]) && (!AtLeast1[1]) ){ if (!internal[1]) return -2; // trajectory outside outer Hexagon. if( Rr > mindist[1]) return -2; if( !internal[0]) return 1; // trajectory contained between inner and outer Hexagon. if( mindist[0] >= Rr) return -1; // trajectory contained in inner Hexagon. return 1; // trajectory contained between inner and outer Hexagon. } else if (AtLeast1[0] && AtLeast1[1] ){ // continuation of if( (!AtLeast1[0]) && ... return 0; } else if (AtLeast1[0]){ return 2; } return 3; } //---------- end of function PndTracking::IntersectionsWithClosedPolygon //----------begin of function PndTracking::IntersectionsWithGapSemicircle UShort_t PndTracking::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 ) { UShort_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 PndTracking::IntersectionsWithGapSemicircle //----------start function PndTracking::IntersectionsWithOpenPolygon UShort_t PndTracking::IntersectionsWithOpenPolygon( //-------- inputs Double_t Oxx, // Track parameter Double_t Oyy, // Track parameter Double_t Rr, // Track parameter UShort_t nSides, // input, n. of Sides of open Polygon. Double_t *a, // coefficient of formula : aX + bY + c = 0 defining Double_t *b, // the Polygon sides. Double_t *c, Double_t *Side_x, // X,Y coordinate of the Sides vertices (in sequence, following Double_t *Side_y, // the Polygon along. //-------- outputs Double_t *XintersectionList, // XintersectionList Double_t *YintersectionList // YintersectionList. ) { // this methods returns the n. of intersections. UShort_t i, is, j, nIntersections, Nintersections; Double_t mindist, distance, //------------------- // a,b,c == coefficients of the implicit equations of the six sides of the Hexagon // centered at 0 : a*x + b*y +c =0; see Gianluigi's logbook on page 277; // the coefficient c has to be multiplied by Erre. tempX[2], tempY[2]; //----------------------- nIntersections=0; for(is=0; is 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 PndTracking::IsInTargetPipe //----------star of function PndTracking::IsInternal bool PndTracking::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 PndTracking::IsInternal //---------------------- begin function PndTracking::LoadPndTrack_TrackCand void PndTracking::LoadPndTrack_TrackCand( bool *keepit, bool *SttSZfit, UShort_t nTotalCandidates, Short_t *Charge, UInt_t nSttTrackCand, Double_t *FI0, Double_t *KAPPA, Double_t info[][7], Double_t SchosenSkew[][MAXSTTHITS], Double_t ZchosenSkew[][MAXSTTHITS] ) { //------- 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 UInt_t ipinco, j, k, l, ncand ; Double_t ddd, dis, Distance, Oxx, Oyy, Ptras, Pxini, Pyini, Pzini, px, py, qop, x, y, Posiz1[3], versor[2] ; TVector3 ErrMomentum, ErrPosition, Momentum, Position ; 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); pTrckCand->setMcTrackId( -1 ); 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 InfoXYZParal( 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 // InfoXYZParal 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); FairTrackParP first( Position, Momentum, ErrPosition, ErrMomentum, Charge[ncand], Position, TVector3(py/Ptras, -px/Ptras, 0.), // first vector defining the plane TVector3(Pzini*px/ddd,Pzini*py/ddd,-Ptras*Ptras/ddd) //second vector defining the plane ); // the last hit k = nTrackCandHit[ncand]-1; if (ListTrackCandHitType[ncand][k] == 0){ // Mvd Pixel Posiz1[0] = XMvdPixel[ ListTrackCandHit[ncand][k] ]; Posiz1[1] = YMvdPixel[ ListTrackCandHit[ncand][k] ]; Posiz1[2] = ZMvdPixel[ ListTrackCandHit[ncand][k] ]; ErrPosition.SetX(sigmaXMvdPixel[ ListTrackCandHit[ncand][k] ]/sqrt(12.)); ErrPosition.SetY(sigmaXMvdPixel[ ListTrackCandHit[ncand][k] ]/sqrt(12.)); ErrPosition.SetZ(sigmaXMvdPixel[ ListTrackCandHit[ncand][k] ]/sqrt(12.)); } else if (ListTrackCandHitType[ncand][k] == 1){ // Mvd Strip Posiz1[0] = XMvdStrip[ ListTrackCandHit[ncand][k] ]; Posiz1[1] = YMvdStrip[ ListTrackCandHit[ncand][k] ]; Posiz1[2] = ZMvdStrip[ ListTrackCandHit[ncand][k] ]; ErrPosition.SetX(sigmaXMvdStrip[ ListTrackCandHit[ncand][k] ]/sqrt(12.)); ErrPosition.SetY(sigmaXMvdStrip[ ListTrackCandHit[ncand][k] ]/sqrt(12.)); ErrPosition.SetZ(sigmaXMvdStrip[ ListTrackCandHit[ncand][k] ]/sqrt(12.)); } else if( ListTrackCandHitType[ncand][k] == 2 ){ // it is a parallel straw hit InfoXYZParal ( info, ListTrackCandHit[ncand][k], Ox[ncand], Oy[ncand], R[ncand], KAPPA[ncand], FI0[ncand], Charge[ncand], Posiz1 ); 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][k] == 3 ){ // it is a skew straw hit Posiz1[0] = Ox[ncand]+R[ncand]*cos(SchosenSkew[ncand][ ListTrackCandHit[ncand][k] ]); Posiz1[1] = Oy[ncand]+R[ncand]*sin(SchosenSkew[ncand][ ListTrackCandHit[ncand][k] ]); Posiz1[2] = ZchosenSkew[ncand][ ListTrackCandHit[ncand][k] ]; 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); // the plane of this FairTrackParP better is perpendicular to // the momentum direction ddd = Ptras*sqrt(Ptras*Ptras+Pzini*Pzini); FairTrackParP last( Position, Momentum, ErrPosition, ErrMomentum, Charge[ncand], Position, TVector3(py/Ptras, -px/Ptras, 0.), // first vector defining the plane TVector3(Pzini*px/ddd,Pzini*py/ddd,-Ptras*Ptras/ddd) //second vector defining the plane ); // loading actually the PndTrack PndTrack *pTrck = new((*fSttMvdPndTrackArray)[ipinco]) PndTrack(first,last,*pTrckCand); pTrck->SetRefIndex(ipinco); pTrck->SetFlag(0); ipinco++; } // end of for(ncand=0, ipinco = 0; ncand< nTotalCandidates; ncand++) } //----------end function PndTracking::LoadPndTrack_TrackCand //-------------------------------- begin function PndTracking::MakeInclusionListStt void PndTracking::MakeInclusionListStt( UInt_t nSttHit, Double_t info[][7] ) { // it needs to be initialized for each event ! size_t len = sizeof(InclusionListStt); memset (InclusionListStt,true,len); len = sizeof(SingleHitListStt); memset (SingleHitListStt,true,len); // fill the inclusion list for Stt, include only first hit for those straws with // multiple hits. for(int i=0; i< nSttHit-1; i++){ if( !InclusionListStt[ i ] ) continue; for(int 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 ) { SingleHitListStt[j]=InclusionListStt[j]= false ; } } // end of for(j=i+1; j< Nhits;; j++) } // end of for(i=0; i< Nhits-1; i++) } //-------------------------------- end of function PndTracking::MakeInclusionListStt //--------------------- begin function PndTracking::MatchMvdHitsToSttTracks void PndTracking::MatchMvdHitsToSttTracks( Double_t delta, UShort_t nSttTrackCand, Double_t *FI0, Double_t *Fifirst, UShort_t *nPixelHitsinTrack, // output UShort_t ListPixelHitsinTrack[][MAXMVDPIXELHITS], // output UShort_t *nStripHitsinTrack, // output UShort_t ListStripHitsinTrack[][MAXMVDSTRIPHITS] // output ) { bool specialcase; UShort_t i,j; Double_t angle, anglemax, anglemin; for(i=0; i=3 ) cout<<"MatchMvdHitsToSttTracks, i= "<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 PndTracking : appena prima arbitration, IVOLTE = "<< IVOLTE<<", Stt track cand = "<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){ 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]+ (MAXMVDPIXELHITS+MAXMVDSTRIPHITS)*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]<(MAXMVDPIXELHITS+MAXMVDSTRIPHITS)*10){//Pixel. ListTrackCandHit[ncand][i] = ListHits[i]; ListTrackCandHitType[ncand][i] = 0; } else { // Strip hits. ListTrackCandHit[ncand][i] = ListHits[i]- (MAXMVDPIXELHITS+MAXMVDSTRIPHITS)*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 ListParallelHits such that the first hits in the // list are alway those closer to the (0,0). // ordering of the hits aaa = atan2( oY, oX); // 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((b1-3.*b1&&aaa<-b1)){//use U as ordering variable; //[case 1 or 3 Gianluigi's Logbook page 285]. for (j = 0; j< nParallelHits; j++){ U[j]=info[ ListParallelHits[j] ][0]/( info[ ListParallelHits[j] ][0]* info[ ListParallelHits[j] ][0]+ info[ ListParallelHits[j] ][1]* info[ ListParallelHits[j] ][1]); } for (j = 0; j< nParallelHits; j++){ tempo[j]= (Int_t) ListParallelHits[j]; } Merge_Sort( nParallelHits, U, tempo); for (j = 0; j< nParallelHits; j++){ ListParallelHits[j] = (UShort_t) tempo[j]; } 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< nParallelHits; j++){ V[j]=info[ ListParallelHits[j] ][1]/( info[ ListParallelHits[j] ][0]* info[ ListParallelHits[j] ][0]+ info[ ListParallelHits[j] ][1]* info[ ListParallelHits[j] ][1]); } for (j = 0; j< nParallelHits; j++){ tempo[j]= (Int_t) ListParallelHits[j]; } Merge_Sort( nParallelHits, V, tempo); for (j = 0; j< nParallelHits; j++){ ListParallelHits[j] = (UShort_t) tempo[j]; } if((aaa<=-3.*b1 || aaa>=3.*b1)){ // case #2; if( Charge == -1){ // inverting the order of the hits. for(i=0;ib1&& .... // FI initial value (at 0,0 vertex) in the Helix reference frame *Fi_initial_helix_referenceframe = atan2(-oY,-oX) ;// this is in order to be coherent // with the calculatation of Fi, which is atan2(oY,oX). // atan2 is defined in [-PI,PI) if ( *Fi_initial_helix_referenceframe <0.) *Fi_initial_helix_referenceframe += 2.*PI; // FI of the last parallel hit in the Helix reference frame *Fi_final_helix_referenceframe = atan2( info[ ListParallelHits[nParallelHits-1] ][1]-oY, info[ ListParallelHits[nParallelHits-1] ][0]-oX ); if ( *Fi_final_helix_referenceframe <0.) *Fi_final_helix_referenceframe += 2.*PI; if ( Charge > 0 ) { if( *Fi_final_helix_referenceframe> *Fi_initial_helix_referenceframe) *Fi_final_helix_referenceframe -= 2.*PI; if( *Fi_final_helix_referenceframe> *Fi_initial_helix_referenceframe) *Fi_final_helix_referenceframe = *Fi_initial_helix_referenceframe; } else { if( *Fi_final_helix_referenceframe< *Fi_initial_helix_referenceframe) *Fi_final_helix_referenceframe += 2.*PI; if( *Fi_final_helix_referenceframe< *Fi_initial_helix_referenceframe) *Fi_final_helix_referenceframe = *Fi_initial_helix_referenceframe; } return; } //----------end of function PndTracking::OrderingParallel //----------begin of function PndTracking::OrderingR_Loading_ListTrackCandHit void PndTracking::OrderingR_Loading_ListTrackCandHit( bool *keepit, UShort_t ncand, Double_t info[][7] ) { UShort_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]; UShort_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 UShort_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 PndTracking::OrderingR_Loading_ListTrackCandHit //----------begin of function PndTracking::OrderingSttSkewandSttParallel void PndTracking::OrderingSttSkewandSttParallel( Double_t oX, Double_t oY, Double_t Rr, UShort_t nSkewhit, UShort_t *ListSkewHits, Double_t *SList, // this is rekated to the skew hits. IMPORTANT : // the index must be the ORIGINAL skew hit number, // therefore SList[Infoskew[ListSkewHits[*]]]. Short_t Charge, UShort_t nParHits, UShort_t *ListParHits, Double_t *U, Double_t *V, UShort_t *BigList // this is the final ordered Parallel+Skew list; // already in NATIVE hit number. ) { UShort_t i, j, tmp[nSkewhit+nParHits], tmpList[nSkewhit]; Int_t index[nSkewhit+nParHits]; Double_t aaa, b1, sign, aux[nSkewhit+nParHits]; // here there is the ordering of the hits, under the assumption that the circumference // in XY goes through (0,0). // Moreover, the code before is supposed to have selected trajectories in XY with (Ox,Oy) // farther from (0,0) by > 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 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. // ordering of the hits aaa = atan2( oY, oX); // 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)){ // case #1 or #3;see Gianluigi's Logbook page 285. if( (aaa>b1&&aaa<3.*b1 && Charge == -1)||( aaa>-3.*b1&&aaa<-b1 && Charge == 1) ) { // for speeding up the ordering taking advantage // that the parallel hits were earlier ordered and // apply the trick of multiplying by -1. sign=-1.; } else { // normal calculation sign=1.; } for (j = 0 ; j< nParHits; j++){ aux[j] = sign*U[j]; } for (j = 0; j< nSkewhit; j++){ // this is U in conformal space aux[j+nParHits]=sign*(oX + Rr*cos(SList[ListSkewHits[j]]))/ (oX*oX+oY*oY+Rr*Rr + 2.*Rr* (oX*cos(SList[ListSkewHits[j]]) +oY*sin(SList[ListSkewHits[j]]))); } } else { // use V as ordering variable // [case 2 and 4 Gianluigi's Logbook page 285]. if( ((aaa<=-3.*b1|| aaa>=3.*b1) && Charge == -1) || ( -b1 <= aaa && aaa <= b1 && Charge == 1) ){ sign=-1.; } else { sign=1.; } for (j = 0 ; j< nParHits; j++){ aux[j] = sign*V[j]; } for (j = 0; j< nSkewhit; j++){ // this is V in conformal space. aux[j+nParHits]=sign*(oY + Rr*sin(SList[ListSkewHits[j]]))/ (oX*oX+oY*oY+Rr*Rr + 2.*Rr* (oX*cos(SList[ListSkewHits[j]]) +oY*sin(SList[ListSkewHits[j]]))); } } // end of if((aaa>b1&& .... for (j = 0 ; j< nParHits; j++){ BigList[j]=ListParHits[j]; index[j] = j; } for (j = 0; j< nSkewhit; j++){ BigList[j+nParHits]=ListSkewHits[j]; index[j+nParHits] = j+nParHits; } Merge_Sort( nSkewhit+nParHits, aux, index); for(i=0, j=0;i= nParHits ){ tmpList[j] = ListSkewHits[index[i]-nParHits]; j++; } } for(i=0;i 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 PndTracking::OrderingUsingConformal //----------begin of function PndTracking::Ordering_Loading_ListTrackCandHit void PndTracking::Ordering_Loading_ListTrackCandHit( bool *keepit, UShort_t FirstCandidate, UShort_t LastCandidate, Double_t info[][7], Double_t Trajectory_Start[][2], Short_t *CHARGE, Double_t SchosenSkew[][MAXSTTHITS] ) { UShort_t ncand; for(ncand=FirstCandidate; ncand< LastCandidate; ncand++){ // for small radius trajectory better the ordering with conformal. if( R[ncand]< RSTRAWDETECTORMAX/2.){ OrderingConformal_Loading_ListTrackCandHit( keepit, ncand, info, 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 PndTracking::Ordering_Loading_ListTrackCandHit //------------------ begin function PndTracking::RefitMvdStt void PndTracking::RefitMvdStt( UShort_t nCandHit, UShort_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; UShort_t i, iparallel; UShort_t MAXIMUMHITSINFIT = 20; if(MAXIMUMHITSINFIT>MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK) MAXIMUMHITSINFIT = MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK; Short_t exitstatus; Double_t dist2, mindis, emme, factor, gamma, qu, ErrorStraw = 0.03, ErrorMvd = 0.01, Xconformal[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK], Yconformal[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK], DriftRadiusconformal[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK], ErrorDriftRadiusconformal[MAXSTTHITSINTRACK+MAXMVDPIXELHITSINTRACK+MAXMVDSTRIPHITSINTRACK]; *status= false; factor=3.; mindis=0.5; // trajectory_vertex[0]=trajectory_vertex[1]=0.; for(i=0, iparallel=0; i2){ cout<<"PndTracking::RefitMvdStt, n. Hits (Mvd+Stt || ) = "<= MAXTRACKSPEREVENT ( = " <2.*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 PndTracking::SeparateInnerOuterParallel //----------begin of function PndTracking::stampafinale void PndTracking::stampafinale( UInt_t nTotalCandidates, bool * keepit ) { int l; cout<<"Evt. n. "<0){ cout<<"\tscitil hit n. "<At(i); cout<<"\tpoint n. "<GetX()<< ", Ymc "<GetY()<<", Zmc "<GetZ() <<", appartiene a MC Track "<< pSttMC->GetTrackID()<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()<GetBranchId(fMvdPixelBranch) = "<< FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)<GetBranchId(fMvdStripBranch) = "<< FairRootManager::Instance()->GetBranchId(fMvdStripBranch)<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()< 0){ if(fi > FI0 ) FI0 += 2.*PI; // Posiz[2] = (FI0-fi)/KAPPA; } else { if(fi < FI0 ) fi += 2.*PI; } Posiz[2] = (fi-FI0)/KAPPA; return; } //---------- end function PndTracking::SttInfoXYZParal //------------------------- begin of function PndTracking::SttMatchedSpurious void PndTracking::SttMatchedSpurious( bool *keepit, UShort_t ntotalHits, Double_t info[][7], UShort_t nTracksFoundSoFar, // those found by PR UShort_t nHitsinTrack[], // n. hits PARALLELI+SKEW, from PR UShort_t ListHitsinTrack[][MAXSTTHITSINTRACK], // from PR UShort_t *nSkewHitsinTrack, // n. hits skew, from PR UShort_t ListSkewHitsinTrack[][MAXSTTHITSINTRACK], // from PR UShort_t *nParalCommon, UShort_t ParalCommonList[][MAXSTTHITSINTRACK], UShort_t *nSpuriParinTrack, UShort_t ParSpuriList[][MAXSTTHITSINTRACK], UShort_t *nSkewCommon, UShort_t SkewCommonList[][MAXSTTHITSINTRACK], UShort_t *nSpuriSkewinTrack, UShort_t SkewSpuriList[][MAXSTTHITSINTRACK], UShort_t *nHitsInMCTrack, UShort_t *nSkewHitsInMCTrack, UShort_t *nMCParalAlone, UShort_t MCParalAloneList[][MAXSTTHITSINTRACK], UShort_t *nMCSkewAlone, UShort_t MCSkewAloneList[][MAXSTTHITSINTRACK], Short_t *daTrackFoundaTrackMC ) { bool flaggo; UShort_t i, jexp, exphit, iHit, enne[MAXTRACKSPEREVENT][MAXSTTHITS]; 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; exphit1) { 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){ 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 PndTracking::SttParalCleanup //----------begin of function PndTracking::SttSkewCleanup bool PndTracking::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, UShort_t nHits, UShort_t *Listofhits, Double_t *S, Double_t info[][7], Double_t RminStrawSkew, Double_t RmaxStrawSkew, Double_t cut, // cut distance (in cm). UShort_t maxnum // max number allowed of failures to pass the cut. ) { bool ConsiderLastHit; Short_t flagSttL, flagSttR, flagOutStt; UShort_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[MAXSTTHITS+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) cout<<"\n\nevt "<1)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)cout<<"in SttSkewCleanup : hit preso!"<1)cout<<"in SttSkewCleanup : hit skew prima di purga = " <1)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){ 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)cout<<"in SttSkewCleanup, reject this track because ibad = "<< ibad <<" and it is > maxnum [="< 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 PndTracking::TrackCleanup //----------begin of function PndTrackFinderReal::TrkAssociatedParallelHitsToHelixQuater UShort_t PndTracking::TrkAssociatedParallelHitsToHelixQuater( bool *ExclusionList, Double_t m, Double_t q, Short_t Status, UShort_t nHitsinTrack, UShort_t *ListHitsinTrack, UInt_t NhitsParallel, Double_t Oxx, Double_t Oyy, Double_t Rr, Double_t info[][7], Double_t infoparalConformal[][5], UShort_t *RConformalIndex, UShort_t *FiConformalIndex, UShort_t nBoxConformal[][NFIDIVCONFORMAL], UShort_t HitsinBoxConformal[][NRDIVCONFORMAL][NFIDIVCONFORMAL], UShort_t *auxListHitsinTrack ) { bool passamin, passamax, Unselected[MAXSTTHITS]; // this has to be at least nSttHits large. Short_t i, i2, j, k, l, l2, l3, itemp, kstart, kend, iFi0,FFimin, FFimax; UShort_t Nextra=8, nFi, Fi, nR, nAssociatedHits, nHit_original; Double_t maxFi, minFi, dist, xx, yy, aaa, angle, r, erre1, erre2, Rin, Rout, Fi0, ddd, fi1, fi2, dx, dy, distance, NTIMES=1.5; // number of Straw radia allowed in association nAssociatedHits=0; size_t len = sizeof(InclusionListStt); memset (Unselected,true,len); // find the range in Fi spanned by the candidate track FFimin = 10000; FFimax = 0; for(j=0; j FFimax ) FFimax = FiConformalIndex[i]; } if( FFimax > 3.*NFIDIVCONFORMAL/4. && FFimin < NFIDIVCONFORMAL/4.) { FFimin = 10000; FFimax = 0; for(j=0; j FFimax ) FFimax = Fi; } } // finding the boundaries in the Conformal plane. The basic assumption is that the range // in Fi is much less that 180 degrees. FFimin -= (Short_t) NFIDIVCONFORMAL/Nextra; FFimax += (Short_t) NFIDIVCONFORMAL/Nextra; if( FFimax - FFimin > NFIDIVCONFORMAL/2 ) { cout<<"something fishy is going on in PndSttTrkAssociatedParallelHitsToHelixQuater!" <<"Range in Fi (rad) is "<<(FFimax - FFimin)*2.*PI/NFIDIVCONFORMAL< 1.e-10 ) { passamax=false; passamin=false; for(itemp=FFimin; itemp<=FFimax;itemp++){ i=itemp; if( i< 0 ) { i += NFIDIVCONFORMAL; } else if (i>=NFIDIVCONFORMAL){ i -= NFIDIVCONFORMAL*( i/NFIDIVCONFORMAL ); } angle = (i+0.5)*2.*PI/NFIDIVCONFORMAL; aaa = cos(angle); if( fabs(cos(angle)) <1.e-10) continue; r = -q/aaa; if(r< radiaConf[0] || r>= 1./RSTRAWDETECTORMIN) continue; for(j=nRdivConformalEffective-1; j>=0;j--){ if( r>= radiaConf[j] ){ nR = j; break; } } for(l=-DELTAnR; l= nRdivConformalEffective ) continue; for( k=0;k NTIMES*STRAWRADIUS ) continue; //------------------- xx=infoparalConformal[ HitsinBoxConformal[k][l2][i] ][0]; dist = fabs( xx +q ); if( AcceptHitsConformal( dist, infoparalConformal[HitsinBoxConformal[k][l2][i]][2], infoparalConformal[ HitsinBoxConformal[k][l2][i] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[k][l2][i]; nAssociatedHits++; } } // end of for( k=0;k=NFIDIVCONFORMAL) i2 -= NFIDIVCONFORMAL; for(l=-2; l<3;l++){ l3 = nR+l; if( l3<0 || l3 >= nRdivConformalEffective ) continue; for( k=0;k NTIMES*STRAWRADIUS ) continue; //------------------- xx=infoparalConformal[ HitsinBoxConformal[k][l3][i2] ][0]; dist = fabs( xx +q ); // if(dist < 3.*infoparalConformal[ HitsinBoxConformal[k][l3][i2] ][2]){ if( AcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[k][l3][i2] ][2], infoparalConformal[ HitsinBoxConformal[k][l3][i2] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[k][l3][i2]; nAssociatedHits++; } } // end of for( k=0;k= nRdivConformalEffective ) continue; for( k=0;k NTIMES*STRAWRADIUS ) continue; //------------------- xx=infoparalConformal[ HitsinBoxConformal[k][l3][i2] ][0]; dist = fabs( xx +q ); // if(dist < 3.*infoparalConformal[ HitsinBoxConformal[k][l3][i2] ][2]){ if( AcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[k][l3][i2] ][2], infoparalConformal[ HitsinBoxConformal[k][l3][i2] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[k][l3][i2]; nAssociatedHits++; } } // end of for( k=0;k x=0 if( FFimax > NRDIVCONFORMAL/4 && Fimin < NRDIVCONFORMAL/4 ) { iFi0 = (Short_t) (NRDIVCONFORMAL/4 ); } else if ( FFimax > 3*NRDIVCONFORMAL/4 && Fimin < 3*NRDIVCONFORMAL/4 ){ iFi0 = (Short_t) (3*NRDIVCONFORMAL/4 ); } else { cout <<"From PndSttTrackFinderReal::PndSttTrkAssociatedParallelHitsToHelixQuater :" <<" inconsistency, 0 associated hits to this track candidate\n"; return 0; } for(itemp=iFi0-5; itemp<=iFi0+5;itemp++){ i=itemp; if( i< 0 ) { i += NFIDIVCONFORMAL; } else if (i>=NFIDIVCONFORMAL){ i -= NFIDIVCONFORMAL*( i/NFIDIVCONFORMAL ); // i -= NFIDIVCONFORMAL; } for(l=0; lNTIMES*STRAWRADIUS ) continue; //------------------- xx=infoparalConformal[ HitsinBoxConformal[k][l][i] ][0]; dist = fabs( xx ); if( AcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[k][l][i] ][2], infoparalConformal[ HitsinBoxConformal[k][l][i] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[k][l][i]; nAssociatedHits++; } } // end of for( k=0;k 1.e-10 ) } else if( fabs(q)> 1.e-10) { // second part of if( Status ==99), in this case y = m*x +q Fi0 = atan(m); // Fi0 belongs to (-PI/2, PI/2] . aaa = atan2(q, -m*q); if(Fi0<0.) { Fi0 += PI; if (Fi0 <0. ) Fi0 =0.;// this is between 0. and PI. if (Fi0 >PI ) Fi0 =PI;// this is between 0. and PI. }; ddd= fabs(q)/sqrt(1.+m*m); for(itemp=FFimin; itemp<=FFimax;itemp++){ i=itemp; if( i< 0 ) { i += NFIDIVCONFORMAL; } else if (i>=NFIDIVCONFORMAL){ i -= NFIDIVCONFORMAL*( i/NFIDIVCONFORMAL ); } // erre1 is the distance from origin of point of intersection of the straight // line of equation y= m*x+q with line of equation y = x*tan(fi1); // when erre1 is < 0 it means the intersection is on the opposite side of the // versor defined by [cos(fi1); sin(fi1)]. // Here we are working in the conformal plane U,V. fi1 = i*2.*(PI/NFIDIVCONFORMAL); if( fabs(sin(fi1)-m*cos(fi1))>1.e-10) { erre1 = q/(sin(fi1)-m*cos(fi1)); } else { erre1 = 99999999999.; } fi2 = (i+1)*2.*(PI/NFIDIVCONFORMAL); if( fabs(sin(fi2)-m*cos(fi2))>1.e-10) { erre2 = q/(sin(fi2)-m*cos(fi2)); } else { erre2 = 99999999999.; } for(j=0; j Rout ){ continue; } } else if(fabs(erre1) < 1.e-10){ if( Fi0 > fi2 || Fi0 < fi1){ continue; } } else if ( erre1 0. ) { continue; } } else if (erre1> Rout && erre2 > Rout && !( fi1<= aaa && aaa<=fi2 && ddd<=Rout ) ) { continue; } for(l=itemp-2; l<=itemp+2; l++){ if( l< 0 ) { l2 = l+NFIDIVCONFORMAL; } else if (l>=NFIDIVCONFORMAL){ l2 = l- NFIDIVCONFORMAL*( l/NFIDIVCONFORMAL ); } else { l2 = l; } if( j-1<0) { kstart=0; } else { kstart = j-1; } if ( j+1 >= NRDIVCONFORMAL ) { kend = NRDIVCONFORMAL; } else { kend = j+2; } for(k=kstart;k NTIMES*STRAWRADIUS ) continue; //------------------- xx=infoparalConformal[ HitsinBoxConformal[l3][k][l2] ][0]; yy=infoparalConformal[ HitsinBoxConformal[l3][k][l2] ][1]; dist = fabs( -yy+ m*xx +q )/sqrt(m*m+1.); if( AcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[l3][k][l2] ][2], infoparalConformal[ HitsinBoxConformal[l3][k][l2] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[l3][k][l2]; Unselected[HitsinBoxConformal[l3][k][l2]]= false; nAssociatedHits++; } } // end of for( l3=0;l3=NFIDIVCONFORMAL){ i -= NFIDIVCONFORMAL*( i/NFIDIVCONFORMAL ); } for(l=0; l NTIMES*STRAWRADIUS ) continue; //------------------- xx=infoparalConformal[ HitsinBoxConformal[k][l][i] ][0]; yy=infoparalConformal[ HitsinBoxConformal[k][l][i] ][1]; dist = fabs( m*xx-yy )/sqrt( m*m+1.); if( AcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[k][l][i] ][2], infoparalConformal[ HitsinBoxConformal[k][l][i] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[k][l][i]; nAssociatedHits++; } } // end of for( k=0;k=3) { cout<<"from TrkAssociatedParallelHitsToHelixQuater, before exiting; nAssociatedHits = " < NTIMES*STRAWRADIUS ) continue; if(angleFi_up) continue; auxListHitsinTrack[nAssociatedHits]= i; nAssociatedHits++; } // end for(i=0; i=2){ cout<<"\n\n---------------------------------------------\n"; cout<<" 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. "<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 PndTracking::WriteMacroAllHitsRestanti //----------start of function PndTracking::WriteMacroParallelHitsGeneral void PndTracking::WriteMacroParallelHitsGeneral( Int_t Nhits, Double_t info[][7], UShort_t nTracksFoundSoFar, bool *keepit, 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]; } 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); //------------------------------------------------------------------------------------------------------------ if(!doMcComparison) return; //---------- parallel straws Macro now con anche le tracce MC 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; 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]-Oyy, POINTS1[j]-Oxx) ; // 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(Oxx,Oyy,Rr,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]-Oxx ; // x component Radial vector of cylinder of trajectory Ry = POINTS1[1+j]-Oyy ; // y direction Radial vector of cylinder of trajectory aaa = sqrt(Rx*Rx+Ry*Ry); SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ; SkewInclWithRespectToS /= Rr; 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]/Rr; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oyy, POINTS1[j]-Oxx) ; // 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 ]-Oyy, XMvdPixel[ ii ]-Oxx); 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 ]-Oyy, XMvdStrip[ ii ]-Oxx); 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 ]-Oyy, XMvdPixel[ ii ]-Oxx); 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 ]-Oyy, XMvdStrip[ ii ]-Oxx); 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,Rr*(Smin-.1*deltaS),zmax+0.1*deltaz,Rr*(Smax+0.1*deltaS)); 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.01*deltaz,Rr*(Smin+0.05*deltaS),zmax+0.01*deltaz,Rr*(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,Rr*(Smin-0.01*deltaS),zmin+0.05*deltaz,Rr*(Smax+0.01*deltaS), Rr*(Smin-0.01*deltaS),Rr*(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]/Rr; if( distance >= info[i][4] + Aellipsis1) continue; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oyy, POINTS1[j]-Oxx) ; // 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],Rr*fi1,Aellipsis1,Rr*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(Oxx,Oyy,Rr,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]-Oxx ; // x component Radial vector of cylinder of trajectory Ry = POINTS1[1+j]-Oyy ; // y direction Radial vector of cylinder of trajectory aaa = sqrt(Rx*Rx+Ry*Ry); SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ; SkewInclWithRespectToS /= Rr; 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]/Rr; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oyy, POINTS1[j]-Oxx) ; // 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],Rr*fi1,Aellipsis1,Rr*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 ]-Oyy, XMvdPixel[ ii ]-Oxx); 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],Rr*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],Rr*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 ]-Oyy, XMvdStrip[ ii ]-Oxx); 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],Rr*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],Rr*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 ]-Oyy, XMvdPixel[ ii ]-Oxx); 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],Rr*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 ]-Oyy, XMvdStrip[ ii ]-Oxx); 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],Rr*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, Rr*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 r_r, Dd, Fifi, Kakka, o_x, o_y, Cx, Cy, Px, Py, carica ; PndMCTrack* pMC; pMC = (PndMCTrack*) fMCTrackArray->At(imc); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track o_x = pMC->GetStartVertex().X(); // X of starting point track o_y = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); r_r = 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 = o_x + Py*1000./(BFIELD*CVEL*carica); Cy = o_y - 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,Rr* 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 PndTracking::WriteMacroSkewAssociatedHitswithMC //--start of function PndTracking::WriteMacroSttParallelAssociatedHitsandMvdwithMC void PndTracking::WriteMacroSttParallelAssociatedHitsandMvdwithMC( Double_t Oxx,Double_t Oyy,Double_t Rr, Double_t primoangolo, Double_t ultimoangolo, UShort_t Nhits, UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], Double_t info[][7], UShort_t iTrack, Int_t iNome, Short_t daSttTrackaMCTrack, UShort_t nParalCommon[MAXTRACKSPEREVENT], UShort_t ParalCommonList[MAXMCTRACKS][MAXSTTHITSINTRACK], UShort_t nSpuriParinTrack[MAXTRACKSPEREVENT], UShort_t ParSpuriList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], UShort_t nMCParalAlone[MAXTRACKSPEREVENT], UShort_t MCParalAloneList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], UShort_t nMvdPixelHitsAssociatedToSttTra, UShort_t ListPixelHitsinTrack[MAXTRACKSPEREVENT][MAXMVDPIXELHITSINTRACK], UShort_t nMvdStripHitsAssociatedToSttTra, UShort_t ListStripHitsinTrack[MAXTRACKSPEREVENT][MAXMVDSTRIPHITSINTRACK], UShort_t nMvdPixelCommon, UShort_t *MvdPixelCommonList, UShort_t nMvdPixelSpuriinTrack, UShort_t *MvdPixelSpuriList, UShort_t nMCMvdPixelAlone, UShort_t *MCMvdPixelAloneList, UShort_t nMvdStripCommon, UShort_t *MvdStripCommonList, UShort_t nMvdStripSpuriinTrack, UShort_t *MvdStripSpuriList, UShort_t nMCMvdStripAlone, UShort_t *MCMvdStripAloneList, UShort_t nSkewHitsinTrack, UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], Double_t *SchosenSkew, UShort_t nSkewCommon[MAXTRACKSPEREVENT], UShort_t SkewCommonList[MAXTRACKSPEREVENT][MAXSTTHITSINTRACK], UShort_t nMCSkewAlone[MAXTRACKSPEREVENT], UShort_t MCSkewAloneList[MAXMCTRACKS][MAXSTTHITSINTRACK] ) { Int_t i, j, i1, ii, index, Kincl, nlow, nup, STATUS; Double_t xmin , xmax, ymin, ymax, dx, dy, diff, d1, d2, delta, deltax, deltay, deltaz, deltaS, factor, zmin, zmax, Smin, Smax, S1, S2, z1, z2, y1, y2,x1,x2, 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, Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; //---------- parallel straws Macro now char nome[300], nome2[300]; sprintf(nome,"MacroSttMvdXYwithMCEvent%dT%d", IVOLTE,iNome); sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); fprintf(MACRO,"{\n"); xmin=1.e20; xmax=-1.e20; ymin=1.e20; ymax=-1.e20; //---- Scitil hits. for( ii=0; ii< nSciTilHitsinTrack[iTrack]; ii++) { i = ListSciTilHitsinTrack[iTrack][0] ; 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]; } //------------- 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=Oxx+Rr*cos(SchosenSkew[ii]); bbb=Oyy+Rr*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" ,Oxx,Oyy,Rr,Rr,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=Oxx+Rr*cos(SchosenSkew[ii]); bbb=Oyy+Rr*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( 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]; 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 r_r, Dd, Fifi, o_x, o_y, Cx, Cy, Px, Py, carica ; PndMCTrack* pMC; im=daSttTrackaMCTrack; pMC = (PndMCTrack*) fMCTrackArray->At(im); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track o_x = pMC->GetStartVertex().X(); // X of starting point track o_y = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); r_r = 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 = o_x + Py*1000./(BFIELD*CVEL*carica); Cy = o_y - 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,r_r,r_r,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 PndTracking::WriteMacroSttParallelAssociatedHitsandMvdwithMC ClassImp(PndTracking)