// ------------------------------------------------------------------------- // ----- PndSttTrackFinderReal source file ----- // ----- Created 28/03/06 by V. Friese ----- // ------------------------------------------------------------------------- #include "glpk.h" // Pnd includes #include "PndSttTrackFinderReal.h" #include "PndSttHit.h" #include "PndSttPoint.h" #include "PndSttHelixHit.h" #include "PndTrackCand.h" #include "PndDetectorList.h" #include "PndSttTube.h" #include #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() { // fVerbose = 1; 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; RminStrawSkewArea = RStrawDetectorMin*2./1.732051 + 18.*StrawRadius ; // delimitation of the skew area RmaxStrawSkewArea = RminStrawSkewArea + 8.*1.732051 *StrawRadius ; } // ------------------------------------------------------------------------- // ----- Standard constructor ------------------------------------------ PndSttTrackFinderReal::PndSttTrackFinderReal(int verbose) { istampa = verbose; 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; RminStrawSkewArea = RStrawDetectorMin*2./1.732051 + 18.*StrawRadius ; // delimitation of the skew area RmaxStrawSkewArea = RminStrawSkewArea + 8.*1.732051 *StrawRadius ; } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndSttTrackFinderReal::~PndSttTrackFinderReal() { } // ------------------------------------------------------------------------- // ----- Public method Init -------------------------------------------- void PndSttTrackFinderReal::Init() { UShort_t i, j, k; Double_t r1, r2, A , rConfMin, rConfMax, tempRadiaConf[nRdivConformal]; MINIMUMHITSPERTRACK=3; // --------------------------- opening files for special purposes N_INTENDED=0; if(istampa >=2 ){ //---- 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 >=2) // ------------------------- fHelixHitProduction = true; IVOLTE=0; ntimes = 0; // Get and check FairRootManager FairRootManager* ioman = FairRootManager::Instance(); // get the MCTrack array fMCTrackArray = (TClonesArray*) ioman->ActivateBranch("MCTrack"); // ------------------------------------------------------------------- // -------------------- initializations for the Kalman with Genfit later /* fSttHitArray=(TClonesArray*) ioman->GetObject("STTHit"); if(fSttHitArray==0){ Error("PndSttKalmanTask2::Init","stt pattern recognition; hit-array not found!"); return ; } // Build hit factory ----------------------------- _theRecoHitFactory = new GFRecoHitFactory(); _theRecoHitFactory->addProducer(3, new GFRecoHitProducer(fSttHitArray)); // */ // -------------------- end initializations for the Kalman with Genfit later // hx = new TH1F("hx", "Associated Z", 400, -200., 200.); 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*RStrawDetectorMax - r1)/nRdivConformal; A = (RStrawDetectorMax - r1)/nRdivConformal; if ( nRdivConformal > 1 ) { for(i = 1; i< nRdivConformal ; i++){ r2 = r1 + A; // tempRadiaConf[nRdivConformal-i] = 1./sqrt(r2); tempRadiaConf[nRdivConformal-i] = 1./r2; r1=r2; } } /* for(i = 1; i< nRdivConformal ; i++){ cout<<"from Init : temporary CONFORMAL radius n. "<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 // ----------------- start of method DoFind Int_t PndSttTrackFinderReal::DoFind(TClonesArray* trackArray, TClonesArray* helixHitArray) { // list of hit numbers of those falling in this cell Short_t Charge[MAXTRACKSPEREVENT], Status[MAXTRACKSPEREVENT], daTrackFoundaTrackMC[MAXTRACKSPEREVENT]; // daMCTrackaTrackFound[MAXMCTRACKS]; UShort_t emme, auxIndex[nmaxHits], OLDinfoparal[nmaxHits], enne[MAXTRACKSPEREVENT][nmaxHits]; UShort_t istep,inclination_type; 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>=2 && IVOLTE <= nmassimo) cout<<"\nda PndSttTrackFinderReal : evento (a partire da 1) n. "<GetEntriesFast(); // num. tracce/evento if(istampa>=2) cout<<"evento n. "<= nmaxHits ) { cout<<"from PndSttTrackFinderReal : # Stt hits = "<= nmaxHits = " <=2 && (nMCTracks != N_INTENDED) ){ cout<<"Gianluigi : n. MC tracks = "< > hitMap; for(Int_t j=0; j= 2 ) cout<<"Gianluigi : da PndSttTrackFinderReal::DoFind : Nhits="<GetRefIndex(); if (ptIndex < 0) continue; // fake or background hit pMCpt = GetPointFromCollections(iHit); // <== FairMCPoint 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(); info[iHit][6]= pMCpt->GetTrackID(); 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 ----------------------------- veritaMC[iHit][0]= pMCpt->GetX(); veritaMC[iHit][1]= pMCpt->GetY(); veritaMC[iHit][2]= pMCpt->GetZ(); FromHitToMCTrack[iHit] = (UShort_t) ( info[iHit][6] + 0.001 ); /* if(info[iHit][5]==1. ){ // associazione hits paralleli FromMCTrackToHit[ FromHitToMCTrack[iHit] ][ nHitsInMCTrack[ FromHitToMCTrack[iHit] ] ] = iHit; nHitsInMCTrack[ FromHitToMCTrack[iHit] ]++ ; } else { // associazione hits skew FromMCTrackToSkewHit[ FromHitToMCTrack[iHit] ][ nSkewHitsInMCTrack[ FromHitToMCTrack[iHit] ] ] = iHit; nSkewHitsInMCTrack[ FromHitToMCTrack[iHit] ]++ ; } // end of if(info[iHit][5]==1. ) */ //--------------- inizio stampaggi, stampe di controllo if (istampa >= 2 && IVOLTE<= nmassimo) { cout <<"iHit "<< iHit << endl; cout <<" hit X, Y, Z space position " << veritaMC[iHit][0] << " " << veritaMC[iHit][1] << " " << veritaMC[iHit][2]<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<<" "<GetTrackID()<= //-------- fine stampaggi } // end of for (Int_t iHit = 0; //--------------- inizio stampaggi // if (istampa >= 2 && IVOLTE<= nmassimo) { if (istampa >= 1 ) { 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 HitsinBoxConformal[nRdivConformal][nFidivConformal][nmaxHits]; // first index -> radial divisions, 2nd index -> azimuthal divisions; 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, ipinco=0; Double_t angle, max1, HoughR, HoughD, HoughFi, HoughKAPPA, HoughFI0, 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], Fi_allowedforskew_low[MAXTRACKSPEREVENT], Fi_allowedforskew_up[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[nmaxHits], Z[nmaxHits], temporeS[nmaxHits], temporeZ[nmaxHits], ZDrift[nmaxHits], ZErrorafterTilt[nmaxHits], temporeZDrift[nmaxHits], temporeZErrorafterTilt[nmaxHits], Posiz[3], 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++) /* // first the parallel straws for(i=0; i< Minclinations[0]-1; i++){ if( ExclusionList[ infoparal[i] ] ){ for(j=i+1; j< Minclinations[0]; j++){ if( fabs(info[ infoparal[i] ][0] - info[ infoparal[j] ][0])<1.e-20 && fabs(info[ infoparal[i] ][1] - info[ infoparal[j] ][1])<1.e-20 ) { ExclusionList[ infoparal[j] ]= false ; ExclusionListbis[ infoparal[j] ]= false ; } } // end of for(j=i+1; j< Minclinations[0]; j++) } // end of if( ExclusionList[ infoparal[i] ] ) } // end of for(i=0; i< Minclinations[0]-1; i++) */ // then the skew straws /* for(i=0; i< NSkewhits-1; i++){ if( ExclusionListSkew[ infoskew[i] ] ){ for(j=i+1; j< NSkewhits; j++){ if( fabs(info[ infoskew[i] ][0] - info[ infoskew[j] ][0])<1.e-20 && fabs(info[ infoskew[i] ][1] - info[ infoskew[j] ][1])<1.e-20 ) { ExclusionListSkew[ infoskew[j] ]= false ; ExclusionListSkewbis[ infoskew[j] ]= false ; } } // end of for(j=1; j< NSkewhits; j++) } // end of if( ExclusionList[ infoskew[i] ] ) } // end of for(i=0; i< Minclinations[0]-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( nHitsinTrack[nTracksFoundSoFar] < MINIMUMHITSPERTRACK) { // cout<<" nHitsinTrack = "<= 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); if( Naux >= MINIMUMOUTERHITSPERTRACK && Naux > 0.7 * Nouter ) break; if( Naux >= MINIMUMOUTERHITSPERTRACK) { // further collection of hits in the inner region but this time strictly connected to the outer ones // first the list of non outer hits for(j=Nouter;j= MINIMUMOUTERHITSPERTRACK) //--------------------------- // finding the rotation angle for best utilization of the MILP procedure for(j=0, rotationcos=0., rotationsin=0.; j= MAXTRACKSPEREVENT ){ cout<<"from PndSttTrackFinderReal : # n. Tracks found so far = "<= MAXTRACKSPEREVENT ( = "<=0 ){ nSkewHitsinTrack[i] = NNN; for(j=0;j=0 ) nSkewHitsinTrack[i]=TemporarynSkewHitsinTrack; for(j=0;j=0 ) // --------- 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] ; } // --------- // ------------------------------------------------------------- fine: ; //---------------------------------------------- some printouts if(istampa>=2 && IVOLTE<= nmassimo) { cout<<" Evento n. "<=2) { cout<<"da TrackFinder Real : nTracksFoundSoFar "<=2) { cout<<"Traccia n. "<0 && nTracksFoundSoFar > 0 ){ /* AssociateFoundTrackstoMC( info, nTracksFoundSoFar, nHitsinTrack, ListHitsinTrack, nSkewHitsinTrack, ListSkewHitsinTrack, daTrackFoundaTrackMC ); } */ AssociateFoundTrackstoMCbis( info, nTracksFoundSoFar, nHitsinTrack, ListHitsinTrack, nSkewHitsinTrack, ListSkewHitsinTrack, daTrackFoundaTrackMC ); } for(jexp=0; jexp=2 ) fprintf(HANDLE, "\n Evento %d NTotaleTracceMC %d ------\n",IVOLTE, nMCTracks); for (ii=0; ii=2 ;ii++){ 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-20 ){ fprintf(HANDLE, " TracciaMC %d; sua FoundTrack associata (n. %d) NONsoddisfaRequisitiMinimi perche' KAPPA troppo piccolo; 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, Px, Py ; 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) //-------------------- 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 Double_t Ptras,Pxini,Pyini,Pzini,dista, qop; if(istampa>2) { cout<<" da TrackFinder Real : ancora nTracksFoundSoFar "<= 2){ ii=daTrackFoundaTrackMC[i]; } dista=sqrt( Ox[i]*Ox[i]+Oy[i]*Oy[i] ); Ptras = R[i]*0.003*BFIELD; Pxini = -Charge[i]*Ptras*Oy[i]/dista; Pyini = Charge[i]*Ptras*Ox[i]/dista; TVector3 posSeed(0.,0.,0.); // direction in starting point if( GoodSkewFit[i] ) { if(fabs(KAPPA[i])>1.e-20 ){ Pzini = -Charge[i]*0.003*BFIELD/KAPPA[i]; } else { Pzini = 999999.; } new((*trackArray)[ipinco]) PndTrackCand; PndTrackCand *pTrckCand = (PndTrackCand*) trackArray->At(ipinco); ipinco++; TVector3 dirSeed(Pxini, Pyini, Pzini); // momentum direction in starting point qop = Charge[i]/dirSeed.Mag(); dirSeed.SetMag(1.); if(istampa > 2){ cout<<" da PndSttTrackFinderReal; caricato PndTrackCand n. "<At(ipinco); ipinco++; TVector3 dirSeed(Pxini, Pyini, Pzini); // momentum direction in starting point qop = Charge[i]/Ptras; // as if Pz=0 if(istampa > 2){ cout<<" da PndSttTrackFinderReal; no Kappa information from fit, caricato PndTrackCand n. "<=3) cout<<"da PndSttTrackFinderReal: DoFind, paralleli, infoparal[ ListHitsinTrack[i][j] ] = "<< infoparal[ ListHitsinTrack[i][j] ]<< ", X = "<< Posiz[0]<< ", Y = "<< Posiz[1]<< ", Z = "<< Posiz[2]<=3) cout<<"DoFind, skew, infoskew[ ListSkewHitsinTrack[i][j] ] = "<< infoskew[ ListSkewHitsinTrack[i][j] ]<<", ListSkewHitsinTrack[i][j] = "<< ListSkewHitsinTrack[i][j]<< ", X = "<< Posiz[0]<< ", Y = "<< Posiz[1]<< ", Z = "<< Posiz[2]< 1.e-20 ) { pos.SetXYZ(Xpos_for_LHeTrack[iHit], Ypos_for_LHeTrack[iHit], Zpos_for_LHeTrack[iHit]); } TVector3 posErr(0, 0, 0); // CHECK put the errors PndSttHelixHit * helixhit = new(clref[size]) PndSttHelixHit(ListPointer_to_Hit[iHit]->GetDetectorID(), ListPointer_to_Hit[iHit]->GetTubeID(), iHit, ListPointer_to_Hit[iHit]->GetRefIndex(), pos, posErr, ListPointer_to_Hit[iHit]->GetIsochrone(), ListPointer_to_Hit[iHit]->GetIsochroneError(), 0.0); } // end of for(iHit=0; iHit0 && nMCTracks=3) cout<<"DoFind, paralleli, n. hits (original notation) = "<< ParalCommonList[i][j] << ", X = "<< Posiz[0]<< ", Y = "<< Posiz[1]<< ", Z = "<< Posiz[2]<=2){ 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=3) cout<<"DoFind, skew, n. hits (original notation) = "<< SkewCommonList[i][j] << ", X = "<< Posiz[0]<< ", Y = "<< Posiz[1]<< ", Z = "<< Posiz[2]<=2){ 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 ) //--------------- end of printouts of comparison with MC for judging algorithm performance //--------------------------------------------------------------------------------------------------------- //------------------------------------- macro di display delle skew if(iplotta && IVOLTE <= nmassimo&& nMCTracks= 0){ WriteMacroSkewAssociatedHitswithMC( KAPPA[i],FI0[i], HoughD, HoughFi, HoughR, info, Nincl,Minclinations,inclination, i,0, nSkewHitsinTrack[i], ListSkewHitsinTrack, nSkewCommon[i], SkewCommonList, daTrackFoundaTrackMC[i], nMCSkewAlone, MCSkewAloneList ); } /* plottamentiSkewconMassimo( imaxima,jmaxima, Nremaining2, RemainingKAPPA, RemainingFI0, KAPPAlow, KAPPAup, FI0low, FI0up, ResultAssociatedHits,HoughR, HoughD, HoughFi); for (i=0; i< ResultAssociatedHits.NAssociatedHits ; i++){ for(j=0; j< ResultAssociatedHits.mAmbiguities[i] ; j++){ hx->Fill(ResultAssociatedHits.AssociatedHitsCoordinates[2][j][i]); } } */ } // 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 nAdmittedRadia*StrawRadius ) continue; if(distance > R ) distance *= -1.; if( KAPPA == 0.) { ResultAssociatedHits.AssociatedHitsCoordinates[0][0][ResultAssociatedHits.NAssociatedHits]= info[i][0]+info[i][3]*dx/distance; // X of the hit ResultAssociatedHits.AssociatedHitsCoordinates[1][0][ResultAssociatedHits.NAssociatedHits]= info[i][1]+info[i][3]*dy/distance; // Y of the hit ResultAssociatedHits.AssociatedHitsCoordinates[2][0][ResultAssociatedHits.NAssociatedHits] = -1.e20; ResultAssociatedHits.mAmbiguities[ResultAssociatedHits.NAssociatedHits]=1; ResultAssociatedHits.AssociatedWireDirection[0][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][0]; ResultAssociatedHits.AssociatedWireDirection[1][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][1]; ResultAssociatedHits.AssociatedWireDirection[2][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][2]; ResultAssociatedHits.AssociatedWireCenter[0][ResultAssociatedHits.NAssociatedHits]=info[i][0]; ResultAssociatedHits.AssociatedWireCenter[1][ResultAssociatedHits.NAssociatedHits]=info[i][1]; ResultAssociatedHits.AssociatedWireCenter[2][ResultAssociatedHits.NAssociatedHits]=info[i][2]; ResultAssociatedHits.AssociatedWireDriftRadius[ResultAssociatedHits.NAssociatedHits]=info[i][3]; ResultAssociatedHits.Hitnumber[ResultAssociatedHits.NAssociatedHits]=i; // temporaneo! ResultAssociatedHits.NAssociatedHits++; ResultAssociatedHits.NAssociatedParallelHits++; continue; } else { diff = FI0-angle; major = 0.5*(KAPPA*(info[i][2]+info[i][4])+diff)/PI ; minor = 0.5*(KAPPA*(info[i][2]-info[i][4])+diff)/PI ; if ( major < minor) { aaa=minor; minor=major; major=aaa; } minor <= 0. ? nlow = (int) minor : nlow = ( (int) minor) + 1; major < 0. ? nup = ((int) major)-1 : nup = (int) major; if( nlow> nup) continue; if( nup-nlow+1 > nmaxAmbiguities ) continue; step = fabs(2.*PI/KAPPA) ; aaa = -diff/KAPPA ; bbb = info[i][0]+info[i][3]*dx/distance; // X of the hit ccc = info[i][1]+info[i][3]*dy/distance; // Y of the hit for(ii=nlow; ii<=nup; ii++){ ResultAssociatedHits.AssociatedHitsCoordinates[0][ii-nlow][ResultAssociatedHits.NAssociatedHits]= bbb; // X of the hit ResultAssociatedHits.AssociatedHitsCoordinates[1][ii-nlow][ResultAssociatedHits.NAssociatedHits]= ccc ; // Y of the hit ResultAssociatedHits.AssociatedHitsCoordinates[2][ii-nlow][ResultAssociatedHits.NAssociatedHits] = aaa + ii*step; // Z of the hit } ResultAssociatedHits.mAmbiguities[ResultAssociatedHits.NAssociatedHits] = nup - nlow + 1; ResultAssociatedHits.AssociatedWireDirection[0][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][0]; ResultAssociatedHits.AssociatedWireDirection[1][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][1]; ResultAssociatedHits.AssociatedWireDirection[2][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][2]; ResultAssociatedHits.AssociatedWireCenter[0][ResultAssociatedHits.NAssociatedHits]=info[i][0]; ResultAssociatedHits.AssociatedWireCenter[1][ResultAssociatedHits.NAssociatedHits]=info[i][1]; ResultAssociatedHits.AssociatedWireCenter[2][ResultAssociatedHits.NAssociatedHits]=info[i][2]; ResultAssociatedHits.AssociatedWireDriftRadius[ResultAssociatedHits.NAssociatedHits]=info[i][3]; ResultAssociatedHits.Hitnumber[ResultAssociatedHits.NAssociatedHits]=i; // temporaneo! ResultAssociatedHits.NAssociatedHits++; ResultAssociatedHits.NAssociatedParallelHits++; } // end if( KAPPA == 0.) } else { //------------------------------ starts here association of skew straw hits // calculation of the intersection points between skew straw axis and trajectory cylinder 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]; 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) ); if( distance >= info[i][4] ) continue; 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; // 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 ) continue; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; // now calculate the distance (along the Z direction) from the edge of the ellipsis and the predicted Helix // trajectory if(fabs(KAPPA) > 1.e-10) { z1 = POINTS1[j+2] + Aellipsis1*Tiltdirection1[0]; y1 = fi1 + Aellipsis1*Tiltdirection1[1] ; zpos1 = (y1 - FI0)/KAPPA; d1 = fabs(zpos1 - z1); z2 = POINTS1[j+2] - Aellipsis1*Tiltdirection1[0]; y2 = fi1 - Aellipsis1*Tiltdirection1[1] ; zpos2 = (y2 - FI0)/KAPPA; d2 = fabs(zpos2 - z2); if ( d2>d1 ) { distance = d1; if( distance < nAdmittedRadia * Aellipsis1 ) { ResultAssociatedHits.AssociatedHitsCoordinates[0][0][ResultAssociatedHits.NAssociatedHits]= R*cos(y1) + Ox; // X of the hit ResultAssociatedHits.AssociatedHitsCoordinates[1][0][ResultAssociatedHits.NAssociatedHits]= R*sin(y1) + Oy; // Y of the hit ResultAssociatedHits.AssociatedHitsCoordinates[2][0][ResultAssociatedHits.NAssociatedHits] = z1; ResultAssociatedHits.mAmbiguities[ResultAssociatedHits.NAssociatedHits]=1; ResultAssociatedHits.AssociatedWireDirection[0][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][0]; ResultAssociatedHits.AssociatedWireDirection[1][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][1]; ResultAssociatedHits.AssociatedWireDirection[2][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][2]; ResultAssociatedHits.AssociatedWireCenter[0][ResultAssociatedHits.NAssociatedHits]=info[i][0]; ResultAssociatedHits.AssociatedWireCenter[1][ResultAssociatedHits.NAssociatedHits]=info[i][1]; ResultAssociatedHits.AssociatedWireCenter[2][ResultAssociatedHits.NAssociatedHits]=info[i][2]; ResultAssociatedHits.AssociatedWireDriftRadius[ResultAssociatedHits.NAssociatedHits]=info[i][3]; ResultAssociatedHits.Hitnumber[ResultAssociatedHits.NAssociatedHits]=-i; // temporaneo! ResultAssociatedHits.NAssociatedHits++; ResultAssociatedHits.NAssociatedSkewHits++; } } else { distance = d2 ; if( distance < nAdmittedRadia * Aellipsis1 ) { ResultAssociatedHits.AssociatedHitsCoordinates[0][0][ResultAssociatedHits.NAssociatedHits]= R*cos(y2) + Ox; // X of the hit ResultAssociatedHits.AssociatedHitsCoordinates[1][0][ResultAssociatedHits.NAssociatedHits]= R*sin(y2) + Oy; // Y of the hit ResultAssociatedHits.AssociatedHitsCoordinates[2][0][ResultAssociatedHits.NAssociatedHits] = z2; ResultAssociatedHits.mAmbiguities[ResultAssociatedHits.NAssociatedHits]=1; ResultAssociatedHits.AssociatedWireDirection[0][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][0]; ResultAssociatedHits.AssociatedWireDirection[1][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][1]; ResultAssociatedHits.AssociatedWireDirection[2][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][2]; ResultAssociatedHits.AssociatedWireCenter[0][ResultAssociatedHits.NAssociatedHits]=info[i][0]; ResultAssociatedHits.AssociatedWireCenter[1][ResultAssociatedHits.NAssociatedHits]=info[i][1]; ResultAssociatedHits.AssociatedWireCenter[2][ResultAssociatedHits.NAssociatedHits]=info[i][2]; ResultAssociatedHits.AssociatedWireDriftRadius[ResultAssociatedHits.NAssociatedHits]=info[i][3]; ResultAssociatedHits.Hitnumber[ResultAssociatedHits.NAssociatedHits]=-i; // temporaneo! ResultAssociatedHits.NAssociatedHits++; ResultAssociatedHits.NAssociatedSkewHits++; } } if(istampa>= 3 && IVOLTE<= nmassimo) { double dista; dista = distance; if(z1 z1 && zpos2 > z1 && zpos1 < z2 && zpos2 < z2 ) dista = -distance; }else{ if( zpos1 > z2 && zpos2 > z2 && zpos1 < z1 && zpos2 < z1 ) dista = -distance;} ; cout<<"distanza minima (con segno) trovata skew-Helix = "< FI0 ? factor = -1. : factor = 1.; ResultAssociatedHits.AssociatedHitsCoordinates[0][0][ResultAssociatedHits.NAssociatedHits]= R*cos(y1+factor*Bellipsis1) + Ox; // X of the hit ResultAssociatedHits.AssociatedHitsCoordinates[1][0][ResultAssociatedHits.NAssociatedHits]= R*sin(y1+factor*Bellipsis1) + Oy; // Y of the hit } else { y2 > FI0 ? factor = -1. : factor = 1.; ResultAssociatedHits.AssociatedHitsCoordinates[0][0][ResultAssociatedHits.NAssociatedHits]= R*cos(y2+factor*Bellipsis1) + Ox; // X of the hit ResultAssociatedHits.AssociatedHitsCoordinates[1][0][ResultAssociatedHits.NAssociatedHits]= R*sin(y2+factor*Bellipsis1) + Oy; // Y of the hit } ResultAssociatedHits.AssociatedHitsCoordinates[2][0][ResultAssociatedHits.NAssociatedHits] = (z1+z2)/2.; ResultAssociatedHits.mAmbiguities[ResultAssociatedHits.NAssociatedHits]=1; ResultAssociatedHits.AssociatedWireDirection[0][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][0]; ResultAssociatedHits.AssociatedWireDirection[1][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][1]; ResultAssociatedHits.AssociatedWireDirection[2][ResultAssociatedHits.NAssociatedHits]=inclination[Kincl][2]; ResultAssociatedHits.AssociatedWireCenter[0][ResultAssociatedHits.NAssociatedHits]=info[i][0]; ResultAssociatedHits.AssociatedWireCenter[1][ResultAssociatedHits.NAssociatedHits]=info[i][1]; ResultAssociatedHits.AssociatedWireCenter[2][ResultAssociatedHits.NAssociatedHits]=info[i][2]; ResultAssociatedHits.AssociatedWireDriftRadius[ResultAssociatedHits.NAssociatedHits]=info[i][3]; ResultAssociatedHits.Hitnumber[ResultAssociatedHits.NAssociatedHits]=-i; // temporaneo! ResultAssociatedHits.NAssociatedHits++; ResultAssociatedHits.NAssociatedSkewHits++; } } // end of if(fabs(KAPPA) > 1.e-10) } // end of for( ii=0; ii<2; ii++) } // end of if( info[i][5] == 1 ) else } // end of for( i=1; i< Nhits; i++) return ResultAssociatedHits; } //----------end of function PndSttTrackFinderReal::PndSttTrkAssociatedHitsToHelix //----------begin function PndSttTrackFinderReal::PndSttTrkAssociatedParallelHitsToHelix UShort_t PndSttTrackFinderReal::PndSttTrkAssociatedParallelHitsToHelix( Double_t Ox, Double_t Oy, Double_t R, Int_t Nhits, Double_t info[][7], UShort_t *auxListHitsinTrack // this is the output ) { Int_t i; UShort_t Nassociatedhits; Double_t dx, dy,distance; // Ox = (D+R)*cos(Fi); // Oy = (D+R)*sin(Fi); // association of hits in the XY plane Nassociatedhits=0; for( i=0; i< Nhits; i++) { // Kincl = (int) info[i][5] - 1; // -------------------- association of the hits from parallel straws // if( info[i][5] == 1. ) { // parallel straws dx = -Ox+info[ infoparal[i] ][0]; dy = -Oy+info[ infoparal[i] ][1]; distance = sqrt(dx*dx+dy*dy); //cout<<"nuov, R "< 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]KAPPAlow && RemainingKAPPA[ii]FI0low && RemainingFI0[ii]= 0){ hXdiffparallel.Fill(ResultAssociatedHits.AssociatedHitsCoordinates[0][j][i]-veritaMC[ResultAssociatedHits.Hitnumber[i]][0]); hYdiffparallel.Fill(ResultAssociatedHits.AssociatedHitsCoordinates[1][j][i]-veritaMC[ResultAssociatedHits.Hitnumber[i]][1]); hZdiffparallel.Fill(ResultAssociatedHits.AssociatedHitsCoordinates[2][j][i]-veritaMC[ResultAssociatedHits.Hitnumber[i]][2]); } else { hXdiffskew.Fill(ResultAssociatedHits.AssociatedHitsCoordinates[0][j][i]-veritaMC[-ResultAssociatedHits.Hitnumber[i]][0]); hYdiffskew.Fill(ResultAssociatedHits.AssociatedHitsCoordinates[1][j][i]-veritaMC[-ResultAssociatedHits.Hitnumber[i]][1]); hZdiffskew.Fill(ResultAssociatedHits.AssociatedHitsCoordinates[2][j][i]-veritaMC[-ResultAssociatedHits.Hitnumber[i]][2]); } } } hfile.Write(nome); hfile.Close(); } //----------end of function PndSttTrackFinderReal::plottamentiSkewconMassimo void PndSttTrackFinderReal::findmaximaDFiR( UShort_t BoxDFiR[nbinD][nbinFi][nbinR], Int_t MINIMUMCOUNTS, Int_t * NumberofMaximaDFiR, Int_t MaximaIndexesDFiR[][3], Int_t * STATUS) { Int_t i, j, iD, iFi, iR, icount, ntotClusters, NinCluster, nClusterElementsFound, nRemai, nRemainingElements, max, nElementsinCluster[MAXElementsOverThresholdinHough]; UShort_t found[MAXElementsOverThresholdinHough][MAXElementsOverThresholdinHough][3], auxDFiRIndex[MAXElementsOverThresholdinHough][3], Remai[MAXElementsOverThresholdinHough][3], RemainingElements[MAXElementsOverThresholdinHough][3], Cluster[MAXElementsOverThresholdinHough][3], ClusterElementsFound[MAXElementsOverThresholdinHough][3]; if(istampa >= 3 && IVOLTE<= nmassimo) { cout<<"da findmaximaDFiR : MINIMUMCOUNTS = "<= "<=3 && IVOLTE<= nmassimo) { cout<<"da findmaximaDFiR : numero di celle superiori al MINIMUM CUT = "<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; iGetOutFile(); TFile* file = FairRootManager::Instance()->GetOutFile(); file->cd(); file->mkdir("PndSttTrackFinderReal"); file->cd("PndSttTrackFinderReal"); // hx->Write(); // delete hx; } //----------end of function PndSttTrackFinderReal::WriteHistograms //----------start of function PndSttTrackFinderReal::WriteMacroParallelHitsGeneral void PndSttTrackFinderReal::WriteMacroParallelHitsGeneral( 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, 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, rrr, angle, minor, major, distance, Rx, Ry, LL, Aellipsis1, Bellipsis1,fi1, fmin, fmax, offset, step, SkewInclWithRespectToS, zpos, zpos1, zpos2, Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; char nome[300], nome2[300]; //---------- parallel straws Macro now sprintf(nome,"MacroGeneralParallelHitsEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); xmin=1.e20; xmax=-1.e20; ymin=1.e20; ymax=-1.e20; for( i=0; i< Nhits; i++) { 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); } } //------------------------------- 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 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); //------------------------------------------------------------------------------------------------------------ // 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( 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 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); } } } fprintf(MACRO,"}\n"); fclose(MACRO); return ; } //----------end of function PndSttTrackFinderReal::WriteMacroParallelHitsGeneralConformalwithMC //----------start of function PndSttTrackFinderReal::WriteMacroParallelHitsGeneralConformalwithMCspecial void PndSttTrackFinderReal::WriteMacroParallelHitsConformalwithMCspecial( Int_t Nhits, // UShort_t iExclude, Double_t auxinfoparalConformal[][5], UShort_t nTracksFoundSoFar, Double_t * ALFA, Double_t * BETA, Double_t * GAMMA, Short_t Status, Double_t *trajectory_vertex ) { Int_t i, j, i1, ii, index, Kincl, nlow, nup, STATUS; Double_t xmin , xmax, ymin, ymax, xl, xu, yl, yu, rrr, cx, cy, Ox[nmaxHits], Oy[nmaxHits], Radius[nmaxHits], alfa, beta, 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) // sprintf(nome,"MacroParallelHitsConformewithMCspecialEvent%dTrack%d", IVOLTE,nTracksFoundSoFar); 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( i== iExclude) continue; // centro sfera in sistema conforme Ox[i] = auxinfoparalConformal[i][0] ; Oy[i] = auxinfoparalConformal[i][1] ; Radius[i] = auxinfoparalConformal[i][2]; 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); 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 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 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, 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 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 fprintf(MACRO,"}\n"); fclose(MACRO); return ; } //----------end of function PndSttTrackFinderReal::WriteMacroParallelAssociatedHitswithMC //----------start of function PndSttTrackFinderReal::WriteMacroSkewAssociatedHits void PndSttTrackFinderReal::WriteMacroSkewAssociatedHits( 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 nMaxima, 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, "MacroParTrack%dSkewTrack%dSkewHitsEvent%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( 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) ); if( distance >= info[i][4] ) continue; 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; // 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 ( 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;++) nohits: ; fprintf(MACRO,"}\n"); fclose(MACRO); } //----------end of function PndSttTrackFinderReal::WriteMacroSkewAssociatedHits //----------start of function PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithMC void PndSttTrackFinderReal::WriteMacroSkewAssociatedHitswithMC( 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 nMaxima, 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, "MacroParTrack%dSkewTrack%dSkewHitswithMCEvent%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( 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) ); if( distance >= info[i][4] ) continue; 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; // 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( 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) ); if( distance >= info[i][4] ) continue; 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; // 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 ( 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 ) 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 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 ( 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;++) 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; } } HitsinBoxConformal[iR][iFi][ nBoxConformal[iR][iFi] ]=(UShort_t) i; nBoxConformal[iR][iFi]++; RConformalIndex[ infoparal[i] ] = iR; FiConformalIndex[ infoparal[i] ] = iFi; } // end of for(i = 0; i< Nparal ; i++) return; } //----------end of function PndSttTrackFinderReal::PndSttBoxConformalFilling //----------begin of function PndSttTrackFinderReal::PndStt_Merge_Sort void PndSttTrackFinderReal::PndStt_Merge_Sort(UShort_t n_ele, Double_t *array, UShort_t *ind) { UShort_t nr, nl, middle, i, ind_left[n_ele], ind_right[n_ele]; Double_t left[n_ele], right[n_ele], result[n_ele]; if( n_ele <= 1) return; middle = n_ele/2 ; for(i=0; i right[0]) { PndStt_Merge(middle, left,ind_left, n_ele-middle, right, ind_right, array, ind); } else { // do the appending for(i=0; i 0 && nr >0){ if( left[nl_curr] <= right[nr_curr]){ result[i] = left[nl_curr]; ind[i] = ind_left[nl_curr]; nl--; nl_curr++; } else { result[i] = right [nr_curr]; ind[i] = ind_right [nr_curr]; nr--; nr_curr++; } i++; } //-------------------- if( nl ==0) { for(j=0; j 0 && i < nHitsinTrack) { nRcell = RConformalIndex[ infoparal[ ListHitsinTrack[i] ] ]; nFicell = FiConformalIndex[ infoparal[ ListHitsinTrack[i] ] ]; //--------------- if (nRcell - NRCELLDISTANCE < 0 ) { nRmin = 0; } else { nRmin = nRcell - NRCELLDISTANCE; } if (nRcell + NRCELLDISTANCE >= nRdivConformalEffective ) { nRmax = nRdivConformalEffective-1; } else { nRmax = nRcell + NRCELLDISTANCE; } for( iR= nRmin ; iR<= nRmax ; iR++){ for( iFi2 = nFicell - NFiCELLDISTANCE ; iFi2<= nFicell + NFiCELLDISTANCE ; iFi2++){ if ( iFi2 < 0 ) { iFi = nFidivConformal + iFi2; } else if ( iFi2 >= nFidivConformal) { iFi = iFi2 - nFidivConformal; } else { iFi = iFi2; } for (j = 0; j< nBoxConformal[iR][iFi]; j++){ if( ExclusionList[ infoparal[ HitsinBoxConformal[iR][iFi][j] ] ] && TemporaryExclusionList[ infoparal[ HitsinBoxConformal[iR][iFi][j] ] ]) { ListHitsinTrack[nHitsinTrack]=HitsinBoxConformal[iR][iFi][j] ; // hit number in the PARALLEL straws scheme //if(IVOLTE==4) cout<<"nFi cell di hit accettato "< 0 && i < nHitsinTrack) // ordering the hits by INCREASING CONFORMAL RADIUS (decreasing space radius) /* for (j = 0; j< nHitsinTrack; j++){ auxIndex[j]=ListHitsinTrack[j]; auxRvalues[j]= info[ infoparal[ ListHitsinTrack[j] ] ][0]* info[ infoparal[ ListHitsinTrack[j] ] ][0]+ info[ infoparal[ ListHitsinTrack[j] ] ][1]* info[ infoparal[ ListHitsinTrack[j] ] ][1]; } PndStt_Merge_Sort( nHitsinTrack, auxRvalues, auxIndex); for (j = 0; j< nHitsinTrack; j++){ ListHitsinTrack[nHitsinTrack-1-j]=auxIndex[j]; } */ return nHitsinTrack; } //----------end of function PndSttTrackFinderReal::PndSttFindTrackPatterninBoxConformal //----------begin of function PndSttTrackFinderReal::PndSttFindTrackPatterninBoxConformalSpecial Short_t PndSttTrackFinderReal::PndSttFindTrackPatterninBoxConformalSpecial( UShort_t NRCELLDISTANCE, UShort_t NFiCELLDISTANCE, UShort_t Nparal, UShort_t NparallelToSearch, UShort_t iSeed, UShort_t *ListHitsinTrackinWhichToSearch, Double_t info[][7], bool ExclusionList[nmaxHits], UShort_t RConformalIndex[nmaxHits], UShort_t FiConformalIndex[nmaxHits], UShort_t nBoxConformal[nRdivConformal][nFidivConformal], UShort_t HitsinBoxConformal[nRdivConformal][nFidivConformal][nmaxHits], UShort_t *OutputListHitsinTrack ) { bool TemporaryExclusionList[nmaxHits]; UShort_t i, i2, j, iR, iFi, nRmin, nRmax, nRemainingHits, nRcell, nFicell, nHitsinTrack, Remaining[nmaxHits], auxIndex[nmaxHits]; Short_t iFi2; Double_t auxRvalues[nmaxHits]; // iSeed is the hit number in the PARALLEL number scheme //cout<<" Nparal "< 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[iR][iFi][j] ] ] && TemporaryExclusionList[ infoparal[ HitsinBoxConformal[iR][iFi][j] ] ]) { OutputListHitsinTrack[nHitsinTrack]=HitsinBoxConformal[iR][iFi][j] ; // hit number in the PARALLEL straws scheme nHitsinTrack++; TemporaryExclusionList[ infoparal[ HitsinBoxConformal[iR][iFi][j] ] ]= 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]=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) cout<<"Stampa dal fitter, prima della correzione displacement *qu, *emme "<<*qu<<", "<<*emme<=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::PndSttTrkAssociatedParallelHitsToHelixBis UShort_t PndSttTrackFinderReal::PndSttTrkAssociatedParallelHitsToHelixBis( Double_t m, Double_t q, Short_t Status, UShort_t nHitsinTrack, UShort_t *ListHitsinTrack, Int_t NhitsParallel, Double_t infoparalConformal[][5], UShort_t *RConformalIndex, UShort_t *FiConformalIndex, UShort_t nBoxConformal[nRdivConformal][nFidivConformal], UShort_t HitsinBoxConformal[nRdivConformal][nFidivConformal][nmaxHits], UShort_t *auxListHitsinTrack ) { bool passamin,passamax; Short_t i, i2, j, k, l, l2, l3, itemp, iFi0, FFimin,FFimax; UShort_t Nextra=8, nFi, Fi, nR, nAssociatedHits; Double_t maxFi, minFi, dist, xx, yy, aaa, angle, r; nAssociatedHits=0; // find the range in Fi spanned by the candidate track FFimin = 10000; FFimax = 0; for(j=0; j FFimax ) FFimax = FiConformalIndex[i]; } if(istampa>=3 && IVOLTE <= nmassimo) { cout<<" evento n. "< 3.*nFidivConformal/4. && FFimin < nFidivConformal/4.) { FFimin = 10000; FFimax = 0; for(j=0; j FFimax ) FFimax = Fi; } } if(istampa>=3 && IVOLTE <= nmassimo) { cout<<" evento n. "< nFidivConformal/2 ) { cout<<"something fishy is going on in PndSttTrkAssociatedParallelHitsToHelixBis!" <<"Range in Fi (rad) is "<<(FFimax - FFimin)*2.*PI/nFidivConformal<=3 && IVOLTE <= nmassimo) { cout<<" evento n. "< 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=nFidivConformal) i2 -= nFidivConformal; for(l=-2; l<3;l++){ l3 = nR+l; if( l3<0 || l3 >= nRdivConformalEffective ) continue; for( k=0;k= nRdivConformalEffective ) continue; for( k=0;k x=0 if( FFimax > nRdivConformal/4 && FFimin < nRdivConformal/4 ) { iFi0 = (Short_t) (nRdivConformal/4 ); } else if ( FFimax > 3*nRdivConformal/4 && FFimin < 3*nRdivConformal/4 ){ iFi0 = (Short_t) (3*nRdivConformal/4 ); } else { cout <<"From PndSttTrackFinderReal::PndSttTrkAssociatedParallelHitsToHelixBis :" <<" 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; l 1.e-10 ) } else if( fabs(q)> 1.e-10) { // second part of if( Status ==99), in this case y = m*x +q passamax=false; passamin=false; for(itemp=FFimin; itemp<=FFimax;itemp++){ i=itemp; if( i< 0 ) { i += nFidivConformal; } else if (i>=nFidivConformal){ i -= nFidivConformal*( i/nFidivConformal ); } angle = (i+0.5)*2.*PI/nFidivConformal; aaa = sin(angle) - m*cos(angle); if( fabs(aaa) < 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=-2; l<3;l++){ l2 = nR+l; if( l2<0 || l2 >= nRdivConformalEffective ) continue; 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= nRdivConformalEffective ) continue; for( k=0;k=nFidivConformal){ i -= nFidivConformal*( i/nFidivConformal ); // i -= nFidivConformal; } for(l=0; l FFimax ) FFimax = FiConformalIndex[i]; } if(istampa>=3 && IVOLTE <= nmassimo) { cout<<" evento n. "< 3.*nFidivConformal/4. && FFimin < nFidivConformal/4.) { FFimin = 10000; FFimax = 0; for(j=0; j FFimax ) FFimax = Fi; } } if(istampa>=3 && IVOLTE <= nmassimo) { cout<<" evento n. "< nFidivConformal/2 ) { cout<<"something fishy is going on in PndSttTrkAssociatedParallelHitsToHelixTris!" <<"Range in Fi (rad) is "<<(FFimax - FFimin)*2.*PI/nFidivConformal<=3 && IVOLTE <= nmassimo) { cout<<" evento n. "< 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=nFidivConformal) i2 -= nFidivConformal; for(l=-2; l<3;l++){ l3 = nR+l; if( l3<0 || l3 >= nRdivConformalEffective ) continue; for( k=0;k= nRdivConformalEffective ) continue; 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::PndSttTrkAssociatedParallelHitsToHelixTris :" <<" 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; l 1.e-10 ) } else if( fabs(q)> 1.e-10) { // second part of if( Status ==99), in this case y = m*x +q Fi0 = atan2(q, -m*q); if(Fi0<0.) { Fi0 += PI; if (Fi0 <0. ) Fi0 =0.; }; ddd= fabs(q)/sqrt(1.+m*m); // cout<<" Fi0 "<=nFidivConformal){ i -= nFidivConformal*( i/nFidivConformal ); } 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.; } // cout<<" iFi = "< Rout ){ continue; } } else if(fabs(erre1) < 1.e-10){ if( Fi0 > fi2 || Fi0 < fi1) continue; } else if ( erre1 Rout && erre2 > Rout && !( fi1<=Fi0 && Fi0<=fi2 && ddd<=Rout ) ) { continue; } // cout<<"accettato ! iFi = "<=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=nFidivConformal){ i -= nFidivConformal*( i/nFidivConformal ); // i -= nFidivConformal; } for(l=0; l 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[l2][i][k] ][0]; dist = fabs( xx +q ); if( PndSttAcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[l2][i][k] ][2], infoparalConformal[ HitsinBoxConformal[l2][i][k] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[l2][i][k]; 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[l3][i2][k] ][0]; dist = fabs( xx +q ); // if(dist < 3.*infoparalConformal[ HitsinBoxConformal[l3][i2][k] ][2]){ if( PndSttAcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[l3][i2][k] ][2], infoparalConformal[ HitsinBoxConformal[l3][i2][k] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[l3][i2][k]; nAssociatedHits++; } } // end of for( k=0;k= nRdivConformalEffective ) continue; for( k=0;k NTIMES*StrawRadius ) continue; //------------------- xx=infoparalConformal[ HitsinBoxConformal[l3][i2][k] ][0]; dist = fabs( xx +q ); // if(dist < 3.*infoparalConformal[ HitsinBoxConformal[l3][i2][k] ][2]){ if( PndSttAcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[l3][i2][k] ][2], infoparalConformal[ HitsinBoxConformal[l3][i2][k] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[l3][i2][k]; 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[l][i][k] ][0]; dist = fabs( xx ); if( PndSttAcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[l][i][k] ][2], infoparalConformal[ HitsinBoxConformal[l][i][k] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[l][i][k]; 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 = atan2(q, -m*q); if(Fi0<0.) { Fi0 += PI; if (Fi0 <0. ) Fi0 =0.; }; 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 ); } 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 Rout && erre2 > Rout && !( fi1<=Fi0 && Fi0<=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[k][l2][l3] ][0]; yy=infoparalConformal[ HitsinBoxConformal[k][l2][l3] ][1]; dist = fabs( -yy+ m*xx +q )/sqrt(m*m+1.); if( PndSttAcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[k][l2][l3] ][2], infoparalConformal[ HitsinBoxConformal[k][l2][l3] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[k][l2][l3]; Unselected[HitsinBoxConformal[k][l2][l3]]= 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[l][i][k] ][0]; yy=infoparalConformal[ HitsinBoxConformal[l][i][k] ][1]; dist = fabs( m*xx-yy )/sqrt( m*m+1.); // if(dist < 3.*infoparalConformal[ HitsinBoxConformal[l][i][k] ][2]){ if( PndSttAcceptHitsConformal( dist, infoparalConformal[ HitsinBoxConformal[l][i][k] ][2], infoparalConformal[ HitsinBoxConformal[l][i][k] ][4] ) ) { auxListHitsinTrack[nAssociatedHits]= HitsinBoxConformal[l][i][k]; nAssociatedHits++; } } // end of for( k=0;k=3) { cout<<"from PndSttTrkAssociatedParallelHitsToHelixQuater, before exiting; nAssociatedHits = " <= info[i][4] ) continue; 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; // 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; } // 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( istampa>=3 && IVOLTE <= nmassimo){ cout<<" Zlast1 "< 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; Double_t allowed_distance = 4.*StrawRadius/sin(3.*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=3 && IVOLTE <= nmassimo) cout<<"From AssociateBetterAfterFitSkewHitsToXYTrack , hit skew n. "<=3 && IVOLTE <= nmassimo) cout<<"From AssociateBetterAfterFitSkewHitsToXYTrack : segno "<<-1+sign*2< zmax ){ tempZ[sign]=fmod( tempZ[sign]-zmax, deltaz) - deltaz; } else if (tempZ[sign]=3 && IVOLTE <= nmassimo) cout<<"From AssociateBetterAfterFitSkewHitsToXYTrack , hit skew n. "<=3 && IVOLTE <= nmassimo) cout<<"From AssociateBetterAfterFitSkewHitsToXYTrack , hit skew n. "<< SkewList[i][0] <<", tempZ "< bbb ) { tempore[NAssociated]=SkewList[i][0]; temporeS[NAssociated]=S[i]; temporeZ[NAssociated]=Z[i]; temporeZDrift[NAssociated]=ZDrift[i]; temporeZErrorafterTilt[NAssociated]=ZErrorafterTilt[i]; NAssociated++; if(istampa==2 && IVOLTE <= nmassimo) cout<<"From AssociateBetterAfterFitSkewHitsToXYTrack , hit skew n. "< 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 { // in this case the equation is 0 = x+*qu . 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]; } PndStt_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]; } } PndStt_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 ) { UShort_t i,j,flag; Double_t old, auxFivalues[nmaxHits], auxRvalues[nmaxHits]; // here there is the ordering of the hits // ordering of the parallel hits 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]; } PndStt_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; } } // finding the charge of the track if( auxFivalues[0] > auxFivalues[nParallelHits-1]) { *Charge = 1; } else { *Charge = -1; } // FI initial value (at 0,0 vertex) in the Helix reference frame *Fi_initial_helix_referenceframe = atan2(oY,oX) + PI; // 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 = auxFivalues[nParallelHits-1]; return; } //----------end of function PndSttTrackFinderReal::PndSttOrderingParallel //----------begin of function PndSttTrackFinderReal::PndSttOrdering void PndSttTrackFinderReal::PndSttOrdering( Double_t oX, Double_t oY, 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 *nTotal, UShort_t *BigList, Short_t * Charge ) { *nTotal = nParallelHits+nSkewHits; UShort_t i,j,flag, aux[*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]; } PndStt_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; } } // finding the charge of the track 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] ]]; } } PndStt_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= "<= 0. ){ // this is the most common case if( Charge < 0.){ teta1=asin(0.5*RStrawDetectorMin/R); teta2=asin(0.5*RStrawDetectorMax/R); tetavertex=atan2( -oX, oY); teta1 +=tetavertex; teta2 +=tetavertex; } else { teta2=asin(0.5*RStrawDetectorMin/R); teta1=asin(0.5*RStrawDetectorMax/R); tetavertex=atan2( oX, -oY); teta1 = tetavertex - teta1; teta2 = tetavertex - teta2; } } else if ( R + tmp1 > RStrawDetectorMin ){ if( Charge < 0.){ teta1=asin(0.5*RStrawDetectorMin/R); teta2= PI - teta1; tetavertex=atan2( -oX, oY); teta1 +=tetavertex; teta2 +=tetavertex; } else { teta2=asin(0.5*RStrawDetectorMin/R); teta1= PI - teta2; tetavertex=atan2( oX, -oY); teta1 = tetavertex - teta1; teta2 = tetavertex - teta2; } } else { // case when the trajectory is too small teta1=-99999.; } // end of if(R + tmp1 - RStrawDetectorMax >= 0. ) if(teta1>-99998) { // add safety margin teta1 -= 2.*StrawRadius/RStrawDetectorMin; teta2 += 2.*StrawRadius/RStrawDetectorMin; if(teta1<0.) { teta1 += 2.*PI; teta2 += 2.*PI; } //------- if(teta1 > 2.*PI) { teta1=fmod(teta1,2.*PI); teta2=fmod(teta2,2.*PI); } *Fi_low_limit=teta1; *Fi_up_limit=teta2; } // end of if(teta1>-99998) //------------ end calculation the maximum fi and minimum fi spanned by this track return; } //---------- end of function PndSttTrackFinderReal::PndSttFindingParallelTrackAngularRange //----------start function PndSttTrackFinderReal::PndSttFindingAllowedAngularRangeforSkew void PndSttTrackFinderReal::PndSttFindingAllowedAngularRangeforSkew( Double_t oX, Double_t oY, Double_t R, Short_t Charge, Double_t *Fi_allowedforskew_low, Double_t *Fi_allowedforskew_up ) { UShort_t i, index, is; Double_t alpha, beta, delta, fi_vertex, c1, c2, maximum, maximum2, minimum, minimum2, x1, x2, y, fi[12], a[] = { 0., -sqrt(3.), sqrt(3.) , 0. , -sqrt(3.), sqrt(3.) } , b[] = { 1., 2. , -2., -1., -2. , 2. }, Erre[] = {RminStrawSkewArea, RmaxStrawSkewArea }, // this is the distance from (0,0) of the side delimiting the Skew area exhagon_side_xlow[6][2] , exhagon_side_xup[6][2] ; // cout<<" Rmin = "<= R) continue; // skew hits not possible. for(i=0;i<2;i++){ c1 = alpha+ 2.*a[is]*b[is]*Erre[i]+ a[is]*beta; delta = c1*c1 - 4.*c2*b[is]*Erre[i]*(b[is]*Erre[i]+beta); if(delta<0.) delta=0.; // delta<0 could never happen in principle, only caused by rounding x1 = 0.5*(-c1 - sqrt(delta))/c2; x2 = 0.5*(-c1 + sqrt(delta))/c2; // cout<<" is = "< 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 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 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) ); if( distance >= info[i][4] ) continue; 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; // 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 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; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if( !inclusionExp[jexp]) continue; massimo = -1; 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( 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 inclusionMC[nTracksFoundSoFar][nmaxHits], inclusionExp[nTracksFoundSoFar]; UShort_t ntoMCtrack[nTracksFoundSoFar], toMCtracklist[nTracksFoundSoFar][nmaxHits], toMCtrackfrequency[nTracksFoundSoFar][nmaxHits]; UShort_t i, j, enne, jtemp,jexp; Short_t itemp, massimo; for(i=0; i -1){ itemp=-1; massimo = -1; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if( !inclusionExp[jexp]) continue; for(i=0; i< ntoMCtrack[jexp]; i++){ if( !inclusionMC[jexp][i]) continue; if( toMCtrackfrequency[jexp][i]>massimo){ massimo=toMCtrackfrequency[jexp][i]; itemp = toMCtracklist[jexp][i]; jtemp = jexp; } } } if( itemp>-1 ){ daTrackFoundaTrackMC[jtemp]=itemp; inclusionExp[jtemp]=false; for(jexp=0; jexp -1) return; } //----------end of function 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; } Posiz[0] = info[infopar][0] + info[infopar][3]*vers[0]/norm; Posiz[1] = info[infopar][1] + info[infopar][3]*vers[1]/norm; if( fabs(KAPPA)<1.e-20 ){ Posiz[2] = -888888888.; return; } fi = atan2(-vers[1],-vers[0]); if(fi<0.) fi += 2.*PI; if ( Charge > 0){ if(fi > FI0 ) FI0 += 2.*PI; Posiz[2] = (FI0-fi)/KAPPA; } else { if(fi < FI0 ) fi += 2.*PI; Posiz[2] = (fi-FI0)/KAPPA; } return; } //----------end of function PndSttTrackFinderReal::PndSttInfoXYZParal //---------- begin of function PndSttTrackFinderReal::PndSttInfoXYZSkew void PndSttTrackFinderReal::PndSttInfoXYZSkew ( // Double_t info[][7], // UShort_t infosk, Double_t Z, // Z coordinate 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 ) { //if(istampa >=3 && IVOLTE <= nmassimo){ if(istampa>=3){ cout<<" stampa da PndSttInfoXYZSkew "<<", Z hit = "<