// ------------------------------------------------------------------------- // ----- PndSttTrackFinderReal source file ----- // ----- by Gianluigi Boca ----- // ------------------------------------------------------------------------- #include "glpk.h" // Pnd includes #include "PndSttTrackFinderReal.h" #include "PndSttHit.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" //#include "GFKalman.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() { iplotta = false; istampa = 0; doMcComparison = false; // fVerbose = 1; Fimin=0.; Fimax=2.*PI; FI0min = 0.; FI0max = 2.*PI; 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; sprintf(fSttBranch,"STTHit"); } // ------------------------------------------------------------------------- // ----- Standard constructor ------------------------------------------ PndSttTrackFinderReal::PndSttTrackFinderReal(int verbose) { istampa = verbose; iplotta = false; doMcComparison = false; MINIMUMOUTERHITSPERTRACK=5; Fimin=0.; Fimax=2.*PI; FI0min = 0.; FI0max = 2.*PI; 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; sprintf(fSttBranch,"STTHit"); } // ------------------------------------------------------------------------- // ----- Second constructor ------------------------------------------ PndSttTrackFinderReal::PndSttTrackFinderReal(int istamp, bool iplott, bool imc) { istampa= istamp; iplotta = iplott; doMcComparison = imc; MINIMUMOUTERHITSPERTRACK=5; Fimin=0.; Fimax=2.*PI; FI0min = 0.; FI0max = 2.*PI; 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; sprintf(fSttBranch,"STTHit"); } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndSttTrackFinderReal::~PndSttTrackFinderReal() { } // ------------------------------------------------------------------------- //----------begin of function PndSttTrackFinderReal::WriteHistograms void PndSttTrackFinderReal::WriteHistograms(){ // TFile* file = FairRootManager::Instance()->GetOutFile(); TFile* file = FairRootManager::Instance()->GetOutFile(); file->cd(); file->mkdir("PndSttTrackFinderReal"); file->cd("PndSttTrackFinderReal"); hdist->Write(); hdistgoodlast->Write(); hdistbadlast->Write(); delete hdist; delete hdistgoodlast; delete hdistbadlast; } //----------end of function PndSttTrackFinderReal::WriteHistograms // ----- Public method Init -------------------------------------------- void PndSttTrackFinderReal::Init() { UShort_t i, j, k; Double_t r1, r2, A , tempRadiaConf[nRdivConformal]; MINIMUMHITSPERTRACK=3; // --------------------------- opening files for special purposes N_INTENDED=0; if(istampa >=1 ){ //---- fetch the n. of tracks MC that were intended to be generated HANDLE = fopen("n_intended_tracks.txt","r"); fscanf(HANDLE,"%d",&N_INTENDED); fclose(HANDLE); //---- apertura file con info su Found tracce su cui si fa Helix fit dopo HANDLE2 = fopen("info_da_PndTrackFinderReal.txt","w"); // ---- open filehandle per statistica sugli hits etc. HANDLE = fopen("statistiche.txt","w"); // --------------- 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; ntimes = 0; // Get and check FairRootManager FairRootManager* ioman = FairRootManager::Instance(); // get the MCTrack array fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack"); // ------------------------------------------------------------------- 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.); if (!ioman) { cout << "-E- PndSttTrackFinderReal::Init: " << "RootManager not instantised, return!" << endl; return; } // 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; r1=r2; } } // now take into account the zone of the skew straws, which is 'empty' as far as the parallel straws is concerned for(i = 1; i< nRdivConformal ; i++){ radiaConf[i] = tempRadiaConf[i]; } 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* trackCandArray, TClonesArray* helixHitArray) {} // CHECK da cancellare // ----------------- start of method DoFind Int_t PndSttTrackFinderReal::DoFind(TClonesArray* trackCandArray, TClonesArray *trackArray, TClonesArray* helixHitArray) { // list of hit numbers of those falling in this cell Short_t Charge[MAXTRACKSPEREVENT], Status[MAXTRACKSPEREVENT], daTrackFoundaTrackMC[MAXTRACKSPEREVENT], enne[MAXTRACKSPEREVENT][nmaxHits]; // daMCTrackaTrackFound[MAXMCTRACKS]; UShort_t emme, auxIndex[nmaxHits], OLDinfoparal[nmaxHits]; UShort_t istep,inclination_type, exitstatus; Double_t aaa, ddd, delta, deltabis, deltaZ, mindis, distanza, fi_hit, ap1, ap2, ap3, carica, cross1, cross2, cross3, 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]; //------------------------------------ modifiche Gianluigi, 9-7-08 IVOLTE++; if(istampa>=1) cout<<"\nda PndSttTrackFinderReal : evento (a partire da 0) n. "<GetEntriesFast(); // num. tracce/evento if(istampa>=1) cout<<"Gianluigi, evento n. "<= nmaxHits ) { cout<<"from PndSttTrackFinderReal : # Stt hits (|| + //) = "<= nmaxHits = " <=1 && (nMCTracks != N_INTENDED) ){ cout<<"Evento n. "< > hitMap; for(Int_t j=0; j= 1 ) cout<<"Gianluigi : da PndSttTrackFinderReal::DoFind : Nhits="<GetRefIndex(); if (ptIndex >= 0) { pMCpt = GetPointFromCollections(iHit); // <== FairMCPoint } // if (ptIndex < 0) continue; // fake or background hit // if (!pMCpt){ // cout<<"from PndSttTrackFinderReal : # MC points pointer missing, return!\n"; // continue; // } // 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 { for (Int_t 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]++; goto jumpout; } } 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]++; jumpout: ; } NSkewhits = 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 ="< 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 UShort_t iExclude, imc,jexp,mchit,exphit, nTracksFoundSoFar, NN, NNN, Nouter, Naux, Nbaux, nTotalHits[MAXTRACKSPEREVENT], iParHit, SeedParallelNumber, OLDnHitsinTrack, TemporarynSkewHitsinTrack, BigList[MAXTRACKSPEREVENT][nmaxHits], TemporarySkewList[nmaxHits][2], nHitsinTrack[MAXTRACKSPEREVENT], nSkewHitsinTrack[MAXTRACKSPEREVENT], tempore[nmaxHits], auxListHitsinTrack[nmaxHits], ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHits], ListHitsinTrackinWhichToSearch[nmaxHits], OLDListHitsinTrack[nmaxHits], OutputListHitsinTrack[nmaxHits], OutputList2HitsinTrack[nmaxHits], ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxHits], 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 // UShort_t HitsinBoxConformal[Nhits][nRdivConformal][nFidivConformal]; UShort_t HitsinBoxConformal[MAXHITSINCELL][nRdivConformal][nFidivConformal]; // first index -> radial divisions, 2nd index -> azimuthal divisions; Short_t flagStt; Int_t status; Int_t i, j,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 ; 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, rotationangle, rotationsin, rotationcos,ptotal, Fi_low_limit[MAXTRACKSPEREVENT], Fi_up_limit[MAXTRACKSPEREVENT], Fi_initial_helix_referenceframe[MAXTRACKSPEREVENT], Fi_final_helix_referenceframe[MAXTRACKSPEREVENT], m[MAXTRACKSPEREVENT], q[MAXTRACKSPEREVENT], ALFA[MAXTRACKSPEREVENT], BETA[MAXTRACKSPEREVENT], GAMMA[MAXTRACKSPEREVENT], R[MAXTRACKSPEREVENT], Ox[MAXTRACKSPEREVENT], Oy[MAXTRACKSPEREVENT], KAPPA[MAXTRACKSPEREVENT], FI0[MAXTRACKSPEREVENT], trajectory_vertex[3], infoparalConformal[nmaxHits][5], auxinfoparalConformal[nmaxHits][5], S[2*nmaxHits], 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]; for(i=0; i< Minclinations[0]; i++){ ExclusionList[ infoparal[i] ]= true ; ExclusionListbis[ infoparal[i] ]= true ; } for(i=0; i< NSkewhits; i++){ ExclusionListSkew[ infoskew[i] ]= true ; ExclusionListSkewbis[ infoskew[i] ]= true ; } //----------------------------------- exclusion of straws with multiple hits // first the parallel straws for(i=0; i< Nhits-1; i++){ if( info[i][5]==1.){ if( ExclusionList[ i ] ){ for(j=i+1; j< Nhits; j++){ if(ExclusionList[ j ] && info[j][5]==1. && fabs(info[i][0] - info[j][0])<1.e-20 && fabs(info[i][1] - info[j][1])<1.e-20 ) { ExclusionList[j]= false ; ExclusionListbis[j]= false ; } } // end of for(j=i+1; j< Nhits;; j++) } // end of if( ExclusionList[ i ] ) } else { // continuation of if( info[i][5]==1.) if( ExclusionListSkew[ i ] ){ for(j=i+1; j< Nhits; j++){ if(ExclusionListSkew[ j ] && info[j][5] != 1. && fabs(info[i][0] - info[j][0])<1.e-20 && fabs(info[i][1] - info[j][1])<1.e-20 ) { ExclusionListSkew[j]= false ; ExclusionListSkewbis[j]= false ; } } // end of for(j=i+1; j< Nhits;; j++) } // end of if( ExclusionList[ i ] ) } // // end of if( info[i][5]==1.) } // end of for(i=0; i< Nhits-1; i++) //----------------------------------- end of exclusion of straws with multiple hits trajectory_vertex[0]=trajectory_vertex[1]=trajectory_vertex[2]=0.; PndSttFromXYtoConformal(trajectory_vertex,info, Minclinations[0],infoparalConformal, &status); PndSttBoxConformalFilling( ExclusionList, infoparalConformal, Minclinations[0], nBoxConformal, HitsinBoxConformal, RConformalIndex, FiConformalIndex); // start the track finding procedure //----- loop over the parallel hits nTracksFoundSoFar=0; // # tracks found // begins the first iteration with more severe cuts on the # hits in track candidate for(iParHit=0; iParHit MAXTRACKSPEREVENT) continue; if( ! ExclusionList[ infoparal[iParHit] ] ) continue; nHitsinTrack[nTracksFoundSoFar] = PndSttFindTrackPatterninBoxConformal( 1, // distance in R cells allowed 2, // distance in Fi cells allowed Minclinations[0], iParHit, info, ExclusionList, RConformalIndex, FiConformalIndex, nBoxConformal, HitsinBoxConformal, &ListHitsinTrack[nTracksFoundSoFar][0] ); if(istampa>=3){ cout<<" \tSttTrackFinderReal : track found n. "<nmaxHitsInTrack) { continue; } //----------------------- // find among the ListHitsinTrack if there are at least a minimum # of hits belonging to the outer part // of the STT system for(j=0, Nouter =0; j= MINIMUMOUTERHITSPERTRACK) { for(i=0; i< Nouter;i++){ ListHitsinTrackinWhichToSearch[i]=ListHitsinTrack[nTracksFoundSoFar][i]; } for(i=0; i< Nouter;i++){ SeedParallelNumber = ListHitsinTrackinWhichToSearch[i]; Naux = PndSttFindTrackPatterninBoxConformalSpecial( 3, // NRCELLDISTANCE 1, // NFiCELLDISTANCE Minclinations[0], Nouter, SeedParallelNumber, ListHitsinTrackinWhichToSearch, info, ExclusionList, RConformalIndex, FiConformalIndex, nBoxConformal, HitsinBoxConformal, OutputListHitsinTrack); // cout<<" Naux = "<= MAXTRACKSPEREVENT ){ cout<<"from PndSttTrackFinderReal : # n. Tracks found so far = "<= MAXTRACKSPEREVENT ( = "<1.e10) { // necessary to load here the Sfinal vector anyway. for(j=0;j=0 ){ nSkewHitsinTrack[i] = NNN; for(j=0;j=0 ) nSkewHitsinTrack[i]=TemporarynSkewHitsinTrack; for(j=0;j=0 ) if( nHitsinTrack[i]+nSkewHitsinTrack[i] < MINIMUMHITSPERTRACK || nHitsinTrack[i]+nSkewHitsinTrack[i]>nmaxHitsInTrack) { keepit[i]=false; continue; } // --------- numbering according to the ORIGINAL hit number for(i1=0; i1< nSkewHitsinTrack[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<<"\tprima di Skew cleanup, IVOLTE = " <=2)cout<<"\thit skew (nativo) n. "<=2)cout<<"\tdopo di Skew cleanup, IVOLTE = " <=2)cout<<"\thit skew (nativo) n. "<0 && nTracksFoundSoFar > 0 && doMcComparison ) { AssociateFoundTrackstoMCbis( keepit, info, nTracksFoundSoFar, nHitsinTrack, ListHitsinTrack, nSkewHitsinTrack, ListSkewHitsinTrack, daTrackFoundaTrackMC ); for(jexp=0; jexp=1 ) fprintf(HANDLE, "\n Evento %d NTotaleTracceMC %d ------\n",IVOLTE, nMCTracks); 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], nSkewHitsinTrack[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; } fprintf(HANDLE, " TracciaMC %d ParHitsMC %d ParMecc %d ParMeccSpuri %d SkewHitsMC %d SkewMecc %d SkewMeccSpuri %d\n", i, nHitsInMCTrack[ii], nParalCommon[ii], nSpuriParinTrack[ii], nSkewHitsInMCTrack[ii], nSkewCommon[ii], nSpuriSkewinTrack[ii] ) ; fprintf(HANDLE, " e corrisponde a track found n. %d\n", ii ); fprintf(HANDLE, " AVENDO %d hits paralleli e %d hits skew non mecciati dalla corrisponde track found\n" , nMCParalAlone[ii],nMCSkewAlone[ii] ); HoughFi = atan2(Oy[ii],Ox[ii]); if(HoughFi<0.) HoughFi += 2.*PI; Double_t Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy ; 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 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 if( istampa>=2){ int NParghost=0, NParhitsghost=0,icc; 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); pTrckCand->setMcTrackId( daTrackFoundaTrackMC[i] ); for(j=0; j< nTotalHits[i]; j++){ pTrckCand->AddHit( FairRootManager::Instance()->GetBranchId(fSttBranch), (Int_t) BigList[i][j] , 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(1., 0., 0.), TVector3(0., 1., 0.) // ); 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<<"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,TypeConf,ALFA,BETA,GAMMA ); //----------------------------------- end macro for display if(doMcComparison) { WriteMacroParallelHitsGeneralConformalwithMC( keepit, Nhits, info, Nincl, Minclinations, inclination,nTracksFoundSoFar,TypeConf,ALFA,BETA,GAMMA ); for(i=0, ii=-1; i= 0){ WriteMacroSkewAssociatedHitswithMC( GoodSkewFit[i], KAPPA[i],FI0[i], HoughD, HoughFi, HoughR, info, Nincl,Minclinations,inclination, i, ii, nSkewHitsinTrack[i], ListSkewHitsinTrack, nSkewCommon[i], SkewCommonList, daTrackFoundaTrackMC[i], nMCSkewAlone, MCSkewAloneList ); } } } // 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 : "<= semilengthSkew1 ){ BAD1[i] = true ; } else { BAD1[i]= false ; } } //---------------- if ( BAD1[0] && BAD1[1] ) { *STATUS = -4; return Result; } //---------------- calculateintersections(Ox,Oy,R,C0x2,C0y2,C0z2,r2,vx2,vy2,vz2, STATUS,POINTS2); if(*STATUS < 0 ) { return Result;} //------------------------- for( i=0; i<2; i++){ j=3*i; distance2[i] = sqrt( (POINTS2[j]-C0x2)*(POINTS2[j]-C0x2) + (POINTS2[1+j]-C0y2)*(POINTS2[1+j]-C0y2) + (POINTS2[2+j]-C0z2)*(POINTS2[2+j]-C0z2) ); if( distance2[i] >= semilengthSkew2 ){ BAD2[i] = true ; } else { BAD2[i]= false ; } } //------------------------- if ( BAD2[0] && BAD2[1] ) { *STATUS = -5; return Result; } //--------------------------------------------------------------------------------------------------------------------- for(k1=0; k1<2;k1++){ // if( BAD1[k1] ) { cout<<"per k1 = "< 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 = r1*aaa/LL; if(istampa >= 3 && IVOLTE<= nmassimo){ cout<<"Lunghezza asse maggiore ellisse1 = "<Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); fprintf(MACRO,"TEllipse* E1 = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\n",POINTS1[i1+2],SSS1,Aellipsis1,Bellipsis1,rotation1); fprintf(MACRO,"TEllipse * E2 = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\n",POINTS2[i2+2],SSS2,Aellipsis2,Bellipsis2,rotation2); fprintf(MACRO,"TLine* L = new TLine(%f,%f,%f,%f);\n", xmin, Result.KAPPA[enne+1][i]*xmin+Result.FI0[enne+1][i], xmax, Result.KAPPA[enne+1][i]*xmax+Result.FI0[enne+1][i] ); fprintf(MACRO,"TGaxis *Assey = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n",xmin+(xmax-xmin)*0.1, ymin+(ymax-ymin)*0.1,xmin+(xmax-xmin)*0.1,ymax-(ymax-ymin)*0.1, ymin+(ymax-ymin)*0.1,ymax-(ymax-ymin)*0.1); fprintf(MACRO,"Assey->Draw();\n"); fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n",xmin+(xmax-xmin)*0.1, ymin+(ymax-ymin)*0.1,xmax-(xmax-xmin)*0.1,ymin+(ymax-ymin)*0.1, xmin+(xmax-xmin)*0.1,xmax-(xmax-xmin)*0.1); fprintf(MACRO,"Assex->Draw();\n"); if( ymin<2.*PI && ymax > 2.*PI) { fprintf(MACRO,"TLine* TWOPI = new TLine(%f,%f,%f,%f);\nTWOPI->SetLineColor(2);\nTWOPI->Draw();\n", xmin,2.*PI,xmax,2.*PI); } if( ymin<0. && ymax > 0.) { fprintf(MACRO,"TLine* BASE = new TLine(%f,%f,%f,%f);\nBASE->SetLineColor(2);\nBASE->Draw();\n", xmin,0.,xmax,0.); } if( ymin<-2.*PI && ymax > -2.*PI) { fprintf(MACRO,"TLine* NEG = new TLine(%f,%f,%f,%f);\nNEG->SetLineColor(2);\nNEG->Draw();\n", xmin,-2.*PI,xmax,-2.*PI); } fprintf(MACRO,"TLine* Teorica = new TLine(%f,%f,%f,%f);\n", xmin, xmin/R + 3.*PI/2., xmax, xmax/R + 3.*PI/2. ); fprintf(MACRO,"Teorica->SetLineColor(4);\n"); fprintf(MACRO,"E1->Draw();\nE2->Draw();\nL->Draw();\nTeorica->Draw();\n}\n"); fclose(MACRO); } // end of if(ntimes < 3) } // end of for(i=0; i 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( UShort_t vec1[2], // input int nListElements, UShort_t List[][2], // input int & nClusterElementsFound, UShort_t ClusterElementsFound[][2], // output int & nRemainingElements, UShort_t RemainingElements[][2] // output ) { int i; UShort_t vec2[2]; nClusterElementsFound=0; nRemainingElements=0; for(i=0; i xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } } if( xmin > 0. ) xmin = 0.; if( xmax < 0.) xmax = 0.; if( ymin > 0. ) ymin = 0.; if( ymax < 0.) ymax = 0.; deltax = xmax-xmin; deltay = ymax - ymin; if( deltax > deltay) { ymin -= 0.5*(deltax-deltay); ymax = ymin+ deltax; delta = deltax; } else { xmin -= 0.5*(deltay-deltax); xmax = xmin+ deltay; delta= deltay; } xmax = xmax + delta*0.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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( 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 = 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.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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( 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); dopo: ; //------------------------------------------------------------------------------------------------------------ // 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,"MacroGeneralParallelHitsConformeEvent%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; 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.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); // 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 for(i=0; iSetLineColor(4);\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->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( 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); } } //------------------ 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( 0.25*ALFA[i]*ALFA[i]+0.25*BETA[i]*BETA[i]-GAMMA[i])/fabs(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); } } } } // ------------------------------ 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], UShort_t nTracksFoundSoFar, bool *TypeConf, Double_t *ALFA, Double_t *BETA, Double_t *GAMMA ) { 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]; //---------- 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,"MacroGeneralParallelHitsConformewithMCEvent%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; 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.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); // 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 for(i=0; iSetLineColor(4);\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->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( 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); } } //------------------ 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.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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, UShort_t Nhits, UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHits], // UShort_t infoparal[nmaxHits], Double_t info[][7], Int_t Nincl, Int_t Minclinations[], Double_t inclination[][3], UShort_t imaxima, Int_t sequencial ) { 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 "< 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.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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"); 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( 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); } 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, UShort_t Nhits, UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHits], Double_t info[][7], UShort_t ifoundtrack, Int_t sequentialNTrack, UShort_t nParalCommon[MAXTRACKSPEREVENT], UShort_t ParalCommonList[MAXTRACKSPEREVENT][nmaxHits], UShort_t nSpuriParinTrack[MAXTRACKSPEREVENT], UShort_t ParSpuriList[MAXTRACKSPEREVENT][nmaxHits], UShort_t nMCParalAlone[MAXTRACKSPEREVENT], UShort_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; // Ox = (D+R)*cos(Fi); // Oy = (D+R)*sin(Fi); //cout<<"da MacroTrackparalleletc. Ox, Oy "< 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.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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"); 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( 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); } //--------------------------- 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) goto carica0 ; 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 ( pMC ) } // end of if( TrackFoundaTrackMC > -1){ //----------- fine parte del MC carica0: ; 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, UShort_t nSkewHitsinTrack, UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxHits] // UShort_t nSkewCommon, // UShort_t SkewCommonList[MAXTRACKSPEREVENT][nmaxHits] ) { 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; for( iii=0; iii< nSkewHitsinTrack; iii++) { i = infoskew[ ListSkewHitsinTrack[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; // 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); // ------ se lo hit e' spurio marcalo in rosso /* for( i1=0; i1SetLineColor(2);\n",index); fuori: ; */ index++; } // end of for( ii=0; ii<2; ii++) } // end of for( i=1; i< Nhits; i++) if(index==0) goto nohits ; if( zmax < zmin ) goto nohits ; if( Smax < Smin ) goto nohits; 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(!goodskewfit)goto niente ; if ( fabs(KAPPA) > 1.e-10 && fabs(KAPPA) < 1.e10 ) { fmin = KAPPA*zmin + FI0; fmax = KAPPA*zmax + FI0; } else { 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;++) niente: ; nohits: ; 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, UShort_t nSkewHitsinTrack, UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxHits], UShort_t nSkewCommon, UShort_t SkewCommonList[MAXTRACKSPEREVENT][nmaxHits], UShort_t daTrackFoundaTrackMC, UShort_t nMCSkewAlone[MAXTRACKSPEREVENT], UShort_t MCSkewAloneList[MAXTRACKSPEREVENT][nmaxHits] ) { 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]; 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; //---------------------------- for( iii=0; iii< nSkewHitsinTrack; iii++) { i = infoskew[ ListSkewHitsinTrack[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; // 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 for( i1=0; i1SetLineColor(2);\n",index); fuori: ; index++; } // end of for( ii=0; ii<2; ii++) } // end of for( iii=0; iii< nSkewHitsinTrack; iii++) //------------------------------- //------ aggiungo in blu eventuali punti della traccia MC che sono non mecciati for( iii=0; iii< nMCSkewAlone[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; // 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); // ------ marca lo hit in blu fprintf(MACRO,"E%d->SetLineColor(4);\n",index); 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==0) goto nohits ; if( zmax < zmin ) goto nohits ; if( Smax < Smin ) goto nohits; 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(!goodskewfit) goto nulla ; if ( fabs(KAPPA) > 1.e-10 && fabs(KAPPA) < 1.e10) { fmin = KAPPA*zmin + FI0; fmax = KAPPA*zmax + FI0; } else { cout<<"PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithMC, 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;++) nulla: ; //--------------------------------------------- qui ci aggiungo le traccie MC // for(imc=0; imcAt(imc); if ( ! pMC ) goto nohits ; 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) goto carica0 ; 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(); FI0 = fmod(Fifi+ PI, 2.*PI); if ( fabs(KAPPA) > 1.e-10 && fabs(KAPPA) < 1.e10) { fmin = KAPPA*zmin + FI0; fmax = KAPPA*zmax + FI0; } else { cout<<"PndSttTrackFinderReal::WriteMacroSkewAssociatedHits, this track found by PR not plotted" <<"\n\t because KAPPA = "<=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;++) nix: ; carica0: ; nohits: ; 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["< 0 && i < nHitsinTrack) { nRcell = RConformalIndex[ infoparal[ OutputListHitsinTrack[i] ] ]; nFicell = FiConformalIndex[ infoparal[ OutputListHitsinTrack[i] ] ]; //cout<<"From PatterninBoxConformal HIT collezionatore ora ha nRcell e nFicell "<= nRdivConformalEffective ) { nRmax = nRdivConformalEffective-1; } else { nRmax = nRcell + NRCELLDISTANCE; } //cout <<" di conseguenza ora nRmin e nRmax sono : "<= nFidivConformal) { iFi = iFi2 - nFidivConformal; } else { iFi = iFi2; } for (j = 0; j< nBoxConformal[iR][iFi]; j++){ if( ExclusionList[ infoparal[ HitsinBoxConformal[j][iR][iFi] ] ] && TemporaryExclusionList[ infoparal[ HitsinBoxConformal[j][iR][iFi] ] ]) { OutputListHitsinTrack[nHitsinTrack]=HitsinBoxConformal[j][iR][iFi] ; // hit number in the PARALLEL straws scheme nHitsinTrack++; TemporaryExclusionList[ infoparal[ HitsinBoxConformal[j][iR][iFi] ] ]= false; nRemainingHits--; } } } } //---------------- i++; } // end while ( nRemainingHits > 0 && i < nHitsinTrack) return nHitsinTrack; } //----------end of function PndSttTrackFinderReal::PndSttFindTrackPatterninBoxConformalSpecial //----------begin of function PndSttTrackFinderReal::PndSttFindTrackStrictCollection Short_t PndSttTrackFinderReal::PndSttFindTrackStrictCollection( UShort_t NFiCELLDISTANCE, UShort_t iSeed, // seed track (parallel notation) as fa as the Fi angle is concerned UShort_t NParallelToSearch, // n. of hits to search in ListHitsinTrackinWhichToSearch UShort_t *ListHitsinTrackinWhichToSearch, bool ExclusionList[nmaxHits], UShort_t FiConformalIndex[nmaxHits], UShort_t *OutputListHitsinTrack ) { UShort_t i, j, iR, iFi, iFiseed, nHitsinTrack; Double_t auxRvalues[nmaxHits]; if(istampa>=3 && IVOLTE <= nmassimo) { cout<<"Da strictcollection, n. elementi delle search list "<=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( ExclusionList[ infoparal[ ListHitsinTrackinWhichToSearch[i] ] ] ) } // end of for(i=0; i=3 && IVOLTE <= nmassimo) { cout<<"Da strictcollection, n. elementi aggiunti "< 1.e-10) { Delta[i] = 3.*auxinfoparalConformal[ i ][2]; // 3 times the Drift Radius } else { Delta[i] = 3.*StrawRadius; } } // sprintf(nome,"GeneralParallelHitsConformeTraccia%dEvent%d.mcs",(nTracksFoundSoFar), IVOLTE); // MACRO = fopen(nome,"w"); // fprintf(MACRO,"NAME FIT\n"); //----------------- write the ROWS section /* fprintf(MACRO,"ROWS\n"); fprintf(MACRO," N OBJECT\n"); for(i=0, ii=0 ; i< NpointsInFit ; i++) { ii++; fprintf(MACRO," L Ap%d\n L Bp%d\n L Cp%d\n L Dp%d\n",i,i,i,i); fprintf(MACRO," L Am%d\n L Bm%d\n L Cm%d\n L Dm%d\n G LAMBDA%d\n",i,i,i,i,i); } */ //-------- // nameRows[0]="OBJECT"; sprintf(&(auxnameRows[0][0]),"OBJECT",i); nameRows[0]=&auxnameRows[0][0]; typeRows[0]=GLP_FR; for(i=0 ; i< NpointsInFit ; i++) { ii=9*i; typeRows[1+ii]=GLP_UP;typeRows[2+ii]=GLP_UP;typeRows[3+ii]=GLP_UP;typeRows[4+ii]=GLP_UP; typeRows[5+ii]=GLP_UP;typeRows[6+ii]=GLP_UP;typeRows[7+ii]=GLP_UP;typeRows[8+ii]=GLP_UP; typeRows[9+ii]=GLP_LO; sprintf(&(auxnameRows[1+ii][0]),"Ap%d",i); nameRows[1+ii]=&auxnameRows[1+ii][0]; sprintf(&(auxnameRows[2+ii][0]),"Bp%d",i); nameRows[2+ii]=&auxnameRows[2+ii][0]; sprintf(&(auxnameRows[3+ii][0]),"Cp%d",i); nameRows[3+ii]=&auxnameRows[3+ii][0]; sprintf(&(auxnameRows[4+ii][0]),"Dp%d",i); nameRows[4+ii]=&auxnameRows[4+ii][0]; sprintf(&(auxnameRows[5+ii][0]),"Am%d",i); nameRows[5+ii]=&auxnameRows[5+ii][0]; sprintf(&(auxnameRows[6+ii][0]),"Bm%d",i); nameRows[6+ii]=&auxnameRows[6+ii][0]; sprintf(&(auxnameRows[7+ii][0]),"Cm%d",i); nameRows[7+ii]=&auxnameRows[7+ii][0]; sprintf(&(auxnameRows[8+ii][0]),"Dm%d",i); nameRows[8+ii]=&auxnameRows[8+ii][0]; sprintf(&(auxnameRows[9+ii][0]),"LAMBDA%d",i); nameRows[9+ii]=&auxnameRows[9+ii][0]; } //----------------- write the COLUMNS section // fprintf(MACRO,"COLUMNS\n"); // Column variable m1 for(i=0, ii=0 ; i< NpointsInFit ; i++) { ii++; // fprintf(MACRO," m1 Ap%d %g Am%d %g\n m1 Bp%d %g Bm%d %g\n", // i,Ox[i],i,Ox[i],i,-Ox[i],i,-Ox[i]); Coefficients[i*4]= Ox[i]; Coefficients[i*4+1]= Ox[i]; Coefficients[i*4+2]= -Ox[i]; Coefficients[i*4+3]= -Ox[i]; } // Column variable m2 for(i=0; i< NpointsInFit ; i++) { // fprintf(MACRO," m2 Ap%d %g Am%d %g\n m2 Bp%d %g Bm%d %g\n", // i,-Ox[i],i,-Ox[i],i,Ox[i],i,Ox[i]); Coefficients[NStructRows+i*4]= -Ox[i]; Coefficients[NStructRows+i*4+1]= -Ox[i]; Coefficients[NStructRows+i*4+2]= Ox[i]; Coefficients[NStructRows+i*4+3]= Ox[i]; } // Column variable q1 for(i=0 ; i< NpointsInFit ; i++) { // fprintf(MACRO," q1 Ap%d 1. Am%d 1.\n q1 Bp%d -1. Bm%d -1.\n", // i,i,i,i); Coefficients[2*NStructRows+i*4]= 1.; Coefficients[2*NStructRows+i*4+1]= 1.; Coefficients[2*NStructRows+i*4+2]= -1.; Coefficients[2*NStructRows+i*4+3]= -1.; } // Column variable q2 for(i=0 ; i< NpointsInFit ; i++) { // fprintf(MACRO," q2 Ap%d -1. Am%d -1.\n q2 Bp%d 1. Bm%d 1.\n", // i,i,i,i); Coefficients[3*NStructRows+i*4]= -1.; Coefficients[3*NStructRows+i*4+1]= -1.; Coefficients[3*NStructRows+i*4+2]= 1.; Coefficients[3*NStructRows+i*4+3]= 1.; } // Column variable lambdap(i) for(i=0 ; i< NpointsInFit ; i++) { // fprintf(MACRO," lamp%d Ap%d %g Bp%d %g\n lamp%d Cp%d %g Dp%d %g\n lamp%d LAMBDA%d 1.\n", // i,i,-M,i,-M, i , i,-M, i, M, i,i); Coefficients[(4+i)*NStructRows+0]= -M; Coefficients[(4+i)*NStructRows+1]= -M; Coefficients[(4+i)*NStructRows+2]= -M; Coefficients[(4+i)*NStructRows+3]= M; Coefficients[(4+i)*NStructRows+4]= 1.; } // Column variable lambdam(i) for(i=0 ; i< NpointsInFit ; i++) { // fprintf(MACRO," lamm%d Am%d %g Bm%d %g\n lamm%d Cm%d %g Dm%d %g\n lamm%d LAMBDA%d 1.\n", // i,i,-M,i,-M, i,i , -M, i, M, i,i); Coefficients[(4+i+NpointsInFit)*NStructRows+0]= -M; Coefficients[(4+i+NpointsInFit)*NStructRows+1]= -M; Coefficients[(4+i+NpointsInFit)*NStructRows+2]= -M; Coefficients[(4+i+NpointsInFit)*NStructRows+3]= M; Coefficients[(4+i+NpointsInFit)*NStructRows+4]= 1.; } // Column variable sigmap(i) for(i=0; i< NpointsInFit ; i++) { // fprintf(MACRO," sigmap%d OBJECT %g Ap%d -1.\n sigmap%d Bp%d -1. Cp%d 1.\n sigmap%d Dp%d -1.\n", // i,1./Delta[i],i,i,i,i,i,i); Coefficients[(4+i+2*NpointsInFit)*NStructRows+0]= 1./Delta[i]; Coefficients[(4+i+2*NpointsInFit)*NStructRows+1]= -1.; Coefficients[(4+i+2*NpointsInFit)*NStructRows+2]= -1.; Coefficients[(4+i+2*NpointsInFit)*NStructRows+3]= 1.; Coefficients[(4+i+2*NpointsInFit)*NStructRows+4]= -1.; } // Column variable sigmam(i) for(i=0 ; i< NpointsInFit ; i++) { // fprintf(MACRO," sigmam%d OBJECT %g Am%d -1.\n sigmam%d Bm%d -1. Cm%d 1.\n sigmam%d Dm%d -1.\n", // i,1./Delta[i],i,i,i,i,i,i); Coefficients[(4+i+3*NpointsInFit)*NStructRows+0]= 1./Delta[i]; Coefficients[(4+i+3*NpointsInFit)*NStructRows+1]= -1.; Coefficients[(4+i+3*NpointsInFit)*NStructRows+2]= -1.; Coefficients[(4+i+3*NpointsInFit)*NStructRows+3]= 1.; Coefficients[(4+i+3*NpointsInFit)*NStructRows+4]= -1.; } // Column variable DUMMY /* for(i=0, ii=0 ; i< NpointsInFit ; i++) { fprintf(MACRO," DUMMY Ap%d 1. Am%d 1.\n DUMMY Bp%d 1. Bm%d 1.\n DUMMY Cp%d 1. Cm%d 1\n", i,i,i,i,i,i); fprintf(MACRO," DUMMY Dp%d 1. Dm%d 1.\n", i,i); } */ for(i=0 ; i< NStructRows ; i++) { Coefficients[(4+4*NpointsInFit)*NStructRows+i]= 1.; } //-------------------- sprintf(&auxStructVarName[0][0],"m1",i); StructVarName[0] = &auxStructVarName[0][0]; // StructVarName[0]="m1"; NRowsInWhichStructVarArePresent[0]= 4*NpointsInFit; sprintf(&auxStructVarName[1][0],"m2",i); StructVarName[1] = &auxStructVarName[1][0]; // StructVarName[1]="m2"; NRowsInWhichStructVarArePresent[1]= 4*NpointsInFit; sprintf(&auxStructVarName[2][0],"q1",i); StructVarName[2] = &auxStructVarName[2][0]; // StructVarName[2]="q1"; NRowsInWhichStructVarArePresent[2]= 4*NpointsInFit; sprintf(&auxStructVarName[3][0],"q2",i); StructVarName[3] = &auxStructVarName[3][0]; // StructVarName[3]="q2"; NRowsInWhichStructVarArePresent[3]= 4*NpointsInFit; for(i=0; i< NpointsInFit ; i++) { sprintf(&auxStructVarName[3+i+1][0],"lamp%d",i); StructVarName[4+i] = &auxStructVarName[4+i][0]; NRowsInWhichStructVarArePresent[4+i]= 5; sprintf(&auxStructVarName[4+NpointsInFit+i][0],"lamm%d",i); StructVarName[4+NpointsInFit+i] = &auxStructVarName[4+NpointsInFit+i][0]; NRowsInWhichStructVarArePresent[4+NpointsInFit+i]= 5; sprintf(&auxStructVarName[4+2*NpointsInFit+i][0],"sigmap%d",i); StructVarName[4+2*NpointsInFit+i] = &auxStructVarName[4+2*NpointsInFit+i][0]; NRowsInWhichStructVarArePresent[4+2*NpointsInFit+i]= 5; sprintf(&auxStructVarName[4+3*NpointsInFit+i][0],"sigmam%d",i); StructVarName[4+3*NpointsInFit+i] = &auxStructVarName[4+3*NpointsInFit+i][0]; NRowsInWhichStructVarArePresent[4+3*NpointsInFit+i]= 5; } sprintf(&auxStructVarName[4+4*NpointsInFit][0],"DUMMY",i); StructVarName[4+4*NpointsInFit] = &auxStructVarName[4+4*NpointsInFit][0]; // StructVarName[4+4*NpointsInFit]="DUMMY"; NRowsInWhichStructVarArePresent[4+4*NpointsInFit]= NStructRows; // for m1, m2, q1, q2 for(i=0; i< 4; i++){ for(ii=0; ii< NpointsInFit;ii++){ sprintf(&aux[i*NStructRows+ii*4][0],"Ap%d",ii); NameRowsInWhichStructVarArePresent[i*NStructRows+ii*4]=&aux[i*NStructRows+ii*4][0]; sprintf(&aux[i*NStructRows+ii*4+1][0],"Am%d",ii); NameRowsInWhichStructVarArePresent[i*NStructRows+ii*4+1]=&aux[i*NStructRows+ii*4+1][0]; sprintf(&aux[i*NStructRows+ii*4+2][0],"Bp%d",ii); NameRowsInWhichStructVarArePresent[i*NStructRows+ii*4+2]=&aux[i*NStructRows+ii*4+2][0]; sprintf(&aux[i*NStructRows+ii*4+3][0],"Bm%d",ii); NameRowsInWhichStructVarArePresent[i*NStructRows+ii*4+3]=&aux[i*NStructRows+ii*4+3][0]; } } // now for the lamp* variables for(i=0; i< NpointsInFit;i++){ sprintf(&aux[(i+4)*NStructRows+0][0],"Ap%d",i); NameRowsInWhichStructVarArePresent[(i+4)*NStructRows+0]= &aux[(i+4)*NStructRows+0][0]; sprintf(&aux[(i+4)*NStructRows+1][0],"Bp%d",i); NameRowsInWhichStructVarArePresent[(i+4)*NStructRows+1]= &aux[(i+4)*NStructRows+1][0]; sprintf(&aux[(i+4)*NStructRows+2][0],"Cp%d",i); NameRowsInWhichStructVarArePresent[(i+4)*NStructRows+2]= &aux[(i+4)*NStructRows+2][0]; sprintf(&aux[(i+4)*NStructRows+3][0],"Dp%d",i); NameRowsInWhichStructVarArePresent[(i+4)*NStructRows+3]= &aux[(i+4)*NStructRows+3][0]; sprintf(&aux[(i+4)*NStructRows+4][0],"LAMBDA%d",i); NameRowsInWhichStructVarArePresent[(i+4)*NStructRows+4]= &aux[(i+4)*NStructRows+4][0]; } // now for the lamm* variables for(i=0; i< NpointsInFit;i++){ sprintf(&aux[(i+4+NpointsInFit)*NStructRows+0][0],"Am%d",i); NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit)*NStructRows+0]= &aux[(i+4+NpointsInFit)*NStructRows+0][0]; sprintf(&aux[(i+4+NpointsInFit)*NStructRows+1][0],"Bm%d",i); NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit)*NStructRows+1]= &aux[(i+4+NpointsInFit)*NStructRows+1][0]; sprintf(&aux[(i+4+NpointsInFit)*NStructRows+2][0],"Cm%d",i); NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit)*NStructRows+2]= &aux[(i+4+NpointsInFit)*NStructRows+2][0]; sprintf(&aux[(i+4+NpointsInFit)*NStructRows+3][0],"Dm%d",i); NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit)*NStructRows+3]= &aux[(i+4+NpointsInFit)*NStructRows+3][0]; sprintf(&aux[(i+4+NpointsInFit)*NStructRows+4][0],"LAMBDA%d",i); NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit)*NStructRows+4]= &aux[(i+4+NpointsInFit)*NStructRows+4][0]; } // now for the sigmap* variables for(i=0; i< NpointsInFit;i++){ sprintf(&aux[(i+4+2*NpointsInFit)*NStructRows+0][0],"OBJECT",i); NameRowsInWhichStructVarArePresent[(i+4+2*NpointsInFit)*NStructRows+0]= &aux[(i+4+2*NpointsInFit)*NStructRows+0][0]; sprintf(&aux[(i+4+2*NpointsInFit)*NStructRows+1][0],"Ap%d",i); NameRowsInWhichStructVarArePresent[(i+4+2*NpointsInFit)*NStructRows+1]= &aux[(i+4+2*NpointsInFit)*NStructRows+1][0]; sprintf(&aux[(i+4+2*NpointsInFit)*NStructRows+2][0],"Bp%d",i); NameRowsInWhichStructVarArePresent[(i+4+2*NpointsInFit)*NStructRows+2]= &aux[(i+4+2*NpointsInFit)*NStructRows+2][0]; sprintf(&aux[(i+4+2*NpointsInFit)*NStructRows+3][0],"Cp%d",i); NameRowsInWhichStructVarArePresent[(i+4+2*NpointsInFit)*NStructRows+3]= &aux[(i+4+2*NpointsInFit)*NStructRows+3][0]; sprintf(&aux[(i+4+2*NpointsInFit)*NStructRows+4][0],"Dp%d",i); NameRowsInWhichStructVarArePresent[(i+4+2*NpointsInFit)*NStructRows+4]= &aux[(i+4+2*NpointsInFit)*NStructRows+4][0]; } // now for the sigmam* variables for(i=0; i< NpointsInFit;i++){ sprintf(&aux[(i+4+3*NpointsInFit)*NStructRows+0][0],"OBJECT",i); NameRowsInWhichStructVarArePresent[(i+4+3*NpointsInFit)*NStructRows+0]= &aux[(i+4+3*NpointsInFit)*NStructRows+0][0]; sprintf(&aux[(i+4+3*NpointsInFit)*NStructRows+1][0],"Am%d",i); NameRowsInWhichStructVarArePresent[(i+4+3*NpointsInFit)*NStructRows+1]= &aux[(i+4+3*NpointsInFit)*NStructRows+1][0]; sprintf(&aux[(i+4+3*NpointsInFit)*NStructRows+2][0],"Bm%d",i); NameRowsInWhichStructVarArePresent[(i+4+3*NpointsInFit)*NStructRows+2]= &aux[(i+4+3*NpointsInFit)*NStructRows+2][0]; sprintf(&aux[(i+4+3*NpointsInFit)*NStructRows+3][0],"Cm%d",i); NameRowsInWhichStructVarArePresent[(i+4+3*NpointsInFit)*NStructRows+3]= &aux[(i+4+3*NpointsInFit)*NStructRows+3][0]; sprintf(&aux[(i+4+3*NpointsInFit)*NStructRows+4][0],"Dm%d",i); NameRowsInWhichStructVarArePresent[(i+4+3*NpointsInFit)*NStructRows+4]= &aux[(i+4+3*NpointsInFit)*NStructRows+4][0]; } // now for the DUMMY variable for(i=0; i< NpointsInFit;i++){ sprintf(&aux[(4+4*NpointsInFit)*NStructRows +8*i][0],"Ap%d",i); NameRowsInWhichStructVarArePresent[(4+4*NpointsInFit)*NStructRows+i*8 ]= &aux[(4+4*NpointsInFit)*NStructRows +8*i][0]; sprintf(&aux[(4+4*NpointsInFit)*NStructRows+1+8*i][0],"Am%d",i); NameRowsInWhichStructVarArePresent[(4+4*NpointsInFit)*NStructRows+1+8*i]= &aux[(4+4*NpointsInFit)*NStructRows+1+8*i][0]; sprintf(&aux[(4+4*NpointsInFit)*NStructRows+2+8*i][0],"Bp%d",i); NameRowsInWhichStructVarArePresent[(4+4*NpointsInFit)*NStructRows+2+8*i]= &aux[(4+4*NpointsInFit)*NStructRows+2+8*i][0]; sprintf(&aux[(4+4*NpointsInFit)*NStructRows+3+8*i][0],"Bm%d",i); NameRowsInWhichStructVarArePresent[(4+4*NpointsInFit)*NStructRows+3+8*i]= &aux[(4+4*NpointsInFit)*NStructRows+3+8*i][0]; sprintf(&aux[(4+4*NpointsInFit)*NStructRows+4+8*i][0],"Cp%d",i); NameRowsInWhichStructVarArePresent[(4+4*NpointsInFit)*NStructRows+4+8*i]= &aux[(4+4*NpointsInFit)*NStructRows+4+8*i][0]; sprintf(&aux[(4+4*NpointsInFit)*NStructRows+5+8*i][0],"Cm%d",i); NameRowsInWhichStructVarArePresent[(4+4*NpointsInFit)*NStructRows+5+8*i]= &aux[(4+4*NpointsInFit)*NStructRows+5+8*i][0]; sprintf(&aux[(4+4*NpointsInFit)*NStructRows+6+8*i][0],"Dp%d",i); NameRowsInWhichStructVarArePresent[(4+4*NpointsInFit)*NStructRows+6+8*i]= &aux[(4+4*NpointsInFit)*NStructRows+6+8*i][0]; sprintf(&aux[(4+4*NpointsInFit)*NStructRows+7+8*i][0],"Dm%d",i); NameRowsInWhichStructVarArePresent[(4+4*NpointsInFit)*NStructRows+7+8*i]= &aux[(4+4*NpointsInFit)*NStructRows+7+8*i][0]; } //----------------- write the RHS section // fprintf(MACRO,"RHS\n"); for(i=0 ; i< NpointsInFit ; i++) { // fprintf(MACRO," BOUND Ap%d %g Bp%d %g\n BOUND Cp%d %g Dp%d %g\n", // i, Oy[i]+auxinfoparalConformal[ i ][2]+2.*M,i, // -Oy[i]-auxinfoparalConformal[ i ][2]+2.*M,i, // Delta[i]+2.*M,i,M-Delta[i]+2.*M); ValueB[i*9] = Oy[i]+auxinfoparalConformal[ i ][2]+2.*M; ValueB[i*9+1]= -Oy[i]-auxinfoparalConformal[ i ][2]+2.*M; ValueB[i*9+2]= Delta[i]+2.*M; ValueB[i*9+3]= M-Delta[i]+2.*M; /* fprintf(MACRO," BOUND Am%d %g Bm%d %g\n BOUND Cm%d %g Dm%d %g\n", i, Oy[i]-auxinfoparalConformal[ i ][2]+2.*M,i, -Oy[i]+auxinfoparalConformal[ i ][2]+2.*M,i, Delta[i]+2.*M,i,M-Delta[i]+2.*M); fprintf(MACRO," BOUND LAMBDA%d 1.\n",i); */ ValueB[i*9+4]= Oy[i]-auxinfoparalConformal[ i ][2]+2.*M; ValueB[i*9+5]= -Oy[i]+auxinfoparalConformal[ i ][2]+2.*M; ValueB[i*9+6]= Delta[i]+2.*M; ValueB[i*9+7]= M-Delta[i]+2.*M; ValueB[i*9+8]= 1.; } //----------------- write the RANGES section // fprintf(MACRO,"RANGES\n"); for(i=0 ; i< NpointsInFit ; i++) { // fprintf(MACRO," RANGE LAMBDA%d 1.\n",i); //--- ValueRanges[i]=1.; sprintf(&auxNameRanges[i][0],"LAMBDA%d",i); NameRanges[i]=&auxNameRanges[i][0]; } //----------------- write the BOUNDS section // fprintf(MACRO,"BOUNDS\n"); for(i=0 ; i< NpointsInFit ; i++) { // fprintf(MACRO," BV Bounds lamp%d\n", i); sprintf(&auxTypeofBound[i][0],"BV"); TypeofBound[i]= &auxTypeofBound[i][0]; // TypeofBound[i]="BV"; sprintf(&auxBoundStructVarName[i][0],"lamp%d",i); BoundStructVarName[i]=&auxBoundStructVarName[i][0]; BoundValue[i]=0.; } for(i=0 ; i< NpointsInFit ; i++) { // fprintf(MACRO," BV Bounds lamm%d\n", i); sprintf(&auxTypeofBound[i+NpointsInFit][0],"BV"); TypeofBound[i+NpointsInFit]= &auxTypeofBound[i+NpointsInFit][0]; // TypeofBound[i+NpointsInFit]="BV"; sprintf(&auxBoundStructVarName[i+NpointsInFit][0],"lamm%d",i); BoundStructVarName[i+NpointsInFit]=&auxBoundStructVarName[i+NpointsInFit][0]; BoundValue[i+NpointsInFit]=0.; } // fprintf(MACRO," FX Bounds DUMMY %g\n",2.*M); sprintf(&auxTypeofBound[2*NpointsInFit][0],"FX"); TypeofBound[2*NpointsInFit]= &auxTypeofBound[2*NpointsInFit][0]; // TypeofBound[2*NpointsInFit]="FX"; sprintf(&auxTypeofBound[2*NpointsInFit][0],"FX"); TypeofBound[2*NpointsInFit]= &auxTypeofBound[2*NpointsInFit][0]; sprintf(&auxBoundStructVarName[2*NpointsInFit][0],"DUMMY"); BoundStructVarName[2*NpointsInFit]=&auxBoundStructVarName[2*NpointsInFit][0]; // BoundStructVarName[2*NpointsInFit]="DUMMY"; BoundValue[2*NpointsInFit]=2.; //----- // fprintf(MACRO,"ENDATA\n"); // fclose(MACRO); //----------------------- funzioni chiamate direttamente /* cout<<"n. punti nel fit "<=3 && IVOLTE <= nmassimo){ sprintf(stringa, "/home/boca/panda/glpk/glpk-4.39/examples/glpsol --min -o soluztrack%dEvent%dstep%d GeneralParallelHitsConformeTraccia%dEvent%d.mcs", nTracksFoundSoFar,IVOLTE,INTERO,nTracksFoundSoFar,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", nTracksFoundSoFar, IVOLTE,INTERO,nTracksFoundSoFar,IVOLTE,INTERO); } /* system(stringa); sprintf(stringa2,"grep -e ' m1 ' -e ' m2 ' -e ' q1 ' -e ' q2 ' soluztrack%dEvent%dstep%d | awk '{print $3}' > temp", nTracksFoundSoFar,IVOLTE,INTERO); system(stringa2); FILE * RESULT = fopen("temp","r"); fscanf(RESULT,"%f",&m1_result); fscanf(RESULT,"%f",&m2_result); fscanf(RESULT,"%f",&q1_result); fscanf(RESULT,"%f",&q2_result); fclose(RESULT); */ m1_result = final_values[0]; m2_result = final_values[1]; q1_result = final_values[2]; q2_result = final_values[3]; /* q_result = q1_result - q2_result; m_result = m1_result ; A = A1_result - A2_result; if( fabs( cose -sine*m_result) > 1.e-12) { *emme = ( m_result*cose + sine ) / (cose -sine*m_result); *qu = q_result/(cose -sine*m_result); Status=1; } else { *emme = m_result*cose + sine ; *qu = q_result; Status=99; } */ //------------------------ transformation of the result in terms of ALFA, BETA, GAMMA *qu = q1_result - q2_result; // *qu = q1_result; // *emme = m1_result ; *emme = m1_result-m2_result ; if(istampa>=3 && IVOLTE <= nmassimo) { cout<<"Stampa dal fitter, risultato : qu, emme non ri-ruotati "<<*qu<<", "<<*emme< 1.e-10) { // trajectory is a circle in XY space ALFA[nTracksFoundSoFar] = *emme/(*qu); BETA[nTracksFoundSoFar] = -1./(*qu); TypeConf[nTracksFoundSoFar]=true; // now take into account the rotation and correct; the only affected quantities are ALFA and BETA alfetta = ALFA[nTracksFoundSoFar]; ALFA[nTracksFoundSoFar] = ALFA[nTracksFoundSoFar]*cose - BETA[nTracksFoundSoFar]*sine; BETA[nTracksFoundSoFar] = alfetta*sine + BETA[nTracksFoundSoFar]*cose; } else if(fabs(*emme)> 1.e-10) { // trajectory is a straight line in XY space of equation y= m*x // the rotation first angle = atan(*emme) + rotationangle; if( fabs(cos(angle)) > 1.e-10 ) { ALFA[nTracksFoundSoFar] = 999999.; BETA[nTracksFoundSoFar] = -ALFA[nTracksFoundSoFar]/tan(angle); } else { // in this case the equation is y = 0. ALFA[nTracksFoundSoFar] = 999999.; BETA[nTracksFoundSoFar] = 0.; TypeConf[nTracksFoundSoFar]=false; } } else { // in this case also the equation in XY plane is y = 0. ALFA[nTracksFoundSoFar] = 999999.; BETA[nTracksFoundSoFar] = 0.; TypeConf[nTracksFoundSoFar]=false; } // now take into account the displacement and correct GAMMA[nTracksFoundSoFar] += (trajectory_vertex[0]*trajectory_vertex[0]+ trajectory_vertex[1]*trajectory_vertex[1] -ALFA[nTracksFoundSoFar]*trajectory_vertex[0]-BETA[nTracksFoundSoFar]*trajectory_vertex[1]); ALFA[nTracksFoundSoFar] -= 2.*trajectory_vertex[0]; BETA[nTracksFoundSoFar] -= 2.*trajectory_vertex[1]; //------------------------ end of transformation of the result in terms of ALFA, BETA, GAMMA //-------- end of taking into account the traslation that was performed and undoing that // taking into account the rotation that was performed and calculate emme and qu in the normal conformal plane if(fabs(cose-*emme*sine)> 1.e-10) { *qu=*qu/(cose-*emme*sine); *emme=(*emme*cose+sine)/(cose-*emme*sine); return 1; } else { // in this case the equation is 0 = x+*qu . if(fabs(sine+*emme*cose) < 1.e-10) { cout<<" From PndSttFitHelixCylinder, situation impossible in principle! Returning -1" < Delta[i]) Delta[i] =aux; } sprintf(nome,"GeneralParallelHitsConformeTraccia%dEvent%d.mcs",(nTracksFoundSoFar), IVOLTE); MACRO = fopen(nome,"w"); fprintf(MACRO,"NAME FIT\n"); //----------------- write the ROWS section fprintf(MACRO,"ROWS\n"); fprintf(MACRO," N OBJECT\n"); for(i=0, ii=0 ; i< nHitsinTrack && ii=3 && IVOLTE <= nmassimo){ sprintf(stringa, "/home/boca/panda/glpk/glpk-4.39/examples/glpsol --min -o soluztrack%dEvent%dstep%d GeneralParallelHitsConformeTraccia%dEvent%d.mcs", nTracksFoundSoFar,IVOLTE,INTERO,nTracksFoundSoFar,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", nTracksFoundSoFar, IVOLTE,INTERO,nTracksFoundSoFar,IVOLTE,INTERO); } system(stringa); // sprintf(stringa2,"grep -e ' m1 ' -e ' q1 ' -e ' q2 ' -e ' A1 ' -e ' A2 ' soluztrack%dEvent%dstep%d | awk '{print $3}' > temp", // (nTracksFoundSoFar),IVOLTE,INTERO); sprintf(stringa2,"grep -e ' m1 ' -e ' q1 ' -e ' q2 ' soluztrack%dEvent%dstep%d | awk '{print $3}' > temp", nTracksFoundSoFar,IVOLTE,INTERO); // sprintf(stringa2,"grep -e ' m1 ' -e ' q1 ' soluztrack%dEvent%dstep%d | awk '{print $3}' > temp", // nTracksFoundSoFar,IVOLTE,INTERO); system(stringa2); FILE * RESULT = fopen("temp","r"); // fscanf(RESULT,"%f",&A1_result); // fscanf(RESULT,"%f",&A2_result); fscanf(RESULT,"%f",&m1_result); // fscanf(RESULT,"%f",&m2_result); fscanf(RESULT,"%f",&q1_result); fscanf(RESULT,"%f",&q2_result); fclose(RESULT); // system ("rm -f temp"); // system ("rm -f soluz* "); // system("rm -f GeneralParallelHitsConformeTraccia*"); /* q_result = q1_result - q2_result; m_result = m1_result ; A = A1_result - A2_result; if( fabs( cose -sine*m_result) > 1.e-12) { emme = ( m_result*cose + sine ) / (cose -sine*m_result); qu = q_result/(cose -sine*m_result); Status=1; } else { emme = m_result*cose + sine ; qu = q_result; Status=99; } */ //------------------------ transformation of the result in terms of ALFA, BETA, GAMMA qu = q1_result - q2_result; // qu = q1_result; emme = m1_result ; //if(istampa) cout<<"Stampa dal fitter, prima della correzione displacement qu, emme "<=3 && IVOLTE <= nmassimo) { cout<<"Stampa dal fitter, risultato : qu, emme non ri-ruotati "< 1.e-10) { // trajectory is a circle in XY space ALFA[nTracksFoundSoFar] = emme/qu; BETA[nTracksFoundSoFar] = -1./qu; TypeConf[nTracksFoundSoFar]=true; // now take into account the rotation and correct; the only affected quantities are ALFA and BETA alfetta = ALFA[nTracksFoundSoFar]; ALFA[nTracksFoundSoFar] = ALFA[nTracksFoundSoFar]*cose - BETA[nTracksFoundSoFar]*sine; BETA[nTracksFoundSoFar] = alfetta*sine + BETA[nTracksFoundSoFar]*cose; } else if(fabs(emme)> 1.e-10) { // trajectory is a straight line in XY space of equatio y= m*x // the rotation first angle = atan(emme) + rotationangle; if( fabs(cos(angle)) > 1.e-10 ) { ALFA[nTracksFoundSoFar] = 999999.; BETA[nTracksFoundSoFar] = -ALFA[nTracksFoundSoFar]/tan(angle); } else { // in this case the equation is y = 0. ALFA[nTracksFoundSoFar] = 999999.; BETA[nTracksFoundSoFar] = 0.; TypeConf[nTracksFoundSoFar]=false; } } else { // in this case also the equation in XY plane is y = 0. ALFA[nTracksFoundSoFar] = 999999.; BETA[nTracksFoundSoFar] = 0.; TypeConf[nTracksFoundSoFar]=false; } // now take into account the displacement and correct GAMMA[nTracksFoundSoFar] += (trajectory_vertex[0]*trajectory_vertex[0]+ trajectory_vertex[1]*trajectory_vertex[1] -ALFA[nTracksFoundSoFar]*trajectory_vertex[0]-BETA[nTracksFoundSoFar]*trajectory_vertex[1]); ALFA[nTracksFoundSoFar] -= 2.*trajectory_vertex[0]; BETA[nTracksFoundSoFar] -= 2.*trajectory_vertex[1]; //------------------------ end of transformation of the result in terms of ALFA, BETA, GAMMA // taking into account the traslation that was performed and undoing that /* if( Status == 99) { // in this case in the conformal space equation is m*X + q = 0 GAMMA[nTracksFoundSoFar] = trajectory_vertex[0]*trajectory_vertex[0]+trajectory_vertex[1]*trajectory_vertex[1] -trajectory_vertex[0]*emme/qu; TypeConf[nTracksFoundSoFar]=true; ALFA[nTracksFoundSoFar]= emme/qu - 2.*trajectory_vertex[0] ; BETA[nTracksFoundSoFar]= - 2.*trajectory_vertex[1]; } else if ( fabs(qu) > 1.e-10 ) { GAMMA[nTracksFoundSoFar]= trajectory_vertex[0]*trajectory_vertex[0]+trajectory_vertex[1]*trajectory_vertex[1] -trajectory_vertex[0]*emme/qu + trajectory_vertex[1]/qu; TypeConf[nTracksFoundSoFar]=true; ALFA[nTracksFoundSoFar]= ( emme/qu - 2.*trajectory_vertex[0] ); BETA[nTracksFoundSoFar]= (-1./qu - 2.*trajectory_vertex[1]); } else { // this is the case when q=0 GAMMA[nTracksFoundSoFar] = -trajectory_vertex[0]*emme + trajectory_vertex[1]; TypeConf[nTracksFoundSoFar]=false; ALFA[nTracksFoundSoFar]= emme; BETA[nTracksFoundSoFar]= -1.; GAMMA[nTracksFoundSoFar]= 0.; } */ //-------- end of taking into account the traslation that was performed and undoing that return 1; } //----------end of function PndSttTrackFinderReal::PndSttFitHelixCylinder2 //----------begin of function PndSttTrackFinderReal::PndSttTrkAssociatedParallelHitsToHelixQuater UShort_t PndSttTrackFinderReal::PndSttTrkAssociatedParallelHitsToHelixQuater( bool ExclusionList[nmaxHits], Double_t m, Double_t q, Short_t Status, UShort_t nHitsinTrack, UShort_t *ListHitsinTrack, Int_t NhitsParallel, Double_t Ox, Double_t Oy, Double_t R, Double_t info[][7], Double_t infoparalConformal[][5], UShort_t *RConformalIndex, UShort_t *FiConformalIndex, UShort_t nBoxConformal[nRdivConformal][nFidivConformal], UShort_t HitsinBoxConformal[][nRdivConformal][nFidivConformal], UShort_t *auxListHitsinTrack ) { bool passamin,passamax, Unselected[nmaxHits]; Short_t i, i2, j, k, l, l2, l3, itemp, kstart, kend, iFi0,FFimin, FFimax; UShort_t Nextra=8, nFi, Fi, nR, nAssociatedHits, nHit_original; Double_t maxFi, minFi, dist, xx, yy, aaa, angle, r, erre1, erre2, Rin, Rout, Fi0, ddd, fi1, fi2, dx, dy, distance, NTIMES=1.5; // number of Straw radia allowed in association nAssociatedHits=0; for(i=0; i 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 && Fimin < nRdivConformal/4 ) { iFi0 = (Short_t) (nRdivConformal/4 ); } else if ( FFimax > 3*nRdivConformal/4 && Fimin < 3*nRdivConformal/4 ){ iFi0 = (Short_t) (3*nRdivConformal/4 ); } else { cout <<"From PndSttTrackFinderReal::PndSttTrkAssociatedParallelHitsToHelixQuater :" <<" inconsistency, 0 associated hits to this track candidate\n"; return 0; } for(itemp=iFi0-5; itemp<=iFi0+5;itemp++){ i=itemp; if( i< 0 ) { i += nFidivConformal; } else if (i>=nFidivConformal){ i -= nFidivConformal*( i/nFidivConformal ); // i -= nFidivConformal; } for(l=0; lNTIMES*StrawRadius ) continue; //------------------- xx=infoparalConformal[ HitsinBoxConformal[k][l][i] ][0]; dist = fabs( xx ); if( 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< NSkewhits; iii++) return NAssociated; } //----------end of function PndSttTrackFinderReal::AssociateSkewHitsToXYTrack //----------begin of function PndSttTrackFinderReal::AssociateBetterAfterFitSkewHitsToXYTrack UShort_t PndSttTrackFinderReal::AssociateBetterAfterFitSkewHitsToXYTrack( UShort_t TemporarynSkewHitsinTrack, // input UShort_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 UShort_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 ) { UShort_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] 1.e-10) { *qu=*qu/(cose-*emme*sine); *emme=(*emme*cose+sine)/(cose-*emme*sine); return 1; } else { // in this case the equation is 0 = x+*qu . if(fabs(sine+*emme*cose) < 1.e-10) { cout<<" From PndSttFitSZspace, situation impossible in principle! Returning -1" < 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], UShort_t nParallelHits, UShort_t *ListParallelHits, UShort_t nSkewHits, UShort_t *ListSkewHits, Double_t *S, UShort_t *Infoparal, UShort_t *Infoskew ) { UShort_t i,j,flag, nTotal = nParallelHits+nSkewHits, 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< nSkewHits; j++){ if( S[j] < PI) { auxFiSkewvalues[j] = S[j] + 2.*PI; } else { auxFiSkewvalues[j] = S[j]; } } } else { for (j = 0; j< nSkewHits; j++){ auxFiSkewvalues[j] = S[j]; } } Merge_Sort( nSkewHits, 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], UShort_t nParallelHits, UShort_t *ListParallelHits, UShort_t *Infoparal, Short_t Charge, Double_t *Fi_initial_helix_referenceframe, Double_t *Fi_final_helix_referenceframe, Double_t *U, Double_t *V ) { UShort_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( UShort_t *Infoparal, UShort_t *Infoskew, Double_t oX, Double_t oY, Double_t Rr, UShort_t nSkewHits, UShort_t *ListSkewHits, Double_t *SList, // this is rekated to the skew hits. IMPORTANT : // the index must be the ORIGINAL skew hit number, // therefore SList[Infoskew[ListSkewHits[*]]]. Short_t Charge, UShort_t nParHits, UShort_t *ListParHits, Double_t *U, Double_t *V, UShort_t *BigList // this is the final ordered Parallel+Skew list; // already in NATIVE hit number. ) { UShort_t i,j, index[nSkewHits+nParHits], tmp[nSkewHits+nParHits], tmpList[nSkewHits]; Double_t aaa, b1, sign, aux[nSkewHits+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< nSkewHits; 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< nSkewHits; 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< nSkewHits; j++){ BigList[j+nParHits]=Infoskew[ListSkewHits[j]]; index[j+nParHits] = j+nParHits; } Merge_Sort( nSkewHits+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; } } // finding the charge of the track not necessary; already done prima. /* if( auxFivalues[0] > auxFivalues[nParallelHits-1]) { *Charge = 1; } else { *Charge = -1; } */ // now ordering of the skew hits if(nSkewHits>0){ if( flag==1) { for (j = 0; j< nSkewHits; 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< nSkewHits; j++){ auxFiSkewvalues[j] = S[Infoskew[ ListSkewHits[j] ]]; } } Merge_Sort( nSkewHits, auxFiSkewvalues, ListSkewHits); } // end of if(nSkewHits>0) // merge the parallel and skew hits for(j = 0;j< nParallelHits; j++){ BigListFi[j]=auxFivalues[j]; BigList[j] = Infoparal[ ListParallelHits[j] ] ; // cout<<"from Ordering, Parallel; j= "< Fi_low_limit and // possibly > 2PI. Short_t * status,//it is a vector; =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], UShort_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,"MacrowithRfromMCParallelHitswithMCEvent%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.05; xmin = xmin - delta*0.05; ymax = ymax + delta*0.05; ymin = ymin - delta*0.05; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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( 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, "MacroParTrack%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) goto nohits ; if( zmax < zmin ) goto nohits ; if( Smax < Smin ) goto nohits; 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 -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], UShort_t nTracksFoundSoFar, UShort_t nHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxHits], UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxHits], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool firstime, inclusionMC[nTracksFoundSoFar][nmaxHits], inclusionExp[nTracksFoundSoFar]; UShort_t ntoMCtrack[nTracksFoundSoFar], toMCtrackfrequency[nTracksFoundSoFar][nmaxHits]; UShort_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], UShort_t infopar, Double_t Ox, Double_t Oy, Double_t R, Double_t KAPPA, Double_t FI0, Short_t Charge, Double_t *Posiz ) { 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 "< cut. ) ) { return false; } if( flagInnerStt == 2) return true; //----------------------------------------------------- if(nOuterHits==0){ return true; // at this point of the code the absence of outer parallel hits // is possible (= track with very high Pz). } //------------ // find the entrance and exit of the track in the Outer Parallel Straw region, Left side. // This region is bounded by a Hexagon (inner), a Circle (outer) and it has // the target gap in the middle. flagOuterStt=FindTrackEntranceExitHexagonCircleLeft( Oxx, Oyy, Rr, Charge, Start, ApotemaMinOuterParStraw, RStrawDetMax, GAP, Xcross, Ycross ); //------------ // find the entrance and exit of the track in the Outer Parallel Straw region, Right side. // This region is bounded by a Hexagon (inner), a Circle (outer) and it has // the target gap in the middle. flagOuterStt=FindTrackEntranceExitHexagonCircleRight( Oxx, Oyy, Rr, Charge, Start, ApotemaMinOuterParStraw, RStrawDetMax, GAP, Xcross, Ycross ); //------------ // find the entrance and exit of the track in the Outer Parallel Straw region. // This region is bounded by a Hexagon (inner), a Circle (outer) and it has // the target gap in the middle. flagOuterStt=FindTrackEntranceExitHexagonCircle( Oxx, Oyy, Rr, Charge, Start, ApotemaMinOuterParStraw, RStrawDetMax, Xcross, Ycross ); if (!(flagOuterStt == 0 || flagOuterStt == 2|| flagOuterStt == -1)) return false; if (flagOuterStt == -99 ){ cout<<"PndSttTrackFinderReal::SttParalCleanup,"<< " contraddiction,outer, nIntersections[0]="<< nIntersections[0]<<"<2, returning false!\n"; return false; } //------------------------- // cleanup of the spurious tracks now using the outer parallel straws. if ( BadTrack_ParStt( Oxx, Oyy, Rr, Charge, Xcross, // Xcross[0]=point of entrance; Xcross[1]=point of exit. Ycross, false, // for the Outer Parallel Stt at this stage don't require // continuity of hits at the external border. nOuterHits, ListOuterHits, info, RStrawDetOuterParMin, RStrawDetMax, 1.1*DiameterStrawTube, // cut of proximity between hits. 1 // maximum allowed # consecutive hits with distance > cut. ) ){ return false; } //---------------------------------------------------------------------------- return true; }; //----------end of function PndSttTrackFinderReal::SttParalCleanup //----------begin of function PndSttTrackFinderReal::SttSkewCleanup bool PndSttTrackFinderReal::SttSkewCleanup( Double_t GAP, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], UShort_t nHits, Double_t *auxS, Double_t RminStrawSkew, Double_t RmaxStrawSkew, bool ConsiderLastHit, Double_t cut, UShort_t maxnum ) { Short_t flagStt; UShort_t i, ibad, nIntersections[2]; Double_t FiStart, r, Distance[nHits+1], Xcross[2], Ycross[2], XintersectionList[12][2], // second index =0 --> inner Hexagon, =1 --> outer. YintersectionList[12][2]; // first index : all the possible intersections // (up to 12 intersections). //------------------------------------------ // find the entrance and exit of the track in the Skew Straw region. // This region is bounded by two Hexagons, and it has the target gap // in the middle. So Left is the left looking from downstream. flagStt=FindTrackEntranceExitbiHexagonLeft( GAP, Oxx, Oyy, Rr, Charge, Start, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, // Apotema is the distance of a Hexagonal side from (0,0) Xcross, Ycross ); flagStt=FindTrackEntranceExitbiHexagonRight( GAP, Oxx, Oyy, Rr, Charge, Start, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, // Apotema is the distance of a Hexagonal side from (0,0) Xcross, Ycross ); // find the entrance and exit of the track in the Skew Straw region. // This region is bounded by two Hexagons, and it has the target gap in the middle. flagStt=FindTrackEntranceExitbiHexagon( Oxx, Oyy, Rr, Charge, Start, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, // Apotema is the distance of a Hexagonal side from (0,0) Xcross, Ycross ); if (!(flagStt == 0 || flagStt == 2)) return false; if (flagStt == -99 ){ cout<<"PndSttTrackFinderReal::SttSkewCleanup,"<< " contraddiction,skew, nIntersections[0]="<< nIntersections[0]<<"<2, returning false!\n"; return false; } //------------------------------------------------------------------------- // if there are also Outer Parallel hits, then require continuity of hits // also at the external border of the inner parallel straw section (ConsiderLastHit=true). ibad=0; Distance[0]= sqrt( (Oxx+Rr*cos(auxS[0])-Xcross[0])*(Oxx+Rr*cos(auxS[0])-Xcross[0]) + (Oyy+Rr*sin(auxS[0])-Ycross[0])*(Oyy+Rr*sin(auxS[0])-Ycross[0]) ); if(Distance[0]>cut){ if(Distance[0]>4.*cut){ return false; } ibad++; } for (i=1; icut){ if(Distance[i]>4.*cut){ return false; } ibad++; } } // end of do (i=1; icut ){ if( Distance[nHits]>4.*cut){ //--------stampaggio if(istampa>=3) cout<<"da SttSkewCleanup,Distance[nHits], distanza catastrofica = "<< Distance[nHits]<<", !Eliminare!\n" "\tcon il printout : Xcross[0] = "<cut ) //-------- stampaggi if(istampa>=3) { for(int ic=0;ic maxnum){ return false; } return true; }; //----------end of function PndSttTrackFinderReal::SttSkewCleanup //----------begin of function PndSttTrackFinderReal::BadTrack_ParStt bool PndSttTrackFinderReal::BadTrack_ParStt( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Xcross[2], // Xcross[0]=point of entrance; // Xcross[1]=point of exit. Double_t Ycross[2], bool ConsiderLastHit, UShort_t nHits, UShort_t* ListHits, Double_t info[][7], Double_t RStrawDetectorParMin, Double_t RStrawDetectorParMax, Double_t cut, UShort_t maxnum ) { UShort_t ibad, ihit; Double_t Distance[nHits+1]; ibad=0; Distance[0]= sqrt( (info[ListHits[0]][0]-Xcross[0])*(info[ListHits[0]][0]-Xcross[0])+ (info[ListHits[0]][1]-Ycross[0])*(info[ListHits[0]][1]-Ycross[0]) ); if(istampa>=3)cout<<"from BadTrack_ParStt, Stt || hit n. (original notation) "<< ListHits[0]<<", Distance = "<cut){ if(Distance[0]>4.*cut){ return true; } ibad++; } for(ihit=1; ihit=3)cout<<"from BadTrack_ParStt, Stt || hit n. (original notation) "<< ListHits[ihit]<<", Distance = "<cut){ if(Distance[ihit]>4.*cut){ return true; } ibad++; } } // end of for(ihit=0,NinDet=0 ;ihit=3)cout<<"from BadTrack_ParStt, Stt || hit n. (original notation) "<< ListHits[nHits-1]<<", Distance to boundary = "<cut ){ if( Distance[nHits]>4.*cut){ return true; } ibad++; } } // end of if( Distance[nHits]>cut ) //-------- stampaggi if(istampa>=3) { for(int ic=0;icFill(Distance[ic]); } if(ConsiderLastHit) hdistgoodlast->Fill( Distance[nHits]); else hdistbadlast->Fill( Distance[nHits]); } //--------------------------------------------- if( ibad > maxnum) return true; return false; } //----------end of function PndSttTrackFinderReal::BadTrack_ParStt //----------start function PndSttTrackFinderReal::IntersectionsWithClosedPolygon Short_t PndSttTrackFinderReal::IntersectionsWithClosedPolygon( //-------- inputs Double_t Ox, Double_t Oy, Double_t R, Double_t Rmi, // min Apotema of hexagonal volume intersected by track; Double_t Rma, // max Apotema of hexagonal volume intersected by track; //-------- outputs UShort_t nIntersections[2], Double_t XintersectionList[][2], Double_t YintersectionList[][2] ) { // return integer convention : // -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. // inner Hexagon --> 0 // outer Hexagon --> 1 bool internal[2], AtLeast1[2]; UShort_t i, is, j, Nintersections; Double_t mindist[2], distance, //------------------- // a,b,c == coefficients of the implicit equations of the six sides of the Hexagon // centered at 0 : a*x + b*y +c =0; see Gianluigi's logbook on page 277; // the coefficient c has to be multiplied by Erre. // The first side is a[] = { 1./sqrt(3.) , 1., -1./sqrt(3.), 1./sqrt(3.) , 1. , -1./sqrt(3.) } , b[] = { 1., 0. , 1., 1., 0. , 1. }, c[] = { -2./sqrt(3.) , -1., 2./sqrt(3.), 2./sqrt(3.), 1., -2./sqrt(3.)}, //---------------------- Erre[] = {Rmi, Rma}, // this is the distance from (0,0) // of the Verteces of the Hexagon delimiting the Skew area tempX[2], tempY[2]; // both hexagon_side_xlow and hexagon_side_xup must be multiplied by appropriate Erre; // sides of Hexagon ordered as a, b, c ..... in Gianluigi's logbook on page 280. Double_t hexagon_side_x[] = { 0., 1. , 1., 0., -1., -1., 0. }, hexagon_side_y[] = { 2./sqrt(3.),1./sqrt(3.),-1./sqrt(3.), -2./sqrt(3.),-1./sqrt(3.), 1./sqrt(3.), 2./sqrt(3.)}; //----------------------- // find intersection with the 6 sides of the small exhagon delimiting the skew straws zone // see on page 277 of Gianluigi's logbook. // status =0, at least 1 intersection with inner Hexagon, at least 1 intersection // with the outer Hexagon; =1, track contained completely between the // two Hexagons; = -1 track contained within inner Hexagon; // =-2 track outside outer Hexagon. for(i=0;i<2;i++){ // i=0 --> inner Hexagon, i= 1 --> outer Hexagon. AtLeast1[i] = false; nIntersections[i]=0; internal[i] = true; mindist[i]=999999.; for(is=0; is<6; is++){ if ( IntersectionCircle_Segment(a[is], b[is], c[is]*Erre[i], hexagon_side_x[is]*Erre[i], hexagon_side_x[is+1]*Erre[i], hexagon_side_y[is]*Erre[i], hexagon_side_y[is+1]*Erre[i], Ox, Oy, R, &Nintersections, tempX, tempY, &distance // distance of (Ox,Oy) from line // defined by a*x+b*y+c=0. ) ){ AtLeast1[i]=true; for(j=0;jdistance) mindist[i]=distance; // the definition of 'internal' here is when the given Point // stays at the same side of the origin (0,0) with respect to // the given line of equation a*x+b*y+c=0. internal[i] = internal[i] && IsInternal(Ox, Oy, a[is], b[is], c[is]*Erre[i] ); } // end of for(is=0; is<6; is++) } // end of for(i=0;i<2;i++) if( (!AtLeast1[0]) && (!AtLeast1[1]) ){ if (!internal[1]) return -2; // trajectory outside outer Hexagon. if( R > mindist[1]) return -2; if( !internal[0]) return 1; // trajectory contained between inner and outer Hexagon. if( mindist[0] >= R) return -1; // trajectory contained in inner Hexagon. return 1; // trajectory contained between inner and outer Hexagon. } else if (AtLeast1[0] && AtLeast1[1] ){ // continuation of if( (!AtLeast1[0]) && ... return 0; } else if (AtLeast1[0]){ return 2; } return 3; } //---------- end of function PndSttTrackFinderReal::IntersectionsWithClosedPolygon //----------start function PndSttTrackFinderReal::IntersectionsWithOpenPolygon UShort_t PndSttTrackFinderReal::IntersectionsWithOpenPolygon( //-------- inputs Double_t Ox, // Track parameter Double_t Oy, // Track parameter Double_t R, // Track parameter UShort_t nSides, // input, n. of Sides of open Polygon. Double_t *a, // coefficient of formula : aX + bY + c = 0 defining Double_t *b, // the Polygon sides. Double_t *c, Double_t *Side_x, // X,Y coordinate of the Sides vertices (in sequence, following Double_t *Side_y, // the Polygon along. //-------- outputs Double_t *XintersectionList, // XintersectionList Double_t *YintersectionList // YintersectionList. ) { // this methods returns the n. of intersections. UShort_t i, is, j, nIntersections, Nintersections; Double_t mindist, distance, //------------------- // a,b,c == coefficients of the implicit equations of the six sides of the Hexagon // centered at 0 : a*x + b*y +c =0; see Gianluigi's logbook on page 277; // the coefficient c has to be multiplied by Erre. tempX[2], tempY[2]; //----------------------- nIntersections=0; for(is=0; is track outside outer perimeter; // 0 --> at least 1 intersection with polygon; // 1 --> track contained completely between the two polygons; // inner Hexagon --> 0 // outer Hexagon --> 1 bool internal, AtLeast1; UShort_t i, is, j, Nintersections; Double_t aaa, maxdistq, distance, //------------------- // a,b,c == coefficients of the implicit equations of the 3 sides of the half inner Hexagon // plus 3 sides of the half outer Hexagon plus 2 vertical sides corresponding to the Gap : // a*x + b*y +c =0; the numbering of these 8 sides foolows the convention of Gianluigi's // logbook on page 286. a[] = {-1./sqrt(3.) , 1., 1./sqrt(3.), 1., 1./sqrt(3.), 1., -1./sqrt(3.), 1.}, b[] = {1., 0., 1., 0., 1., 0., 1., 0.}, c[] = {-2.*Ama/sqrt(3.),Ama, 2.*Ama/sqrt(3.),vgap/2.,2.*Ami/sqrt(3.),Ami, -2.*Ami/sqrt(3.),vgap/2.}, //---------------------- tempX[2], tempY[2]; // side_x and side_y are ordered as side1, side2, ..... in Gianluigi's logbook on page 286. Double_t side_x[] = { -vgap/2., -Ama , -Ama, -vgap/2., -vgap/2., -Ami, -Ami, -vgap/2., -vgap/2.}, side_y[] = {(-0.5*vgap+2.*Ama)/sqrt(3.), Ama/sqrt(3.), -Ama/sqrt(3.), -(-0.5*vgap+2.*Ama)/sqrt(3.), -(-0.5*vgap+2.*Ami)/sqrt(3.), -Ami/sqrt(3.), Ami/sqrt(3.), (-0.5*vgap+2.*Ami)/sqrt(3.), (-0.5*vgap+2.*Ama)/sqrt(3.) }; //----------------------- // find intersections (maximum 16) with the 8 sides. AtLeast1 = false; *nIntersections =0; internal = true; maxdistq=-9999.; for(is=0; is<8; is++){ aaa = (side_x[is]-Ox)*(side_x[is]-Ox)+(side_y[is]-Oy)*(side_y[is]-Oy); if(aaa>maxdistq) maxdistq=aaa; aaa = (side_x[is+1]-Ox)*(side_x[is+1]-Ox)+(side_y[is+1]-Oy)*(side_y[is+1]-Oy); if(aaa>maxdistq) maxdistq=aaa; if ( IntersectionCircle_Segment( a[is], b[is], c[is], side_x[is], side_x[is+1], side_y[is], side_y[is+1], Ox, Oy, R, &Nintersections, tempX, tempY, &distance // distance of (Ox,Oy) from line // defined by a*x+b*y+c=0. ) ){ AtLeast1=true; for(j=0;j track outside outer perimeter; // 0 --> at least 1 intersection with polygon; // 1 --> track contained completely between the two polygons; // inner Hexagon --> 0 // outer Hexagon --> 1 bool internal, AtLeast1; UShort_t i, is, j, Nintersections; Double_t aaa, maxdistq, distance, //------------------- // a,b,c == coefficients of the implicit equations of the 3 sides of the half inner Hexagon // plus 3 sides of the half outer Hexagon plus 2 vertical sides corresponding to the Gap : // a*x + b*y +c =0; the numbering of these 8 sides foolows the convention of Gianluigi's // logbook on page 286. a[] = {1./sqrt(3.) , 1., -1./sqrt(3.), 1., -1./sqrt(3.), 1., 1./sqrt(3.), 1.}, b[] = {1., 0., 1., 0., 1., 0., 1., 0.}, c[] = {-2.*Ama/sqrt(3.),-Ama, 2.*Ama/sqrt(3.),-vgap/2.,2.*Ami/sqrt(3.),-Ami, -2.*Ami/sqrt(3.),-vgap/2.}, //---------------------- tempX[2], tempY[2]; // side_x and side_y are ordered as side1, side2, ..... in Gianluigi's logbook on page 286. Double_t side_x[] = { vgap/2., Ama , Ama, vgap/2., vgap/2., Ami, Ami, vgap/2., vgap/2.}, side_y[] = {(-0.5*vgap+2.*Ama)/sqrt(3.), Ama/sqrt(3.), -Ama/sqrt(3.), -(-0.5*vgap+2.*Ama)/sqrt(3.), -(-0.5*vgap+2.*Ami)/sqrt(3.), -Ami/sqrt(3.), Ami/sqrt(3.), (-0.5*vgap+2.*Ami)/sqrt(3.), (-0.5*vgap+2.*Ama)/sqrt(3.) }; //----------------------- // find intersections (maximum 16) with the 8 sides. AtLeast1 = false; *nIntersections =0; internal = true; maxdistq=-9999.; for(is=0; is<8; is++){ aaa = (side_x[is]-Ox)*(side_x[is]-Ox)+(side_y[is]-Oy)*(side_y[is]-Oy); if(aaa>maxdistq) maxdistq=aaa; aaa = (side_x[is+1]-Ox)*(side_x[is+1]-Ox)+(side_y[is+1]-Oy)*(side_y[is+1]-Oy); if(aaa>maxdistq) maxdistq=aaa; if ( IntersectionCircle_Segment( a[is], b[is], c[is], side_x[is], side_x[is+1], side_y[is], side_y[is+1], Ox, Oy, R, &Nintersections, tempX, tempY, &distance // distance of (Ox,Oy) from line // defined by a*x+b*y+c=0. ) ){ AtLeast1=true; for(j=0;j 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 UShort_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 ) { UShort_t nIntersectionsCircle; Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; nIntersectionsCircle=0; aaa = sqrt(Oxx*Oxx+Oyy*Oyy); // preliminary condition for having intersections between trajectory and Circle. if( !( aaa >= Rr + Rma || Rr >= aaa + Rma) && !( aaa + Rr <= Rma || aaa - Rr >= Rma ) ) { // now the calculation FI0 = atan2(-Oyy,-Oxx); cosFi = (aaa*aaa + Rr*Rr - Rma*Rma)/(2.*Rr*aaa); if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.; Fi = acos(cosFi); // (x1, y1) and (x2,y2) are the intersections between the trajectory // and the circle, in the laboratory reference frame. x1 = Oxx+Rr*cos(FI0 - Fi); y1 = Oyy+Rr*sin(FI0 - Fi); theta1 = atan2(y1,x1); // in this way theta1 is between -PI and PI. x2 = Oxx+Rr*cos(FI0 + Fi); y2 = Oyy+Rr*sin(FI0 + Fi); theta2 = atan2(y2,x2); // in this way theta2 is between -PI and PI. // Theta1, Theta2 = angle of the edges of the outer circle + Gap in the laboratory frame. // Theta1, Theta2 are angles between -PI/2 and +PI/2. if(!left){ // Right (looking into the beam) Semicircle. Theta2 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.); Theta1 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.); if( Theta1<= theta1 && theta1 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x1; YintersectionList[nIntersectionsCircle]=y1; nIntersectionsCircle++; } if( Theta1<= theta2 && theta2 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x2; YintersectionList[nIntersectionsCircle]=y2; nIntersectionsCircle++; } } else { // Left (looking into the beam) Semicircle. Theta2 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.); Theta1 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.); if( Theta1<= theta1 || theta1 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x1; YintersectionList[nIntersectionsCircle]=y1; nIntersectionsCircle++; } if( Theta1<= theta2 || theta2 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x2; YintersectionList[nIntersectionsCircle]=y2; nIntersectionsCircle++; } } } // end of if (!( a >= Rr + Rma || ..... //---------------------------- end of calculation intersection with outer circle. return nIntersectionsCircle; } //----------end of function 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, UShort_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 ) { UShort_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. UShort_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. UShort_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 ) { UShort_t i, j; // this method works under the hypothesis that there are at least 2 intersections. if (nIntersections<2) return; if(Charge > 0) { UShort_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)++; }else{ ListInnerHits[ *nInnerHits ] = ListHits[ihit]; (*nInnerHits)++; } } } //----------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; UShort_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; UShort_t nIntersections; Double_t FiStart, XintersectionList[16], // all the possible intersections YintersectionList[16]; // (up to 16 intersections). // The following is the form of the Left (looking from downstream into the beam) biHexagon // geometrical shape considered in this method : // /* /| / | / | / / / / / / / / | | | | | | | | \ \ \ \ \ \ \ \ \ | \ | \| */ // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; flag=IntersectionsWithClosedbiHexagonLeft( vgap, Oxx, Oyy, Rr, ApotemaMin, // Apotema of the inner Hexagon, ApotemaMax,// Apotema of the outer Hexagon. &nIntersections, XintersectionList, // XintersectionList[..][0] --> inner polygon, // XintersectionList[..][1] --> outer polygon. YintersectionList ); // IMPORTANT : // this is true because here it is assumed that the track comes from (0,0,0) // otherwise the code must be changed! if (!(flag == 0)) return flag; if( nIntersections<2) return -1; //------- among all possible intersection find the entrance point (Xcross[0], Ycross[0]) // and the exit point (Xcross[1], Ycross[1]) of this track. //-------- the starting point of the track. FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); // 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; UShort_t nIntersections; Double_t FiStart, XintersectionList[16], // all the possible intersections YintersectionList[16]; // (up to 16 intersections). // The following is the form of the Left (looking from downstream into the beam) biHexagon // geometrical shape considered in this method : // /* |\ | \ | \ \ \ \ \ | | | | / / / / / / | / | / |/ */ // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; flag=IntersectionsWithClosedbiHexagonRight( vgap, Oxx, Oyy, Rr, ApotemaMin, // Apotema of the inner Hexagon, ApotemaMax,// Apotema of the outer Hexagon. &nIntersections, XintersectionList, // 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; UShort_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; //------------------ UShort_t nIntersectionsCircle, nIntersections; Double_t FiStart, XintersectionList[12], YintersectionList[12]; // all the possible intersections (up to 12 intersections). // finding all possible intersections with inner parallel straw region. Double_t Side_x[] = { -GAP/2., -GAP/2. , -ApotemaMin, -ApotemaMin, -GAP/2., -GAP/2. }, Side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-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; //------------------ UShort_t nIntersectionsCircle, nIntersections; Double_t FiStart, XintersectionList[12], YintersectionList[12]; // all the possible intersections (up to 10 intersections). // finding all possible intersections with inner parallel straw region. Double_t Side_x[] = { GAP/2., GAP/2. , ApotemaMin, ApotemaMin, GAP/2., GAP/2. }, Side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-GAP)/sqrt(3.), ApotemaMin/sqrt(3.), -ApotemaMin/sqrt(3.), -(2.*ApotemaMin-GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)}, a[] = {1., 1./sqrt(3.), 1., -1./sqrt(3.), 1.}, b[] = {0., 1., 0., 1., 0.}, c[] = {-GAP/2.,-2.*ApotemaMin/sqrt(3.),-ApotemaMin, 2.*ApotemaMin/sqrt(3.), -GAP/2.}; nIntersections=IntersectionsWithOpenPolygon( Oxx, Oyy, Rr, 5, // n. Sides of open Polygon. a, // coefficient of formula : aX + bY + c = 0 defining the Polygon sides. b, c, Side_x, // X,Y coordinate of the Sides vertices (in sequence, following Side_y, // the Polygon along. XintersectionList, // XintersectionList YintersectionList // YintersectionList. ); //------------------------------------------------------------------------- // finding intersections of trajectory [assumed to originate from (0,0) ] // with outer semicircle, the Left part. nIntersectionsCircle=IntersectionsWithGapSemicircle( Oxx, // input from trajectory Oyy, // input from trajectory Rr, // input from trajectory GAP, // input, vertical gap in XY plane of STT detector. false, // true --> Left semicircle, false --> Right semicircle. Rma, // radius of external Stt containment. &XintersectionList[nIntersections], // output, X list of intersections (maximal 2). &YintersectionList[nIntersections] ); nIntersections += nIntersectionsCircle; //-------- the starting point of the track. if(nIntersections<2) return -1; FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); // 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 ClassImp(PndSttTrackFinderReal)