// ------------------------------------------------------------------------- // ----- PndSttTrackFinderReal source file ----- // ----- by Gianluigi Boca ----- // ------------------------------------------------------------------------- #include "glpk.h" // Pnd includes #include "PndSttTrackFinderReal.h" #include "PndSttHit.h" #include "PndSciTHit.h" #include "PndSttPoint.h" #include "PndSttHelixHit.h" #include "PndTrackCand.h" #include "PndTrack.h" #include "PndDetectorList.h" #include "PndSttTube.h" #include #include "FairTrackParP.h" #include "FairMCPoint.h" #include "FairRootManager.h" // ROOT includes #include "TClonesArray.h" #include "TDatabasePDG.h" #include "TRandom.h" #include "TH1.h" #include "TH2.h" #include "TH3.h" // C++ includes #include #include using std::cout; using std::cin; using std::endl; using std::map; // ----- Default constructor ------------------------------------------- PndSttTrackFinderReal::PndSttTrackFinderReal() { doMcComparison = false; Fimin=0.; Fimax=2.*PI; FI0min = 0.; FI0max = 2.*PI; iplotta = false; istampa = 0; MINIMUMOUTERHITSPERTRACK=5; nSciTilHits=0; stepD=(Dmax-Dmin)/nbinD; stepFi=(Fimax-Fimin)/nbinFi; stepR=(Rmax-Rmin)/nbinR; stepKAPPA=(KAPPAmax-KAPPAmin)/nbinKAPPA; stepFI0=(FI0max-FI0min)/nbinFI0; stepfineKAPPA=2.*DELTA_KAPPA/nbinKAPPA; stepfineFI0=2.*DELTA_FI0/nbinFI0; YesSciTil = false ; Initialization_ClassVariables(); sprintf(fSttBranch,"STTHit"); } // ------------------------------------------------------------------------- // ----- Standard constructor ------------------------------------------ PndSttTrackFinderReal::PndSttTrackFinderReal(int verbose) { doMcComparison = false; Fimin=0.; Fimax=2.*PI; FI0min = 0.; FI0max = 2.*PI; iplotta = false; istampa = verbose; MINIMUMOUTERHITSPERTRACK=5; nSciTilHits=0; stepD=(Dmax-Dmin)/nbinD; stepFi=(Fimax-Fimin)/nbinFi; stepR=(Rmax-Rmin)/nbinR; stepKAPPA=(KAPPAmax-KAPPAmin)/nbinKAPPA; stepFI0=(FI0max-FI0min)/nbinFI0; stepfineKAPPA=2.*DELTA_KAPPA/nbinKAPPA; stepfineFI0=2.*DELTA_FI0/nbinFI0; YesSciTil = false ; Initialization_ClassVariables(); sprintf(fSttBranch,"STTHit"); } // ------------------------------------------------------------------------- // ----- Second constructor ------------------------------------------ PndSttTrackFinderReal::PndSttTrackFinderReal(int istamp, bool iplott, bool imc) { doMcComparison = imc; Fimin=0.; Fimax=2.*PI; FI0min = 0.; FI0max = 2.*PI; iplotta = iplott; istampa= istamp; MINIMUMOUTERHITSPERTRACK=5; nSciTilHits=0; stepD=(Dmax-Dmin)/nbinD; stepFi=(Fimax-Fimin)/nbinFi; stepR=(Rmax-Rmin)/nbinR; stepKAPPA=(KAPPAmax-KAPPAmin)/nbinKAPPA; stepFI0=(FI0max-FI0min)/nbinFI0; stepfineKAPPA=2.*DELTA_KAPPA/nbinKAPPA; stepfineFI0=2.*DELTA_FI0/nbinFI0; YesSciTil = false ; Initialization_ClassVariables(); sprintf(fSttBranch,"STTHit"); } // ------------------------------------------------------------------------- // ----- Third constructor ------------------------------------------ PndSttTrackFinderReal::PndSttTrackFinderReal(int istamp, bool iplott, bool imc, bool doSciTil) { doMcComparison = imc; Fimin=0.; Fimax=2.*PI; FI0min = 0.; FI0max = 2.*PI; iplotta = iplott; istampa= istamp; MINIMUMOUTERHITSPERTRACK=5; nSciTilHits=0; stepD=(Dmax-Dmin)/nbinD; stepFi=(Fimax-Fimin)/nbinFi; stepR=(Rmax-Rmin)/nbinR; stepKAPPA=(KAPPAmax-KAPPAmin)/nbinKAPPA; stepFI0=(FI0max-FI0min)/nbinFI0; stepfineKAPPA=2.*DELTA_KAPPA/nbinKAPPA; stepfineFI0=2.*DELTA_FI0/nbinFI0; YesSciTil = doSciTil ; Initialization_ClassVariables(); sprintf(fSttBranch,"STTHit"); } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndSttTrackFinderReal::~PndSttTrackFinderReal() { } // ------------------------------------------------------------------------- // -------------------- begin PndSttTrackFinderReal::Initialization_ClassVariables void PndSttTrackFinderReal::Initialization_ClassVariables() { // this is only for initializing the Class Variables. char zero; size_t len; // booleans : len = sizeof(InclusionListSciTil); memset (InclusionListSciTil,0,len); len = sizeof(TypeConf); memset (TypeConf,0,len); // int : IVOLTE=-1; // Short_t : nRdivConformalEffective=0; nSciTilHits=0; nSttSkewhit=0; nMCTracks=0; len = sizeof(infoparal); memset(infoparal,0,len); len = sizeof(infoskew); memset(infoskew,0,len); len = sizeof(nHitsInMCTrack); memset(nHitsInMCTrack,0,len); len = sizeof(nSciTilHitsinTrack); memset(nSciTilHitsinTrack,0,len); len = sizeof(nSttSkewhitInMCTrack); memset(nSttSkewhitInMCTrack,0,len); len = sizeof(ListSciTilHitsinTrack); memset(ListSciTilHitsinTrack,0,len); // Double_t : Fimax=0.; FI0min=0.; FI0max=0.; stepD=0.; stepFi=0.; stepR=0.; stepKAPPA=0.; stepFI0=0.; stepfineKAPPA=0.; stepfineFI0=0.; SEMILENGTH_STRAIGHT=0.; ZCENTER_STRAIGHT=0.; len = sizeof(veritaMC); memset (veritaMC,0,len); len = sizeof(ALFA); memset (ALFA,0,len); len = sizeof(BETA); memset (BETA,0,len); len = sizeof(GAMMA); memset (GAMMA,0,len); len = sizeof(radiaConf); memset (radiaConf,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(posizSciTil); memset (posizSciTil,0,len); len = sizeof(S_SciTilHitsinTrack); memset (S_SciTilHitsinTrack,0,len); // pointers : hdist=NULL; hdistgoodlast=NULL, hdistbadlast=NULL; HANDLE=NULL; HANDLE2=NULL; HANDLEXYZ=NULL; PHANDLEX=NULL; PHANDLEY=NULL; PHANDLEZ=NULL; SHANDLEX=NULL; SHANDLEY=NULL; SHANDLEZ=NULL; fMCTrackArray=NULL; fSciTPointArray=NULL; fSciTHitArray=NULL; pMCtr=NULL; fSttHitArray=NULL; } // -------------------- end of PndSttTrackFinderReal::Initialization_ClassVariables //----------begin of function PndSttTrackFinderReal::WriteHistograms void PndSttTrackFinderReal::WriteHistograms(){ TFile* file = FairRootManager::Instance()->GetOutFile(); file->cd(); file->mkdir("PndSttTrackFinderReal"); file->cd("PndSttTrackFinderReal"); if(iplotta){ hdist->Write(); hdistgoodlast->Write(); hdistbadlast->Write(); delete hdist; delete hdistgoodlast; delete hdistbadlast; } } //----------end of function PndSttTrackFinderReal::WriteHistograms // ----- Public method Init -------------------------------------------- void PndSttTrackFinderReal::Init() { Short_t i; Double_t r1, r2, A , tempRadiaConf[nRdivConformal]; MINIMUMHITSPERTRACK=3; // --------------------------- 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("statistichePndTrackFinderReal.txt","w"); // --------------- if(istampa >=3 ) 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) // ------------------------- fHelixHitProduction = true; // IVOLTE=-1; // Get and check FairRootManager FairRootManager* ioman = FairRootManager::Instance(); if (!ioman) { cout << "-E- PndSttTrackFinderReal::Init: " << "RootManager not instantiated, return!" << endl; return; } // get the MCTrack array fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack"); // get the SciTil point array // fSciTPointArray = (TClonesArray*) ioman->GetObject("SciTPoint"); // get the SciTil hit array if(YesSciTil) { fSciTHitArray = (TClonesArray*) ioman->GetObject("SciTHit"); } else { fSciTHitArray = NULL; } // ------------------------------------------------------------------- if(iplotta){ hdist = new TH1F("hdist", "distance (cm)", 20, 0., 10.); hdistgoodlast = new TH1F("hDistanceTrulyInnerLast", "distance (cm)", 20, 0., 10.); hdistbadlast = new TH1F("hDistanceNonLastInner", "distance (cm)", 20, 0., 10.); } // calculate the boundaries of the Box in Conformal Space, see Gianluigi logbook on pag. 210-211 radiaConf[0] = 1./RStrawDetectorMax; r1 = RStrawDetectorMin; A = (RStrawDetectorMax - r1)/nRdivConformal; if ( nRdivConformal > 1 ) { for(i = 1; i< nRdivConformal ; i++){ r2 = r1 + A; // tempRadiaConf[nRdivConformal-i] = 1./r2; radiaConf[nRdivConformal-i] = 1./r2; r1=r2; } } nRdivConformalEffective = nRdivConformal; //-------------------- end of method Init -------------------------------------------- } // -------------------- starting PndSttTrackFinderReal::GetHitFromCollections PndSttHit* PndSttTrackFinderReal::GetHitFromCollections(Int_t hitCounter) { PndSttHit *retval = NULL; Int_t relativeCounter = hitCounter; for (Int_t collectionCounter = 0; collectionCounter < fHitCollectionList.GetEntries(); collectionCounter++) { Int_t size = ((TClonesArray *)fHitCollectionList.At(collectionCounter))->GetEntriesFast(); if (relativeCounter < size) { retval = (PndSttHit*) ((TClonesArray *)fHitCollectionList.At(collectionCounter))->At(relativeCounter); break; } else { relativeCounter -= size; } } return retval; } // -------------------- end of PndSttTrackFinderReal::GetHitFromCollections // -------------------- starting PndSttTrackFinderReal::GetPointFromCollections FairMCPoint* PndSttTrackFinderReal::GetPointFromCollections(Int_t hitCounter) { FairMCPoint *retval = NULL; Int_t relativeCounter = hitCounter; for (Int_t collectionCounter = 0; collectionCounter < fHitCollectionList.GetEntries(); collectionCounter++) { Int_t size = ((TClonesArray *)fHitCollectionList.At(collectionCounter))->GetEntriesFast(); if (relativeCounter < size) { Int_t tmpHit = ((PndSttHit*) ((TClonesArray *)fHitCollectionList.At(collectionCounter))->At(relativeCounter))->GetRefIndex(); retval = (FairMCPoint*) ((TClonesArray *)fPointCollectionList.At(collectionCounter))->At(tmpHit); break; } else { relativeCounter -= size; } } return retval; } // -------------------- end of PndSttTrackFinderReal::GetPointFromCollections Int_t PndSttTrackFinderReal::DoFind( TClonesArray* , TClonesArray* ) {return 0;} // CHECK da cancellare // ----------------- start of method DoFind Int_t PndSttTrackFinderReal::DoFind(TClonesArray* trackCandArray, TClonesArray *trackArray, TClonesArray* helixHitArray) { bool flaggo, outcome; Int_t i, j, iaccept, ii,jj,k,kk,n1,n2,n3, i1,j1,k1,imaxima, jmaxima, iofmax, jofmax, kofmax, index[nmaxHits], integer, Ncirc, nSkew1, nSkew2, n_K, n_FI0, n_R, n_D, n_Fi,STATUS, Nremaining, Nremaining2, NumberofMaximaDFiR, NumberofMaximaKFI0, nAssociatedParallelHits ; Short_t emme, exitstatus, inclination_type, istep, // nSciTilHitsinTrack[MAXTRACKSPEREVENT], auxIndex[nmaxHits], OLDinfoparal[nmaxHits]; Short_t Charge[MAXTRACKSPEREVENT], Status[MAXTRACKSPEREVENT], daTrackFoundaTrackMC[MAXTRACKSPEREVENT], enne[MAXTRACKSPEREVENT][nmaxHits]; Double_t aaa, ddd, delta, deltabis, deltaZ, mindis, distanza, fi_hit, ap1, ap2, ap3, carica, cross1, cross2, cross3, dummy, lowlimit[MAXTRACKSPEREVENT], uplimit[MAXTRACKSPEREVENT], info[nmaxHits][7], WDX, WDY, WDZ, auxRvalues[nmaxHits], inclination[nmaxinclinationversors][3]; PndSttHit* ListPointer_to_Hit[nmaxHits]; inclination[0][0]=inclination[0][1]=0., inclination[0][2]=1.; Int_t Nincl = 1, Ninclinate, iHit; Int_t Minclinations[nmaxinclinationversors]; bool InclusionList[nmaxHits], // list of parallel hits already assigned to a track or with multiple hits. InclusionListSkew[nmaxHits], // list of skew hits with multiple hits InclusionListbis[nmaxHits], // list of || hits ONLY with multiple hits InclusionListSkewbis[nmaxHits], // list of skew hits ONLY with multiple hits (duplicate of // InclusionListSkew) // TypeConf[MAXTRACKSPEREVENT], // if TypeConf[]=false --> the track is a line in the Conformal space, // if TypeConf[]=true it is a crf; TypeConfSkew[MAXTRACKSPEREVENT], GoodTrack[MAXTRACKSPEREVENT]; // yes = track found passes minimal requirements Short_t iParHit, nFicell, nRcell; Short_t exphit, iExclude, imc, jexp, mchit, nTracksFoundSoFar, NNN, Nouter, nTotalHits[MAXTRACKSPEREVENT], TemporarynSttSkewhitinTrack, BigList[MAXTRACKSPEREVENT][nmaxHitsInTrack], TemporarySkewList[nmaxHits][2], nHitsinTrack[MAXTRACKSPEREVENT], nSttSkewhitinTrack[MAXTRACKSPEREVENT], tempore[nmaxHits], ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHitsInTrack], ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxHitsInTrack], nMCParalAlone[MAXTRACKSPEREVENT], nMCSkewAlone[MAXTRACKSPEREVENT], MCParalAloneList[MAXTRACKSPEREVENT][nmaxHits], MCSkewAloneList[MAXTRACKSPEREVENT][nmaxHits], nParalCommon[MAXTRACKSPEREVENT], ParalCommonList[MAXTRACKSPEREVENT][nmaxHits], nSpuriParinTrack[MAXTRACKSPEREVENT], ParSpuriList[MAXTRACKSPEREVENT][nmaxHits], nSkewCommon[MAXTRACKSPEREVENT], SkewCommonList[MAXTRACKSPEREVENT][nmaxHits], nSpuriSkewinTrack[MAXTRACKSPEREVENT], SkewSpuriList[MAXTRACKSPEREVENT][nmaxHits], RConformalIndex[nmaxHits], // given a Hit number it gives its radial box number FiConformalIndex[nmaxHits], // given a Hit number it gives its azimuthal box number nBoxConformal[nRdivConformal][nFidivConformal]; // first index -> radial divisions, 2nd index -> azimuthal divisions; n. of // hits falling in this cell Short_t HitsinBoxConformal[MAXHITSINCELL][nRdivConformal][nFidivConformal]; // first index -> radial divisions, 2nd index -> azimuthal divisions; Int_t status; int iimax,jjmax, ncount; Double_t angle, Distance, max1, HoughR, HoughD, HoughFi, HoughKAPPA, HoughFI0, Px, Py, tempR, tempD,tempFi,tempKAPPA,tempFI0, Rlow, Rup, Dlow, Dup,Filow,Fiup, KAPPAlow, KAPPAup, FI0low, FI0up, gamma; Double_t D, Fi, ptotal, Fi_low_limit[MAXTRACKSPEREVENT], Fi_up_limit[MAXTRACKSPEREVENT], Fi_initial_helix_referenceframe[MAXTRACKSPEREVENT], Fi_final_helix_referenceframe[MAXTRACKSPEREVENT], R[MAXTRACKSPEREVENT], Ox[MAXTRACKSPEREVENT], Oy[MAXTRACKSPEREVENT], KAPPA[MAXTRACKSPEREVENT], FI0[MAXTRACKSPEREVENT], trajectory_vertex[3], infoparalConformal[nmaxHits][5], S[2*nmaxHits], tmpErrorZDrift[nmaxHitsInTrack], tmpS[nmaxHitsInTrack], tmpZ[nmaxHitsInTrack], tmpZDrift[nmaxHitsInTrack], Z[2*nmaxHits], temporeS[nmaxHits], temporeZ[nmaxHits], ZDrift[2*nmaxHits], ZErrorafterTilt[2*nmaxHits], temporeZDrift[nmaxHits], temporeZErrorafterTilt[nmaxHits], Posiz[3], Posiz1[3], Posiz2[3], versor[3], // U[MAXTRACKSPEREVENT][nmaxHits], // V[MAXTRACKSPEREVENT][nmaxHits], Xpos_for_LHeTrack[nmaxHits], Ypos_for_LHeTrack[nmaxHits], Zpos_for_LHeTrack[nmaxHits], Sfinal[MAXTRACKSPEREVENT][nmaxHits], Zfinal[MAXTRACKSPEREVENT][nmaxHits], ZDriftfinal[MAXTRACKSPEREVENT][nmaxHits], ZErrorafterTiltfinal[MAXTRACKSPEREVENT][nmaxHits]; Float_t RemainingR[nmaxHits], RemainingFi[nmaxHits], RemainingD[nmaxHits], RemainingCX[nmaxHits], RemainingCY[nmaxHits], RemainingKAPPA[nmaxHits], RemainingFI0[nmaxHits]; bool GoodSkewFit[nmaxHits]; bool temporflag[8], IFLAG; CalculatedCircles Result; CalculatedHelix HelixResult; // AssociatedHitsToHelix ResultAssociatedHits ; char nomef[100], nome[30],titolo[100]; //------------------------------------ modifiche Gianluigi, 9-7-08 IVOLTE++; if(istampa>0) cout<<"\nEntering in PndSttTrackFinderReal : evt (starting from 0) n. "<GetEntriesFast(); // num. tracce/evento if (nMCTracks ==0){ cout<<"da PndSttTrackFinderReal : N. di MC truth tracks = 0, return -1!\n" < MAXMCTRACKS){ cout<<"da PndSttTrackFinderReal : N. di MC truth tracks = "< MAXMCTRACKS = "<=1) cout<<"Gianluigi, evento n. "< nmaxHits ) { cout<<"da PndSttTrackFinderReal : N. di Stt Hits = "< nmaxHits (="<GetEntriesFast(); if(istampa>0) cout<<"da PndSttTrackFinderReal, 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 PndSttTrackFinderReal, evt "<At(j); posiz = pPndSciTHit->GetPosition(); if(istampa>0)cout<<"da PndSttTrackFinderReal, evt. "<0){ //-----------stampe. if(istampa>0){ cout<<"da PndSttTrackFinderReal, dopo purga di SciTil; evt. "< > hitMap; for(j=0; j= 1 ) cout<<"Gianluigi : da PndSttTrackFinderReal::DoFind : Nhits="<GetRefIndex(); if (ptIndex >= 0) { pMCpt = GetPointFromCollections(iHit); // <== FairMCPoint } // tubeID CHECK added Int_t tubeID = pMhit->GetTubeID(); PndSttTube *tube = (PndSttTube*) fTubeArray->At(tubeID); // real hit center of tube coordinates TVector3 center = tube->GetPosition(); // drift radius Double_t dradius = pMhit->GetIsochrone(); // wire direction TVector3 wiredirection = tube->GetWireDirection(); // "real" MC coordinates (in + out)/2. // TVector3 mcpoint; // pMCpt->Position(mcpoint); if(wiredirection.Z() >=0.) { WDX = wiredirection.X(); WDY = wiredirection.Y(); WDZ = wiredirection.Z(); } else { WDX = -wiredirection.X(); WDY = -wiredirection.Y(); WDZ = -wiredirection.Z(); } // stampe di controllo info[iHit][0]= tube->GetPosition().X(); info[iHit][1]= tube->GetPosition().Y(); info[iHit][2]= tube->GetPosition().Z(); info[iHit][3]= dradius; info[iHit][4]= tube->GetHalfLength(); if (ptIndex >= 0) { info[iHit][6]= pMCpt->GetTrackID(); }else{ info[iHit][6]= -20.; } if( fabs( WDX )< 0.00001 && fabs( WDY )< 0.00001 ){ info[iHit][5]= 1.; infoparal[Minclinations[0]]= iHit ; Minclinations[0]++; ZCENTER_STRAIGHT = info[iHit][2]; // this works because just few lines below there is the SEMILENGTH_STRAIGHT = info[iHit][4]; // requirement that Minclinations[0] > 2 (= at least 3 parallel straws) } else { flaggo=true; for (i=2; i<=Nincl;i++) { if (fabs( WDX-inclination[i-1][0] )< 0.00001 && fabs( WDY -inclination[i-1][1])< 0.00001 && fabs( WDZ -inclination[i-1][2])< 0.00001 ){ info[iHit][5]= i; infoskew[Ninclinate]= iHit; Ninclinate++; Minclinations[i-1]++; flaggo=false; break; } } if(flaggo){ Nincl++; inclination[Nincl-1][0]=(Double_t) WDX; inclination[Nincl-1][1]=(Double_t) WDY; inclination[Nincl-1][2]=(Double_t) WDZ; info[iHit][5]= Nincl; infoskew[Ninclinate]= iHit; Ninclinate++; Minclinations[Nincl-1]++; } // end if(flaggo) } nSttSkewhit = Ninclinate; // calcoli validi solo per il MC ----------------------------- if (ptIndex >= 0) { veritaMC[iHit][0]= pMCpt->GetX(); veritaMC[iHit][1]= pMCpt->GetY(); veritaMC[iHit][2]= pMCpt->GetZ(); } // FromHitToMCTrack[iHit] = (Short_t) ( info[iHit][6] + 0.001 ); //--------------- inizio stampaggi, stampe di controllo if (istampa >= 2 && IVOLTE<= nmassimo) { cout <<"iHit "<< iHit << endl; if (ptIndex < 0) { cout<<"from PndSttTrackFinderReal...this hit must be noise (RefIndex = "<GetTrackID()<GetPosition().X() << " " << tube->GetPosition().Y() << " " << tube->GetPosition().Z() << "; R = "<GetPosition().X()*tube->GetPosition().X()+tube->GetPosition().Y()*tube->GetPosition().Y())<< endl; cout <<" wire direction, X, Y, Z (Z direction set always positive)" << WDX<<" "<= //-------- fine stampaggi } // end of for (Int_t iHit = 0; //--------------- inizio stampaggi // if (istampa >= 2 && IVOLTE<= nmassimo) { if (istampa >= 2 ) { cout<<"Gianluigi : da PndSttTrackFinderReal::DoFind : Nhits totali =" < MAXTRACKSPEREVENT) continue; nRcell = -1; // because SciTil hit is outside of the Stt system. Fi = atan2(posizSciTil[i][1],posizSciTil[i][0]) ; if ( Fi < 0. ) Fi += 2.*PI; nFicell = (Short_t) (0.5*nFidivConformal*Fi/PI); if(nFicell > nFidivConformal ) { nFicell = nFidivConformal; } else if (nFicell<0) { nFicell = 0; } outcome = FindTrackInXYProjection( -i-1, // seed hit; it is negative for SciTil Hits. nRcell, nFicell, Minclinations, info, InclusionList, RConformalIndex, FiConformalIndex, nBoxConformal, HitsinBoxConformal, nTracksFoundSoFar, nHitsinTrack, ListHitsinTrack, trajectory_vertex, infoparalConformal, posizSciTil[i][0], posizSciTil[i][1], &S_SciTilHitsinTrack[nTracksFoundSoFar][0], Ox, Oy, R, Fi_low_limit, Fi_up_limit, Fi_initial_helix_referenceframe, Fi_final_helix_referenceframe, Charge, &U[nTracksFoundSoFar][0], &V[nTracksFoundSoFar][0] ); if(!outcome){ continue; } for(j=0; j0){ cout<<"PndTrackFinderReal, evt. "<= MAXTRACKSPEREVENT ){ cout<<"from PndSttTrackFinderReal : # n. Tracks found so far = "<= MAXTRACKSPEREVENT ( = "< MAXTRACKSPEREVENT) continue; if( ! InclusionList[ infoparal[iParHit] ] ) continue; nRcell = RConformalIndex[ infoparal[iParHit] ]; nFicell = FiConformalIndex[infoparal[iParHit]]; outcome = FindTrackInXYProjection( iParHit, // seed hit; it is positive for STT Hits. nRcell, nFicell, Minclinations, info, InclusionList, RConformalIndex, FiConformalIndex, nBoxConformal, HitsinBoxConformal, nTracksFoundSoFar, nHitsinTrack, ListHitsinTrack, 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 &dummy, Ox, Oy, R, Fi_low_limit, Fi_up_limit, Fi_initial_helix_referenceframe, Fi_final_helix_referenceframe, Charge, &U[nTracksFoundSoFar][0], &V[nTracksFoundSoFar][0] ); if(!outcome) continue; // -------- here the track and its hits were found, filling the Inclusion list for(j=0; j0){ cout<<"PndTrackFinderReal, evt. "<= MAXTRACKSPEREVENT ){ cout<<"from PndSttTrackFinderReal : # n. Tracks found so far = "<= MAXTRACKSPEREVENT ( = "<=2){ for(j=0;j nmaxHitsInTrack ) { if(nmaxHitsInTrack-nHitsinTrack[i]>0)nSttSkewhitinTrack[i]=nmaxHitsInTrack-nHitsinTrack[i]; else nSttSkewhitinTrack[i]=0; } // here there are up to 2 SciTil hits in track. if(nSciTilHitsinTrack[i]>0){ for(j=0;j SciTil hit. tmpErrorZDrift, Fi_initial_helix_referenceframe[i], // this is an input; NHITSINFIT, // maximum n. hits in fit. &KAPPA[i] ); // Status[i] is negative (-99) when m = 0. and (-100) as result when the fit with glpk // failed. if(Status[i] < 0 || fabs(KAPPA[i])>1.e10) { // necessary to load here the Sfinal vector anyway. for(j=0;j=0 ){ nSttSkewhitinTrack[i] = NNN; // limit the total # of hits in track to nmaxHitsInTrack. if( nSttSkewhitinTrack[i]+nHitsinTrack[i] > nmaxHitsInTrack ) { if(nmaxHitsInTrack-nHitsinTrack[i]>0)nSttSkewhitinTrack[i]=nmaxHitsInTrack-nHitsinTrack[i]; else nSttSkewhitinTrack[i]=0; } for(j=0;j=0 ) nSttSkewhitinTrack[i]=TemporarynSttSkewhitinTrack; // limit the total # of hits in track to nmaxHitsInTrack. if( nSttSkewhitinTrack[i]+nHitsinTrack[i] > nmaxHitsInTrack ) { if(nmaxHitsInTrack-nHitsinTrack[i]>0)nSttSkewhitinTrack[i]=nmaxHitsInTrack-nHitsinTrack[i]; else nSttSkewhitinTrack[i]=0; } for(j=0;j=0 ) if( nHitsinTrack[i]+nSttSkewhitinTrack[i] < MINIMUMHITSPERTRACK || nHitsinTrack[i]+nSttSkewhitinTrack[i]>nmaxHitsInTrack) { keepit[i]=false; continue; } // --------- numbering according to the ORIGINAL hit number for(i1=0; i1< nSttSkewhitinTrack[i]; i1++){ Sfinal[i][ infoskew[ ListSkewHitsinTrack[i][i1] ] ] = S[i1] ; Zfinal[i][ infoskew[ ListSkewHitsinTrack[i][i1] ] ] = Z[i1] ; ZDriftfinal[i][ infoskew[ ListSkewHitsinTrack[i][i1] ] ] = ZDrift[i1] ; ZErrorafterTiltfinal[i][ infoskew[ ListSkewHitsinTrack[i][i1] ] ] = ZErrorafterTilt[i1] ; } // --------- // ------------------------------------------------------------- } // end of for(i=0; i=2) { cout<<"Evt. n. "<=2){cout<<"\tPndSttTrackFinderReal, fine di procedura, IVOLTE = " <0 && nTracksFoundSoFar > 0 && doMcComparison ) { AssociateFoundTrackstoMCbis( keepit, info, nTracksFoundSoFar, nHitsinTrack, ListHitsinTrack, nSttSkewhitinTrack, ListSkewHitsinTrack, daTrackFoundaTrackMC ); for(jexp=0; jexpAt(i); if ( ! pMCtr ) continue; Double_t 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 PndSttTrackFinderReal, evento (cominciando da 0) n. "<=1 ){ fprintf(HANDLE, "\n Evento %d NTotaleTracceMC %d ------\n", IVOLTE, nMCTracksaccettabili); int ibene=0; if(nMCTracksaccettabili>0){ for(ii=0; ii0) fprintf(HANDLE,"\tn. volte almeno 1 traccia MC accettabile e' ricostruita %d\n" ,ibene); }// end of if(istampa>=1 ) //---------------- for (ii=0; ii=1 ;ii++){ if(!keepit[ii]) continue; fprintf(HANDLE,"----------------------------------------------------------\n"); i=daTrackFoundaTrackMC[ii]; if( i <0 ) { fprintf(HANDLE, " No TracciaMC associated to found track n. %d in pattern recognition, with %d Hits ||, %d skew hits, %f Radius \n " ,ii,nHitsinTrack[ii], nSttSkewhitinTrack[ii], R[ii] ); continue; } if( ! ( GoodSkewFit[ii] ) ) { fprintf(HANDLE, " TracciaMC %d; sua FoundTrack associata (n. %d) NONsoddisfaRequisitiMinimi perche' in Z-S KAPPA e' risultata || asse S\n",i,ii); continue; } if(fabs(KAPPA[ii])<1.e-10 ){ fprintf(HANDLE, " TracciaMC %d; sua FoundTrack associata (n. %d) NONsoddisfaRequisitiMinimi perche' KAPPA troppo piccolo; KAPPA = %g\n", i,ii,KAPPA[ii]); continue; // } else if (fabs(KAPPA[ii])>1.e-10 ){ // fprintf(HANDLE, //" TracciaMC %d; sua FoundTrack associata (n. %d) NONsoddisfaRequisitiMinimi perche' KAPPA troppo grande; KAPPA = %g\n", // i,ii,KAPPA[ii]); // continue; } Double_t dista=sqrt( Ox[ii]*Ox[ii]+Oy[ii]*Oy[ii] ); if(fabs(dista)<1.e-20 ){ fprintf(HANDLE, " TracciaMC %d; sua FoundTrack associata (n. %d) NONsoddisfaRequisitiMinimi perche' centro Helix Cilinder trovato dista solo %g da (0,0)\n", i,ii,dista); continue; } pMCtr = (PndMCTrack*) fMCTrackArray->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 Px = pMCtr->GetMomentum().X(); Py = pMCtr->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); Fifi = atan2(Cy, Cx); // MC truth Fifi angle of circle of Helix trajectory if(Fifi<0.) Fifi += 2.*PI; 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, HoughFi, Kakka, KAPPA[ ii ], fmod(Fifi+ PI, 2.*PI), // FI0 da MC truth FI0[ ii ] ); } // end of for (ii=0; ii=2 ;ii++) //--------------ghosts // fa il conto delle ghost solo sugli eventi che hanno almeno 1 traccia MC accettabile. int NParghost=0, NParhitsghost=0,icc; if( istampa>=1 && nMCTracksaccettabili>0 ){ for(icc=0; icc=2) } // end of if( nMCTracks >0 && nTracksFoundSoFar > 0 && doMcComparison ) //-------------------- fine della sezione sul confronto tra MC truth e tracce trovate //---------------------------------------------------------------------------------- // loading the hits found and associated to a track in a PndTrackCand class; // a class per each track int ipinco=0, ipanco=0; Int_t flag; Double_t Ptras,Pxini,Pyini,Pzini,dista, qop; PndTrackCand *pTrckCand; PndTrack* pTrck; TVector3 Momentum,ErrMomentum,Position,ErrPosition; for(i=0; iAt(ipinco); if(fabs(KAPPA[i])>1.e-20 ){ Pzini = -Charge[i]*0.003*BFIELD/KAPPA[i]; // PMAX is the maximum Pz tolerable; it is set at 100 for now. if(fabs(Pzini) < PMAX) flag =0 ; else flag=-1; TVector3 dirSeed(Pxini, Pyini, Pzini); // momentum direction in starting point qop = Charge[i]/dirSeed.Mag(); dirSeed.SetMag(1.); pTrckCand->setTrackSeed(posSeed, dirSeed, qop); if(doMcComparison){ pTrckCand->setMcTrackId( daTrackFoundaTrackMC[i] ); } else { pTrckCand->setMcTrackId(-1); } for(j=0; j< nTotalHits[i]; j++){ pTrckCand->AddHit( FairRootManager::Instance()->GetBranchId(fSttBranch), (Int_t) BigList[i][j] , j); } // add the SciTil hit(s). if(nSciTilHitsinTrack[i]>0){ for(j=0; j< nSciTilHitsinTrack[i]; j++){ pTrckCand->AddHit( 1001,// my personal hit type for SciTil. (Int_t) ListSciTilHitsinTrack[i][j] , nTotalHits[i]+j); } } //-------- do relevant calculation for this track and load the PndTrack class //---- first hit if(fabs(info[ BigList[i][0] ][5] -1.)< 1.e-10 ) {// it is a parallel straw PndSttInfoXYZParal( info, BigList[i][0], Ox[i], Oy[i], R[i], KAPPA[i], FI0[i], Charge[i], Posiz1 ); } else { // it is a skew straw PndSttInfoXYZSkew( Zfinal[i][ BigList[i][0] ], // Z coordinate of selected Skew hit ZDriftfinal[i][ BigList[i][0] ], // drift distance IN Z DIRECTION only, of Skew hit Sfinal[i][ BigList[i][0] ], Ox[i], Oy[i], R[i], KAPPA[i], FI0[i], Charge[i], Posiz1 ); } // if all X, Y, Z positions were found for the first hit, go on if( Posiz1[0] > -777777776. ) { //---- last hit if(fabs(info[ BigList[i][nTotalHits[i]-1] ][5] -1.)< 1.e-10 ) {// it is a parallel straw PndSttInfoXYZParal ( info, BigList[i][nTotalHits[i]-1], Ox[i], Oy[i], R[i], KAPPA[i], FI0[i], Charge[i], Posiz2 ); } else { // it is a skew straw PndSttInfoXYZSkew ( Zfinal[i][ BigList[i][nTotalHits[i]-1] ], // Z coordinate of selected Skew hit ZDriftfinal[i][ BigList[i][nTotalHits[i]-1] ], // drift distance IN Z DIRECTION only, of Skew hit Sfinal[i][ BigList[i][nTotalHits[i]-1] ], Ox[i], Oy[i], R[i], KAPPA[i], FI0[i], Charge[i], Posiz2 ); } if( Posiz2[0] > -777777776. ) { // load in FairTrackParP first the relevant quantities // of the first hit Position.SetX( Posiz1[0] ); Position.SetY( Posiz1[1] ); Position.SetZ( Posiz1[2] ); ErrPosition.SetX(0.02); // 200 microns ErrPosition.SetY(0.02); // 200 microns ErrPosition.SetZ(1.); // 1 cm // calculate Px and Py versor[0] = Ox[i]-Posiz1[0]; versor[1] = Oy[i]-Posiz1[1]; Distance = sqrt(versor[0]*versor[0]+versor[1]*versor[1]); // I already know from PndSttInfoXYZ... that versor is not zero. versor[0] /= Distance; versor[1] /= Distance; Px = -Charge[i]*Ptras*versor[1]; Py = Charge[i]*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[i], Position, TVector3(Py/Ptras, -Px/Ptras, 0.), TVector3(Pzini*Px/ddd,Pzini*Py/ddd, -Ptras*Ptras/ddd) ); // load in FairTrackParP first the relevant quantities // of the last hit Position.SetX( Posiz2[0] ); Position.SetY( Posiz2[1] ); Position.SetZ( Posiz2[2] ); ErrPosition.SetX(0.02); // 200 microns ErrPosition.SetY(0.02); // 200 microns ErrPosition.SetZ(1.); // 1 cm // calculate Px and Py versor[0] = Ox[i]-Posiz2[0]; versor[1] = Oy[i]-Posiz2[1]; Distance = sqrt(versor[0]*versor[0]+versor[1]*versor[1]); // I already know from PndSttInfoXYZ... that versor is not zero. versor[0] /= Distance; versor[1] /= Distance; Px = -Charge[i]*Ptras*versor[1]; Py = Charge[i]*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 FairTrackParP last( Position, Momentum, ErrPosition, ErrMomentum, Charge[i], Position, TVector3(Py/Ptras, -Px/Ptras, 0.), TVector3(Pzini*Px/ddd,Pzini*Py/ddd, -Ptras*Ptras/ddd) ); //-------- load PndTrack object pTrck = new((*trackArray)[ipanco]) PndTrack(first, last, *pTrckCand); pTrck->SetRefIndex(ipanco); pTrck->SetFlag(flag); // at this point flag = 0 for Pz-777777776.) //case of something wrong with last hit FairTrackParP first( TVector3(-99999., -99999., -99999.), // dummy Position TVector3(-99999., -99999., -99999.), // dummy Momentum TVector3(-99999., -99999., -99999.), // dummy ErrPosition TVector3(-99999., -99999., -99999.), // dummy ErrMomentum 0, // dummy Charge TVector3(-99999., -99999., -99999.), // dummy Position again TVector3(1., 0., 0.), // dummy direction versor TVector3(0., 1., 0.) // dummy direction versor ); FairTrackParP last( TVector3(-99999., -99999., -99999.), // dummy Position TVector3(-99999., -99999., -99999.), // dummy Momentum TVector3(-99999., -99999., -99999.), // dummy ErrPosition TVector3(-99999., -99999., -99999.), // dummy ErrMomentum 0, // dummy Charge TVector3(-99999., -99999., -99999.), // dummy Position again TVector3(1., 0., 0.), // dummy direction versor TVector3(0., 1., 0.) // dummy direction versor ); //-------- load PndTrack object pTrck = new((*trackArray)[ipanco]) PndTrack(first, last, *pTrckCand); pTrck->SetRefIndex(ipanco); pTrck->SetFlag(-1); ipanco++; } // end of if( Posiz2[0]>-777777776. ) } else { // continuation of if( Posiz1[0] > -777777776.) //case of something wrong with first hit FairTrackParP first( TVector3(-99999., -99999., -99999.), // dummy Position TVector3(-99999., -99999., -99999.), // dummy Momentum TVector3(-99999., -99999., -99999.), // dummy ErrPosition TVector3(-99999., -99999., -99999.), // dummy ErrMomentum 0, // dummy Charge TVector3(-99999., -99999., -99999.), // dummy Position again TVector3(1., 0., 0.), // dummy direction versor TVector3(0., 1., 0.) // dummy direction versor ); FairTrackParP last( TVector3(-99999., -99999., -99999.), // dummy Position TVector3(-99999., -99999., -99999.), // dummy Momentum TVector3(-99999., -99999., -99999.), // dummy ErrPosition TVector3(-99999., -99999., -99999.), // dummy ErrMomentum 0, // dummy Charge TVector3(-99999., -99999., -99999.), // dummy Position again TVector3(1., 0., 0.), // dummy direction versor TVector3(0., 1., 0.) // dummy direction versor ); //-------- load PndTrack object pTrck = new((*trackArray)[ipanco]) PndTrack(first, last, *pTrckCand); pTrck->SetRefIndex(ipanco); pTrck->SetFlag(-1); ipanco++; } // end of if( Posiz1[0] > -777777776.) } else { // continuation of if(fabs(KAPPA[i])>1.e-20 ) Pzini = 999999.; TVector3 dirSeed(Pxini, Pyini, Pzini); // momentum direction in starting point qop = Charge[i]/dirSeed.Mag(); dirSeed.SetMag(1.); pTrckCand->setTrackSeed(posSeed, dirSeed, qop); pTrckCand->setMcTrackId( daTrackFoundaTrackMC[i] ); for(j=0; j< nTotalHits[i]; j++){ pTrckCand->AddHit( // FairRootManager::Instance()->GetBranchId("STTHit"), FairRootManager::Instance()->GetBranchId(fSttBranch), (Int_t) BigList[i][j] , j); } // load dummy quantities in PndTrack FairTrackParP first( TVector3(-99999., -99999., -99999.), // dummy Position TVector3(-99999., -99999., -99999.), // dummy Momentum TVector3(-99999., -99999., -99999.), // dummy ErrPosition TVector3(-99999., -99999., -99999.), // dummy ErrMomentum 0, // dummy Charge TVector3(-99999., -99999., -99999.), // dummy Position again TVector3(1., 0., 0.), // dummy direction versor TVector3(0., 1., 0.) // dummy direction versor ); FairTrackParP last( TVector3(-99999., -99999., -99999.), // dummy Position TVector3(-99999., -99999., -99999.), // dummy Momentum TVector3(-99999., -99999., -99999.), // dummy ErrPosition TVector3(-99999., -99999., -99999.), // dummy ErrMomentum 0, // dummy Charge TVector3(-99999., -99999., -99999.), // dummy Position again TVector3(1., 0., 0.), // dummy direction versor TVector3(0., 1., 0.) // dummy direction versor ); //-------- load PndTrack object pTrck = new((*trackArray)[ipanco]) PndTrack(first, last, *pTrckCand); pTrck->SetRefIndex(ipanco); pTrck->SetFlag(-1); ipanco++; } // end of if(fabs(KAPPA[i])>1.e-20 ) } else { // continuation of if( GoodSkewFit[i] ) // case in which there is no // skew hits in this track. // This fact is signalled by Pzini = 9999.; new((*trackCandArray)[ipinco]) PndTrackCand; pTrckCand = (PndTrackCand*) trackCandArray->At(ipinco); TVector3 dirSeed(Pxini, Pyini, Pzini); // momentum direction in starting point qop = Charge[i]/Ptras; // as if Pz=0 // dirSeed.SetMag(1.); pTrckCand->setTrackSeed(posSeed, dirSeed, qop); pTrckCand->setMcTrackId( daTrackFoundaTrackMC[i] ); for(j=0; j< nTotalHits[i]; j++){ pTrckCand->AddHit( // FairRootManager::Instance()->GetBranchId("STTHit"), FairRootManager::Instance()->GetBranchId(fSttBranch), (Int_t) BigList[i][j] , j); } // load dummy quantities in PndTrack FairTrackParP first( TVector3(-99999., -99999., -99999.), // dummy Position TVector3(-99999., -99999., -99999.), // dummy Momentum TVector3(-99999., -99999., -99999.), // dummy ErrPosition TVector3(-99999., -99999., -99999.), // dummy ErrMomentum 0, // dummy Charge TVector3(-99999., -99999., -99999.), // dummy Position again TVector3(1., 0., 0.), // dummy direction versor TVector3(0., 1., 0.) // dummy direction versor ); FairTrackParP last( TVector3(-99999., -99999., -99999.), // dummy Position TVector3(-99999., -99999., -99999.), // dummy Momentum TVector3(-99999., -99999., -99999.), // dummy ErrPosition TVector3(-99999., -99999., -99999.), // dummy ErrMomentum 0, // dummy Charge TVector3(-99999., -99999., -99999.), // dummy Position again TVector3(1., 0., 0.), // dummy direction versor TVector3(0., 1., 0.) // dummy direction versor ); //-------- load PndTrack object pTrck = new((*trackArray)[ipanco]) PndTrack(first, last, *pTrckCand); pTrck->SetRefIndex(ipanco); pTrck->SetFlag(-1); ipanco++; } // end of if( GoodSkewFit[i] ) //--- now increment the : number of PndTrackCand counter ipinco++; } // end of for(i=0; i0 && nTracksFoundSoFar > 0 && doMcComparison ) { for(i=0; i=2) cout<<"PndTrackFinderReal::DoFind, paralleli, n. hits (original notation) = "<< ParalCommonList[i][j] << ", X = "<< Posiz[0]<< ", Y = "<< Posiz[1]<< ", Z = "<< Posiz[2]<=2){ if( info[ParalCommonList[i][j]][6]<0. ){ fprintf(PHANDLEX,"Noise hit\n"); fprintf(PHANDLEY,"Noise hit\n"); fprintf(PHANDLEZ,"Noise hit\n"); } else { fprintf(PHANDLEX,"%g\n", veritaMC[ ParalCommonList[i][j] ] [0] - Posiz[0]); fprintf(PHANDLEY,"%g\n", veritaMC[ ParalCommonList[i][j] ] [1] - Posiz[1]); fprintf(PHANDLEZ,"%g\n", veritaMC[ ParalCommonList[i][j] ] [2] - Posiz[2]); } } } // end of for( j=0; j=2) cout<<"DoFind, skew, n. hits (original notation) = "<< SkewCommonList[i][j] << ", X = "<< Posiz[0]<< ", Y = "<< Posiz[1]<< ", Z = "<< Posiz[2]<=2){ if( info[SkewCommonList[i][j]][6]<0. ){ fprintf(SHANDLEX,"Noise hit\n"); fprintf(SHANDLEY,"Noise hit\n"); fprintf(SHANDLEZ,"Noise hit\n"); } else { fprintf(SHANDLEX,"%g\n", veritaMC[ SkewCommonList[i][j] ] [0] - Posiz[0]); fprintf(SHANDLEY,"%g\n", veritaMC[ SkewCommonList[i][j] ] [1] - Posiz[1]); fprintf(SHANDLEZ,"%g\n", veritaMC[ SkewCommonList[i][j] ] [2] - Posiz[2]); } } } // end of for( j=0; j0 && nTracksFoundSoFar > 0 && doMcComparison ) //--------------- end of printouts of comparison with MC for judging algorithm performance //--------------------------------------------------------------------------------------------------------- //------------------------------------- macro di display. if(iplotta && IVOLTE <= nmassimo){ //-------------------------------------------------- macro for display //------ //----------------fine stampa WriteMacroParallelHitsGeneral( keepit, Nhits, info, Nincl, Minclinations, inclination, nTracksFoundSoFar ); //----------------------------------- end macro for display if(doMcComparison) { WriteMacroParallelHitsGeneralConformalwithMC( keepit, Nhits, info, Nincl, Minclinations, inclination, nTracksFoundSoFar ); } for(i=0, ii=-1; i0){ for(j=0;j= 0){ WriteMacroSkewAssociatedHitswithMC( GoodSkewFit[i], KAPPA[i],FI0[i], HoughD, HoughFi, HoughR, info, Nincl,Minclinations,inclination, i, ii, nSttSkewhitinTrack[i], ListSkewHitsinTrack, nSkewCommon[i], SkewCommonList, daTrackFoundaTrackMC[i], nMCSkewAlone, MCSkewAloneList, nSciTilHitsinTrack[i], ESSE, ZETA ); } } } // end of if (iplotta) } // end of for(i=0; i0.){ S = d/a; T = -c/a; BBB = y1; CCC = (S-T*Radius)*(S-T*Radius)-Radius*Radius-2.*x1*(S-T*Radius)-2.*i1*Radius*r1+x1*x1+y1*y1-r1*r1; DELTA = BBB*BBB-CCC; if(DELTA<0.) { continue; } else if (DELTA ==0.){ Nsol++; R[Nsol-1] = Radius ; CX[Nsol-1] = S-T*Radius; CY[Nsol-1] = -BBB; } else{ Nsol+=2; R[Nsol-2] = Radius ; CX[Nsol-2] = S-T*Radius; CY[Nsol-2] = -BBB+sqrt(DELTA); R[Nsol-1] = Radius ; CX[Nsol-1] = S-T*Radius; CY[Nsol-1] = -BBB-sqrt(DELTA); } } // end if ( Radius >0.) } else if(bp*c-b*cp != 0.) { Radius = (-bp*d+b*dp)/(bp*c-b*cp) ; if( Radius > 0. && b != 0. ){ S = (d+c*Radius)/b; T = a/b; AAA = 1+T*T; BBB = -(T*S+x1-T*y1); CCC = S*S-Radius*Radius-2.*y1*S-2.*i1*r1*Radius+x1*x1+y1*y1-r1*r1; DELTA= BBB*BBB-AAA*CCC; //cout<<"AAA ="< 0.){ Nsol++; R[Nsol-1] = solution; CX[Nsol-1] = A+B*R[Nsol-1]; CY[Nsol-1] = C+D*R[Nsol-1]; } solution=(-BBB-radix)/AAA; //cout<<"-BBB[j] radix AAA[j] : "<<-BBB <<"; "<< radix <<"; "<< AAA<<"; solution=- :"< 0.){ Nsol++; R[Nsol-1] = solution; CX[Nsol-1] = A+B*R[Nsol-1]; CY[Nsol-1] = C+D*R[Nsol-1]; } } // end of if(AAA == 0.) } // end for(i3 } // end for(i2 } // end for(i1 } // end if ( aaa == 0.) //cout<<"----------------------------------------\n"; // cout<<"n. soluzioni : "< 2.*STRAWRADIUS ) continue; auxListHitsinTrack[Nassociatedhits]=i; Nassociatedhits++; // } // end of if( info[i][5] == 1 ) else } // end of for( i=1; i< Nhits; i++) return Nassociatedhits; } //----------end of function PndSttTrackFinderReal::PndSttTrkAssociatedParallelHitsToHelix //----------start of function PndSttTrackFinderReal::plottamentiParalleleGenerali void PndSttTrackFinderReal::plottamentiParalleleGenerali( Int_t Nremaining, Float_t * RemainingR, Float_t * RemainingD, Float_t * RemainingFi, Float_t * RemainingCX, Float_t * RemainingCY, bool *Goodflag ) { int itemp, i, j, ii, jj; Double_t D, Fi; char nome[300],titolo[300]; sprintf(nome,"HoughGeneralPlotsEvent%d.root",IVOLTE); TFile hfile(nome,"RECREATE", "STT pattern recognition"); //---------- sprintf(titolo,"Cxofcircle"); TH1F hCX(titolo,titolo,nbinCX,CXmin,CXmax); //---------- sprintf(titolo,"Cyofcircle"); TH1F hCY(titolo,titolo,nbinCY,CYmin,CYmax); //----- sprintf(titolo,"Radiusofcircle"); TH1F hR(titolo,titolo,nbinR,Rmin,Rmax); //----- sprintf(titolo,"CxabscissavsCy"); TH2F hCX_CY(titolo,titolo,nbinCX,CXmin,CXmax,nbinCY,CYmin,CYmax); //----- sprintf(titolo,"CxabscissavsRadius"); TH2F hCX_R(titolo,titolo,nbinCX,CXmin,CXmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"CyabscissavsRadius"); TH2F hCY_R(titolo,titolo,nbinCY,CYmin,CYmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"CxabscissavsCyordinatevsRadius"); TH3F hCX_CY_R(titolo,titolo,nbinCX,CXmin,CXmax,nbinCY,CYmin,CYmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"Distanceofclosestapproachofcircle"); TH1F hD(titolo,titolo,nbinD,Dmin,Dmax); //---------- sprintf(titolo,"Firadofcircle"); TH1F hFi(titolo,titolo,nbinFi,Fimin,Fimax); //----- sprintf(titolo,"DabscissavsFi"); TH2F hD_Fi(titolo,titolo,nbinD,Dmin,Dmax,nbinFi,Fimin,Fimax); //----- sprintf(titolo,"DabscissavsRadius"); TH2F hD_R(titolo,titolo,nbinD,Dmin,Dmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"FiabscissavsRadius"); TH2F hFi_R(titolo,titolo,nbinFi,Fimin,Fimax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"DabscissavsFiordinatevsRadius"); TH3F hD_Fi_R(titolo,titolo,nbinD,Dmin,Dmax,nbinFi,Fimin,Fimax,nbinR,Rmin,Rmax); //----- //--------------------------------------------------- // gli stessi plots ma per le soluzioni che sono le piu' vicine alla verita' MC //---------- sprintf(titolo,"Cxofcircle_truth"); TH1F hCXverita(titolo,titolo,nbinCX,CXmin,CXmax); //---------- sprintf(titolo,"Cyofcircle_truth"); TH1F hCYverita(titolo,titolo,nbinCY,CYmin,CYmax); //----- sprintf(titolo,"Radiusofcircle_truth"); TH1F hRverita(titolo,titolo,nbinR,Rmin,Rmax); //----- sprintf(titolo,"CxabscissavsCy_truth"); TH2F hCX_CYverita(titolo,titolo,nbinCX,CXmin,CXmax,nbinCY,CYmin,CYmax); //----- sprintf(titolo,"CxabscissavsRadius_truth"); TH2F hCX_Rverita(titolo,titolo,nbinCX,CXmin,CXmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"CyabscissavsRadius_truth"); TH2F hCY_Rverita(titolo,titolo,nbinCY,CYmin,CYmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"CxabscissavsCyordinatevsRadius_truth"); TH3F hCX_CY_Rverita(titolo,titolo,nbinCX,CXmin,CXmax,nbinCY,CYmin,CYmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"Distanceofclosestapproachofcircle_truth"); TH1F hDverita(titolo,titolo,nbinD,Dmin,Dmax); //---------- sprintf(titolo,"Firadofcircle_truth"); TH1F hFiverita(titolo,titolo,nbinFi,Fimin,Fimax); //----- sprintf(titolo,"DabscissavsFi_truth"); TH2F hD_Fiverita(titolo,titolo,nbinD,Dmin,Dmax,nbinFi,Fimin,Fimax); //----- sprintf(titolo,"DabscissavsRadius_truth"); TH2F hD_Rverita(titolo,titolo,nbinD,Dmin,Dmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"FiabscissavsRadius_truth"); TH2F hFi_Rverita(titolo,titolo,nbinFi,Fimin,Fimax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"DabscissavsFiordinatevsRadius_truth"); TH3F hD_Fi_Rverita(titolo,titolo,nbinD,Dmin,Dmax,nbinFi,Fimin,Fimax,nbinR,Rmin,Rmax); //----- //--------------------------------------------------- // gli stessi plots ma per le soluzioni che NON sono le piu' vicine alla verita' MC //---------- sprintf(titolo,"Cxofcircle_nontruth"); TH1F hCXnonverita(titolo,titolo,nbinCX,CXmin,CXmax); //---------- sprintf(titolo,"Cyofcircle_nontruth"); TH1F hCYnonverita(titolo,titolo,nbinCY,CYmin,CYmax); //----- sprintf(titolo,"Radiusofcircle_nontruth"); TH1F hRnonverita(titolo,titolo,nbinR,Rmin,Rmax); //----- sprintf(titolo,"CxabscissavsCy_nontruth"); TH2F hCX_CYnonverita(titolo,titolo,nbinCX,CXmin,CXmax,nbinCY,CYmin,CYmax); //----- sprintf(titolo,"CxabscissavsRadius_nontruth"); TH2F hCX_Rnonverita(titolo,titolo,nbinCX,CXmin,CXmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"CyabscissavsRadius_nontruth"); TH2F hCY_Rnonverita(titolo,titolo,nbinCY,CYmin,CYmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"CxabscissavsCyordinatevsRadius_nontruth"); TH3F hCX_CY_Rnonverita(titolo,titolo,nbinCX,CXmin,CXmax,nbinCY,CYmin,CYmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"Distanceofclosestapproachofcircle_nontruth"); TH1F hDnonverita(titolo,titolo,nbinD,Dmin,Dmax); //---------- sprintf(titolo,"Firadofcircle_nontruth"); TH1F hFinonverita(titolo,titolo,nbinFi,Fimin,Fimax); //----- sprintf(titolo,"DabscissavsFi_nontruth"); TH2F hD_Finonverita(titolo,titolo,nbinD,Dmin,Dmax,nbinFi,Fimin,Fimax); //----- sprintf(titolo,"DabscissavsRadius_nontruth"); TH2F hD_Rnonverita(titolo,titolo,nbinD,Dmin,Dmax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"FiabscissavsRadius_nontruth"); TH2F hFi_Rnonverita(titolo,titolo,nbinFi,Fimin,Fimax,nbinR,Rmin,Rmax); //----- sprintf(titolo,"DabscissavsFiordinatevsRadius_nontruth"); TH3F hD_Fi_Rnonverita(titolo,titolo,nbinD,Dmin,Dmax,nbinFi,Fimin,Fimax,nbinR,Rmin,Rmax); //----- //--------------------------------------------------- for(itemp = 0; itempRlow && RemainingR[itemp]Dlow && RemainingD[itemp]Filow && RemainingFi[itemp]3 || vec1[i]-vec2[i]<-3) return false; } return true; } //----------end of function PndSttTrackFinderReal::iscontiguous void PndSttTrackFinderReal::clustering2( Short_t vec1[2], // input int nListElements, Short_t List[][2], // input int & nClusterElementsFound, Short_t ClusterElementsFound[][2], // output int & nRemainingElements, Short_t RemainingElements[][2] // output ) { int i; Short_t vec2[2]; nClusterElementsFound=0; nRemainingElements=0; for(i=0; i xmax) xmax = posizSciTil[i][0]; if (posizSciTil[i][1] < ymin) ymin = posizSciTil[i][1]; if (posizSciTil[i][1] > ymax) ymax = posizSciTil[i][1]; } //------ for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws 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]; } } 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); disegnaAssiXY(MACRO,xmin,xmax,ymin,ymax); //---- disegna gli Scitil. for( i=0; i< nSciTilHits; i++) { disegnaSciTilHit( MACRO, i, posizSciTil[i][0], posizSciTil[i][1], 0 ); } //------------------------ for( i=0; i< Nhits; i++) { 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); } } //------------------------------- plotting all the tracks found for(i=0; iSetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", i,aaa,bbb,rrr,rrr,i,i,i); } else { if(fabs(BETA[i]) < 1.e-10){ // fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,- GAMMA[i]/ALFA[i],ymin,- GAMMA[i]/ALFA[i],ymax); fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,0.,ymin,0.,ymax); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } else { // yl = -xmin*ALFA[i]/BETA[i] - GAMMA[i]/BETA[i]; // yu = -xmax*ALFA[i]/BETA[i] - GAMMA[i]/BETA[i]; yl = -xmin*ALFA[i]/BETA[i]; yu = -xmax*ALFA[i]/BETA[i]; fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } } if(istampa>= 3 && IVOLTE <= nmassimo) { cout<<"stampa da WriteMacroParallelHitsGeneral , for trackfoundsofar n. "< xmax) xmax = posizSciTil[i][0]; if (posizSciTil[i][1] < ymin) ymin = posizSciTil[i][1]; if (posizSciTil[i][1] > ymax) ymax = posizSciTil[i][1]; } //------ for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws 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]; } } 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); disegnaAssiXY(MACRO,xmin,xmax,ymin,ymax); //---- disegna gli Scitil. for( i=0; i< nSciTilHits; i++) { /* fprintf(MACRO,"TMarker* SciT%d = new TMarker(%f,%f,%d);\n", i,posizSciTil[i][0],posizSciTil[i][1],30); fprintf(MACRO,"SciT%d->SetMarkerSize(1.5);\n",i); fprintf(MACRO,"SciT%d->SetMarkerColor(1);\nSciT%d->Draw();\n" ,i,i); */ disegnaSciTilHit( MACRO, i, posizSciTil[i][0], posizSciTil[i][1], 0 ); } //------------------------ for( i=0; i< Nhits; i++) { 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); } } //--------------------------- ora le tracce MC Int_t icode; Double_t Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Px, Py, carica ; PndMCTrack* pMC; for(int im=0;imAt(im); if ( ! pMC ) continue; icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if( fabs(carica)<0.1) continue; Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); fprintf(MACRO,"TEllipse* MC%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nMC%d->SetFillStyle(0);\nMC%d->SetLineColor(3);\nMC%d->Draw();\n", im,Cx,Cy,Rr,Rr,im,im,im); } // end of for(int im=0;imSetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", i,aaa,bbb,rrr,rrr,i,i,i); } else { if(fabs(BETA[i]) < 1.e-10){ fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,- GAMMA[i]/ALFA[i],ymin,- GAMMA[i]/ALFA[i],ymax); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } else { yl = -xmin*ALFA[i]/BETA[i] - GAMMA[i]/BETA[i]; yu = -xmax*ALFA[i]/BETA[i] - GAMMA[i]/BETA[i]; fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } } } // ----------- fprintf(MACRO,"}\n"); fclose(MACRO); } // end of if(doMcComparison) //------------------------------------------------------------------------------------------------------------ // ora riplotto tutto nello spazio conforme usando la trasformazione u= x/(x**2+y**2) e v = y/(x**2+y**2) // // aggiungo anche la grigliatura dello spazio conforme usata per la box del pattern recognition sprintf(nome,"MacroSttParConformeEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); xmin=1.e20; xmax=-1.e20; ymin=1.e20; ymax=-1.e20; //--- SciTil info Double_t USciTil[nmaxSciTilHits], VSciTil[nmaxSciTilHits]; for( i=0; i< nSciTilHits; i++) { Double_t erre = posizSciTil[i][0]*posizSciTil[i][0]+ posizSciTil[i][1]*posizSciTil[i][1]; USciTil[i] = posizSciTil[i][0]/erre; VSciTil[i] = posizSciTil[i][1]/erre; if (USciTil[i] < xmin) xmin = USciTil[i]; if (USciTil[i] > xmax) xmax = USciTil[i]; if (VSciTil[i] < ymin) ymin = VSciTil[i]; if (VSciTil[i] > ymax) ymax = VSciTil[i]; } //------ for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws // centro sfera in sistema conforme gamma = info[i][0]*info[i][0] + info[i][1]*info[i][1] - info[i][3]*info[i][3]; Ox[i] = info[i][0] / gamma; Oy[i] = info[i][1] / gamma; Radius[i] = info[i][3]/gamma; if (Ox[i]-Radius[i] < xmin) xmin = Ox[i]-Radius[i]; if (Ox[i]+Radius[i] > xmax) xmax = Ox[i]+Radius[i]; if (Oy[i]-Radius[i] < ymin) ymin = Oy[i]-Radius[i]; if (Oy[i]+Radius[i] > ymax) ymax = Oy[i]+Radius[i]; } } 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); // ora la grigliatura con i cerchi fprintf(MACRO,"TEllipse* Griglia%d = new TEllipse(0.,0.,%f,%f,0.,360.);\nGriglia%d->SetLineColor(4);\nGriglia%d->Draw();\n", nRdivConformalEffective,1./RStrawDetectorMin,1./RStrawDetectorMin,nRdivConformalEffective,nRdivConformalEffective); for( i=nRdivConformalEffective-1; i>=0 ; i--) { fprintf(MACRO,"TEllipse* Griglia%d = new TEllipse(0.,0.,%f,%f,0.,360.);\nGriglia%d->SetLineColor(4);\nGriglia%d->Draw();\n", i,radiaConf[i],radiaConf[i],i,i); } //--------------------- // ora la grigliatura con i segmenti blu per delimitare la zona degli Stt. for(i=0; iSetLineColor(4);\nSeg%d->Draw();\n", i,x1,y1,x2,y2,i,i); } //--------------------- // ora la grigliatura con i segmenti magenta per comprendere la zona degli SciTil. for(i=0; iSetLineColor(6);\nSeg%d->Draw();\n", i,x1,y1,x2,y2,i,i); } //--------------------- fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n",xmin,0.,xmax,0.,xmin,xmax); fprintf(MACRO,"Assex->SetTitle(\"U \");\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(\"V \");\n"); fprintf(MACRO,"Assey->SetTitleOffset(1.5);\n"); fprintf(MACRO,"Assey->Draw();\n"); // plot degli Hits Stt. for( i=0; i< Nhits; i++) { 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,Ox[i],Oy[i],Radius[i],Radius[i],i,i); } } // plot degli Hit SciTil for( i=0; i< nSciTilHits; i++) { disegnaSciTilHit( MACRO, i, posizSciTil[i][0], posizSciTil[i][1], 2 // 2--> disegno SciTil in conforme. ); /* fprintf(MACRO,"TMarker* SciT%d = new TMarker(%f,%f,%d);\n", i,USciTil[i],VSciTil[i],30); fprintf(MACRO,"SciT%d->SetMarkerSize(1.1);\n",i); fprintf(MACRO,"SciT%d->SetMarkerColor(1);\nSciT%d->Draw();\n" ,i,i); */ } //------------------ plotting all the tracks found for(i=0; iSetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } } else { yl = -xmin*ALFA[i]/BETA[i] - 1./BETA[i]; yu = -xmax*ALFA[i]/BETA[i] - 1./BETA[i]; fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } } else { // in this case TypConf = 0 --> straight line in XY if( fabs(GAMMA[i]) > 1.e-10) { aaa = -0.5*ALFA[i]/GAMMA[i]; bbb = -0.5*BETA[i]/GAMMA[i]; rrr = sqrt( aaa*aaa+bbb*bbb); fprintf(MACRO,"TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,0.,360.);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", i,aaa,bbb,rrr,rrr,i,i,i); } else { if(fabs(BETA[i])>1.e-10) { yl = -xmin*ALFA[i]/BETA[i] ; yu = -xmax*ALFA[i]/BETA[i] ; fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } else { fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,0.,ymin,0.,ymax); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } } } } // ------------------------------ fprintf(MACRO,"}\n"); fclose(MACRO); return ; } //----------end of function PndSttTrackFinderReal::WriteMacroParallelHitsGeneral //----------start of function PndSttTrackFinderReal::WriteMacroParallelHitsGeneralConformalwithMC void PndSttTrackFinderReal::WriteMacroParallelHitsGeneralConformalwithMC( bool * keepit, Int_t Nhits, Double_t info[][7], Int_t Nincl, Int_t Minclinations[], Double_t inclination[][3], Short_t nTracksFoundSoFar ) { Int_t i, j, i1, ii, index, Kincl, nlow, nup, STATUS; Double_t xmin , xmax, ymin, ymax, xl, xu, yl, yu, rrr, Ox[nmaxHits], Oy[nmaxHits], Radius[nmaxHits], 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, angle, minor, major, distance, Rx, Ry, LL, Aellipsis1, Bellipsis1,fi1, fmin, fmax, offset, step, SkewInclWithRespectToS, zpos, zpos1, zpos2, zl[200],zu[200]; Double_t USciTil[nmaxSciTilHits], VSciTil[nmaxSciTilHits]; //---------- parallel straws Macro now char nome[300], nome2[300]; FILE * MACRO; // ora riplotto tutto nello spazio conforme usando la trasformazione u= x/(x**2+y**2) e v = y/(x**2+y**2) // // aggiungo anche la grigliatura dello spazio conforme usata per la box del pattern recognition sprintf(nome,"MacroSttParConformewithMCEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); xmin=1.e20; xmax=-1.e20; ymin=1.e20; ymax=-1.e20; //--- SciTil info for( i=0; i< nSciTilHits; i++) { Double_t erre = posizSciTil[i][0]*posizSciTil[i][0]+ posizSciTil[i][1]*posizSciTil[i][1]; USciTil[i] = posizSciTil[i][0]/erre; VSciTil[i] = posizSciTil[i][1]/erre; if (USciTil[i] < xmin) xmin = USciTil[i]; if (USciTil[i] > xmax) xmax = USciTil[i]; if (VSciTil[i] < ymin) ymin = VSciTil[i]; if (VSciTil[i] > ymax) ymax = VSciTil[i]; } //------ for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws // centro sfera in sistema conforme gamma = info[i][0]*info[i][0] + info[i][1]*info[i][1] - info[i][3]*info[i][3]; Ox[i] = info[i][0] / gamma; Oy[i] = info[i][1] / gamma; Radius[i] = info[i][3]/gamma; if (Ox[i]-Radius[i] < xmin) xmin = Ox[i]-Radius[i]; if (Ox[i]+Radius[i] > xmax) xmax = Ox[i]+Radius[i]; if (Oy[i]-Radius[i] < ymin) ymin = Oy[i]-Radius[i]; if (Oy[i]+Radius[i] > ymax) ymax = Oy[i]+Radius[i]; } } 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); // ora la grigliatura con i cerchi fprintf(MACRO,"TEllipse* Griglia%d = new TEllipse(0.,0.,%f,%f,0.,360.);\nGriglia%d->SetLineColor(4);\nGriglia%d->Draw();\n", nRdivConformalEffective,1./RStrawDetectorMin,1./RStrawDetectorMin,nRdivConformalEffective,nRdivConformalEffective); for( i=nRdivConformalEffective-1; i>=0 ; i--) { fprintf(MACRO,"TEllipse* Griglia%d = new TEllipse(0.,0.,%f,%f,0.,360.);\nGriglia%d->SetLineColor(4);\nGriglia%d->Draw();\n", i,radiaConf[i],radiaConf[i],i,i); } //--------------------- // ora la grigliatura con i segmenti blu per delimitare la zona degli Stt. for(i=0; iSetLineColor(4);\nSeg%d->Draw();\n", i,x1,y1,x2,y2,i,i); } //--------------------- // ora la grigliatura con i segmenti magenta per comprendere la zona degli SciTil. for(i=0; iSetLineColor(6);\nSeg%d->Draw();\n", i,x1,y1,x2,y2,i,i); } //--------------------- fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n",xmin,0.,xmax,0.,xmin,xmax); fprintf(MACRO,"Assex->SetTitle(\"U \");\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(\"V \");\n"); fprintf(MACRO,"Assey->SetTitleOffset(1.5);\n"); fprintf(MACRO,"Assey->Draw();\n"); // plot degli Hits Stt for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nE%d->SetLineWidth(2);\nE%d->SetFillStyle(0);\nE%d->Draw();\n", i,Ox[i],Oy[i],Radius[i],Radius[i],i,i,i); } } // plot degli Hit SciTil for( i=0; i< nSciTilHits; i++) { disegnaSciTilHit( MACRO, i, posizSciTil[i][0], posizSciTil[i][1], 2 // 2--> disegno SciTil in conforme. ); /* fprintf(MACRO,"TMarker* SciT%d = new TMarker(%f,%f,%d);\n", i,USciTil[i],VSciTil[i],30); fprintf(MACRO,"SciT%d->SetMarkerSize(1.1);\n",i); fprintf(MACRO,"SciT%d->SetMarkerColor(1);\nSciT%d->Draw();\n" ,i,i); */ } //------------------ plotting all the tracks found for(i=0; i 1.e-10) { aaa = -0.5*ALFA[i]/GAMMA[i]; bbb = -0.5*BETA[i]/GAMMA[i]; rrr = sqrt( aaa*aaa+bbb*bbb-1./GAMMA[i]); if( fabs(rrr/GAMMA[i]) < 30.) { fprintf(MACRO, "TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,0.,360.);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", i,aaa,bbb,rrr,rrr,i,i,i); } else{ yl = -xmin*ALFA[i]/BETA[i] - 1./BETA[i]; yu = -xmax*ALFA[i]/BETA[i] - 1./BETA[i]; fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } } else { if(fabs(BETA[i]) < 1.e-10){ if(fabs(ALFA[i])<1.e-10) { continue; } else { fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,-1./ALFA[i],ymin,- 1./ALFA[i],ymax); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } } else { yl = -xmin*ALFA[i]/BETA[i] - 1./BETA[i]; yu = -xmax*ALFA[i]/BETA[i] - 1./BETA[i]; fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } } } else { // in this case TypConf = 0 --> straight line in XY if( fabs(GAMMA[i]) > 1.e-10) { aaa = -0.5*ALFA[i]/GAMMA[i]; bbb = -0.5*BETA[i]/GAMMA[i]; rrr = sqrt( aaa*aaa+bbb*bbb); fprintf(MACRO,"TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,0.,360.);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", i,aaa,bbb,rrr,rrr,i,i,i); } else { if(fabs(BETA[i])>1.e-10) { yl = -xmin*ALFA[i]/BETA[i] ; yu = -xmax*ALFA[i]/BETA[i] ; fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } else { fprintf(MACRO,"TLine* ris%d = new TLine(%f,%f,%f,%f);\n",i,0.,ymin,0.,ymax); fprintf(MACRO,"ris%d->SetLineColor(2);\n",i); fprintf(MACRO,"ris%d->Draw();\n",i); } } } } // ------------------------------ // plotting all the tracks MC generated Int_t icode; Double_t Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Px, Py, carica ; PndMCTrack* pMC; for(i=0;iAt(i); if ( ! pMC ) continue; icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if (fabs(carica)<0.1 ) continue; Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); gamma = -Rr*Rr + Cx*Cx+Cy*Cy; if(fabs(gamma)< 0.001) { if(Cy != 0.) { yl = xmin*(-Cx/Cy) + 0.5/Cy; yu = xmax*(-Cx/Cy) + 0.5/Cy; xl = xmin; xu = xmax; } else { yl = ymin; yu = ymax; xu=xl = 0.5/Cx; } fprintf(MACRO,"TLine* MCris%d = new TLine(%f,%f,%f,%f);\n",i,xl,yl,xu,yu); fprintf(MACRO,"MCris%d->SetLineStyle(2);\n",i); fprintf(MACRO,"MCris%d->SetLineColor(3);\n",i); fprintf(MACRO,"MCris%d->SetLineWidth(1);\n",i); fprintf(MACRO,"MCris%d->Draw();\n",i); } else { if(fabs(Rr/gamma) > 1.) { if(fabs(Cy)>0.001 ) { yl = -xmin*Cx/Cy+0.5/Cy; yu = -xmax*Cx/Cy+0.5/Cy; fprintf(MACRO,"TLine* MCline%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"MCline%d->SetLineColor(3);\n",i); fprintf(MACRO,"MCline%d->Draw();\n",i); } else { fprintf(MACRO,"TLine* MCline%d = new TLine(%f,%f,%f,%f);\n",i,2.*CxMC[i],ymin,2.*CxMC[i],ymax); fprintf(MACRO,"MCline%d->SetLineColor(2);\n",i); fprintf(MACRO,"MCline%d->Draw();\n",i); } } else { fprintf(MACRO, "TEllipse* MCcerchio%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nMCcerchio%d->SetLineColor(3);\n", i,Cx/gamma,Cy/gamma,Rr/fabs(gamma),Rr/fabs(gamma),i); fprintf(MACRO,"MCcerchio%d->SetFillStyle(0);\nMCcerchio%d->SetLineStyle(2);\nMCcerchio%d->SetLineWidth(1);\nMCcerchio%d->Draw();\n", i,i,i,i); } } } // end of for(i=0;i xmax) xmax = Ox[i]+Radius[i]; if (Oy[i]-Radius[i] < ymin) ymin = Oy[i]-Radius[i]; if (Oy[i]+Radius[i] > ymax) ymax = Oy[i]+Radius[i]; } 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); fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n",xmin,0.,xmax,0.,xmin,xmax); 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->Draw();\n"); for( i=0; i< Nhits; i++) { // if( i== iExclude) continue; fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nE%d->SetLineWidth(2);\nE%d->SetFillStyle(0);\nE%d->Draw();\n", i,Ox[i],Oy[i],Radius[i],Radius[i],i,i,i); } //------------------ plotting the track found i = nTracksFoundSoFar; if( Status != 99){ rrr = sqrt( 0.25*ALFA[nTracksFoundSoFar]*ALFA[nTracksFoundSoFar] + 0.25*BETA[nTracksFoundSoFar]*BETA[nTracksFoundSoFar] -GAMMA[nTracksFoundSoFar]); alfa = ALFA[nTracksFoundSoFar]+2.*trajectory_vertex[0]; beta = BETA[nTracksFoundSoFar]+2.*trajectory_vertex[1]; cx = -0.5*alfa; cy = -0.5*beta; gamma = cx*cx+cy*cy-rrr*rrr; if(fabs(gamma) > 1.e-8){ fprintf(MACRO, "TEllipse* FoundTrack%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nFoundTrack%d->SetLineColor(2);\n", i,cx/gamma,cy/gamma,rrr/fabs(gamma),rrr/fabs(gamma),i); fprintf(MACRO,"FoundTrack%d->SetFillStyle(0);\nFoundTrack%d->SetLineStyle(2);\nFoundTrack%d->SetLineWidth(1);\nFoundTrack%d->Draw();\n", i,i,i,i); } else { if( fabs(beta) > 1.e-8){ yl = -xmin*alfa/beta-1./beta ; yu = -xmax*alfa/beta-1./beta ; fprintf(MACRO,"TLine* FoundTrack%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"FoundTrack%d->SetLineColor(2);\n",i); fprintf(MACRO,"FoundTrack%d->Draw();\n",i); } else { fprintf(MACRO,"TLine* FoundTrack%d = new TLine(%f,%f,%f,%f);\n",i,-1./alfa,ymin,-1./alfa,ymax); fprintf(MACRO,"FoundTrack%d->SetLineColor(2);\n",i); fprintf(MACRO,"FoundTrack%d->Draw();\n",i); } } } else { alfa = ALFA[nTracksFoundSoFar]; beta = BETA[nTracksFoundSoFar]; gamma = GAMMA[nTracksFoundSoFar]+ALFA[nTracksFoundSoFar]*trajectory_vertex[0]+BETA[nTracksFoundSoFar]*trajectory_vertex[1]; if( fabs(beta) > 1.e-8){ yl = -xmin*alfa/beta-1./beta ; yu = -xmax*alfa/beta-1./beta ; fprintf(MACRO,"TLine* FoundTrack%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"FoundTrack%d->SetLineColor(2);\n",i); fprintf(MACRO,"FoundTrack%d->Draw();\n",i); } else { fprintf(MACRO,"TLine* FoundTrack%d = new TLine(%f,%f,%f,%f);\n",i,-1./alfa,ymin,-1./alfa,ymax); fprintf(MACRO,"FoundTrack%d->SetLineColor(2);\n",i); fprintf(MACRO,"FoundTrack%d->Draw();\n",i); } } // ------------------------------ // plotting all the tracks MC generated //cout<<" TRASLAZIONE IN "<At(i); if ( ! pMC ) continue; icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if(fabs(carica)<0.1) continue; Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); // for(i=0; iSetLineStyle(2);\n",i); fprintf(MACRO,"MCris%d->SetLineColor(3);\n",i); fprintf(MACRO,"MCris%d->SetLineWidth(1);\n",i); fprintf(MACRO,"MCris%d->Draw();\n",i); } else { if(fabs(Rr/gamma) > 1.) { if(fabs(Cy)>0.001 ) { yl = -xmin*Cx/Cy+0.5/Cy; yu = -xmax*Cx/Cy+0.5/Cy; fprintf(MACRO,"TLine* MCline%d = new TLine(%f,%f,%f,%f);\n",i,xmin,yl,xmax,yu); fprintf(MACRO,"MCline%d->SetLineColor(3);\n",i); fprintf(MACRO,"MCline%d->Draw();\n",i); } else { fprintf(MACRO,"TLine* MCline%d = new TLine(%f,%f,%f,%f);\n",i,2.*CxMC[i],ymin,2.*CxMC[i],ymax); fprintf(MACRO,"MCline%d->SetLineColor(2);\n",i); fprintf(MACRO,"MCline%d->Draw();\n",i); } } else { fprintf(MACRO, "TEllipse* MCcerchio%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nMCcerchio%d->SetLineColor(3);\n", i,cx/gamma,cy/gamma,Rr/fabs(gamma),Rr/fabs(gamma),i); fprintf(MACRO, "MCcerchio%d->SetFillStyle(0);\nMCcerchio%d->SetLineStyle(2);\nMCcerchio%d->SetLineWidth(1);\nMCcerchio%d->Draw();\n", i,i,i,i); } } } fprintf(MACRO,"}\n"); fclose(MACRO); return ; } //----------end of function PndSttTrackFinderReal::WriteMacroParallelHitsGeneralConformalwithMCspecial //----------start of function PndSttTrackFinderReal::WriteMacroParallelAssociatedHits void PndSttTrackFinderReal::WriteMacroParallelAssociatedHits( Double_t Ox, Double_t Oy, Double_t R, Short_t Nhits, Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHitsInTrack], Double_t info[][7], Int_t Nincl, Int_t Minclinations[], Double_t inclination[][3], Short_t imaxima, Int_t sequencial, Short_t nscitilhitsintrack, Short_t *listscitilhitsintrack ) { 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, 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]; // Ox = (D+R)*cos(Fi); // Oy = (D+R)*sin(Fi); //cout<<"da MacroTrackparalleletc. Ox, Oy "<0) { for(j=0;j 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 = infoparal[ ListHitsinTrack[imaxima][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]; } 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); fprintf(MACRO,"TEllipse* TC = new TEllipse(%f,%f,%f,%f,0.,360.);\n",Ox,Oy,R,R); fprintf(MACRO,"TC->SetLineColor(2);\nTC->SetFillStyle(0);\nTC->Draw();\n"); disegnaAssiXY(MACRO,xmin,xmax,ymin,ymax); for( ii=0; ii< Nhits; ii++) { i = infoparal[ ListHitsinTrack[imaxima][ii] ] ; 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); } //---- disegna gli Scitil. if(nscitilhitsintrack>0) { for(j=0;jSetMarkerSize(1.5);\n",i); fprintf(MACRO,"SciT%d->SetMarkerColor(1);\nSciT%d->Draw();\n" ,i,i); */ disegnaSciTilHit( MACRO, i, posizSciTil[i][0], posizSciTil[i][1], 0 ); } } //------------------------ fprintf(MACRO,"}\n"); fclose(MACRO); return ; } //----------end of function PndSttTrackFinderReal::WriteMacroParallelAssociatedHits //----------start of function PndSttTrackFinderReal::WriteMacroParallelAssociatedHitswithMC void PndSttTrackFinderReal::WriteMacroParallelAssociatedHitswithMC( Double_t Ox, Double_t Oy, Double_t R, Short_t TrackFoundaTrackMC, Short_t Nhits, Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHitsInTrack], Double_t info[][7], Short_t ifoundtrack, Int_t sequentialNTrack, Short_t nscitilhitsintrack, Short_t *listscitilhitsintrack, Short_t nParalCommon[MAXTRACKSPEREVENT], Short_t ParalCommonList[MAXTRACKSPEREVENT][nmaxHits], Short_t nSpuriParinTrack[MAXTRACKSPEREVENT], Short_t ParSpuriList[MAXTRACKSPEREVENT][nmaxHits], Short_t nMCParalAlone[MAXTRACKSPEREVENT], Short_t MCParalAloneList[MAXTRACKSPEREVENT][nmaxHits] ) { Int_t i, j, i1, ii, index, imaxima, 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, 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]; imaxima=ifoundtrack; //---------- parallel straws Macro now char nome[300], nome2[300]; sprintf(nome,"MacroSttParwithMCEvent%dT%d", IVOLTE,sequentialNTrack); 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; //---- Scitil hits. if(nscitilhitsintrack>0) { for(j=0;j 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 = infoparal[ ListHitsinTrack[imaxima][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]; } 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); fprintf(MACRO,"TEllipse* TC = new TEllipse(%f,%f,%f,%f,0.,360.);\n",Ox,Oy,R,R); fprintf(MACRO,"TC->SetLineColor(2);\nTC->SetFillStyle(0);\nTC->Draw();\n"); disegnaAssiXY(MACRO,xmin,xmax,ymin,ymax); for( ii=0; ii< nParalCommon[ifoundtrack]; ii++) { i = ParalCommonList[ifoundtrack][ii] ; fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nE%d->SetFillStyle(0);", i,info[i][0],info[i][1],info[i][3],info[i][3],i); fprintf(MACRO,"E%d->Draw();\n",i); } for( ii=0; ii< nSpuriParinTrack[ifoundtrack]; ii++) { i = ParSpuriList[ifoundtrack][ii] ; fprintf(MACRO, "TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nE%d->SetFillStyle(0);", i,info[i][0],info[i][1],info[i][3],info[i][3],i); fprintf(MACRO,"E%d->SetLineColor(2);\n",i); fprintf(MACRO,"E%d->Draw();\n",i); } for( ii=0; ii< nMCParalAlone[ifoundtrack]; ii++) { i = MCParalAloneList[ifoundtrack][ii] ; fprintf(MACRO, "TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nE%d->SetFillStyle(0);", i,info[i][0],info[i][1],info[i][3],info[i][3],i); fprintf(MACRO,"E%d->SetLineColor(4);\n",i); fprintf(MACRO,"E%d->Draw();\n",i); } //---- disegna gli Scitil. if(nscitilhitsintrack>0) { for(j=0;jSetMarkerSize(1.5);\n",i); fprintf(MACRO,"SciT%d->SetMarkerColor(1);\nSciT%d->Draw();\n" ,i,i); */ disegnaSciTilHit( MACRO, i, posizSciTil[i][0], posizSciTil[i][1], 0 ); } } //------------------------ //--------------------------- ora le tracce MC Int_t icode, im; Double_t Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Px, Py, carica ; PndMCTrack* pMC; if( TrackFoundaTrackMC > -1){ im=TrackFoundaTrackMC; pMC = (PndMCTrack*) fMCTrackArray->At(im); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if( fabs(carica)>=0.1){ Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); fprintf(MACRO,"TEllipse* MC%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nMC%d->SetFillStyle(0);\nMC%d->SetLineColor(3);\nMC%d->Draw();\n", im,Cx,Cy,Rr,Rr,im,im,im); } // end of if( fabs(carica)>=0.1) } // end of if ( pMC ) } // end of if( TrackFoundaTrackMC > -1){ //----------- fine parte del MC fprintf(MACRO,"}\n"); fclose(MACRO); return ; } //----------end of function PndSttTrackFinderReal::WriteMacroParallelAssociatedHitswithMC //----------start of function PndSttTrackFinderReal::WriteMacroSkewAssociatedHits void PndSttTrackFinderReal::WriteMacroSkewAssociatedHits( bool goodskewfit, Double_t KAPPA, Double_t FI0, Double_t D, Double_t Fi, Double_t R, Double_t info[][7], Int_t Nincl, Int_t Minclinations[], Double_t inclination[][3], Int_t imaxima, Int_t sequentialNTrack, Short_t nSttSkewhitinTrack, Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxHitsInTrack], Short_t nscitilhits, Double_t *ESSE, Double_t *ZETA ) { Int_t i, j, i1, ii, iii, index, Kincl, nlow, nup, STATUS, imc, Nmin, Nmax; Double_t xmin , xmax, ymin, ymax, Ox, Oy, dx, dy, diff, d1, d2, delta, deltax, deltay, deltaz, deltaS, factor, zmin, zmax, Smin, Smax, S1, S2, z1, z2, y1, y2, vx1, vy1, vz1, C0x1, C0y1, C0z1, aaa, bbb, ccc, angle, minor, major, distance, Rx, Ry, LL, Aellipsis1, Bellipsis1,fi1, fmin, fmax, offset, step, SkewInclWithRespectToS, zpos, zpos1, zpos2, Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; //------------------- skew straws hits Macro now char nome2[300],nome[300]; FILE *MACRO; sprintf(nome,"MacroSttSkewEvent%dT%d",IVOLTE,sequentialNTrack); sprintf(nome2, "%s.C",nome); MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); //KAPPA = 1./166.67 ; FI0 = 1.5*PI; Smin=zmin = 1.e10; Smax=zmax = -zmin; index=0; // prima lo hits SciTil for(i=0;iSmax ) Smax=ESSE[i]; if( ESSE[i]zmax ) zmax=ZETA[i]; if( ZETA[i] 1.e-10){ Tiltdirection1[0] = vz1/bbb; Tiltdirection1[1] = SkewInclWithRespectToS/bbb; } else { Tiltdirection1[0] = 1.; Tiltdirection1[1] = 0.; } LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; if( distance >= info[i][4] + Aellipsis1) continue; // checks that the projected ellipsis doesn't go out the boundaries of both the skew straw and the trajectory cylinder //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; Double_t rotation1 = 180.*atan2(Tiltdirection1[1],Tiltdirection1[0])/PI; fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\nE%d->SetFillStyle(0);\n", index,POINTS1[j+2],fi1,Aellipsis1,Bellipsis1,rotation1,index); // ------ se lo hit e' spurio marcalo in rosso index++; } // end of for( ii=0; ii<2; ii++) } // end of for( i=1; i< Nhits; i++) if(!(index+nscitilhits==0 || zmax < zmin || Smax < Smin )) { aaa = Smax-Smin; Smin -= aaa*1.; Smax += aaa*1.; aaa = zmax-zmin; zmin -= aaa*0.05; zmax += aaa*0.05; if(Smax > 2.*PI) Smax = 2.*PI; if( Smin < 0.) Smin = 0.; // Smin=0.; // Smax=2.*PI; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",zmin,Smin,zmax,Smax); for( ii=0; ii< index; ii++) { fprintf(MACRO,"E%d->Draw();\n",ii); } deltaz = zmax-zmin; deltaS = Smax-Smin; fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", zmin+0.05*deltaz,Smin+0.05*deltaS,zmax-0.05*deltaz,Smin+0.05*deltaS,zmin+0.05*deltaz,zmax-0.05*deltaz); fprintf(MACRO,"Assex->Draw();\n"); fprintf(MACRO,"TGaxis *Assey = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", zmin+0.05*deltaz,Smin+0.05*deltaS,zmin+0.05*deltaz,Smax-0.05*deltaS,Smin+0.05*deltaS,Smax-0.05*deltaS); fprintf(MACRO,"Assey->Draw();\n"); //------------ // plot di eventuali hits SciTil; for(i=0;iDraw();\n",i); disegnaSciTilHit( MACRO, i, ZETA[i], ESSE[i], 1 ); } //------------ // -------------------------------- // plot della traccia trovata dal finder if(goodskewfit) { // if ( fabs(KAPPA) > 1.e-10 && fabs(KAPPA) < 1.e10 ) { if ( fabs(KAPPA) <= 1.e-10 || fabs(KAPPA) >= 1.e10 ) { cout<<"PndSttTrackFinderReal::WriteMacroSkewAssociatedHits," <<" this track found by PR not plotted" <<"\n\t because 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* FOUND%d = new TLine(%f,%f,%f,%f);\nFOUND%d->SetLineColor(2);\nFOUND%d->Draw();\n", i-Nmin,z1,0.,z2, 2.*PI,i-Nmin,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) } // end of if ( fabs(KAPPA) <= 1.e-10 || fabs(KAPPA) >= 1.e10 ) } // end of if(goodskewfit) } // end of if(!(index+nscitilhits==0||zmax < zmin||Smax < Smin)) fprintf(MACRO,"}\n"); fclose(MACRO); } //----------end of function PndSttTrackFinderReal::WriteMacroSkewAssociatedHits //----------start of function PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithMC void PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithMC( bool goodskewfit, Double_t KAPPA, Double_t FI0, Double_t D, Double_t Fi, Double_t R, Double_t info[][7], Int_t Nincl, Int_t Minclinations[], Double_t inclination[][3], Int_t imaxima, Int_t sequentialNTrack, Short_t nSttSkewhitinTrack, Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxHitsInTrack], Short_t nSkewCommon, Short_t SkewCommonList[MAXTRACKSPEREVENT][nmaxHits], Short_t daTrackFoundaTrackMC, Short_t nMCSkewAlone[MAXTRACKSPEREVENT], Short_t MCSkewAloneList[MAXTRACKSPEREVENT][nmaxHits], Short_t nscitilhits, Double_t *ESSE, Double_t *ZETA ) { bool flaggo; Int_t i, j, i1, ii, iii, index, Kincl, nlow, nup, STATUS, imc, Nmin, Nmax; Short_t Lista[nmaxHits]; Double_t xmin , xmax, ymin, ymax, Ox, Oy, dx, dy, diff, d1, d2, delta, deltax, deltay, deltaz, deltaS, factor, zmin, zmax, Smin, Smax, S1, S2, z1, z2, y1, y2, vx1, vy1, vz1, C0x1, C0y1, C0z1, aaa, bbb, ccc, angle, minor, major, distance, Rx, Ry, LL, Aellipsis1, Bellipsis1,fi1, fmin, fmax, offset, step, SkewInclWithRespectToS, zpos, zpos1, zpos2, Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; TDatabasePDG *fdbPDG; TParticlePDG *fParticle; //------------------- skew straws hits Macro now char nome2[300],nome[300]; FILE *MACRO; sprintf(nome, "MacroSttSkewwithMCEvent%dT%d", IVOLTE,sequentialNTrack); sprintf(nome2, "%s.C",nome); MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); //KAPPA = 1./166.67 ; FI0 = 1.5*PI; Smin=zmin = 1.e10; Smax=zmax = -zmin; index=0; // prima lo (gli) hits SciTil for(i=0;iSmax ) Smax=ESSE[i]; if( ESSE[i]zmax ) zmax=ZETA[i]; if( ZETA[i] 1.e-10){ Tiltdirection1[0] = vz1/bbb; Tiltdirection1[1] = SkewInclWithRespectToS/bbb; } else { Tiltdirection1[0] = 1.; Tiltdirection1[1] = 0.; } LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; if( distance >= info[i][4] + Aellipsis1) continue; // checks that the projected ellipsis doesn't go out the boundaries of both the skew straw and the trajectory cylinder //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; Lista[index] = i+ii*10000; // do la possibilita' di plottare 2 hits skew che vengono dalla // stessa skew straw. Double_t rotation1 = 180.*atan2(Tiltdirection1[1],Tiltdirection1[0])/PI; fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\nE%d->SetFillStyle(0);\n", i+ii*10000,POINTS1[j+2],fi1,Aellipsis1,Bellipsis1,rotation1,i+ii*10000); // index,POINTS1[j+2],fi1,Aellipsis1,Bellipsis1,rotation1,index); // ------ se lo hit e' spurio marcalo in rosso flaggo=true; for( i1=0; i1SetLineColor(2);\n",i+ii*10000); index++; } // end of for( ii=0; ii<2; ii++) } // end of for( iii=0; iii< nSttSkewhitinTrack; iii++) //------------------------------- //------ aggiungo in blu eventuali punti della traccia MC che sono non mecciati for( iii=0; iii< nMCSkewAlone[imaxima]; iii++) { i = MCSkewAloneList[imaxima][iii]; Kincl = (int) info[i][5] - 1; aaa = sqrt(inclination[Kincl][0]*inclination[Kincl][0]+inclination[Kincl][1]*inclination[Kincl][1]+ inclination[Kincl][2]*inclination[Kincl][2]); vx1 = inclination[Kincl][0]/aaa; vy1 = inclination[Kincl][1]/aaa; vz1 = inclination[Kincl][2]/aaa; C0x1 = info[i][0]; C0y1 = info[i][1]; C0z1 = info[i][2]; Ox = (R+D)*cos(Fi); Oy = (R+D)*sin(Fi); calculateintersections(Ox,Oy,R,C0x1,C0y1,C0z1,info[i][3], vx1,vy1,vz1, &STATUS,POINTS1); if(STATUS < 0 ) continue ; for( ii=0; ii<2; ii++){ j=3*ii; distance = sqrt( (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) ); Rx = POINTS1[j]-Ox ; // x component Radial vector of cylinder of trajectory Ry = POINTS1[1+j]-Oy ; // y direction Radial vector of cylinder of trajectory aaa = sqrt(Rx*Rx+Ry*Ry); SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ; SkewInclWithRespectToS /= R; bbb = sqrt( SkewInclWithRespectToS*SkewInclWithRespectToS + vz1*vz1); // the tilt direction of this ellipse is (1,0) when major axis along Z direction if( bbb > 1.e-10){ Tiltdirection1[0] = vz1/bbb; Tiltdirection1[1] = SkewInclWithRespectToS/bbb; } else { Tiltdirection1[0] = 1.; Tiltdirection1[1] = 0.; } LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; if( distance >= info[i][4] + Aellipsis1) continue; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; Double_t rotation1 = 180.*atan2(Tiltdirection1[1],Tiltdirection1[0])/PI; fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\nE%d->SetFillStyle(0);\n", i+ii*10000,POINTS1[j+2],fi1,Aellipsis1,Bellipsis1,rotation1,i+ii*10000); // ------ marca lo hit in blu fprintf(MACRO,"E%d->SetLineColor(4);\n",i+ii*10000); Lista[index] = i+ii*10000; // do la possibilita' di plottare 2 hits skew che vengono dalla // stessa skew straw. index++; } // 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 if( !(index+nscitilhits==0|| zmax < zmin || Smax < Smin ) ){ aaa = Smax-Smin; Smin -= aaa*1.; Smax += aaa*1.; aaa = zmax-zmin; zmin -= aaa*0.05; zmax += aaa*0.05; if(Smax > 2.*PI) Smax = 2.*PI; if( Smin < 0.) Smin = 0.; // Smin=0.; // Smax=2.*PI; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",zmin,Smin,zmax,Smax); for( ii=0; ii< index; ii++) { // fprintf(MACRO,"E%d->Draw();\n",ii); fprintf(MACRO,"E%d->Draw();\n",Lista[ii]); } deltaz = zmax-zmin; deltaS = Smax-Smin; fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", zmin+0.05*deltaz,Smin+0.05*deltaS,zmax-0.05*deltaz,Smin+0.05*deltaS,zmin+0.05*deltaz,zmax-0.05*deltaz); fprintf(MACRO,"Assex->Draw();\n"); fprintf(MACRO,"TGaxis *Assey = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", zmin+0.05*deltaz,Smin+0.05*deltaS,zmin+0.05*deltaz,Smax-0.05*deltaS,Smin+0.05*deltaS,Smax-0.05*deltaS); fprintf(MACRO,"Assey->Draw();\n"); //------------ // plot di eventuali hits SciTil; for(i=0;iDraw();\n",i); disegnaSciTilHit( MACRO, i, ZETA[i], ESSE[i], 1 ); } //------------ // -------------------------------- // plot della traccia trovata dal finder if(goodskewfit) { if ( fabs(KAPPA) > 1.e-10 && fabs(KAPPA) < 1.e10) { fmin = KAPPA*zmin + FI0; fmax = KAPPA*zmax + FI0; 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* FOUND%d = new TLine(%f,%f,%f,%f);\nFOUND%d->SetLineColor(2);\nFOUND%d->Draw();\n", i-Nmin,z1,0.,z2, 2.*PI,i-Nmin,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) }else{// continuation of if ( fabs(KAPPA) > 1.e-10 && fabs(KAPPA) < 1.e10) cout<<"PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithMC," <<" this track found by PR not plotted" <<"\n\t because KAPPA = "<At(imc); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla fdbPDG= TDatabasePDG::Instance(); fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if( fabs(carica)>=0.1) { Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); Fifi = atan2(Cy, Cx); // MC truth Fifi angle of circle of Helix trajectory if(Fifi<0.) Fifi += 2.*PI; if( fabs( pMC->GetMomentum().Z() )< 1.e-20) KAPPA = 99999999.; else KAPPA = -carica*0.001*BFIELD*CVEL/pMC->GetMomentum().Z(); FI0 = fmod(Fifi+ PI, 2.*PI); if ( fabs(KAPPA) > 1.e-10 && fabs(KAPPA) < 1.e10) { fmin = KAPPA*zmin + FI0; fmax = 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, 2.*PI,imc,i-Nmin,imc,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) } else { // end of if ( fabs(KAPPA) > 1.e-10 && fabs(KAPPA) < 1.e10) cout<<"PndSttTrackFinderReal::WriteMacroSkewAssociatedHits, this track found by PR not plotted" <<"\n\t because KAPPA = "<=0.1) } // end of if ( pMC ) } // end of if( !(index+nscitilhits==0|| zmax < zmin || Smax < Smin ) ) fprintf(MACRO,"}\n"); fclose(MACRO); } //----------end of function PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithMC //----------begin of function PndSttTrackFinderReal::PndSttfromXYtoConformal void PndSttTrackFinderReal::PndSttFromXYtoConformal( Double_t trajectory_vertex[3], Double_t info[][7], Int_t Nparal, Double_t infoparalConformal[][5], Int_t *status ) { // do the transformation in the conformal space : u= x/(x**2+y**2), v= y/(x**2+y**2) for each hit from parallel // straws; also the equidrift radius changes. // Double_t gamma, x, y, r; for(int i=0; i nFidivConformal ) { iFi = nFidivConformal; } else if (iFi<0) { iFi = 0; } Double_t RRR = sqrt(infoparalConformal[i][0]*infoparalConformal[i][0]+ infoparalConformal[i][1]*infoparalConformal[i][1]); // cout<<"from Fillng : hit parallel n. "<0; j--){ if( RRR> radiaConf[j] ){ iR = j; break; } } if( nBoxConformal[iR][iFi] >= MAXHITSINCELL ){ cout<<"Warning from PndSttTrackFinderReal::PndSttBoxConformalFilling :" <<"\n\tcontent in nBoxConformal["<=3 && IVOLTE <= nmassimo)cout<<"Da strictcollection, iFiseed "< iFiseed if( -iFiseed + iFi <= NFiCELLDISTANCE ) { OutputListHitsinTrack[nHitsinTrack]=ListHitsinTrackinWhichToSearch[i]; nHitsinTrack++; } else { if( -iFi + nFidivConformal + iFiseed<= NFiCELLDISTANCE ) { OutputListHitsinTrack[nHitsinTrack]=ListHitsinTrackinWhichToSearch[i]; nHitsinTrack++; } } } // end of if( iFi == iFiseed ) } // end of if( InclusionList[ infoparal[ ListHitsinTrackinWhichToSearch[i] ] ] ) } // end of for(i=0; i=3 && IVOLTE <= nmassimo) { cout<<"Da strictcollection, n. elementi aggiunti "< 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 ); // 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( PndSttAcceptHitsConformal( 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( PndSttAcceptHitsConformal( 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( PndSttAcceptHitsConformal( 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 && FFimin < nRdivConformal/4 ) { iFi0 = (Short_t) (nRdivConformal/4 ); } else if ( FFimax > 3*nRdivConformal/4 && FFimin < 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( PndSttAcceptHitsConformal( 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( PndSttAcceptHitsConformal( 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 ); // 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( PndSttAcceptHitsConformal( 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 PndSttTrkAssociatedParallelHitsToHelixQuater, before exiting; nAssociatedHits = " < NTIMES*STRAWRADIUS ) continue; if(angleFi_up) continue; auxListHitsinTrack[nAssociatedHits]= i; nAssociatedHits++; } // end for(i=0; i 1.e-10){ Tiltdirection1[0] = vz1/bbb; Tiltdirection1[1] = SkewInclWithRespectToS/bbb; } else { Tiltdirection1[0] = 1.; Tiltdirection1[1] = 0.; } LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; if( distance >= info[i][4] + Aellipsis1 ) continue; // checks that the projected ellipsis doesn't go out the boundaries of both the skew straw and the trajectory cylinder // if( // fabs(POINTS1[j+2]-info[i][2]) > SEMILENGTH_STRAIGHT- Aellipsis1 || // distance + bbb > info[i][4] // the ellipsis goes out of the boundaries of the skew straw // ) { // if( istampa) cout<<"the ellipsis goes out of the boundaries of the skew straw, hit n. "< Fi_up_limit) continue; // } else { // if( Sprime > Fi_up_limit) continue; // } 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; } // second check, to see if this hits are in the allowed Fi range in the Helix reference frame /* if ( S[NAssociated] < Fi_allowedforskew_low) { if (S[NAssociated] + 2.*PI > Fi_allowedforskew_up) continue; } else { if (S[NAssociated] > Fi_allowedforskew_up) continue; } */ //--------------------------- end check Z[NAssociated] = POINTS1[j+2]; ZDrift[NAssociated] = Aellipsis1*Tiltdirection1[0]; ZErrorafterTilt[NAssociated] = StrawDriftError*aaa*Tiltdirection1[0]/LL; SkewList[NAssociated][0] = iii; // n. skew hit in skew 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< nSttSkewhit; iii++) return NAssociated; } //----------end of function PndSttTrackFinderReal::AssociateSkewHitsToXYTrack //----------begin of function PndSttTrackFinderReal::AssociateBetterAfterFitSkewHitsToXYTrack Short_t PndSttTrackFinderReal::AssociateBetterAfterFitSkewHitsToXYTrack( Short_t TemporarynSttSkewhitinTrack, // input Short_t SkewList[nmaxHits][2], // input, list of selected skew hits (in skew numbering) Double_t *S, // input, S coordinate of selected Skew hit Double_t *Z, // input, Z coordinate of selected Skew hit Double_t *ZDrift, // input, drift distance IN Z DIRECTION only, of selected Skew hit Double_t *ZErrorafterTilt, // input, Radius taking into account the tilt, IN Z DIRECTION only, of selected Skew hit Double_t KAPPA, // input, KAPPA result of fit Double_t FI0, // input, FI0 result of fit Short_t *tempore, // output result, associated skew hits Double_t *temporeS, // output, associated skew hit S Double_t *temporeZ, // output, associated skew hits Z Double_t *temporeZDrift, // output, associated skew hit Z drift Double_t *temporeZErrorafterTilt, // output, associated skew hits Z error after tilt Int_t *STATUS // output ) { Short_t NAssociated; Short_t sign; Int_t i, j, i1, ii, iii, Kincl, nlow, nup; Double_t bbb, tempZ[2], zmin, zmax, deltaz, zdist[2], zdist1, zdist2; Double_t allowed_distance = 4.*STRAWRADIUS/sin(SKEWinclination_DEGREES*PI/180.); if(fabs(KAPPA)<1.e-20) { *STATUS=-1; return 0; } NAssociated=0; if(KAPPA>0) { 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]=3){ cout<<"from FitHelixCylinder,prima di rotazione, Evento "<=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) { *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 PndSttTrackFinderReal::FitSZspace //----------begin of function PndSttTrackFinderReal::PndSttFitSZspacebis Short_t PndSttTrackFinderReal::PndSttFitSZspacebis( Short_t nSttSkewhitinTrack, Double_t *S, Double_t *Z, Double_t *DriftRadius, Double_t FInot, Short_t NMAX, Double_t *emme ) { // definition of variables for the glpsol solver // ROWS (for read_rows function) // Short_t NpointsInFit = nSttSkewhitinTrack-NMAX <0 ? nSttSkewhitinTrack : NMAX; int nRows= NpointsInFit*9 +1; int typeRows[nRows]; char * nameRows[nRows]; char auxnameRows[nRows][20]; //------- end ROWS information //--------begin COLUMNS information int NStructVar=5+NpointsInFit*4; // number of structural variables int NStructRows = 8*NpointsInFit ; // maximum number of ROWS in which a structural variable can be found double final_values[NStructVar]; int NRowsInWhichStructVarArePresent[NStructVar]; char *StructVarName[NStructVar]; char auxStructVarName[NStructVar][20]; char *NameRowsInWhichStructVarArePresent[NStructVar*NStructRows]; char aux[NStructVar*NStructRows][20]; // double Coefficients[NStructVar][NStructRows]; double Coefficients[NStructVar*NStructRows]; //--------end COLUMNS information //--------begin RHS information double ValueB[9*NpointsInFit]; //--------end RHS information //--------begin RANGES information int nRanges = NpointsInFit; double ValueRanges[nRanges]; char *NameRanges[nRanges]; char auxNameRanges[nRanges][20]; //--------end RANGES information //--------start BOUNDS information int nBounds=2*NpointsInFit+3; // int nBounds=2*NpointsInFit+1; double BoundValue[nBounds]; char *BoundStructVarName[nBounds]; char auxBoundStructVarName[nBounds][20]; char *TypeofBound[nBounds]; char auxTypeofBound[nBounds][20]; //--------end BOUNDS information //---------------------------------------------------- Double_t M = 50., m_result, q_result, A, alfetta, angle, offsety, Ox[nmaxHits], Oy[nmaxHits], Delta[nmaxHits]; Short_t i, ii; Short_t Status; char nome[300], stringa[300], stringa2[300]; // FILE * MACRO ; float m1_result,m2_result, q1_result,q2_result, A1_result, A2_result; // -- // rotation of 90 degrees for(i=0;i 1.e-10) { *emme=-1./(*emme); return 1; } else { return -99; } } //----------end of function PndSttTrackFinderReal::PndSttFitSZspacebis //----------begin of function PndSttTrackFinderReal::PndSttFitwithKalman /* void PndSttTrackFinderReal::PndSttFitwithKalman( Double_t oX, Double_t oY, Double_t Pxini, Double_t Pyini, Double_t Pzini, Double_t Ptras, Double_t info[][7], Short_t nParallelHits, Short_t *ListParallelHits, Short_t nSttSkewhit, Short_t *ListSkewHits, Double_t *S, Short_t *Infoparal, Short_t *Infoskew ) { Short_t i,j,flag, nTotal = nParallelHits+nSttSkewhit, BigList[nTotal]; Double_t old, auxFivalues[nmaxHits], auxFiSkewvalues[nmaxHits], auxRvalues[nmaxHits], BigListFi[nTotal]; // here there is the ordering of the hits // ordering of the parallel hits first for (j = 0; j< nParallelHits; j++){ auxRvalues[j]= info[ infoparal[ ListParallelHits[j] ] ][0]* info[ infoparal[ ListParallelHits[j] ] ][0]+ info[ infoparal[ ListParallelHits[j] ] ][1]* info[ infoparal[ ListParallelHits[j] ] ][1]; } Merge_Sort( nParallelHits, auxRvalues, ListParallelHits); for (j = 0; j< nParallelHits; j++){ auxFivalues[j] = atan2( info[ infoparal[ ListParallelHits[j] ] ][1]-oY, info[ infoparal[ ListParallelHits[j] ] ][0]-oX); if( auxFivalues[j] < 0. ) auxFivalues[j] += 2.*PI; } // fixing possible discontinuity between fi<2*PI and fi>0. for (old=auxFivalues[0],flag=0, j = 1; j< nParallelHits; j++){ if( fabs(old - auxFivalues[j]) > PI ) { flag=1; break; } else { old=auxFivalues[j]; } } if( flag==1) { for (j = 0; j< nParallelHits; j++){ if( auxFivalues[j] < PI) auxFivalues[j] += 2.*PI; } } // now ordering of the skew hits if( flag==1) { for (j = 0; j< nSttSkewhit; j++){ if( S[j] < PI) { auxFiSkewvalues[j] = S[j] + 2.*PI; } else { auxFiSkewvalues[j] = S[j]; } } } else { for (j = 0; j< nSttSkewhit; j++){ auxFiSkewvalues[j] = S[j]; } } Merge_Sort( nSttSkewhit, auxFiSkewvalues, ListSkewHits); // merge the parallel and skew hits for(j = 0;j< nParallelHits; j++){ BigListFi[j]=auxFivalues[j]; BigList [j] = Infoparal[ ListParallelHits[i] ] ; } for(j=0; j auxFivalues[nParallelHits]) { for(j=0; jGetMomentum(); Double_t momerr; // = 4% of mc mom // momerr = 0.04 * mcStartMom.Mag(); momerr = 0.04 * Ptras; pterr = momerr; plerr = momerr; TVector3 StartMomErr = TVector3(pterr, pterr, plerr); // Int_t pdg = mctrack->GetPdgCode(); Int_t pdg = 13; TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(pdg); Double_t fCharge= fParticle->Charge()/3.; TVector3 u(0.,1.,0.); TVector3 v(0.,0.,1.); GFDetPlane pl(StartPos,u,v); GFAbsTrackRep* rep = 0; GeaneTrackRep *grep = new GeaneTrackRep(fPro,pl,StartMom,StartPosErr,StartMomErr,fCharge,pdg); grep->setPropDir(1); // propagate in flight direction! rep=grep; GFTrack* trk = new GFTrack(rep); GFTrackCand *cand = new GFTrackCand(); int detId; for(int iPoint = 0; iPoint < nTotal; iPoint++) { PndSttHit * currenthit = (PndSttHit*) GetHitFromCollections(BigList[iPoint]); if(!currenthit) continue; detId = currenthit->GetDetectorID() ; cand->addHit(detId, BigList[iPoint]); } // end of for(int iPoint = 0; iPoint < nTotal; iPoint++) trk->setCandidate(*cand); // here the candidate is copied! //---- now the Kalman trk->addHitVector(_theRecoHitFactory->createMany(trk->getCand())); GFKalman k; k.setLazy(1); k.setNumIterations(1); k.processTrack(trk); //------ estrazione delle info dal Kalman secondo Lia TVector3 dum = trk->getCardinalRep()->getMom(); //------------ delete grep; delete trk; delete cand; return; } */ //----------end of function PndSttTrackFinderReal::PndSttFitwithKalman //----------begin of function PndSttTrackFinderReal::PndSttOrderingParallel void PndSttTrackFinderReal::PndSttOrderingParallel( Double_t oX, Double_t oY, Double_t info[][7], Short_t nParallelHits, Short_t *ListParallelHits, Short_t *Infoparal, Short_t Charge, // input Double_t *Fi_initial_helix_referenceframe, // output Double_t *Fi_final_helix_referenceframe, // output Double_t *U, // output Double_t *V // output ) { Short_t i,j, tmp[nParallelHits]; Double_t aaa, b1, firstR2, lastR2, aux[nParallelHits]; // 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 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((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< nParallelHits; j++){ U[j]=info[ Infoparal[ ListParallelHits[j] ] ][0]/( info[ Infoparal[ ListParallelHits[j] ] ][0]* info[ Infoparal[ ListParallelHits[j] ] ][0]+ info[ Infoparal[ ListParallelHits[j] ] ][1]* info[ Infoparal[ ListParallelHits[j] ] ][1]); } Merge_Sort( nParallelHits, U, ListParallelHits); 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[ Infoparal[ ListParallelHits[j] ] ][1]/( info[ Infoparal[ ListParallelHits[j] ] ][0]* info[ Infoparal[ ListParallelHits[j] ] ][0]+ info[ Infoparal[ ListParallelHits[j] ] ][1]* info[ Infoparal[ ListParallelHits[j] ] ][1]); } Merge_Sort( nParallelHits, V, ListParallelHits); 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[ Infoparal[ ListParallelHits[nParallelHits-1] ] ][1]-oY, info[ Infoparal[ 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 PndSttTrackFinderReal::PndSttOrderingParallel //----------begin of function PndSttTrackFinderReal::PndSttOrderingSkewandParallel void PndSttTrackFinderReal::PndSttOrderingSkewandParallel( Short_t *Infoparal, Short_t *Infoskew, Double_t oX, Double_t oY, Double_t Rr, Short_t nSkewhit, Short_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, Short_t nParHits, Short_t *ListParHits, Double_t *U, Double_t *V, Short_t *BigList // this is the final ordered Parallel+Skew list; // already in NATIVE hit number. ) { Short_t i,j, index[nSkewhit+nParHits], tmp[nSkewhit+nParHits], tmpList[nSkewhit]; 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]; // BigList[j]=Infoparal[ListParHits[j]]; // index[j] = j; } for (j = 0; j< nSkewhit; j++){ // this is U in conformal space aux[j+nParHits]=sign*(oX + Rr*cos(SList[Infoskew[ListSkewHits[j]]]))/ (oX*oX+oY*oY+Rr*Rr + 2.*Rr* (oX*cos(SList[Infoskew[ListSkewHits[j]]]) +oY*sin(SList[Infoskew[ListSkewHits[j]]]))); // BigList[j+nParHits]=Infoskew[ListSkewHits[j]]; // index[j+nParHits] = j+nParHits; } } 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]; // BigList[j]=Infoparal[ListParHits[j]]; // index[j] = j; } for (j = 0; j< nSkewhit; j++){ // this is V in conformal space. aux[j+nParHits]=sign*(oY + Rr*sin(SList[Infoskew[ListSkewHits[j]]]))/ (oX*oX+oY*oY+Rr*Rr + 2.*Rr* (oX*cos(SList[Infoskew[ListSkewHits[j]]]) +oY*sin(SList[Infoskew[ListSkewHits[j]]]))); // BigList[j+nParHits]=Infoskew[ListSkewHits[j]]; // index[j+nParHits] = j+nParHits; } } // end of if((aaa>b1&& .... for (j = 0 ; j< nParHits; j++){ BigList[j]=Infoparal[ListParHits[j]]; index[j] = j; } for (j = 0; j< nSkewhit; j++){ BigList[j+nParHits]=Infoskew[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;i0. for (old=auxFivalues[0],flag=0, j = 1; j< nParallelHits; j++){ if( fabs(old - auxFivalues[j]) > PI ) { flag=1; break; } else { old=auxFivalues[j]; } } if( flag==1) { for (j = 0; j< nParallelHits; j++){ if( auxFivalues[j] < PI) auxFivalues[j] += 2.*PI; } } // now ordering of the skew hits if(nSkewHit>0){ if( flag==1) { for (j = 0; j< nSkewHit; j++){ if( S[Infoskew[ ListSkewHits[j] ]] < PI) { auxFiSkewvalues[j] = S[Infoskew[ ListSkewHits[j] ]] + 2.*PI; } else { auxFiSkewvalues[j] = S[Infoskew[ ListSkewHits[j] ]]; } } } else { for (j = 0; j< nSkewHit; j++){ auxFiSkewvalues[j] = S[Infoskew[ ListSkewHits[j] ]]; } } Merge_Sort( nSkewHit, auxFiSkewvalues, ListSkewHits); } // end of if(nSkewHit>0) // merge the parallel and skew hits for(j = 0;j< nParallelHits; j++){ BigListFi[j]=auxFivalues[j]; BigList[j] = Infoparal[ ListParallelHits[j] ] ; } for(j=0; j0) Merge_Sort(*nTotal, BigListFi, BigList); //--------- end ordering return; } //----------end of function PndSttTrackFinderReal::PndSttOrdering //----------start function PndSttTrackFinderReal::PndSttFindingParallelTrackAngularRange void PndSttTrackFinderReal::PndSttFindingParallelTrackAngularRange( Double_t oX, Double_t oY, Double_t R, Short_t Charge, Double_t *Fi_low_limit, // Fi (in XY Helix frame) lower limit using // the Stt detector minimum/maximum radius // Fi_low_limit is ALWAYS between 0. and 2PI Double_t *Fi_up_limit,// Fi (in XY Helix frame) upper limit using // the Stt detector maximum/minimum radius // Fi_up_limit is ALWAYS > Fi_low_limit and // possibly > 2PI. Short_t * status,// *status =0, all well; =1, // track contained completely between RMin // and RMax; = -1 track contained within RMin; // =-2 track outside RMax. Double_t Rmi, // Rmin of cylindrical volume intersected by track; Double_t Rma // Rmax of cylindrical volume intersected by track; ) { // The hits ALWAYS are contained between Fi_low_limit and Fi_up_limit; // the Point at (0,0) is NEVER contained between Fi_low_limit and Fi_up_limit. // -------------- 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. // working in the hypothesis that the starting point of the track is near (0,0) so that // R_vertex < RStrawDetectorMin 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 + R <= Rmi ) // in this case there might be hits at radius < Rmi. { *status = -1 ;return;} if( a >= R + Rma || R >= a + Rma) // in this case there can be no hits with radius < Rma. { *status = -2;return;} if( a - R >= Rmi ) intersection_inner = false; else intersection_inner = true; if( a + R <= Rma || a - R >= 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 + R*R - Rma*Rma)/(2.*R*a); if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.; Fi = acos(cosFi); } if( intersection_inner ){ cosfi = (a*a + R*R - Rmi*Rmi)/(2.*R*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 PndSttTrackFinderReal::PndSttFindingParallelTrackAngularRange //----------start of function PndSttTrackFinderReal::WriteMacroParallelHitswithRfromMC void PndSttTrackFinderReal::WriteMacroParallelHitswithRfromMC( Int_t Nhits, Double_t info[][7], Short_t nTracksFoundSoFar, // bool *TypeConf, Double_t *Ox, Double_t *Oy, Short_t * daTrackFoundaTrackMC ) { 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 con anche le tracce MC sprintf(nome,"MacroSttwithRfromMCParallelHitswithMCEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); // 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++) { if( info[i][5] == 1 ) { // parallel straws 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]; } } 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); disegnaAssiXY(MACRO,xmin,xmax,ymin,ymax); for( i=0; i< Nhits; i++) { 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); } } for( i=0; i< nMCTracks ; i++) { Int_t icode; Double_t Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Px, Py, carica, KAPPA, FI0mc ; PndMCTrack* pMC; pMC = (PndMCTrack*) fMCTrackArray->At(i); if ( ! pMC ) continue ; icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if( fabs(carica)<0.1) continue; Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); if( fabs( pMC->GetMomentum().Z() )< 1.e-20) KAPPA = 99999999.; else KAPPA = -carica*0.001*BFIELD*CVEL/pMC->GetMomentum().Z(); FI0mc = fmod(Fifi+ PI, 2.*PI); fprintf(MACRO,"TEllipse* MC%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nMC%d->SetFillStyle(0);\nMC%d->SetLineColor(3);\nMC%d->Draw();\n", i,Cx,Cy,Rr,Rr,i,i,i); } //------------------------------- plotting all the tracks found, with R from corresponding MC track for(i=0; iAt(daTrackFoundaTrackMC[i]); if ( ! pMC ) continue ; icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if( fabs(carica)<0.1) continue; Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); if( fabs( pMC->GetMomentum().Z() )< 1.e-20) KAPPA = 99999999.; else KAPPA = -carica*0.001*BFIELD*CVEL/pMC->GetMomentum().Z(); FI0mc = fmod(Fifi+ PI, 2.*PI); rrr = Rr; double angolo = atan2(Oy[i],Ox[i]); aaa = rrr*cos(angolo); bbb = rrr*sin(angolo); fprintf(MACRO,"TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,0.,360.);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", i,aaa,bbb,rrr,rrr,i,i,i); } // ----------- fprintf(MACRO,"}\n"); fclose(MACRO); //------------------------------------------------------------------------------------------------------------ return ; } //----------end of function PndSttTrackFinderReal::WriteMacroParallelHitswithRfromMC //----------start of function PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithRfromMC void PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithRfromMC( Double_t KAPPA,Double_t FI0,Double_t D,Double_t Fi,Double_t R, Int_t Nhits, Double_t info[][7], Int_t Nincl, Int_t Minclinations[], Double_t inclination[][3], Int_t imaxima, Int_t nMaxima ) { Int_t i, j, i1, ii, index, Kincl, nlow, nup, STATUS, imc, Nmin, Nmax; Double_t xmin , xmax, ymin, ymax, Ox, Oy, dx, dy, diff, d1, d2, delta, deltax, deltay, deltaz, deltaS, factor, zmin, zmax, Smin, Smax, S1, S2, z1, z2, y1, y2, vx1, vy1, vz1, C0x1, C0y1, C0z1, aaa, bbb, ccc, angle, minor, major, distance, Rx, Ry, LL, Aellipsis1, Bellipsis1,fi1, fmin, fmax, offset, step, SkewInclWithRespectToS, zpos, zpos1, zpos2, Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; //------------------- skew straws hits Macro now char nome2[300],nome[300]; FILE *MACRO; sprintf(nome, "MacroSttParTrack%dSkewTrack%dSkewHitswithRfromMCEvent%d",imaxima,nMaxima, IVOLTE); sprintf(nome2, "%s.C",nome); MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); //KAPPA = 1./166.67 ; FI0 = 1.5*PI; Smin=zmin = 1.e10; Smax=zmax = -zmin; index=0; for( i=0; i< Nhits; i++) { if( info[i][5] == 1. ) continue; // exclude parallel straws Kincl = (int) info[i][5] - 1; aaa = sqrt(inclination[Kincl][0]*inclination[Kincl][0]+inclination[Kincl][1]*inclination[Kincl][1]+ inclination[Kincl][2]*inclination[Kincl][2]); vx1 = inclination[Kincl][0]/aaa; vy1 = inclination[Kincl][1]/aaa; vz1 = inclination[Kincl][2]/aaa; C0x1 = info[i][0]; C0y1 = info[i][1]; C0z1 = info[i][2]; Ox = (R+D)*cos(Fi); Oy = (R+D)*sin(Fi); calculateintersections(Ox,Oy,R,C0x1,C0y1,C0z1,info[i][3], vx1,vy1,vz1, &STATUS,POINTS1); if(STATUS < 0 ) continue ; for( ii=0; ii<2; ii++){ j=3*ii; distance = sqrt( (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) ); Rx = POINTS1[j]-Ox ; // x component Radial vector of cylinder of trajectory Ry = POINTS1[1+j]-Oy ; // y direction Radial vector of cylinder of trajectory aaa = sqrt(Rx*Rx+Ry*Ry); SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ; SkewInclWithRespectToS /= R; bbb = sqrt( SkewInclWithRespectToS*SkewInclWithRespectToS + vz1*vz1); // the tilt direction of this ellipse is (1,0) when major axis along Z direction if( bbb > 1.e-10){ Tiltdirection1[0] = vz1/bbb; Tiltdirection1[1] = SkewInclWithRespectToS/bbb; } else { Tiltdirection1[0] = 1.; Tiltdirection1[1] = 0.; } LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; if( distance >= info[i][4] + Aellipsis1) continue; // checks that the projected ellipsis doesn't go out the boundaries of both the skew straw and the trajectory cylinder if( fabs(POINTS1[j+2]-ZCENTER_STRAIGHT) > SEMILENGTH_STRAIGHT- Aellipsis1 || distance + bbb > info[i][4] // the ellipsis goes out of the boundaries of the skew straw ) { cout<<"the ellipsis goes out of the boundaries of the skew straw, hit n. "< 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* E%d = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\nE%d->SetFillStyle(0);\n", index,POINTS1[j+2],fi1,Aellipsis1,Bellipsis1,rotation1,index); index++; } // end of for( ii=0; ii<2; ii++) } // end of for( i=1; i< Nhits; i++) if(index!=0 && zmax >= zmin && Smax >= Smin ){ aaa = Smax-Smin; Smin -= aaa*1.; Smax += aaa*1.; aaa = zmax-zmin; zmin -= aaa*0.05; zmax += aaa*0.05; if(Smax > 2.*PI) Smax = 2.*PI; if( Smin < 0.) Smin = 0.; // Smin=0.; // Smax=2.*PI; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",zmin,Smin,zmax,Smax); for( ii=0; ii< index; ii++) { fprintf(MACRO,"E%d->Draw();\n",ii); } deltaz = zmax-zmin; deltaS = Smax-Smin; fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", zmin+0.05*deltaz,Smin+0.05*deltaS,zmax-0.05*deltaz,Smin+0.05*deltaS,zmin+0.05*deltaz,zmax-0.05*deltaz); fprintf(MACRO,"Assex->Draw();\n"); fprintf(MACRO,"TGaxis *Assey = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", zmin+0.05*deltaz,Smin+0.05*deltaS,zmin+0.05*deltaz,Smax-0.05*deltaS,Smin+0.05*deltaS,Smax-0.05*deltaS); fprintf(MACRO,"Assey->Draw();\n"); // -------------------------------- // plot della traccia trovata dal finder 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* FOUND%d = new TLine(%f,%f,%f,%f);\nFOUND%d->SetLineColor(2);\nFOUND%d->Draw();\n", i-Nmin,z1,0.,z2, 2.*PI,i-Nmin,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) //--------------------------------------------- qui ci aggiungo le traccie MC for(imc=0; imcAt(imc); if ( ! pMC ) continue ; icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if( fabs(carica)<0.1) continue; Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); if( fabs( pMC->GetMomentum().Z() )< 1.e-20) KAPPA = 99999999.; else KAPPA = -carica*0.001*BFIELD*CVEL/pMC->GetMomentum().Z(); FI0mc = fmod(Fifi+ PI, 2.*PI); // KAPPA=MCtruthTrkInfo[12][imc]; FI0 = FI0mc; 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, 2.*PI,imc,i-Nmin,imc,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) } // end of for(imc=0; imc= zmin && Smax >= Smin ) fprintf(MACRO,"}\n"); fclose(MACRO); } //----------end of function PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithRfromMC //----------begin of function PndSttTrackFinderReal::AssociateFoundTrackstoMC void PndSttTrackFinderReal::AssociateFoundTrackstoMC( Double_t info[][7], Short_t nTracksFoundSoFar, Short_t *nHitsinTrack, Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHits], Short_t *nSttSkewhitinTrack, Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxHits], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] // Short_t *daMCTrackaTrackFound ) { bool firstime, flaggo, inclusionMC[nTracksFoundSoFar][nmaxHits], inclusionExp[nTracksFoundSoFar]; Short_t ntoMCtrack[nTracksFoundSoFar], toMCtracklist[nTracksFoundSoFar][nmaxHits], toMCtrackfrequency[nTracksFoundSoFar][nmaxHits]; Short_t i, j, jtemp,jexp; Short_t enne, 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; inclusionMC[jtemp][itemp]=false; } } // end while ( itemp > -1) return; } //----------end of function PndSttTrackFinderReal::AssociateFoundTrackstoMC //----------begin of function PndSttTrackFinderReal::AssociateFoundTrackstoMCbis void PndSttTrackFinderReal::AssociateFoundTrackstoMCbis( bool *keepit, Double_t info[][7], Short_t nTracksFoundSoFar, Short_t nHitsinTrack[MAXTRACKSPEREVENT], Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHitsInTrack], Short_t nSttSkewhitinTrack[MAXTRACKSPEREVENT], Short_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxHitsInTrack], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool firstime, flaggo, inclusionMC[nTracksFoundSoFar][nmaxHits], inclusionExp[nTracksFoundSoFar]; Short_t ntoMCtrack[nTracksFoundSoFar], toMCtrackfrequency[nTracksFoundSoFar][nmaxHits]; Short_t i, j, jtemp,jexp; Short_t enne, itemp, massimo, toMCtracklist[nTracksFoundSoFar][nmaxHits]; 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 PndSttTrackFinderReal::AssociateFoundTrackstoMCbis //---------- begin of function PndSttTrackFinderReal::PndSttInfoXYZParal void PndSttTrackFinderReal::PndSttInfoXYZParal( Double_t info[][7], Short_t infopar, Double_t Ox, Double_t Oy, Double_t R, Double_t KAPPA, Double_t FI0, Short_t Charge, Double_t *Posiz // output. ) { // Posiz = position (on the drift radius) of the hit whose 'parallel' scheme // number is infopar. Double_t fi, norm, vers[2]; vers[0] = Ox - info[infopar][0]; vers[1] = Oy - info[infopar][1]; norm = sqrt( vers[0]*vers[0] + vers[1]*vers[1] ); if(norm < 1.e-20) { Posiz[0] = -999999999.; return; } if( fabs( R - fabs( norm - info[infopar][3] ) ) // distance trajectory-drift radius < fabs( R - (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( R - 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-10 ){ Posiz[2] = -888888888.; return; } else if ( fabs(KAPPA) > 1.e10){ Posiz[2] = -777777777.; 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 PndSttTrackFinderReal::PndSttInfoXYZParal //---------- begin of function PndSttTrackFinderReal::PndSttInfoXYZSkew void PndSttTrackFinderReal::PndSttInfoXYZSkew( Double_t Z, // Z coordinate (center wire) of selected Skew hit Double_t ZDrift, // drift distance IN Z DIRECTION only, of selected Skew hit Double_t S, Double_t Ox, Double_t Oy, Double_t R, Double_t KAPPA, Double_t FI0, Short_t Charge, Double_t *Posiz // output ) { Short_t sign; Double_t bbb, tempZ[2], zmin, zmax, deltaz, zdist[2], zdist1, zdist2; Double_t Zline; if(fabs(KAPPA)< 1.e-10) { Posiz[0] = -999999999.; return; } else if (fabs(KAPPA)> 1.e10) { Posiz[0] = -888888888.; return; } Posiz[0] = Ox + R * cos(S); Posiz[1] = Oy + R * sin(S); if( Charge > 0 ) { if( S > FI0 ) { if(istampa>=3){ cout<<"from PndSttInfoXYZSkew : inconsistency, FI0 is not the maximum for this track "< 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 PndSttTrackFinderReal::IntersectionCircle_Segment //----------begin of function PndSttTrackFinderReal::IntersectionsWithGapSemicircle Short_t PndSttTrackFinderReal::IntersectionsWithGapSemicircle( Double_t Oxx, Double_t Oyy, Double_t Rr, Double_t GAP, bool left, // if true --> Left semicircle; false --> Right semicircle. Double_t Rma, Double_t *XintersectionList, Double_t *YintersectionList ) { Short_t nIntersectionsCircle; Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; nIntersectionsCircle=0; aaa = sqrt(Oxx*Oxx+Oyy*Oyy); // preliminary condition for having intersections between trajectory and Circle. if( !( aaa >= Rr + Rma || Rr >= aaa + Rma) && !( aaa + Rr <= Rma || aaa - Rr >= Rma ) ) { // now the calculation FI0 = atan2(-Oyy,-Oxx); cosFi = (aaa*aaa + Rr*Rr - Rma*Rma)/(2.*Rr*aaa); if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.; Fi = acos(cosFi); // (x1, y1) and (x2,y2) are the intersections between the trajectory // and the circle, in the laboratory reference frame. x1 = Oxx+Rr*cos(FI0 - Fi); y1 = Oyy+Rr*sin(FI0 - Fi); theta1 = atan2(y1,x1); // in this way theta1 is between -PI and PI. x2 = Oxx+Rr*cos(FI0 + Fi); y2 = Oyy+Rr*sin(FI0 + Fi); theta2 = atan2(y2,x2); // in this way theta2 is between -PI and PI. // Theta1, Theta2 = angle of the edges of the outer circle + Gap in the laboratory frame. // Theta1, Theta2 are angles between -PI/2 and +PI/2. if(!left){ // Right (looking into the beam) Semicircle. Theta2 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.); Theta1 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.); if( Theta1<= theta1 && theta1 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x1; YintersectionList[nIntersectionsCircle]=y1; nIntersectionsCircle++; } if( Theta1<= theta2 && theta2 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x2; YintersectionList[nIntersectionsCircle]=y2; nIntersectionsCircle++; } } else { // Left (looking into the beam) Semicircle. Theta2 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.); Theta1 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.); if( Theta1<= theta1 || theta1 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x1; YintersectionList[nIntersectionsCircle]=y1; nIntersectionsCircle++; } if( Theta1<= theta2 || theta2 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x2; YintersectionList[nIntersectionsCircle]=y2; nIntersectionsCircle++; } } } // end of if (!( a >= Rr + Rma || ..... //---------------------------- end of calculation intersection with outer circle. return nIntersectionsCircle; } //----------end of function PndSttTrackFinderReal::IntersectionsWithGapSemicircle //----------star of function PndSttTrackFinderReal::IsInternal bool PndSttTrackFinderReal::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 PndSttTrackFinderReal::IsInternal //----------star of function PndSttTrackFinderReal::ChooseEntranceExit void PndSttTrackFinderReal::ChooseEntranceExit( Double_t Oxx, Double_t Oyy, Short_t flag, Short_t Charge, Double_t FiStart, Short_t nIntersections[2], Double_t XintersectionList[][2], // second index =1 -->inner polygon; Double_t YintersectionList[][2], // second index =2 -->outer polygon. Double_t Xcross[2], // output Double_t Ycross[2] // output ) { Short_t i, j; // this method works under the hypothesis that flag=0 or 2 only. if(flag == 0) { if(Charge > 0) { for(j=0;j<2;j++){ // j=0 --> inner polygon; j=1 --> outer polygon. Short_t auxIndex[nIntersections[j]]; Double_t fi[nIntersections[j]]; for( i=0;i FiStart) fi[i] -= 2.*PI; if( fi[i] > FiStart) fi[i] = FiStart; auxIndex[i]=i; } // end of for( i=0;i inner polygon; j=1 --> outer polygon. Short_t auxIndex[nIntersections[j]]; Double_t fi[nIntersections[j]]; for( i=0;i 0) { for( i=0;i FiStart) fi[i] -= 2.*PI; if( fi[i] > FiStart) fi[i] = FiStart; auxIndex[i]=i; } // end of for( i=0;iinner polygon; Double_t *YintersectionList, // second index =2 -->outer polygon. Double_t Xcross[2], // output Double_t Ycross[2] // output ) { Short_t i, j; // this method works under the hypothesis that there are at least 2 intersections. if (nIntersections<2) return; if(Charge > 0) { Short_t auxIndex[nIntersections]; Double_t fi[nIntersections]; for( i=0;i FiStart) fi[i] -= 2.*PI; if( fi[i] > FiStart) fi[i] = FiStart; auxIndex[i]=i; } // end of for( i=0;i2.*ApotemaMaxInnerParStraw/sqrt(3.) ){ // outer Parallel hit. ListOuterHits[ *nOuterHits ] = ListHits[ihit]; (*nOuterHits)++; if(info[ListHits[ihit]][0]<0.){ ListOuterHitsLeft[ *nOuterHitsLeft ] = ListHits[ihit]; (*nOuterHitsLeft)++; } else { ListOuterHitsRight[ *nOuterHitsRight ] = ListHits[ihit]; (*nOuterHitsRight)++; } }else{ ListInnerHits[ *nInnerHits ] = ListHits[ihit]; (*nInnerHits)++; if(info[ListHits[ihit]][0]<0.){ ListInnerHitsLeft[ *nInnerHitsLeft ] = ListHits[ihit]; (*nInnerHitsLeft)++; } else { ListInnerHitsRight[ *nInnerHitsRight ] = ListHits[ihit]; (*nInnerHitsRight)++; } } } } //----------end of function PndSttTrackFinderReal::SeparateInnerOuterParallel //----------begin of function PndSttTrackFinderReal::FindTrackEntranceExitbiHexagon Short_t PndSttTrackFinderReal::FindTrackEntranceExitbiHexagon( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t ApotemaMax, Double_t Xcross[2], Double_t Ycross[2] ) { Short_t flag; Short_t nIntersections[2]; Double_t FiStart, XintersectionList[12][2], // first index =0 --> inner Hexagon, =1 --> outer. YintersectionList[12][2]; // second index : all the possible intersections // (up to 12 intersections). // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -2 --> 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. flag=IntersectionsWithClosedPolygon( Oxx, Oyy, Rr, ApotemaMin, // Rmin of the inner part of parallele straws, ApotemaMax,// max Apotema of the inner part of parallele straws. 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 || flag == 2)) return flag; if (flag == 2 &&nIntersections[0]<2 ){ cout<<"PndSttTrackFinderReal::FindTrackEntranceExitbiHexagon,"<< " contraddiction, nIntersections[0]="<< nIntersections[0]<<"<2, returning -99!\n"; return -99; } //------- 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); // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. ChooseEntranceExit( Oxx, Oyy, flag, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); return flag; } //----------end of function PndSttTrackFinderReal::FindTrackEntranceExitbiHexagon //----------begin of function PndSttTrackFinderReal::FindTrackEntranceExitbiHexagonLeft Short_t PndSttTrackFinderReal::FindTrackEntranceExitbiHexagonLeft( Double_t vgap, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t ApotemaMax, Double_t Xcross[2], Double_t Ycross[2] ) { Short_t flag; Short_t nIntersections; Double_t FiStart, XintersectionList[16], // all the possible intersections YintersectionList[16]; // (up to 16 intersections). // The following is the form of the Left (looking from downstream into the beam) biHexagon // geometrical shape considered in this method : // /* /| / | / | / / / / / / / / | | | | | | | | \ \ \ \ \ \ \ \ \ | \ | \| */ // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; flag=IntersectionsWithClosedbiHexagonLeft( vgap, Oxx, Oyy, Rr, ApotemaMin, // Apotema of the inner Hexagon, ApotemaMax,// Apotema of the outer Hexagon. &nIntersections, XintersectionList, // XintersectionList[..][0] --> inner polygon, // XintersectionList[..][1] --> outer polygon. YintersectionList ); // IMPORTANT : // this is true because here it is assumed that the track comes from (0,0,0) // otherwise the code must be changed! if (!(flag == 0)) return flag; if( nIntersections<2) return -1; //------- among all possible intersection find the entrance point (Xcross[0], Ycross[0]) // and the exit point (Xcross[1], Ycross[1]) of this track. //-------- the starting point of the track. FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); // 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 PndSttTrackFinderReal::FindTrackEntranceExitbiHexagonLeft //----------begin of function PndSttTrackFinderReal::FindTrackEntranceExitbiHexagonRight Short_t PndSttTrackFinderReal::FindTrackEntranceExitbiHexagonRight( Double_t vgap, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t ApotemaMax, Double_t Xcross[2], Double_t Ycross[2] ) { Short_t flag; Short_t nIntersections; Double_t FiStart, XintersectionList[16], // all the possible intersections YintersectionList[16]; // (up to 16 intersections). // The following is the form of the Left (looking from downstream into the beam) biHexagon // geometrical shape considered in this method : // /* |\ | \ | \ \ \ \ \ | | | | / / / / / / | / | / |/ */ // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; flag=IntersectionsWithClosedbiHexagonRight( vgap, Oxx, Oyy, Rr, ApotemaMin, // Apotema of the inner Hexagon, ApotemaMax,// Apotema of the outer Hexagon. &nIntersections, XintersectionList, // 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); // 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 PndSttTrackFinderReal::FindTrackEntranceExitbiHexagonRight //----------begin of function PndSttTrackFinderReal::FindTrackEntranceExitHexagonCircle Short_t PndSttTrackFinderReal::FindTrackEntranceExitHexagonCircle( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t ApotemaMax, Double_t Xcross[2], Double_t Ycross[2] ) { Short_t flag; Short_t nIntersections[2]; Double_t FiStart, XintersectionList[12][2], // second index =0 --> inner Hexagon, =1 --> outer. YintersectionList[12][2]; // first index : all the possible intersections // (up to 12 intersections). // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -2 --> 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. flag=IntersectionsWithClosedPolygon( Oxx, Oyy, Rr, ApotemaMin, // Rmin of the inner part of parallele straws, ApotemaMax,// max Apotema of the inner part of parallele straws. 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 || flag == 2)) return flag; if (flag == 2 &&nIntersections[0]<2 ){ cout<<"PndSttTrackFinderReal::FindTrackEntranceExitbiHexagon,"<< " contraddiction, nIntersections[0]="<< nIntersections[0]<<"<2, returning -99!\n"; return -99; } //------- 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); // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. ChooseEntranceExit( Oxx, Oyy, flag, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); //---------------------- return flag; } //----------end of function PndSttTrackFinderReal::FindTrackEntranceExitHexagonCircle //----------begin of function PndSttTrackFinderReal::FindTrackEntranceExitHexagonCircleLeft Short_t PndSttTrackFinderReal::FindTrackEntranceExitHexagonCircleLeft( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t Rma, // outer radius of the Circle. Double_t GAP, Double_t Xcross[2], Double_t Ycross[2] ) { // This methods finds the intersections between a trajectory coming from (0,0) and // parameters Oxx, Oyy, Rr with the closed geometrical figure (in XY) formed by the STT // external Left semicircle and the Outer STT parallel Left straw semi-Hexagon + Gap for // the pellet target target. // It returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections. Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; //------------------ Short_t nIntersectionsCircle, nIntersections; Double_t FiStart, XintersectionList[12], YintersectionList[12]; // all the possible intersections (up to 12 intersections). // finding all possible intersections with inner parallel straw region. Double_t Side_x[] = { -GAP/2., -GAP/2. , -ApotemaMin, -ApotemaMin, -GAP/2., -GAP/2. }, Side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-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. 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); // 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 PndSttTrackFinderReal::FindTrackEntranceExitHexagonCircleLeft //----------begin of function PndSttTrackFinderReal::FindTrackEntranceExitHexagonCircleRight Short_t PndSttTrackFinderReal::FindTrackEntranceExitHexagonCircleRight( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t Rma, // outer radius of the Circle. Double_t GAP, Double_t Xcross[2], Double_t Ycross[2] ) { // This methods finds the intersections between a trajectory coming from (0,0) and // parameters Oxx, Oyy, Rr with the closed geometrical figure (in XY) formed by the STT // external Left semicircle and the Outer STT parallel Left straw semi-Hexagon + Gap for // the pellet target target. // It returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections. Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; //------------------ Short_t nIntersectionsCircle, nIntersections; Double_t FiStart, XintersectionList[12], YintersectionList[12]; // all the possible intersections (up to 10 intersections). // finding all possible intersections with inner parallel straw region. Double_t Side_x[] = { GAP/2., GAP/2. , ApotemaMin, ApotemaMin, GAP/2., GAP/2. }, Side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-GAP)/sqrt(3.), ApotemaMin/sqrt(3.), -ApotemaMin/sqrt(3.), -(2.*ApotemaMin-GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)}, a[] = {1., 1./sqrt(3.), 1., -1./sqrt(3.), 1.}, b[] = {0., 1., 0., 1., 0.}, c[] = {-GAP/2.,-2.*ApotemaMin/sqrt(3.),-ApotemaMin, 2.*ApotemaMin/sqrt(3.), -GAP/2.}; nIntersections=IntersectionsWithOpenPolygon( Oxx, Oyy, Rr, 5, // n. Sides of open Polygon. a, // coefficient of formula : aX + bY + c = 0 defining the Polygon sides. b, c, Side_x, // X,Y coordinate of the Sides vertices (in sequence, following Side_y, // the Polygon along. XintersectionList, // XintersectionList YintersectionList // YintersectionList. ); //------------------------------------------------------------------------- // finding intersections of trajectory [assumed to originate from (0,0) ] // with outer semicircle, the Left part. nIntersectionsCircle=IntersectionsWithGapSemicircle( Oxx, // input from trajectory Oyy, // input from trajectory Rr, // input from trajectory GAP, // input, vertical gap in XY plane of STT detector. false, // true --> Left semicircle, false --> Right semicircle. Rma, // radius of external Stt containment. &XintersectionList[nIntersections], // output, X list of intersections (maximal 2). &YintersectionList[nIntersections] ); nIntersections += nIntersectionsCircle; //-------- the starting point of the track. if(nIntersections<2) return -1; FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); // 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 PndSttTrackFinderReal::FindTrackEntranceExitHexagonCircleRight //----------begin of function PndSttTrackFinderReal::FindIntersectionsOuterCircle Short_t PndSttTrackFinderReal::FindIntersectionsOuterCircle( Double_t oX, Double_t oY, Double_t R, Double_t Rma, Double_t Xcross[2], Double_t Ycross[2] ) { // return -1 --> non intersection; // return 0 --> 2 intersections. Double_t a, cosFi, Fi, FI0; a = sqrt(oX*oX+oY*oY); // case with no intersections. if( a >= R + Rma || R >= a + Rma || a + R <= Rma) return -1; FI0 = atan2(-oY,-oX); cosFi = (a*a + R*R - Rma*Rma)/(2.*R*a); if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.; Fi = acos(cosFi); Xcross[0] = oX + R*cos(FI0+Fi); Ycross[0] = oY + R*sin(FI0+Fi); Xcross[1] = oX + R*cos(FI0-Fi); Ycross[1] = oY + R*sin(FI0-Fi); return 0; } //----------end of function PndSttTrackFinderReal::FindIntersectionsOuterCircle //----------begin of function PndSttTrackFinderReal::CalculateArcLength Double_t PndSttTrackFinderReal::CalculateArcLength( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t charge, Double_t Xcross[2], // entrance-exit point Double_t Ycross[2] // entrance-exit point ) { Double_t dis, theta1, theta2; theta1 = atan2(Ycross[0]-Oyy, Xcross[0]- Oxx); theta2 = atan2(Ycross[1]-Oyy, Xcross[1]- Oxx); if(charge>0){ // the rotation was clockwise. dis = theta1-theta2; } else { // the rotation was counterclockwise. dis = theta2-theta1; } if(dis<0.) dis += 2.*PI; if(dis<0.) dis =0.; dis *= Rr; return dis; } //----------end of function PndSttTrackFinderReal::CalculateArcLength //----------begin of function PndSttTrackFinderReal::IsInsideCircle bool PndSttTrackFinderReal::IsInsideArc( Double_t Oxx, Double_t Oyy, Short_t Charge, Double_t Xcross[2], Double_t Ycross[2], Double_t f // f should be between 0 and 2PI. ) { Double_t f1, f2; // Xcross[0],Ycross[0] is the point of entrance. f1 = atan2(Ycross[0]-Oyy, Xcross[0]-Oxx); if(f1<0.) f1+= 2.*PI; if(f1<0.) f1= 0.; f2 = atan2(Ycross[1]-Oyy, Xcross[1]-Oxx); if(f2<0.) f2+= 2.*PI; if(f2<0.) f2= 0.; if(Charge<0){ if(f1 > f2 ) f2 +=2.*PI; if(f1 > f2 ) f2 = f1; if( f>f1){ if(f 0. if(f1 < f2 ) f1 +=2.*PI; if(f1 < f2 ) f1 = f2; if( f>f2){ if(f 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 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((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++){ U[j]=XY[j][0]/(XY[j][0]*XY[j][0]+XY[j][1]*XY[j][1]); } 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++){ V[j]=XY[j][1]/(XY[j][0]*XY[j][0]+XY[j][1]*XY[j][1]); } 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 PndSttTrackFinderReal::OrderingUsingConformal //----------begin of function PndSttTrackFinderReal::FindTrackInXYProjection bool PndSttTrackFinderReal::FindTrackInXYProjection( Short_t iHit, // seed hit; if it is negative it is a SciTil hit. Short_t nRcell, // R cell of the seed hit; can be negative when SciTil hit. Short_t nFicell, // Fi cell of the seed hit; Int_t *Minclinations, Double_t info[nmaxHits][7], bool *InclusionList, Short_t *RConformalIndex, Short_t *FiConformalIndex, Short_t nBoxConformal[nRdivConformal][nFidivConformal], Short_t HitsinBoxConformal[MAXHITSINCELL][nRdivConformal][nFidivConformal], Short_t nTracksFoundSoFar, Short_t *nHitsinTrack, Short_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHitsInTrack], Double_t *trajectory_vertex, Double_t infoparalConformal[nmaxHits][5], Double_t posizSciTilx, Double_t posizSciTily, Double_t *S, Double_t *Ox, Double_t *Oy, Double_t *R, Double_t *Fi_low_limit, Double_t *Fi_up_limit, Double_t *Fi_initial_helix_referenceframe, Double_t *Fi_final_helix_referenceframe, Short_t * Charge, Double_t *U, Double_t *V ) { //--------------- Short_t i, j, Naux, Nbaux, NN, Nouter, auxListHitsinTrack[nmaxHits], OutputListHitsinTrack[nmaxHits], OutputList2HitsinTrack[nmaxHits], ListHitsinTrackinWhichToSearch[nmaxHits]; Short_t flagStt, status; Double_t aaa, m, q, rotationangle, rotationcos, rotationsin, auxinfoparalConformal[nmaxHits+1][5]; //---------------- nHitsinTrack[nTracksFoundSoFar] = PndSttFindTrackPatterninBoxConformal( 1, // distance in R cells allowed 2, // distance in Fi cells allowed Minclinations[0], iHit, // seed hit; if it is negative it is a SciTil hit. nRcell, nFicell, info, InclusionList, RConformalIndex, FiConformalIndex, nBoxConformal, HitsinBoxConformal, &ListHitsinTrack[nTracksFoundSoFar][0] ); //---------stampe. if(istampa>0){ cout<<"PndSttTrackFinderReal, evt. "<nmaxHitsInTrack) { 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= MINIMUMOUTERHITSPERTRACK) { for(i=0; i< Nouter;i++){ ListHitsinTrackinWhichToSearch[i]= ListHitsinTrack[nTracksFoundSoFar][i]; } for(i=0; i< Nouter;i++){ Naux = PndSttFindTrackPatterninBoxConformalSpecial( 3, // NRCELLDISTANCE 1, // NFiCELLDISTANCE Minclinations[0], Nouter, ListHitsinTrackinWhichToSearch[i], // seed hit. ListHitsinTrackinWhichToSearch, info, InclusionList, RConformalIndex, FiConformalIndex, nBoxConformal, HitsinBoxConformal, OutputListHitsinTrack); if( Naux >= MINIMUMOUTERHITSPERTRACK && Naux > 0.7 * Nouter ) break; if( Naux >= MINIMUMOUTERHITSPERTRACK) { // 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]<=nmaxHitsInTrack) { for(j=0;j= .... } // end of if( Naux >= MINIMUMOUTERHITSPERTRACK) } // end of for(i=0; i< Nouter;i++) } // end of if( Nouter >= MINIMUMOUTERHITSPERTRACK) if( nHitsinTrack[nTracksFoundSoFar] < MINIMUMHITSPERTRACK || nHitsinTrack[nTracksFoundSoFar]>nmaxHitsInTrack) { 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;jnmaxHitsInTrack) 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. PndSttOrderingParallel( Ox[nTracksFoundSoFar], Oy[nTracksFoundSoFar], info, nHitsinTrack[nTracksFoundSoFar], &ListHitsinTrack[nTracksFoundSoFar][0], infoparal, 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 PndSttFindingParallelTrackAngularRange( 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. //-------------------------------------------------- macro for display for(j=0; j disegna in XY, 1 --> in SZ, altro --> in UV. ) { double x1,x2,y1,y2,L,R,RR; L=DIMENSIONSCITIL/2.; if(tipo==0){ // SciTil disegnate in XY. R = sqrt(posx*posx+posy*posy); x1 = posx + posy*L/R; x2 = posx - posy*L/R; y1 = posy - posx*L/R; y2 = posy + posx*L/R; } 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; R = sqrt(RR); x1 = posx + posy*L/R; x1 /= RR; x2 = posx - posy*L/R; x2 /= RR; y1 = posy - posx*L/R; y1 /= RR; y2 = posy + posx*L/R; 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); /* fprintf(MACRO,"TMarker* SciT%d = new TMarker(%f,%f,%d);\n", ScitilHit,posx,posy,30); fprintf(MACRO,"SciT%d->SetMarkerSize(1.5);\n",ScitilHit); fprintf(MACRO,"SciT%d->SetMarkerColor(1);\nSciT%d->Draw();\n" ,ScitilHit,ScitilHit); */ } //----------end of function PndSttTrackFinderReal::disegnaSciTilHit //----------begin of function PndSttTrackFinderReal::disegnaassiXY void PndSttTrackFinderReal::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 PndSttTrackFinderReal::disegnaassiXY ClassImp(PndSttTrackFinderReal)