// ------------------------------------------------------------------------- // ----- 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_t verbose) { fVerbose = 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]; // --------------------------- opening files for special purposes /* 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"); // --------------- 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) */ // ------------------------- 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!" << 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) { UShort_t auxIndex[nmaxHits], OLDinfoparal[nmaxHits]; UShort_t istep,inclination_type; Double_t aaa, ddd, delta, deltabis, deltaZ, mindis, distanza, fi_hit, ap1, ap2, ap3, cross1, cross2, cross3, lowlimit[MAXMCTRACKS], uplimit[MAXMCTRACKS], 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<20) cout<<"da PndSttTrackFinderReal : evento n. "<GetEntriesFast(); // num. tracce/evento if(istampa>=2 && (nMCTracks != N_INTENDED) ){ cout<<"Gianluigi : n. MC tracks = "<MAXMCTRACKS){ return -10; } for (Int_t iMCTrack = 0; iMCTrack < nMCTracks; iMCTrack++) { nHitsInMCTrack[iMCTrack] = 0; // initialization nSkewHitsInMCTrack[iMCTrack] = 0; // initialization pMCtr = (PndMCTrack*) fMCTrackArray->At(iMCTrack); if ( ! pMCtr ) continue; Double_t R, D, Fi, Ox, Oy, Cx, Cy, Px, Py ; Int_t icode; icode = pMCtr->GetPdgCode() ; // PDG code of track Ox = pMCtr->GetStartVertex().X(); // X of starting point track Oy = pMCtr->GetStartVertex().Y(); // Y of starting point track Px = pMCtr->GetMomentum().X(); Py = pMCtr->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); R = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla if( icode > 0 ) { Cx = Ox -Py*1000./(BFIELD*CVEL); // MC truth X of center of circle of Helix trajectory Cy = Oy + Px*1000./(BFIELD*CVEL); // MC truth Y of center of circle of Helix trajectory } else { Cx = Ox + Py*1000./(BFIELD*CVEL); // MC truth X of center of circle of Helix trajectory Cy = Oy - Px*1000./(BFIELD*CVEL); // MC truth Y of center of circle of Helix trajectory } Fi = atan2(Cy, Cx); // MC truth Fi angle of circle of Helix trajectory if(Fi<0.) Fi += 2.*PI; D = sqrt( Cx*Cx+Cy*Cy) - R; if(istampa>=3){ Double_t gomma = Cx*Cx+Cy*Cy -R*R; cout<<"verita MC, traccia n. "<GetStartVertex().Z(); // Z of starting point track MCtruthTrkInfo[3][iMCTrack] = Px; // Px at starting point of track MCtruthTrkInfo[4][iMCTrack] = Py; // Py at starting point of track MCtruthTrkInfo[5][iMCTrack] = pMCtr->GetMomentum().Z(); // Pz at starting point of track MCtruthTrkInfo[6][iMCTrack] = D; // D of Helix of track projected in XY plane MCtruthTrkInfo[7][iMCTrack] = Fi; // Fi of Helix of track projected in XY plane MCtruthTrkInfo[8][iMCTrack] = R; // R of Helix of track projected in XY plane MCtruthTrkInfo[9][iMCTrack] = Cx; // Cx of Helix of track projected in XY plane MCtruthTrkInfo[10][iMCTrack] = Cy; // Cy of Helix of track projected in XY plane MCtruthTrkInfo[11][iMCTrack] = icode ; // PDG code of track TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); MCtruthTrkInfo[14][iMCTrack] = fParticle->Charge()/3. ; // charge of track MCtruthTrkInfo[12][iMCTrack] = -MCtruthTrkInfo[14][iMCTrack]*0.001*BFIELD*CVEL/ MCtruthTrkInfo[5][iMCTrack] ; // KAPPA of Helix of track in cm*-1 // fabs(MCtruthTrkInfo[5][iMCTrack]) ; // KAPPA of Helix of track in cm*-1 // double tempoang = atan2(-MCtruthTrkInfo[10][iMCTrack]+MCtruthTrkInfo[1][iMCTrack], // -MCtruthTrkInfo[9][iMCTrack]+MCtruthTrkInfo[0][iMCTrack]); // if(tempoang<0.) tempoang+= 2.*PI; MCtruthTrkInfo[13][iMCTrack] = fmod(Fi+ PI, 2.*PI) ; // FI0 of Helix of track } // end of for (Int_t iMCTrack = 0; iMCTrack < nMCTracks; iMCTrack++) // Initialise control counters Int_t nNoTrack = 0; Int_t nNoSttPoint = 0; Int_t nNoSttHit = 0; // Create pointers to hit and SttPoint PndSttHit* pMhit = NULL; FairMCPoint* pMCpt = NULL; // PndSttTrack* pTrck[MAXTRACKSPEREVENT] ; // this is for the hits found by pattern recognition // Number of STT hits Int_t Nhits = 0; for (Int_t hitListCounter = 0; hitListCounter < fHitCollectionList.GetEntries(); hitListCounter++) { Nhits += ((TClonesArray *)fHitCollectionList.At(hitListCounter))->GetEntriesFast(); } if(Nhits > nmaxHits ) { return -10; } // Declare some variables outside the loops Int_t trackIndex = 0; // STTTrack index // Create STL map from MCtrack index to number of valid SttHits map > hitMap; for(Int_t j=0; j= 1 ) cout<<"Gianluigi : da PndSttTrackFinderReal::DoFind : Nhits="<GetRefIndex(); if (ptIndex < 0) continue; // fake or background hit pMCpt = GetPointFromCollections(iHit); // <== FairMCPoint if (!pMCpt) 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<20) { 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<20) { if (istampa >= 2 ) { cout<<"Gianluigi : da PndSttTrackFinderReal::DoFind : Nhits totali ="<= 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=0 ){ if (NNN < MINIMUMHITSPERTRACK) { continue; } GoodSkewFit[i]= true; nSkewHitsinTrack[i] = NNN; for(j=0;j=0 ) // --------- only for the comparison with the MC truth; numbering according to the ORIGINAL hit number for(i1=0; i1< nSkewHitsinTrack[i]; i1++){ Sfinal[i][ infoskew[ ListSkewHitsinTrack[i][i1] ] ] = S[i1] ; Zfinal[i][ infoskew[ ListSkewHitsinTrack[i][i1] ] ] = Z[i1] ; ZDriftfinal[i][ infoskew[ ListSkewHitsinTrack[i][i1] ] ] = ZDrift[i1] ; ZErrorafterTiltfinal[i][ infoskew[ ListSkewHitsinTrack[i][i1] ] ] = ZErrorafterTilt[i1] ; } // --------- end of only for the comparison with the MC truth // ------------------------------------------------------------- // ordering the parallel and skew hits; determining the charge of this track PndSttOrdering( Ox[i], Oy[i], info, nHitsinTrack[i], &ListHitsinTrack[i][0], nSkewHitsinTrack[i], &ListSkewHitsinTrack[i][0], S, infoparal, infoskew, &nTotalHits[i], &BigList[i][0], &Charge[i] ); //---------------------------------------------- some printouts if(istampa>=2 && IVOLTE<20) { cout<<" Evento n. "<0 && nTracksFoundSoFar > 0 ){ AssociateFoundTrackstoMC( nTracksFoundSoFar, nHitsinTrack, ListHitsinTrack, nSkewHitsinTrack, ListSkewHitsinTrack, daTrackFoundaTrackMC, daMCTrackaTrackFound ); } //--------------------------------------------------------------------------------------------------------- //-------------------- inizio della sezione sul confronto tra MC truth e tracce trovate // prima gli hits paralleli --------------------- for(jexp=0; jexp=2 ) fprintf(HANDLE, "\n Evento %d NTotaleTracceMC %d ------\n",IVOLTE, nMCTracks); for (i=0; i=2 ;i++){ fprintf(HANDLE,"----------------------------------------------------------\n"); ii=daMCTrackaTrackFound[i]; if( ii <0 ) { fprintf(HANDLE," TracciaMC %d has not a corresponding found track in pattern recognition\n",i); continue; } if( ! ( nHitsInMCTrack[i]>= MINIMUMHITSPERTRACK && nSkewHitsInMCTrack[i] >= MINIMUMSKEWHITSPERTRACK && R_MC[i] > RStrawDetectorMin ) ) { fprintf(HANDLE," TracciaMC %d NONsoddisfaRequisitiMinimi with %d Hits paral. in MC track, %d skew hits in MC track, %f Radius in MC track ",i, nHitsInMCTrack[i], nSkewHitsInMCTrack[i], R_MC[i] ); fprintf(HANDLE," (MINIMUMHITSPERTRACK = %d, MINIMUMSKEWHITSPERTRACK=%d, minRadius=%f\n", MINIMUMHITSPERTRACK, MINIMUMSKEWHITSPERTRACK, RStrawDetectorMin ); 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 dista solo %g da (0,0)\n", i,ii,dista); continue; } fprintf(HANDLE, " TracciaMC %d ParHitsMC %d ParMecc %d ParMeccSpuri %d SkewHitsMC %d SkewMecc %d SkewMeccSpuri %d\n", i, nHitsInMCTrack[i], nParalCommon[daMCTrackaTrackFound[i]], nSpuriParinTrack[daMCTrackaTrackFound[i]], nSkewHitsInMCTrack[i], nSkewCommon[daMCTrackaTrackFound[i]], nSpuriSkewinTrack[daMCTrackaTrackFound[i]] ) ; fprintf(HANDLE, " e corrisponde a track found n. %d\n", daMCTrackaTrackFound[i] ); fprintf(HANDLE, " AVENDO %d hits paralleli e %d hits skew non mecciati dalla corrisponde track found\n", nMCParalAlone[i],nMCSkewAlone[i] ); HoughFi = atan2(Oy[daMCTrackaTrackFound[i]],Ox[daMCTrackaTrackFound[i]]); if(HoughFi<0.) HoughFi += 2.*PI; fprintf(HANDLE, " R_MC %g R %g Fi_MC %g Fi %g KAPPA_MC %g KAPPA %g FI0_MC %g FI0 %g\n", MCtruthTrkInfo[8][i], R[ daMCTrackaTrackFound[i] ], MCtruthTrkInfo[7][i], HoughFi, MCtruthTrkInfo[12][i], KAPPA[ daMCTrackaTrackFound[i] ], MCtruthTrkInfo[13][i], FI0[ daMCTrackaTrackFound[i] ] ); } // end of for (i=0; i=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 associates to a track in a PndTrackCand class; a class per each track for(i=0; iAt(ipinco); // PndTrackCand *pTrckCand = (PndTrackCand*) (*trackArray)[ipinco]; ipinco++; Double_t inv=R[i]*KAPPA[i]; TVector3 dirSeed(Pxini, Pyini, Pzini); // momentum direction in starting point Double_t qop = Charge[i]/dirSeed.Mag(); dirSeed.SetMag(1.); TVector3 posSeed(0., 0., 0.); // direction in starting point pTrckCand->setTrackSeed(posSeed, dirSeed, qop); pTrckCand->setMcTrackId( daTrackFoundaTrackMC[i] ); for(j=0; j< nTotalHits[i]; j++){ pTrckCand->AddHit(kSttHit, (Int_t) BigList[i][j] , j); } } // end of for(i=0; i=3) cout<<"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 ) { for(i=0; i=2){ fprintf(HANDLE2,"Evento n. %d Found track %d messa in PndTrackCand",IVOLTE, i); fprintf(HANDLE2, " R_MC %g R %g Fi_MC %g Fi %g KAPPA_MC %g KAPPA %g FI0_MC %g FI0 %g\n", MCtruthTrkInfo[8][ii], R[ i ], MCtruthTrkInfo[7][ii], HoughFi, MCtruthTrkInfo[12][ii], KAPPA[ i ], MCtruthTrkInfo[13][ii], // FI0[ i ] fmod(HoughFi+ PI, 2.*PI) ); } // end of if (istampa>=2) for( j=0; j=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 < 20){ for(i=0; i= 0){ WriteMacroSkewAssociatedHitswithMC( KAPPA[i],FI0[i], HoughD, HoughFi, HoughR, info, Nincl,Minclinations,inclination, i,0, nSkewHitsinTrack[i], ListSkewHitsinTrack, nSkewCommon[i], SkewCommonList, daTrackFoundaTrackMC[i], nMCSkewAlone[ daTrackFoundaTrackMC[i] ], 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<20){ 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<20) { 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<20) { cout<<"da findmaximaDFiR : MINIMUMCOUNTS = "<= "<=3 && IVOLTE<20) { 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<20) { 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); } } for( i=0; i< nMCTracks ; i++) { 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,CxMC[i],CyMC[i],R_MC[i],R_MC[i],i,i,i); // cout<<" info MC buone : MCtrack n. "<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){ 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 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(R_MC[i]/gamma) > 1.) { if(fabs(CyMC[i])>0.001 ) { yl = -xmin*CxMC[i]/CyMC[i]+0.5/CyMC[i]; yu = -xmax*CxMC[i]/CyMC[i]+0.5/CyMC[i]; 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,CxMC[i]/gamma,CyMC[i]/gamma,R_MC[i]/fabs(gamma),R_MC[i]/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 "<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(R_MC[i]/gamma) > 1.) { if(fabs(CyMC[i])>0.001 ) { yl = -xmin*CxMC[i]/CyMC[i]+0.5/CyMC[i]; yu = -xmax*CxMC[i]/CyMC[i]+0.5/CyMC[i]; 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,R_MC[i]/fabs(gamma),R_MC[i]/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::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, UShort_t MCSkewAloneList[MAXMCTRACKS][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%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( 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; imc= 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 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; } 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<20) { cout<<"Da strictcollection, n. elementi delle search list "<=3 && IVOLTE<20)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<20) { 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<20){ 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<20) { 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<20){ 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<20) { 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<20) { cout<<" evento n. "< 3.*nFidivConformal/4. && FFimin < nFidivConformal/4.) { FFimin = 10000; FFimax = 0; for(j=0; j FFimax ) FFimax = Fi; } } if(istampa>=3 && IVOLTE<20) { 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<20) { 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<20) { cout<<" evento n. "< 3.*nFidivConformal/4. && FFimin < nFidivConformal/4.) { FFimin = 10000; FFimax = 0; for(j=0; j FFimax ) FFimax = Fi; } } if(istampa>=3 && IVOLTE<20) { 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<20) { 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(istampa>=3 && IVOLTE<20) { cout<<" evento n. "< 3.*nFidivConformal/4. && FFimin < nFidivConformal/4.) { FFimin = 10000; FFimax = 0; for(j=0; j FFimax ) FFimax = Fi; } } if(istampa>=3 && IVOLTE<20) { cout<<" evento n. "< nFidivConformal/2 ) { cout<<"something fishy is going on in PndSttTrkAssociatedParallelHitsToHelixQuater!" <<"Range in Fi (rad) is "<<(FFimax - FFimin)*2.*PI/nFidivConformal<=3 && IVOLTE<20) { 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 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(l=0; l<3;l++) // ------- special cases if((nR == nRdivConformalEffective-1 && passamin && ! passamax) || (nR==0 && passamax && !passamin) ) { // do the last two Fi columns if(nR == nRdivConformalEffective-1) passamax=true; if(nR == 0) passamin=true; for(l2=1;l2<3;l2++){ i2 = i+l2; if(i2>=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(l=-2; l<3;l++) } // end of for(l2=0;l2<2;l2++) return nAssociatedHits; } else if ((nR == nRdivConformalEffective-1 && ! passamin && !passamax) || (nR==0 && !passamax && !passamin)){ if(nR == nRdivConformalEffective-1) passamax=true; if(nR == 0) passamin=true; for(l2=1;l2<3;l2++){ i2 = i-l2; if(i2= 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(l=-2; l<3;l++) } // end of for(l2=0;l2<2;l2++) } // end of if((nR == nRdivConformalEffective-1 && passamin) || (nR==0 && passamax) ) } // end of for(itemp=Fimin; itemp<=FFimax;itemp++) } else { // q=0 --> 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<20){ 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; if(istampa>=3 && IVOLTE<20) cout<<" Zmax "<=3 && IVOLTE<20) cout<<"From AssociateBetterAfterFitSkewHitsToXYTrack , hit skew n. "<=3 && IVOLTE<20) cout<<"From AssociateBetterAfterFitSkewHitsToXYTrack : segno "<<-1+sign*2< zmax ){ tempZ[sign]=fmod( tempZ[sign]-zmax, deltaz) - deltaz; } else if (tempZ[sign]=3 && IVOLTE<20) cout<<"From AssociateBetterAfterFitSkewHitsToXYTrack , hit skew n. "<=3 && IVOLTE<20) 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<20) 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( 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[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++) { 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,CxMC[i],CyMC[i],R_MC[i],R_MC[i],i,i,i); // cout<<" info MC buone : MCtrack n. "<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; imc= 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= 3 && IVOLTE<20) cout<<" evento n. "<