#include "glpk.h" #include "PndSttMvdTracking.h" #include "PndSttHit.h" #include "PndSttTrack.h" #include "PndSttPoint.h" #include "PndSttHelixHit.h" #include "PndSttSingleStraw.h" #include "PndSttTube.h" #include "PndSttMapCreator.h" #include "PndSdsHit.h" #include "PndSdsMCPoint.h" #include "PndTrackCand.h" #include "PndTrackCandHit.h" #include "PndTrack.h" #include "FairRootManager.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "FairTrackParP.h" #include "TGeoManager.h" #include "TClonesArray.h" #include "TGeoVolume.h" #include "TVector3.h" #include "TRandom.h" #include "TH1F.h" #include "TMath.h" #include "TCanvas.h" #include "TGeoTube.h" #include #include using namespace std; // ----- Default constructor ------------------------------------------- PndSttMvdTracking::PndSttMvdTracking() : FairTask("STT Stt-Mvd Tracking") { fPersistence = kTRUE; fVerbose = 0; istampa = 0; iplotta = false; doMcComparison = false; sprintf(fSttBranch,"STTHit"); sprintf(fMvdPixelBranch,"MVDHitsPixel"); sprintf(fMvdStripBranch,"MVDHitsStrip"); } // ------------------------------------------------------------------------- PndSttMvdTracking::PndSttMvdTracking(Int_t verbose) : FairTask("STT Stt-Mvd Tracking") { fPersistence = kTRUE; fVerbose = verbose; istampa = verbose; iplotta = false; doMcComparison = false; sprintf(fSttBranch,"STTHit"); sprintf(fMvdPixelBranch,"MVDHitsPixel"); sprintf(fMvdStripBranch,"MVDHitsStrip"); } // ------------------------------------------------------------------------- PndSttMvdTracking::PndSttMvdTracking(int istamp, bool iplot, bool imc) : FairTask("STT Stt-Mvd Tracking") { fPersistence = kTRUE; istampa = istamp; iplotta = iplot; doMcComparison = imc; sprintf(fSttBranch,"STTHit"); sprintf(fMvdPixelBranch,"MVDHitsPixel"); sprintf(fMvdStripBranch,"MVDHitsStrip"); } // ----- Destructor ---------------------------------------------------- PndSttMvdTracking::~PndSttMvdTracking() { } // ------------------------------------------------------------------------- // ----- Public method Init -------------------------------------------- InitStatus PndSttMvdTracking::Init() { ZCENTER_STRAIGHT = 35.; SEMILENGTH_STRAIGHT = 75.; /* hx = new TH1F("hx", "x: mc - reco", 100, -1, 1); hzresvsslope = new TH2F("hzresvsslope", "z: mc - reco vs slope", 100, -3.5, 3.5, 100, -3., 3.); */ if(iplotta){ hdeltaRPixel = new TH1F("hdeltaRPixel", "distance MC Pixel point from trajectory in XY plane", 100, -1, 1); hdeltaRStrip = new TH1F("hdeltaRStrip", "distance MC Strip point from trajectory in XY plane", 100, -1, 1); hdeltaRPixel2 = new TH1F( "hdeltaRPixel2", "distance MC point from trajectory in XY plane (Pixels)", 100, -10, 10); hdeltaRStrip2 = new TH1F( "hdeltaRStrip2", "distance MC point from trajectory in XY plane (Strips)", 100, -10,10); } IVOLTE=-1; // --------------------------- opening files for special purposes int N_INTENDED=0; if(istampa >=1 ){ //---- fetch the n. of tracks MC that were intended to be generated // HANDLE = fopen("n_intended_tracks.txt","r"); // fscanf(HANDLE,"%d",&N_INTENDED); // fclose(HANDLE); //---- apertura file con info su Found tracce su cui si fa Helix fit dopo HANDLE2 = fopen("info_da_PndTrackFinderReal.txt","w"); // ---- open filehandle per statistica sugli hits etc. HANDLE = fopen("statistiche.txt","w"); // --------------- // HANDLEXYZ = fopen("infoPndTrackFinderRealXYZ.txt","w"); // --------------- open file delle info su deltaX, Y, Z degli hits in comune tra tracce trovate e MC PHANDLEX = fopen("deltaParXmio.txt","w"); PHANDLEY = fopen("deltaParYmio.txt","w"); PHANDLEZ = fopen("deltaParZmio.txt","w"); SHANDLEX = fopen("deltaSkewXmio.txt","w"); SHANDLEY = fopen("deltaSkewYmio.txt","w"); SHANDLEZ = fopen("deltaSkewZmio.txt","w"); } // end of if(istampa >=1) // ------------------------- // Get RootManager FairRootManager* ioman = FairRootManager::Instance(); if ( ! ioman ) { cout << "-E- PndSttMvdTracking::Init: " << "RootManager not instantiated, return!" << endl; return kFATAL; } // ----- maps of STT tubes // CHECK added PndSttMapCreator *mapper = new PndSttMapCreator(fSttParameters); fSttTubeArray = mapper->FillTubeArray(); //---------------------------------------------------- end map // get the MCTrack array fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack"); if ( ! fMCTrackArray) { cout << "-E- PndSttMvdTracking::Init: No MCTrack array, return!" << endl; return kERROR; } // Get SttTrackCand array dal pattern recognition di STT fSttTrackCandArray = (TClonesArray*) ioman->GetObject("STTTrackCand"); if ( ! fSttTrackCandArray) { cout << "-E- PndSttMvdTracking::Init: No SttTrack Cand array, return!" << endl; return kERROR; } // Get input array questi sono i MC point di STT fSttPointArray = (TClonesArray*) ioman->GetObject("STTPoint"); if ( ! fSttPointArray ) { cout << "-W- PndSttHelixHitProducer::Init: " << "No STTPoint array, return!" << endl; return kERROR; } // Get input array hit di STT after digi fSttHitArray = (TClonesArray*) ioman->GetObject(fSttBranch); // fSttHitArray = (TClonesArray*) ioman->GetObject("STTHit"); if ( ! fSttHitArray ) { cout << "-W- PndSttMvdTracking::Init: " << "No STTHit array, return!" << endl; return kERROR; } // Create and register output array for PndTrackCand of Stt+Mvd combined fSttMvdPndTrackCandArray = new TClonesArray("PndTrackCand"); ioman->Register("SttMvdTrackCand","SttMvd",fSttMvdPndTrackCandArray, kTRUE); // Create and register output array for PndTrack of Stt+Mvd combined fSttMvdPndTrackArray = new TClonesArray("PndTrack"); ioman->Register("SttMvdTrack","SttMvd",fSttMvdPndTrackArray, kTRUE); // ------------------------- get the Mvd hits fMvdPixelHitArray = (TClonesArray*) ioman->GetObject(fMvdPixelBranch); // fMvdPixelHitArray = (TClonesArray*) ioman->GetObject("MVDHitsPixel"); if ( !fMvdPixelHitArray){ std::cout << "-W- PndSttMvdTracking::Init: " << "No MVD Pixel hitArray, return!" << std::endl; return kERROR; } fMvdStripHitArray = (TClonesArray*) ioman->GetObject(fMvdStripBranch); // fMvdStripHitArray = (TClonesArray*) ioman->GetObject("MVDHitsStrip"); if ( !fMvdStripHitArray){ std::cout << "-W- PndSttMvdTracking::Init: " << "No MVD Strip hitArray, return!" << std::endl; return kERROR; } // ------------------------- get the Mvd track candidates fMvdTrackCandArray = (TClonesArray*) ioman->GetObject("MVDRiemannTrackCand"); if ( !fMvdTrackCandArray){ std::cout << "-W- PndSttMvdTracking::Init: " << "No MVD TrackCand Array, return!" << std::endl; return kERROR; } cout << "-I- PndSttMvdTracking: Initialization successfull" << endl; // ------------------------- get the Mvd MC points fMvdMCPointArray = (TClonesArray*) ioman->GetObject("MVDPoint"); if ( !fMvdMCPointArray){ std::cout << "-W- PndSttMvdTracking::Init: " << "No MVD MC Point Array, return!" << std::endl; return kERROR; } cout << "-I- PndSttMvdTracking: Initialization successfull" << endl; return kSUCCESS; } // ------------------------------------------------------------------------- // CHECK added void PndSttMvdTracking::SetParContainers() { FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb(); fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar"); } void PndSttMvdTracking::WriteHistograms(){ /* TFile* file = FairRootManager::Instance()->GetOutFile(); file->cd(); file->mkdir("PndSttHelixHit"); file->cd("PndSttHelixHit"); hx->Write(); delete hx; */ TFile* file = FairRootManager::Instance()->GetOutFile(); file->cd(); file->mkdir("PndSttMvdTracking"); file->cd("PndSttMvdTracking"); hdeltaRPixel->Write(); hdeltaRStrip->Write(); hdeltaRPixel2->Write(); hdeltaRStrip2->Write(); delete hdeltaRPixel; delete hdeltaRStrip; delete hdeltaRPixel2; delete hdeltaRStrip2; } // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- // ----- Public method Exec -------------------------------------------- void PndSttMvdTracking::Exec(Option_t* opt) { bool Status, keepit[MAXTRACKSPEREVENT], status[MAXTRACKSPEREVENT], SttSZfit[MAXTRACKSPEREVENT], Mvdhits[MAXTRACKSPEREVENT]; Short_t Candidato, Charge, statusflag[MAXTRACKSPEREVENT], iflag, iHit, i, j, k, resultFitSZagain[MAXTRACKSPEREVENT], tubeID ; Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT], CHARGE[MAXTRACKSPEREVENT]; UShort_t l, kall; UShort_t FromHitToMCTrack[nmaxSttHits], nMCParalAlone[MAXTRACKSPEREVENT], nMCSkewAlone[MAXTRACKSPEREVENT], MCParalAloneList[MAXTRACKSPEREVENT][nmaxSttHits], MCSkewAloneList[MAXTRACKSPEREVENT][nmaxSttHits], nHitsInMCTrack[MAXTRACKSPEREVENT], nSkewHitsInMCTrack[MAXTRACKSPEREVENT], nParalCommon[MAXTRACKSPEREVENT], ParalCommonList[MAXTRACKSPEREVENT][nmaxSttHits], nSpuriParinTrack[MAXTRACKSPEREVENT], ParSpuriList[MAXTRACKSPEREVENT][nmaxSttHits], nSkewCommon[MAXTRACKSPEREVENT], SkewCommonList[MAXTRACKSPEREVENT][nmaxSttHits], nSpuriSkewinTrack[MAXTRACKSPEREVENT], SkewSpuriList[MAXTRACKSPEREVENT][nmaxSttHits], TemporarySkewList[2*nmaxSttHits][2]; UShort_t nTotalCandidates, nRemainingCandidates, nMvdMCPoint, nSttHit, nSttParHit, nSttSkewHit, nSttHelixTrack, nSttMCPoint, nSttTrackCand, TemporarynSkewHitsinTrack, ListAllParHits[nmaxSttHits], ListAllSkewHits[nmaxSttHits], nSttHitsinTrack[MAXTRACKSPEREVENT], ListSttHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], ListSttHitsinTrackType[MAXTRACKSPEREVENT][nmaxSttHitsInTrack] ; UShort_t ipinco, ncand, n, nalone; Int_t nrounds0, nrounds1, ipunto; Double_t Dist, Distance, Fi_final_helix_referenceframe, HoughFi[MAXTRACKSPEREVENT], Ntras, Phi, Ptras, Rad, TanL, Z, ddd, dis, dista, dista0, dista0_1, dista1, dista1_1, delta, emme, gap, highqualitycut, Pxini, px, Pyini, py, Pzini, qop, x, y, s[2], Start[3], versor[2], z[2], zeta0, zeta1, zdrift[2], zerror[2], primoangolo[MAXTRACKSPEREVENT], ultimoangolo[MAXTRACKSPEREVENT], AloneX[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], AloneY[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack]; Double_t ALFA[MAXTRACKSPEREVENT], BETA[MAXTRACKSPEREVENT], GAMMA[MAXTRACKSPEREVENT], KAPPA[MAXTRACKSPEREVENT], WDX[nmaxSttHits], WDY[nmaxSttHits], WDZ[nmaxSttHits], Ox[MAXTRACKSPEREVENT], Oy[MAXTRACKSPEREVENT], R[MAXTRACKSPEREVENT], FI0[MAXTRACKSPEREVENT], Fifirst[MAXTRACKSPEREVENT], Fi_low_limit[MAXTRACKSPEREVENT], Fi_up_limit[MAXTRACKSPEREVENT], Posiz1[3], Posiz2[3], Px[MAXTRACKSPEREVENT], Py[MAXTRACKSPEREVENT], Pz[MAXTRACKSPEREVENT], S[2*nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits], // multiplication by 2 in the ZED[2*nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits], // rather improbable chance that DriftRadius[2*nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits],// all skew hits have double ErrorDriftRadius[2*nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits],// solutions Sbis[2*nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits][2], // multiplication by 2 in the ZEDbis[2*nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits][2], // rather improbable chance that DriftRadiusbis[2*nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits][2],// all skew hits have double ErrorDriftRadiusbis[2*nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits][2],// solutions SchosenPixel[MAXTRACKSPEREVENT][nmaxMvdPixelHits], SchosenStrip[MAXTRACKSPEREVENT][nmaxMvdStripHits], SchosenSkew[MAXTRACKSPEREVENT][nmaxSttHits], // NO multiplication by 2 here because for the // skew hits only one // solution is selected. ZchosenPixel[MAXTRACKSPEREVENT][nmaxMvdPixelHits], ZchosenStrip[MAXTRACKSPEREVENT][nmaxMvdStripHits], ZchosenSkew[MAXTRACKSPEREVENT][nmaxSttHits], ErrorchosenPixel[MAXTRACKSPEREVENT][nmaxMvdPixelHits], ErrorchosenStrip[MAXTRACKSPEREVENT][nmaxMvdStripHits], ErrorchosenSkew[MAXTRACKSPEREVENT][nmaxSttHits], info[nmaxSttHits][7], TemporaryS[2*nmaxSttHits], TemporaryZ[2*nmaxSttHits], TemporaryZDrift[2*nmaxSttHits], TemporaryZErrorafterTilt[2*nmaxSttHits]; TVector3 Momentum,ErrMomentum,Position,ErrPosition; FairMCPoint * pSttMCPoint[nmaxSttHits]; PndTrackCand * pMvdTrackCand, * pSttTrackCand; PndSttTrack * pSttHelixTrack; PndSttTube * pSttTube; PndSttHit * pSttHit; PndSdsHit * pMvdPixelHit, * pMvdStripHit; PndSdsMCPoint * pMvdMCPoint; PndMCTrack* pMCtr; PndTrackCandHit pndtrackcandhit; fSttMvdPndTrackCandArray->Delete(); fSttMvdPndTrackArray->Delete(); IVOLTE++; if(istampa>=1) {cout<<"from PndSttMvdTracking, IVOLTE = "<GetEntriesFast(); // num. tracce/evento if (istampa >= 1 && IVOLTE<20) { cout<<"da PndSttMvdTracking : evento (partendo da 0) N. "<< IVOLTE<< "\n N. di MC truth tracks : "<GetEntriesFast(); if (nSttMCPoint ==0){ cout<<"da PndSttMvdTracking : N. di Stt MC points = 0\n"<nmaxSttHits){ cout<<"da PndSttMvdTracking : N. di Stt MC points = "< nmaxSttHits ("<At(i); } //------------------------------ //---------- fetching the STT hits nSttHit = fSttHitArray->GetEntriesFast(); if (nSttHit ==0){ cout<<"da PndSttMvdTracking : N. di Stt Hits = 0, return!\n"< nmaxSttHits) { cout<<"da PndSttMvdTracking : N. di Stt Hits = "< nmaxSttHits (="<= 1 && IVOLTE<20) { cout<<"da PndSttMvdTracking : evento (partendo da 0) N. "<< IVOLTE<< "\n\tN. totale Hits in STT : "<GetEntriesFast(); nMvdStripHit = fMvdStripHitArray->GetEntriesFast(); if(nMvdPixelHit>nmaxMvdPixelHits){ cout<<"from PndSttMvdTracking, nMvdPixelHit is > maximum allowed (" <nmaxMvdStripHits){ cout<<"from PndSttMvdTracking, nMvdStripHit is > maximum allowed (" <GetEntriesFast(); if(nMvdMCPoint>nmaxMvdPixelHits+nmaxMvdStripHits) { cout<<"from PndSttMvdTracking, nMvdMCPoint = "< the maximum number allowed ("<2&& IVOLTE<20) cout<<"N. MC Points delle Mvd = "<GetEntriesFast(); if (nMvdTrackCand> MAXMVDTRACKSPEREVENT) { cout<<"da PndSttMvdTracking : N. of MvdTrackCand = "<< nMvdTrackCand<<" and it is > MAXMVDTRACKSPEREVENT (="<=3 && IVOLTE<20)cout<<"da PndSttMvdTracking : n. totale di Mvd track cand = " <GetEntriesFast(); if (nSttTrackCand ==0){ cout<<"da PndSttMvdTracking : N. of SttTrackCand = 0, return!\n"< MAXTRACKSPEREVENT) { cout<<"da PndSttMvdTracking : N. of nSttTrackCand = "<< nSttTrackCand <<" and it is > MAXTRACKSPEREVENT (="<=2 && IVOLTE<20){ cout<<"N. totale di PndTrackCand del PR solo = " <At(i); ipunto= pSttHit->GetRefIndex();// right way to extract the corrisponding MC point. tubeID = pSttHit->GetTubeID(); pSttTube = (PndSttTube *) fSttTubeArray->At(tubeID); TVector3 center = pSttTube->GetPosition(); // drift radius Double_t dradius = pSttHit->GetIsochrone(); // wire direction TVector3 wiredirection = pSttTube->GetWireDirection(); if(wiredirection.Z() >=0.) { WDX[i] = wiredirection.X(); WDY[i] = wiredirection.Y(); WDZ[i] = wiredirection.Z(); } else { WDX[i] = -wiredirection.X(); WDY[i] = -wiredirection.Y(); WDZ[i] = -wiredirection.Z(); } info[i][0]= pSttTube->GetPosition().X(); info[i][1]= pSttTube->GetPosition().Y(); info[i][2]= pSttTube->GetPosition().Z(); info[i][3]= dradius; info[i][4]= pSttTube->GetHalfLength(); if(ipunto>=0) { info[i][6]= pSttMCPoint[ipunto]->GetTrackID(); } else { info[i][6]= -10.; } if( fabs( WDX[i] )< 0.00001 && fabs( WDY[i] )< 0.00001 ){ info[i][5]= 1.; ListAllParHits[nSttParHit]=i; nSttParHit++; ZCENTER_STRAIGHT = info[i][2]; // this works because just few lines below there is the SEMILENGTH_STRAIGHT = info[i][4]; // requirement that Minclinations[0] > 2 (= at least 3 parallel straws) } else { info[i][5]= 99.; // to signal that it is a skew straw. ListAllSkewHits[nSttSkewHit]=i; nSttSkewHit++; } //--------------- inizio stampaggi, stampe di controllo if (istampa >= 3) { cout <<"da PndSttMvdTracking, Stt iHit "<< i << " e n. punto MC ottenuto con RefIndex = " <GetPosition().X() << " " << pSttTube->GetPosition().Y() << " " << pSttTube->GetPosition().Z() << "; R = "<GetPosition().X()*pSttTube->GetPosition().X()+ pSttTube->GetPosition().Y()*pSttTube->GetPosition().Y())<< ", suo drift radius = "<GetX() << " " << pSttMCPoint[ipunto]->GetY() << " " << pSttMCPoint[ipunto]->GetZ()<GetPosition().X() << " " << pSttTube->GetPosition().Y() << " " << pSttTube->GetPosition().Z() << "; R = "<GetPosition().X()*pSttTube->GetPosition().X()+ pSttTube->GetPosition().Y()*pSttTube->GetPosition().Y())<< ", suo drift radius = "<GetTrackID()<= //-------- fine stampaggi } // end of for( i= 0; i< nSttHit; i++) // fill the exclusion list for those straws with multiple hits // first the parallel straws for(i=0; i< nSttHit-1; i++){ if( !ExclusionListStt[ i ] ) continue; for(j=i+1; j< nSttHit; j++){ if(ExclusionListStt[ j ] && fabs(info[i][0] - info[j][0])<1.e-20 && fabs(info[i][1] - info[j][1])<1.e-20 ) { ExclusionListStt[j]= false ; } } // end of for(j=i+1; j< Nhits;; j++) } // end of for(i=0; i< Nhits-1; i++) //----------------------------------- end of exclusion of straws with multiple hits //------------------------------------------------------------- fine hits STT // --------------------------------------------- estraggo gli HITS Pixel MVD for( i= 0; i< nMvdPixelHit; i++){ pMvdPixelHit = (PndSdsHit *) fMvdPixelHitArray->At(i); // ipunto = pMvdPixelHit->GetRefIndex(); TVector3 temp = pMvdPixelHit->GetPosition(); XMvdPixel[i] = temp.X(); YMvdPixel[i] = temp.Y(); ZMvdPixel[i] = temp.Z(); sigmaXMvdPixel[i] = pMvdPixelHit->GetDx(); sigmaYMvdPixel[i] = pMvdPixelHit->GetDy(); sigmaZMvdPixel[i] = pMvdPixelHit->GetDz(); refindexMvdPixel[i] = pMvdPixelHit->GetRefIndex(); // cout<<"\tPixel n. "<GetBranchId(fMvdPixelBranch) = "<< FairRootManager::Instance()->GetBranchId(fMvdPixelBranch)<GetBranchId(fMvdStripBranch) = "<< FairRootManager::Instance()->GetBranchId(fMvdStripBranch)<=2 && IVOLTE<20) //------------------------------------------ stampaggi // --------------------------------------------- get MC Points of MVD Short_t MCPointtoMCTrackID[nmaxMvdPixelHits+nmaxMvdStripHits]; Double_t XMvdMCPoint[nmaxMvdPixelHits+nmaxMvdStripHits], YMvdMCPoint[nmaxMvdPixelHits+nmaxMvdStripHits], ZMvdMCPoint[nmaxMvdPixelHits+nmaxMvdStripHits]; for( i= 0; i< nMvdMCPoint; i++){ pMvdMCPoint = (PndSdsMCPoint*) fMvdMCPointArray->At(i); TVector3 position; pMvdMCPoint->Position(position); XMvdMCPoint[i]=position.X(); YMvdMCPoint[i]=position.Y(); ZMvdMCPoint[i]=position.Z(); MCPointtoMCTrackID[i]= pMvdMCPoint->GetTrackID(); } // ---------------------------------------------------------------------- // ------------------------------------------ estraggo le altre info della trackcand del MVD bool inMvdTrackCandPixel[nmaxMvdPixelHits], inMvdTrackCandStrip[nmaxMvdStripHits]; for(j=0; jAt(i); TVector3 dirSeed=pMvdTrackCand->getDirSeed(); TVector3 posSeed=pMvdTrackCand->getPosSeed(); qop = pMvdTrackCand->getQoverPseed(); nHitMvdTrackCand[i] = pMvdTrackCand->GetNHits(); // n. hits in questa track cand UShort_t kPixel,kStrip; for(j=0, k=0, kPixel=0, kStrip=0; jGetSortedHit(j); ListHitMvdTrackCand[i][k] = pndtrackcandhit.GetHitId(); // questo e' il n. Hit nativo che posso usare // per estrarre tutte le info che voglio. Se il n. e' -1 // dopo non lo considero; chiedere a Tobias cos'e' perche' // anche il suo detID e' -1 : ne' Pixel ne' Strip...! ListHitTypeMvdTrackCand[i][k] = pndtrackcandhit.GetDetId(); // questo in realta' e' il Branch dello // Hit che viene usato - stupidamente - per dire che e' // un Pixel. Roba da matti. // Se e' -1 dovrebbe essere noise ma a in quale Pixel // o Strip?? Mistero. // the following case should never happen (in principle), therefore don't // consider that hit. if(ListHitMvdTrackCand[i][k]<0 || ListHitTypeMvdTrackCand[i][k]<0)continue; if( ListHitTypeMvdTrackCand[i][k]== FairRootManager::Instance()->GetBranchId(fMvdPixelBranch) && kPixel < nmaxMvdPixelHitsInTrack){ inMvdTrackCandPixel[ ListHitMvdTrackCand[i][k] ]= true; kPixel++; } else if( ListHitTypeMvdTrackCand[i][k]== FairRootManager::Instance()->GetBranchId(fMvdStripBranch) && kStrip < nmaxMvdStripHitsInTrack){ inMvdTrackCandStrip[ ListHitMvdTrackCand[i][k] ]= true; kStrip++; } else { // this is the case should (in principle) never happen. continue; // ignore this hit. } k++; } // end of for(j=0; j=0.){ ListMvdDSPixelHitNotTrackCand[nMvdDSPixelHitNotTrackCand] = i; nMvdDSPixelHitNotTrackCand++; } else { ListMvdUSPixelHitNotTrackCand[nMvdDSPixelHitNotTrackCand] = i; nMvdUSPixelHitNotTrackCand++; } } } for( i= 0; i< nMvdStripHit ; i++){ if( ! inMvdTrackCandStrip[i] ){ if( ZMvdStrip[i]>=0.){ ListMvdDSStripHitNotTrackCand[nMvdDSStripHitNotTrackCand] = i; nMvdDSStripHitNotTrackCand++; } else { ListMvdUSStripHitNotTrackCand[nMvdUSStripHitNotTrackCand] = i; nMvdUSStripHitNotTrackCand++; } } } if(istampa>2){ for( i= 0; i< nMvdTrackCand ; i++){ cout<<" PndSttMvdTracking, MvdTrackCand n. "<At(i); nSttHitsinTrack[i] = pSttTrackCand->GetNHits(); // n. hits in questa track cand if( nSttHitsinTrack[i]> nmaxSttHits) { cout<<"da PndSttMvdTracking : N. di Stt Hits in trackcand "< nmaxSttHits (="<GetSortedHit(0); x = info[ pndtrackcandhit.GetHitId() ][0]; // this is in the middle of the tube y = info[ pndtrackcandhit.GetHitId() ][1]; // this is in the middle of the tube for(j=0,nSttParHitsinTrack[i]=0,nSttSkewHitsinTrack[i]=0; jGetSortedHit(j); ListSttHitsinTrack[i][j] = pndtrackcandhit.GetHitId(); // # hit of Stt if ( fabs(info[ pndtrackcandhit.GetHitId() ][5]- 1.) < 0.0001) { ListSttHitsinTrackType[i][j] = 2; ListSttParHitsinTrack[i][nSttParHitsinTrack[i]] = pndtrackcandhit.GetHitId(); // # hit of Stt nSttParHitsinTrack[i]++; } else { ListSttHitsinTrackType[i][j] = 3; ListSttSkewHitsinTrack[i][nSttSkewHitsinTrack[i]] = pndtrackcandhit.GetHitId(); // # hit of Stt nSttSkewHitsinTrack[i]++; } } // end of for(j=0; jgetDirSeed(); TVector3 posSeed=pSttTrackCand->getPosSeed(); qop = pSttTrackCand->getQoverPseed(); daTrackFoundaTrackMC[i]=pSttTrackCand->getMcTrackId(); Px[i] = dirSeed.X(); Py[i] = dirSeed.Y(); Pz[i] = dirSeed.Z(); if( dirSeed.Z() < 9990.) { // normal case with skew hits information existing SttSZfit[i]=true; Ntras = sqrt(dirSeed.X()*dirSeed.X()+dirSeed.Y()*dirSeed.Y() ) ; R[i] = Ntras /(0.003*BFIELD*fabs(qop)); if(qop <0.){ Ox[i] = -R[i]*dirSeed.Y()/Ntras; Oy[i] = R[i]*dirSeed.X()/Ntras; KAPPA[i]=0.003*BFIELD*fabs(qop)/dirSeed.Z(); CHARGE[i]=-1; } else { Ox[i] = R[i]*dirSeed.Y()/Ntras; Oy[i] = -R[i]*dirSeed.X()/Ntras; KAPPA[i]=-0.003*BFIELD*fabs(qop)/dirSeed.Z(); CHARGE[i]=1; } } else if (dirSeed.Z()-10000.<0.) { // case without skew straws info present SttSZfit[i]=false; Ptras = sqrt(dirSeed.X()*dirSeed.X()+dirSeed.Y()*dirSeed.Y() ) ; R[i] = Ptras /(0.003*BFIELD); if(qop <0.){ Ox[i] = -R[i]*dirSeed.Y()/Ptras; Oy[i] = R[i]*dirSeed.X()/Ptras; CHARGE[i]=-1; } else { Ox[i] = R[i]*dirSeed.Y()/Ptras; Oy[i] = -R[i]*dirSeed.X()/Ptras; CHARGE[i]=1; } } else { // case in which KAPPA = 0 SttSZfit[i]=true; if(qop <0.){ Ox[i] = -R[i]*dirSeed.Y()/Ntras; Oy[i] = R[i]*dirSeed.X()/Ntras; CHARGE[i]=-1; } else { Ox[i] = R[i]*dirSeed.Y()/Ptras; Oy[i] = -R[i]*dirSeed.X()/Ptras; CHARGE[i]=1; } KAPPA[i]=0.; } Fifirst[i] = atan2( y-Oy[i], x-Ox[i]); FI0[i] = atan2(-Oy[i], -Ox[i]); HoughFi[i] = FI0[i] + PI; if(HoughFi[i]<0.) HoughFi[i]=0.; if( FI0[i] < 0. ) FI0[i]+= 2.*PI; if( FI0[i] < 0. ) FI0[i]=0.; if( Fifirst[i] < 0. ) Fifirst[i]+= 2.*PI;if( Fifirst[i] < 0. ) Fifirst[i]=0.; } // end of for( i= 0; i< nSttTrackCand; i++) //--------------- end fetching of the Stt PndTrackCand from PR of the STT //----------- treat here the special pathological case when there are no Mvd hits at all. //----------- Copy simply the quantities from the STT alone PndTrackCand. if(nMvdPixelHit+nMvdStripHit<=0){ for( ncand= 0; ncand< nSttTrackCand; ncand++){ int nbuoni, nparbuoni, nskewbuoni; for(j=0, nbuoni=0,nparbuoni=0, nskewbuoni=0;j fabs(ddd) ){ dista0=ddd; zeta0 = z[0]-zdrift[0]; } else { zeta0 = z[0]+zdrift[0]; } dista1 = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]+zdrift[1],s[1],&nrounds1); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]-zdrift[1],s[1],&nrounds1); if( fabs(dista1)> fabs(ddd) ){ dista1=ddd; zeta1 = z[1]-zdrift[1]; } else { zeta1 = z[1]+zdrift[1]; } if( fabs(dista1) fabs(ddd) ){ ZchosenSkew[ncand][k]=z[0]-zdrift[0]; } else { ZchosenSkew[ncand][k]=z[0]+zdrift[0]; } ListTrackCandHit[ncand][nbuoni]= ListSttSkewHitsinTrack[ncand][nskewbuoni]= ListSttHitsinTrack[ncand][j]; ListTrackCandHitType[ncand][nbuoni]= ListSttHitsinTrackType[ncand][j]; nskewbuoni++; nbuoni++; } else if( z[1]<999998.){ dista = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand] ,z[1]+zdrift[1],s[1],&nrounds0); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]-zdrift[1],s[1],&nrounds1); SchosenSkew[ncand][k]=s[1]; if( fabs(dista)> fabs(ddd) ){ ZchosenSkew[ncand][k]=z[1]-zdrift[1]; } else { ZchosenSkew[ncand][k]=z[1]+zdrift[1]; } ListTrackCandHit[ncand][nbuoni]= ListSttSkewHitsinTrack[ncand][nskewbuoni]= ListSttHitsinTrack[ncand][j]; ListTrackCandHitType[ncand][nbuoni]= ListSttHitsinTrackType[ncand][j]; nskewbuoni++; nbuoni++; } // end of if( z[0]<999998. && z[1]<999998.) } else { // continuation of if(ListSttHitsinTrackType[ncand][j]==3) // so, here it must be a parallel hit. ListTrackCandHit[ncand][nbuoni]= ListSttParHitsinTrack[ncand][nparbuoni]= ListSttHitsinTrack[ncand][j]; ListTrackCandHitType[ncand][nbuoni]= ListSttHitsinTrackType[ncand][j]; nparbuoni++; nbuoni++; } // end of if(ListSttHitsinTrackType[ncand][j]==3) } // end of for(j=0; jFill(dis); hdeltaRPixel2->Fill(dis); } for( j= 0; j< nMvdStripHit; j++){ dis = sqrt( pow(XMvdStrip[j]-Ox[i],2)+pow(YMvdStrip[j]-Oy[i],2)) - R[i]; hdeltaRStrip->Fill(dis); hdeltaRStrip2->Fill(dis); } } // end of for( i= 0; i< nSttTrackCand; i++) } //-------- fine plottamenti //------------------- //------------------- //------------------- //------------------- //------------------- //------------------- //------------------- start the combined Mvd-Stt PR //------------------- //------------------- //------------------- //------------------- //------------------- //---- find the angular range (in Fi) allowed for the STT hits, with the present Ox,Oy and R // of the track candidates, and for the Mvd hits (FI0 and Fi_low_limit ). for( ncand= 0; ncand< nSttTrackCand; ncand++){ keepit[ncand]=true; PndSttFindingParallelTrackAngularRange( Ox[ncand], Oy[ncand], R[ncand], CHARGE[ncand], &Fi_low_limit[ncand], // Fi (in XY Helix frame) lower limit using // the Stt detector minimum/maximum radius // Fi_low_limit is ALWAYS between 0. and 2PI &Fi_up_limit[ncand], // Fi (in XY Helix frame) upper limit using // the Stt detector maximum/minimum radius // Fi_up_limit is ALWAYS > Fi_low_limit and // possibly > 2PI. &statusflag[ncand],//it is a vector; =0, all well; =1, track contained completely between RMin // and RMax; = -1 track contained within RMin; =-2 track outside RMax. RStrawDetectorMin, RStrawDetectorMax ); if( statusflag[ncand] == -1) { Fi_low_limit[ncand] = -99999.; } else if (statusflag[ncand] == -2){ keepit[ncand] = false; }; } // end of for( ncand= 0; ncand< nSttTrackCand; ncand++) //--------------------- here call to the function that matches Mvd hits with Stt tracks delta=0.5; // parameter of proximity for associating Mvd hits to Stt tracks highqualitycut=0.2; // parameter of proximity for associating Mvd hits to Stt tracks // This method matches the Mvd hits to the found tracks. MatchMvdHitsToSttTracks2( keepit,// input and output. If there are no MVD hits associate, // keepit is set to false. delta, highqualitycut, nSttTrackCand, Ox, Oy, R, FI0, Fi_low_limit, CHARGE, nMvdPixelHitsinTrack, // output ListMvdPixelHitsinTrack, // output nMvdStripHitsinTrack, // output ListMvdStripHitsinTrack // output ); for( i= 0; i< nSttTrackCand; i++){ if( nMvdPixelHitsinTrack[i] > nmaxMvdPixelHitsInTrack){ cout<<"from PndSttMvdTracking, nMvdPixelHitsinTrack["< nmaxMvdPixelHitsInTrack (" < nmaxMvdStripHitsInTrack){ cout<<"from PndSttMvdTracking, nMvdStripHitsinTrack["< nmaxMvdStripHitsInTrack (" <=2&& IVOLTE<20){ cout<<"da PndSttMvdTracking ; n. SttTrackCand totali = "<=0) //----------------- end of section with match Mvd hits with Stt hits //---------------------------------------------------------- // ordering all the hits belonging to the candidate track, by increasing R; // forming the new track with Mvd+Stt hits UShort_t tempmvdindex[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], tempmvdtype[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], auxIndex[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack]; Double_t auxR[nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack]; nTotalCandidates = nSttTrackCand; // nSttTrackCand is already <= MAXTRACKSPEREVENT. for(ncand=0; ncand< nTotalCandidates; ncand++){ if(!keepit[ncand]) continue; nTrackCandHit[ncand] =nSttHitsinTrack[ncand]+ nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; // adding the Mvd hits (Pixel and Strips) for(i=0; i< nMvdPixelHitsinTrack[ncand]; i++){ auxR[i] = XMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]* XMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]+ YMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]* YMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ] ; tempmvdindex[i]=ListMvdPixelHitsinTrack[ncand][i]; tempmvdtype[i]=0; auxIndex[i] = i; } for(i=0; i< nMvdStripHitsinTrack[ncand]; i++){ auxR[i+nMvdPixelHitsinTrack[ncand]] = XMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]* XMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]+ YMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]* YMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ] ; tempmvdindex[i+nMvdPixelHitsinTrack[ncand]]= ListMvdStripHitsinTrack[ncand][i]; tempmvdtype[i+nMvdPixelHitsinTrack[ncand]]=1; auxIndex[i+nMvdPixelHitsinTrack[ncand]]= i+nMvdPixelHitsinTrack[ncand]; } // ordering the Mvd Hits if( nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] >0){ Merge_Sort( nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand], auxR, auxIndex); // constructing the ordered new Track Candidate now for(i=0; i< nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; i++){ ListTrackCandHit[ncand][i] = tempmvdindex[ auxIndex[i] ]; ListTrackCandHitType[ncand][i] = tempmvdtype[ auxIndex[i] ]; } } // end of if( nMvdPixelHitsinTrack[ncand]+ for(i=0; i=2) cout<<"da PndSttMvdTracking, IVOLTE = "<< IVOLTE<<", ncand = " < Pixel; = 1 --> Strip; = 2 Stt parallel; 3 = Stt Skew. // There are no -1 type hits at this point. if(!(ListTrackCandHitType[ncand][i]==1 || ListTrackCandHitType[ncand][i]==0) ) continue; if( ListTrackCandHitType[ncand][i] == 0 ){ d = fabs( sqrt( (Ox[ncand]-XMvdPixel[ListTrackCandHit[ncand][i]]) *(Ox[ncand]-XMvdPixel[ListTrackCandHit[ncand][i]]) +(Oy[ncand]-YMvdPixel[ListTrackCandHit[ncand][i]]) *(Oy[ncand]-YMvdPixel[ListTrackCandHit[ncand][i]]) ) - R[ncand]); if( d < diff ){ diff = d; iexcl=i; trajectory_vertex[0]= XMvdPixel[ListTrackCandHit[ncand][i]]; trajectory_vertex[1]= YMvdPixel[ListTrackCandHit[ncand][i]]; } } else { d = fabs( sqrt( (Ox[ncand]-XMvdStrip[ListTrackCandHit[ncand][i]]) *(Ox[ncand]-XMvdStrip[ListTrackCandHit[ncand][i]]) +(Oy[ncand]-YMvdStrip[ListTrackCandHit[ncand][i]]) *(Oy[ncand]-YMvdStrip[ListTrackCandHit[ncand][i]]) ) - R[ncand]); if( d < diff ){ diff = d; iexcl=i; trajectory_vertex[0]= XMvdStrip[ListTrackCandHit[ncand][i]]; trajectory_vertex[1]= YMvdStrip[ListTrackCandHit[ncand][i]]; } } } // end of for(i=0, diff = 1.e20;i<..... if( diff > 0.5 ) { trajectory_vertex[0]=trajectory_vertex[1]=0.; iexcl=-1; } RefitMvdStt( nTrackCandHit[ncand], //this is input &ListTrackCandHit[ncand][0],//this is both input and output &ListTrackCandHitType[ncand][0],//this is both input and output info, rotationangle, trajectory_vertex, iexcl, &ALFA[ncand], // output of the fit &BETA[ncand], // output of the fit &GAMMA[ncand],// output of the fit &status[ncand] // fit status; true = successful ); if(status[ncand]){ Ox[ncand] = -ALFA[ncand]/2.; Oy[ncand] = -BETA[ncand]/2.; R[ncand] = Ox[ncand]*Ox[ncand]+Oy[ncand]*Oy[ncand]-GAMMA[ncand]; if( R[ncand] > 0. ) R[ncand]=sqrt(R[ncand]) ; else R[ncand]=0.; FI0[ncand] = atan2(-Oy[ncand], -Ox[ncand]); if( FI0[ncand] < 0. ) FI0[ncand]+= 2.*PI; } if(R[ncand]2) { cout<<"PndSttMvdTracking, prima di MatchMvdHitsToSttTracksagain, IVOLTE = "<< IVOLTE<<", ncand "< Fi_low_limit and // possibly > 2PI. // Fi_low_limit and Fi_up_limit are undefined // when statusflag is negative. &statusflag[ncand],//it is a vector; =0, all well; =1, track contained completely between RMin // and RMax; = -1 track contained within RMin; =-2 track outside RMax. RStrawDetectorMin, RStrawDetectorMax ); // Fi_low_limit set to -99999. when it is contained in Mvd // region completely. if( statusflag[ncand] == -1 ) Fi_low_limit[ncand] = -99999.; if( statusflag[ncand] == -2) keepit[ncand] = false; } // end of for(ncand=0; ncand< nTotalCandidates; ncand++) //--------------------- end of refit the Helix in XY plane using Stt + Mvd associated hits //--------------------- here call to the function that matches Mvd hits with the new // circular trajectory in XY found for the second time, after first refit delta=0.5; // parameter of proximity for associating Mvd hits to Stt tracks // highqualitycut=0.3; // parameter of proximity for associating Mvd hits to Stt tracks highqualitycut=0.5; // parameter of proximity for associating Mvd hits to Stt tracks MatchMvdHitsToSttTracksagain( keepit, Mvdhits, delta, highqualitycut, nSttTrackCand, Ox, Oy, R, FI0, Fi_low_limit, // because here we deal with hits in Mvd region CHARGE, nMvdPixelHitsinTrack, // input and output ListMvdPixelHitsinTrack, // input and output nMvdStripHitsinTrack, // input and output ListMvdStripHitsinTrack // input and output ); //--------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // use the risult just obtained from the fit in XY to redo the association of the Skew Straw hits for(ncand=0; ncand< nTotalCandidates; ncand++) { if(!keepit[ncand]) continue; if( ! Mvdhits[ncand]) continue; if( statusflag[ncand] == -1 ) continue; // this is when the XY circle is contained in the // the Mvd region completely; skip the association of // the Skews. TemporarynSkewHitsinTrack = AssociateSkewHitsToXYTrack( ExclusionListStt, nSttSkewHit, ListAllSkewHits, Ox[ncand], // input : X of center of XY plane circle Oy[ncand], // input : Y of center of XY plane circle R[ncand], // input : Radius of XY plane circle info, WDX, WDY, WDZ, Fi_low_limit[ncand], // Fi (in Helix XY frame) lower limit using the Stt detector minimum/maximum radius Fi_up_limit[ncand], // Fi (in Helix XY frame) upper limit using the Stt detector maximum/minimum radius CHARGE[ncand], TemporarySkewList, // output, list of selected skew hits (in skew numbering) TemporaryS, // output, S coordinate of selected Skew hit TemporaryZ, // output, Z coordinate of selected Skew hit (center wire) TemporaryZDrift, // output, drift distance IN Z DIRECTION only, of selected Skew hit TemporaryZErrorafterTilt // output, Radius taking into account the tilt, IN Z DIRECTION only, of selected Skew hit ); nSttSkewHitsinTrack[ncand]=TemporarynSkewHitsinTrack; // it can be also zero! for(j=0;j2) { cout<<"PndSttMvdTracking, dopo di AssociateSkewHitsToXYTrack, IVOLTE = "<< IVOLTE<<", ncand "<1.e-10) { ErrorDriftRadiusbis[kall][0]=ErrorDriftRadius[i]=zerror[0]; } else { ErrorDriftRadiusbis[kall][0]=ErrorDriftRadius[i]=0.5; } i++; ZEDbis[kall][1]=ZED[i]=z[1]; Sbis[kall][1]=S[i] = s[1]; DriftRadiusbis[kall][1]=DriftRadius[i]=zdrift[1]; if( fabs(zdrift[1]) >1.e-10) { ErrorDriftRadiusbis[kall][1]=ErrorDriftRadius[i]=zerror[1]; } else { ErrorDriftRadiusbis[kall][1]=ErrorDriftRadius[i]=0.5; } i++; if(!Mvdhits[ncand]) // in this case calculate also ZchosenSkew and SchosenSkew { dista0 = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[0]+zdrift[0],s[0],&nrounds0); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[0]-zdrift[0],s[0],&nrounds0); if( fabs(dista0)> fabs(ddd) ){ dista0=ddd; zeta0 = z[0]-zdrift[0]; } else { zeta0 = z[0]+zdrift[0]; } dista1 = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]+zdrift[1],s[1],&nrounds1); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]-zdrift[1],s[1],&nrounds1); if( fabs(dista1)> fabs(ddd) ){ dista1=ddd; zeta1 = z[1]-zdrift[1]; } else { zeta1 = z[1]+zdrift[1]; } if( fabs(dista1)1.e-10) { ErrorDriftRadiusbis[kall][0]=ErrorDriftRadius[i]=zerror[0]; } else { ErrorDriftRadiusbis[kall][0]=ErrorDriftRadius[i]=0.5; } ZEDbis[kall][1]=999999.; i++; if(!Mvdhits[ncand]) // in this case calculate also ZchosenSkew and SchosenSkew { dista = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[0]+zdrift[0],s[0],&nrounds0); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[0]-zdrift[0],s[0],&nrounds0); SchosenSkew[ncand][k]=s[0]; if( fabs(dista)> fabs(ddd) ){ ZchosenSkew[ncand][k]=z[0]-zdrift[0]; } else { ZchosenSkew[ncand][k]=z[0]+zdrift[0]; } }// end of if(Mvdhits[ncand]) } else if( z[1]<999998.){ ZEDbis[kall][1]=ZED[i]=z[1]; Sbis[kall][1]=S[i] = s[1]; DriftRadiusbis[kall][1]=DriftRadius[i]=zdrift[1]; if( fabs(zdrift[1]) >1.e-10) { ErrorDriftRadiusbis[kall][1]=ErrorDriftRadius[i]=zerror[1]; } else { ErrorDriftRadiusbis[kall][1]=ErrorDriftRadius[i]=0.5; } ZEDbis[kall][0]=999999.; i++; if(!Mvdhits[ncand]) { // in this case calculate also ZchosenSkew and SchosenSkew dista = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand] ,z[1]+zdrift[1],s[1],&nrounds0); ddd = Dist_SZ(R[ncand],KAPPA[ncand],FI0[ncand], z[1]-zdrift[1],s[1],&nrounds1); SchosenSkew[ncand][k]=s[1]; if( fabs(dista)> fabs(ddd) ){ ZchosenSkew[ncand][k]=z[1]-zdrift[1]; } else { ZchosenSkew[ncand][k]=z[1]+zdrift[1]; } } // end of if(Mvdhits[ncand]) } else { ZEDbis[kall][1]=ZEDbis[kall][0]=999999.; } } // end of for(j=0, i = nMvdPixelHitsinTrack[ncand]+ j=nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; // finding if there are discontinuity at 0 for fi value of the Mvd Hit. // In case of discontinuity at 0, add 2*PI to fi of those hits with fi in the 1st quadrant. // This is necessary because the discontinuities would make the fit // in the SZ plane fail. // In this discontinuity fixing, the value FI0 of the vertex (0,0) is also included. // If there is discontinuity fixing, the values of S[i] AND POSSIBLY Fi_initial_helix_referenceframe[i] // might be modified (+2.*PI) from now on. if(Mvdhits[ncand]){ // in this case there is at least 1 Mvd hits associated to Stt track. if( j ==1){ nSttSkewHitsinTrack[ncand]<5 ? j += nSttSkewHitsinTrack[ncand] : j += 5; } else if (j==2){ // trick when 2 Mvd hit are very discordant double dot = ZED[0]*ZED[1] + (S[0]-FI0[ncand])*(S[1]-FI0[ncand]); double modulo0 = ZED[0]*ZED[0] + (S[0]-FI0[ncand])*(S[0]-FI0[ncand]); double modulo1 = ZED[1]*ZED[1] + (S[1]-FI0[ncand])*(S[1]-FI0[ncand]); if( modulo0 > 1.e-20 && modulo1 > 1.e-20 ){ dot /= (sqrt(modulo1*modulo0)); if(dot < 0.866 ) { nSttSkewHitsinTrack[ncand]<5?j +=nSttSkewHitsinTrack[ncand]:j += 5; } } } FixDiscontinuitiesFiangleinSZplane( j, // n hits in track. S, // S can be modified by +-2*PI if necessary. &FI0[ncand], // this remains unchanged. CHARGE[ncand] // this remains unchanged. ); //--------------------- here do the fit again in the SZ space if there are Mvd hits. // For this, reordering of the Mvd hits is not necessary. resultFitSZagain[ncand] = FitSZspace( j, // n. hits to be fitted S, ZED, DriftRadius, ErrorDriftRadius, FI0[ncand], 20, // maximum number allowed in the fit &emme ); if( resultFitSZagain[ncand]==1){ KAPPA[ncand] = emme; // in case it was a TrackCand from Stt alone, originally with no Skew // hit and therefore without KAPPA information, signal that now the // SZ fit worked (thanks to the Mvd hits) by setting SttSZfit[ncand]=true. if( ncand<= nSttTrackCand ) SttSZfit[ncand]=true; } else { keepit[ncand]=false; } //------------------------------------------- // use the risult just obtained from the fit in SZ to reject the spurious Skew Straw hits // and the Mvd spurious hits; also in this function there is the calculation of the Z position // of the SKEW hits and the MVD hits, for a given track candidate (ie for a given Helix // circle in the XY plane) UShort_t MaxTurns; if(R[ncand] < RStrawDetectorMax/2.){ if(-CHARGE[ncand]*KAPPA[ncand]>0.) // this means Pz>0. MaxTurns=(UShort_t)(0.5*fabs((ZCENTER_STRAIGHT+SEMILENGTH_STRAIGHT) *KAPPA[ncand])/PI); else MaxTurns=(UShort_t) (0.5*fabs((ZCENTER_STRAIGHT-SEMILENGTH_STRAIGHT)*KAPPA[ncand])/PI); } else { MaxTurns=0; } if(keepit[ncand]) EliminateSpuriousSZ( MaxTurns, &nMvdPixelHitsinTrack[ncand], // input and output &ListMvdPixelHitsinTrack[ncand][0],// input and output &nMvdStripHitsinTrack[ncand], // input and output &ListMvdStripHitsinTrack[ncand][0],// input and output &nSttSkewHitsinTrack[ncand], // input and output &ListSttSkewHitsinTrack[ncand][0],// input and output Sbis, // input, position of the central wire ZEDbis, // input, position of the central wire. DriftRadiusbis, // input ErrorDriftRadiusbis, // input &SchosenPixel[ncand][0], // this value from now on &SchosenStrip[ncand][0], // can also be > 2PI or < 2PI when &SchosenSkew[ncand][0], // the particle makes more than 1 turn. &ZchosenPixel[ncand][0], &ZchosenStrip[ncand][0], &ZchosenSkew[ncand][0], &ErrorchosenPixel[ncand][0], &ErrorchosenStrip[ncand][0], &ErrorchosenSkew[ncand][0], KAPPA[ncand], FI0[ncand], R[ncand] ); //------------------------ } // end of if(Mvdhits[ncand]) // First cleanup based on the absence of Mvd hits // reject the candidate if it is NOT contained in the pipe and // therefore it should have at least 1 Mvd hit but it has none. if( (!IsInTargetPipe( Ox[ncand], Oy[ncand], R[ncand], FI0[ncand], KAPPA[ncand], CHARGE[ncand], VERTICALGAP/2.) ) && nMvdStripHitsinTrack[ncand]+nMvdPixelHitsinTrack[ncand]==0) { keepit[ncand]=false; } } // end of for(ncand=0; ncand< nTotalCandidates; ncand++) //--- redo association of parallel Stt straw hits to this track, after better refit. if(istampa>2&& IVOLTE<20){ cout<<"da PndSttMvdTracking, prima di CollectParSttHitsagain ; n. SttTrackCand totali = "<=0) CollectParSttHitsagain( keepit, Mvdhits, info, nSttParHit, ListAllParHits, 0, // starting candidate n. (included). nTotalCandidates, // ending candidate n. (excluded). Ox, Oy, R, KAPPA, FI0, Fi_low_limit, Fi_up_limit, nSttParHitsinTrack, // input and output ListSttParHitsinTrack // input and output ); // ordering all the hits belonging to the candidate track, by increasing R; // from candidate n. 0 to candidate n. nTotalCandidates-1; loading ListTrackCandHit. // the array ordered are : // ListTrackCandHit, ListTrackCandHitType, ListSttParHitsinTrack, ListSttSkewHitsinTrack. Ordering_Loading_ListTrackCandHit(keepit,0,nTotalCandidates,info); //------------- cleanup section. Start[0]=0.; Start[1]=0.; Start[2]=0.; gap = (Double_t) (VERTICALGAP); for(ncand=0, nRemainingCandidates=0; ncand< nTotalCandidates; ncand++){ if(!keepit[ncand]) continue; UShort_t &nHitsPar = nSttParHitsinTrack[ncand]; UShort_t &nHitsSkew = nSttSkewHitsinTrack[ncand]; Double_t auxS[nHitsSkew]; for(i=0;i=2&&IVOLTE<20){ cout<<"\n------------------------------------------------\n" <<"\tIVOLTE = "<=3){ cout<<"\n------------------------------------------------\n" <<"\tIVOLTE = "<= MAXTRACKSPEREVENT ) goto dontdoit ; // if there are already // too many candidates // don't do anything. // load the UsedPixel and UsedStrip vectors. bool UsedPixel[nmaxMvdPixelHits], UsedStrip[nmaxMvdStripHits]; UShort_t List[nmaxMvdPixelHits+nmaxMvdStripHits], ListType[nmaxMvdPixelHits+nmaxMvdStripHits]; Double_t centrex, centrey, radius; for(i=0;iGetBranchId(fMvdPixelBranch)){ if(UsedPixel[ ListHitMvdTrackCand[i][j] ]) continue; AloneX[nalone] = XMvdPixel[ ListHitMvdTrackCand[i][j] ]; AloneY[nalone] = YMvdPixel[ ListHitMvdTrackCand[i][j] ]; List[nalone] = ListHitMvdTrackCand[i][j]; ListType[nalone] = 0; nalone++; }else {// at this point this is a Strip hit; already made sure // earlier in the code that there is no third possibility. // if (ListHitTypeMvdTrackCand[i][j] == // this condition in principle at this // // point is always satisfied. // FairRootManager::Instance()->GetBranchId(fMvdStripBranch)){ if(UsedStrip[ ListHitMvdTrackCand[i][j] ]) continue; AloneX[nalone] = XMvdStrip[ ListHitMvdTrackCand[i][j] ]; AloneY[nalone] = YMvdStrip[ ListHitMvdTrackCand[i][j] ]; List[nalone] = ListHitMvdTrackCand[i][j]; ListType[nalone] = 1; nalone++; } } // end of for(j=0;j 0. ) R[nTotalCandidates]= sqrt(R[nTotalCandidates]) ; else R[nTotalCandidates]=0.; FI0[nTotalCandidates] = atan2(-Oy[nTotalCandidates], -Ox[nTotalCandidates]); if( FI0[nTotalCandidates] < 0. ) FI0[nTotalCandidates]+= 2.*PI; nSttParHitsinTrack[nTotalCandidates]=0; nSttSkewHitsinTrack[nTotalCandidates]=0; nMvdPixelHitsinTrack[nTotalCandidates]=0; nMvdStripHitsinTrack[nTotalCandidates]=0; for(j=0; j0. flag = false; for(i=1;i 1.5*PI ){ flag=true; break; } } if( flag ){ if( esse[0] < PI) Sini = esse[0]+2.*PI; if( esse[nTrackCandHit[ncand]-1] < PI) Slast=esse[nTrackCandHit[ncand]-1]+2.*PI; } else { Sini = esse[0]; Slast= esse[nTrackCandHit[ncand]-1]; } if( Sini > Slast) { CHARGE[ncand] = 1; } else { CHARGE[ncand] = -1; } // it is necessary first to calculate Fi_low_limit and Fi_up_limit. PndSttFindingParallelTrackAngularRange( Ox[ncand], Oy[ncand], R[ncand], CHARGE[ncand], &Fi_low_limit[ncand], // Fi (in XY Helix frame) lower limit using // the Stt detector minimum/maximum radius // Fi_low_limit is ALWAYS between 0. and 2PI &Fi_up_limit[ncand], // Fi (in XY Helix frame) upper limit using // the Stt detector maximum/minimum radius // Fi_up_limit is ALWAYS > Fi_low_limit and // possibly > 2PI. &statusflag[ncand],// vector; =0, all well; =1, track contained completely between RMin // and RMax; = -1 track contained within RMin; =-2 track outside RMax. RStrawDetectorMin, RStrawDetectorMax ); // discard tracks completely outside RMin (since they are formed by Mvd hits). if(statusflag[ncand]== -2 || statusflag[ncand]== 1){ keepit[ncand]=false; continue; } // take care of possible discontinuities at 0. FixDiscontinuitiesFiangleinSZplane( nTrackCandHit[ncand], // n hits in track. S, // S can be modified by +-2*PI if necessary. &FI0[ncand], // this remains unchanged. CHARGE[ncand] // this remains unchanged. ); // ----------------------------- fit in SZ with the Mvd tracks resultFitSZagain[ncand] = FitSZspace( nMvdOnly, // n. hits to be fitted S, ZED, DriftRadius, ErrorDriftRadius, FI0[ncand], 20, // maximum number allowed in the fit &emme ); if( resultFitSZagain[ncand]==1){ KAPPA[ncand] = emme; } else { keepit[ncand]=false; } } // end of for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++) // ----------------------------- end fit in SZ with the Mvd tracks //--- try to associate of parallel Stt straw hits to this new tracks. CollectParSttHitsagain( keepit, Mvdhits, info, nSttParHit, ListAllParHits, nSttTrackCand, // starting candidate n. (included). nTotalCandidates, // ending candidate n. (excluded). Ox, Oy, R, KAPPA, FI0, Fi_low_limit, Fi_up_limit, nSttParHitsinTrack, // input and output ListSttParHitsinTrack // input and output ); for(ncand=nSttTrackCand; ncand< nTotalCandidates; ncand++){ // try to attach skew hits to the new tracks (the result can also be 0). if(!keepit[ncand]) continue; // when statusflag[ncand]<0 the track does not intersect Stt region. if( statusflag[ncand]<0) continue; nSttSkewHitsinTrack[ncand] = AssociateSkewHitsToXYTrack( ExclusionListStt, nSttSkewHit, ListAllSkewHits, Ox[ncand], // input : X of center of XY plane circle Oy[ncand], // input : Y of center of XY plane circle R[ncand], // input : Radius of XY plane circle info, WDX, WDY, WDZ, Fi_low_limit[ncand], // Fi (in Helix XY frame) lower limit using the Stt detector minimum/maximum radius Fi_up_limit[ncand], // Fi (in Helix XY frame) upper limit using the Stt detector maximum/minimum radius CHARGE[ncand], TemporarySkewList, // output, list of selected skew hits (in skew numbering) TemporaryS, // output, S coordinate of selected Skew hit TemporaryZ, // output, Z coordinate of selected Skew hit (center wire) TemporaryZDrift, // output, drift distance IN Z DIRECTION only, of selected Skew hit TemporaryZErrorafterTilt // output, Radius taking into account the tilt, IN Z DIRECTION only, of selected Skew hit ); for(j=0;j=2){ cout<<"Evt. n. "<0 && nTotalCandidates > 0 && doMcComparison ){ int nmid,nn; Double_t X1[MAXTRACKSPEREVENT], Y1[MAXTRACKSPEREVENT], X2[MAXTRACKSPEREVENT], Y2[MAXTRACKSPEREVENT], X3[MAXTRACKSPEREVENT], Y3[MAXTRACKSPEREVENT]; for(i=0; i0 ){ // clockwise. if( FI0[i] < angle ) angle -= 2.*PI; if( FI0[i] < angle ) angle =0.; } else{ // counterclockwise. if( FI0[i] > angle ) angle += 2.*PI; if( FI0[i] > angle ) angle = FI0[i]; } middle = (FI0[i]+angle)/2.; X2[i] = Ox[i] + R[i]*cos( middle ); Y2[i] = Oy[i] + R[i]*sin( middle ); } // end of for(i=0; i=2){ // for(i=0;i=3){ // for(ncand=0; ncand< nSttTrackCand; ncand++){ int nbuone=-1; for(ncand=0; ncand< nTotalCandidates; ncand++){ if(!keepit[ncand]) continue; nbuone++; //----------------- ora la traccia MC corrispondente a questa traccia Stt // cout<<"\n\nda PndSttMvdTracking : la PndTrackCand n. "<-1 ) { Int_t icode, im; Double_t aaa,Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Pxx, Pyy, carica ; PndMCTrack* pMC; im=daTrackFoundaTrackMC[ncand]; pMC = (PndMCTrack*) fMCTrackArray->At(im); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Pxx = pMC->GetMomentum().X(); Pyy = pMC->GetMomentum().Y(); aaa = sqrt( Pxx*Pxx + Pyy*Pyy); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if(fabs(carica)<1.e-5) { cout<<"...associata ad una traccia neutra, assurdo!\n"; continue ;} Cx = Oxx + Pyy*1000./(BFIELD*CVEL*carica); Cy = Oyy - Pxx*1000./(BFIELD*CVEL*carica); cout<<" e' associata alla traccia MC n. "<=2) //-------------fine stampaggi // ora il confronto per il meeting di Groningen //---------- conteggio delle tracce MC accettabili!! int citata; int nMCTracksaccettabili=0; for (i=0;iAt(i); if ( ! pMCtr ) continue; Double_t aaa, carica, Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Pxx, Pyy ; Int_t icode; icode = pMCtr->GetPdgCode() ; // PDG code of track Oxx = pMCtr->GetStartVertex().X(); // X of starting point track Oyy = pMCtr->GetStartVertex().Y(); // Y of starting point track Pxx = pMCtr->GetMomentum().X(); Pyy = pMCtr->GetMomentum().Y(); aaa = sqrt( Pxx*Pxx + Pyy*Pyy); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla if(istampa>2){ TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if(fabs(carica)<1.e-5) continue; Cx = Oxx + Pyy*1000./(BFIELD*CVEL*carica); Cy = Oyy - Pxx*1000./(BFIELD*CVEL*carica); cout<<"da PndSttMvdTracking, evento (cominciando da 0) n. "<=1 ) fprintf(HANDLE, "\n Evento %d NTotaleTracceMC %d ------\n",IVOLTE, nMCTracks); if(istampa>=1 ) { fprintf(HANDLE, "\n Evento %d NTotaleTracceMC %d ------\n",IVOLTE, nMCTracksaccettabili); int ii, ibuone=-1; Double_t HoughFiii; //for (ii=0; ii=3 ;ii++){ for (ii=0; iiAt(i); if ( ! pMCtr ){ fprintf(HANDLE, " MC track n. %d doesn't have pointer to MC Track TClones Array\n", i); continue; } fprintf(HANDLE, " TracciaMC %d ParHitsMC %d ParMecc %d ParMeccSpuri %d SkewHitsMC %d SkewMecc %d SkewMeccSpuri %d\n", i, nHitsInMCTrack[ii], nParalCommon[ii], nSpuriParinTrack[ii], nSkewHitsInMCTrack[ii], nSkewCommon[ii], nSpuriSkewinTrack[ii] ) ; fprintf(HANDLE, " e corrisponde a track found n. %d\n", ibuone ); fprintf(HANDLE, " AVENDO %d hits paralleli e %d hits skew non mecciati dalla corrisponde track found\n" , nMCParalAlone[ii],nMCSkewAlone[ii] ); HoughFiii = atan2(Oy[ii],Ox[ii]); if(HoughFiii<0.) HoughFiii += 2.*PI; Double_t aaa, carica, Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Pxx, Pyy ; Int_t icode; icode = pMCtr->GetPdgCode() ; // PDG code of track Oxx = pMCtr->GetStartVertex().X(); // X of starting point track Oyy = pMCtr->GetStartVertex().Y(); // Y of starting point track Pxx = pMCtr->GetMomentum().X(); Pyy = pMCtr->GetMomentum().Y(); aaa = sqrt( Pxx*Pxx + Pyy*Pyy); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if(fabs(carica)<1.e-5) fprintf(HANDLE," MC track n. %d e' neutra, assurdo!\n",i); Cx = Oxx + Pyy*1000./(BFIELD*CVEL*carica); Cy = Oyy - Pxx*1000./(BFIELD*CVEL*carica); Fifi = atan2(Cy, Cx); // MC truth Fifi angle of circle of Helix trajectory if(Fifi<0.) Fifi += 2.*PI; Double_t Kakka ; if( fabs( pMCtr->GetMomentum().Z() )< 1.e-20) Kakka = 99999999.; else Kakka = -carica*0.001*BFIELD*CVEL/pMCtr->GetMomentum().Z(); fprintf(HANDLE, " R_MC %g R %g Fi_MC %g Fi %g KAPPA_MC %g KAPPA %g FI0_MC %g FI0 %g\n", Rr, R[ ii ], Fifi, HoughFiii, Kakka, KAPPA[ ii ], fmod(Fifi+ PI, 2.*PI), // FI0 da MC truth FI0[ ii ] ); //------------------------ // HoughFi = atan2(Oy[i],Ox[i]); // if(HoughFi<0.) HoughFi += 2.*PI; fprintf(HANDLE2,"Evento n. %d Found track %d messa in PndTrackCand",IVOLTE, ii); fprintf(HANDLE2, " R_MC %g R %g Fi_MC %g Fi %g KAPPA_MC %g KAPPA %g FI0_MC %g FI0 %g\n", Rr, R[ ii ], Fifi, HoughFiii, Kakka, KAPPA[ ii ], fmod(Fifi+ PI, 2.*PI), FI0[ ii ] ); //---------------------------------- } // end of for (ii=0; ii=1) //--------------ghosts if( istampa>=1){ int NParghost=0, NParhitsghost=0,icc; for(icc=0; icc=1) } // end of if( nMCTracks >0 && nSttTrackCand > 0 && doMcComparison ) //------------------ end section with comparison MC Mvd hits - associated hits to a certain track //------- load the new PndTrackCand ; each track has the STT and the Mvd hits associated //------- also load the new PndTrack ; each track has the STT and the Mvd hits associated Double_t Oxx, Oyy; // Int_t iflaggo; for(ncand=0, ipinco = 0; ncand< nTotalCandidates; ncand++){ if(!keepit[ncand]) continue; // case in which there was no Skew hits and no KAPPA info and that // candidate could not be associated to any Mvd hits --> no KAPPA information! if(ncand1.e-20 ){ Pzini = -CHARGE[ncand]*0.003*BFIELD/KAPPA[ncand]; if(fabs(Pzini) > PMAX) continue; } else { continue; } // PndTrackCand Array loading new((*fSttMvdPndTrackCandArray)[ipinco]) PndTrackCand; PndTrackCand *pTrckCand = (PndTrackCand*) fSttMvdPndTrackCandArray->At(ipinco); TVector3 dirSeed(Pxini,Pyini,Pzini); // momentum direction in starting point qop = CHARGE[ncand]/dirSeed.Mag(); dirSeed.SetMag(1.); pTrckCand->setTrackSeed(posSeed, dirSeed, qop); pTrckCand->setMcTrackId( daTrackFoundaTrackMC[ncand] ); // pTrckCand->sorted=true; for(j=0; j< nTrackCandHit[ncand]; j++){ switch (ListTrackCandHitType[ncand][j]){ case 0: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId(fMvdPixelBranch), (Int_t)ListTrackCandHit[ncand][j],j); break; case 1: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId(fMvdStripBranch), (Int_t)ListTrackCandHit[ncand][j],j); break; case 2: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId(fSttBranch), (Int_t)ListTrackCandHit[ncand][j],j); break; case 3: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId(fSttBranch), (Int_t)ListTrackCandHit[ncand][j],j); break; } } // PndTrack Array loading // last = last hit in found track; first = first hit in found track; // pTrckCand = pointer to the corresponding PndTrackCand. // the first hit if (ListTrackCandHitType[ncand][0] == 0){ // Mvd Pixel Posiz1[0] = XMvdPixel[ ListTrackCandHit[ncand][0] ]; Posiz1[1] = YMvdPixel[ ListTrackCandHit[ncand][0] ]; Posiz1[2] = ZMvdPixel[ ListTrackCandHit[ncand][0] ]; ErrPosition.SetX(sigmaXMvdPixel[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); ErrPosition.SetY(sigmaXMvdPixel[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); ErrPosition.SetZ(sigmaXMvdPixel[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); } else if (ListTrackCandHitType[ncand][0] == 1){ // Mvd Strip Posiz1[0] = XMvdStrip[ ListTrackCandHit[ncand][0] ]; Posiz1[1] = YMvdStrip[ ListTrackCandHit[ncand][0] ]; Posiz1[2] = ZMvdStrip[ ListTrackCandHit[ncand][0] ]; ErrPosition.SetX(sigmaXMvdStrip[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); ErrPosition.SetY(sigmaXMvdStrip[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); ErrPosition.SetZ(sigmaXMvdStrip[ ListTrackCandHit[ncand][0] ]/sqrt(12.)); } else if( ListTrackCandHitType[ncand][0] == 2 ){ // it is a parallel straw hit PndSttInfoXYZParal( info, ListTrackCandHit[ncand][0], Ox[ncand], Oy[ncand], R[ncand], KAPPA[ncand], FI0[ncand], CHARGE[ncand], Posiz1 ); // cases in which the calculation of the position failed, see // PndSttInfoXYZParal method. if(Posiz1[2]<-888888887. || Posiz1[0] < -999999998.) continue; ErrPosition.SetX(0.02); // 200 microns ErrPosition.SetY(0.02); // 200 microns ErrPosition.SetZ(1.); // 1 cm } else if ( ListTrackCandHitType[ncand][0] == 3 ){ // it is a skew straw hit Posiz1[0] = Ox[ncand]+R[ncand]*cos(SchosenSkew[ncand][ ListTrackCandHit[ncand][0] ]); Posiz1[1] = Oy[ncand]+R[ncand]*sin(SchosenSkew[ncand][ ListTrackCandHit[ncand][0] ]); Posiz1[2] = ZchosenSkew[ncand][ ListTrackCandHit[ncand][0] ]; ErrPosition.SetX(0.02); // 200 microns ErrPosition.SetY(0.02); // 200 microns ErrPosition.SetZ(1.); // 1 cm } Position.SetX( Posiz1[0] ); Position.SetY( Posiz1[1] ); Position.SetZ( Posiz1[2] ); versor[0] = Ox[ncand]-Posiz1[0]; versor[1] = Oy[ncand]-Posiz1[1]; Distance = sqrt(versor[0]*versor[0]+versor[1]*versor[1]); versor[0] /= Distance; versor[1] /= Distance; px = -CHARGE[ncand]*Ptras*versor[1]; py = CHARGE[ncand]*Ptras*versor[0]; Momentum.SetX(px); Momentum.SetY(py); Momentum.SetZ(Pzini); ErrMomentum.SetX(0.05*Ptras); // set at 5% all the times. ErrMomentum.SetY(0.05*Ptras); // set at 5% all the times. ErrMomentum.SetZ(0.05*Pzini); // set at 5% all the times. // the plane of this FairTrackParP better is perpendicular to // the momentum direction ddd = Ptras*sqrt(Ptras*Ptras+Pzini*Pzini); if(istampa>=3){ cout<<" evento = "<=2) if(iplotta && IVOLTE=2&&IVOLTE<20){ cout<<"\n\n---------------------------------------------\n"; cout<<" evt. n. "<0.){ if( ultimoangolo[i]> primoangolo[i]) ultimoangolo[i]-=2.*PI; primo=ultimoangolo[i]*180./PI; ultimoangolo[i]=primoangolo[i]*180./PI; primoangolo[i]=primo; }else{ if( ultimoangolo[i]0 && doMcComparison) { for( j=0;jGetX(); MCSkewAloneY[ MCSkewAloneList[i][j] ]=pSttMCPoint[MCSkewAloneList[i][j]]->GetY(); } WriteMacroSttParallelAssociatedHitsandMvdwithMC( Ox[i], Oy[i], R[i], primoangolo[i],ultimoangolo[i], nSttParHitsinTrack[i], ListSttParHitsinTrack, info, i, k, // questo si usa solo per il nome della Macro. daTrackFoundaTrackMC[i], nParalCommon, ParalCommonList, nSpuriParinTrack, ParSpuriList, nMCParalAlone, MCParalAloneList, nMvdPixelHitsinTrack[i], ListMvdPixelHitsinTrack, nMvdStripHitsinTrack[i], ListMvdStripHitsinTrack, nMvdPixelCommon[i], &MvdPixelCommonList[i][0], nMvdPixelSpuriinTrack[i], &MvdPixelSpuriList[i][0], nMCMvdPixelAlone[i], &MCMvdPixelAloneList[i][0], nMvdStripCommon[i], &MvdStripCommonList[i][0], nMvdStripSpuriinTrack[i], &MvdStripSpuriList[i][0], nMCMvdStripAlone[i], &MCMvdStripAloneList[i][0], nSttSkewHitsinTrack[i], ListSttSkewHitsinTrack, &SchosenSkew[i][0], nSkewCommon, SkewCommonList, nMCSkewAlone, MCSkewAloneList ); } // end of if( nSttParHitsinTrack[i]+nMvdPixelHitsinTrack[i]+ // nMvdStripHitsinTrack[i]>0 && doMcComparison) if( nSttSkewHitsinTrack[i]+nMvdPixelHitsinTrack[i]+ nMvdStripHitsinTrack[i]>0 && doMcComparison){ WriteMacroSkewAssociatedHitswithMC( KAPPA[i],FI0[i], Ox[i], Oy[i], R[i], CHARGE[i], info, WDX,WDY,WDZ, i, k, // questo si usa solo per il nome della Macro. nSttSkewHitsinTrack[i], ListSttSkewHitsinTrack, nSkewCommon[i], SkewCommonList, daTrackFoundaTrackMC[i], nMCSkewAlone[i], MCSkewAloneList, nMvdPixelHitsinTrack, ListMvdPixelHitsinTrack, nMvdStripHitsinTrack, ListMvdStripHitsinTrack, nMvdPixelCommon[i], &MvdPixelCommonList[i][0], nMvdPixelSpuriinTrack[i], &MvdPixelSpuriList[i][0], nMCMvdPixelAlone[i], &MCMvdPixelAloneList[i][0], nMvdStripCommon[i], &MvdStripCommonList[i][0], nMvdStripSpuriinTrack[i], &MvdStripSpuriList[i][0], nMCMvdStripAlone[i], &MCMvdStripAloneList[i][0] ); } // end of if( nSttSkewHitsinTrack[i]+nMvdPixelHitsinTrack[i]+ // nMvdStripHitsinTrack[i]>0 && doMcComparison) } // end of for( i= 0; i< nSttTrackCand; i++) i=0; WriteMacroParallelHitsGeneral( nSttHit, info, // nSttTrackCand, nTotalCandidates, keepit, Ox,Oy,R, FI0, ultimoangolo, primoangolo ); /* for(int kk=-200; kk<= 150; kk += 50){ float time = (float) kk; WriteMacroParallelHitsGeneralspecial( time, // backgound time nSttHit, info, nTotalCandidates, Ox,Oy,R, FI0, ultimoangolo, primoangolo ); } */ WriteMacroAllHitsRestanti( nSttHit, nSttParHit, nSttSkewHit, info, keepit, nTotalCandidates, nTrackCandHit, ListTrackCandHit, ListTrackCandHitType ); } // end of if(iplotta) //--------------------- fine plottamenti -------------------------------------------- return; } //---------------------- end of PndSttMvdTracking::Exec //-------------------------------- begin of function PndSttMvdTracking::calculateintersections void PndSttMvdTracking::calculateintersections(Double_t Ox,Double_t Oy, Double_t R,Double_t C0x,Double_t C0y, Double_t C0z,Double_t r,Double_t vx,Double_t vy,Double_t vz, Int_t *STATUS, Double_t* POINTS) { //------------------------------------------- Double_t P1x, P1y, P1z, P2x, P2y, P2z; Double_t AAA, DELTA, ax, ay, aaa; /* INPUTS : Ox, Oy = abscissa and ordinate of the center of the circular trajectory of the particle; R = radius of such trajectory; C0x, C0y, Coz = x, y, z coordinates of a point belonging to the axis of the skewed straw; r = radius of equidrift of such skewed straw; vx, vy, vz = versor of the direction along which the skewed straw lies. OUTPUTS : P1x, P1y, P1z = x, y, z coordinates of the point intersection between the particle trajectory circle and the equidrift cylinder of the skewed straw calculated as a function of theta (first solution); P2x, P2y, P2z = x, y, z coordinates of the point intersection between the particle trajectory circle and the equidrift cylinder of the skewed straw calculated as a function of theta (second solution); P1x, P1y, .... , P2z are stored in the array POINTS[0-5]; the status is stored in *STATUS. */ //-------------------------------------------------------------------------- // from the resolving formula on page 23 of Gianluigi's notes, setting r=0. ax = C0x - Ox; ay = C0y - Oy; DELTA = R*R*(vx*vx+vy*vy) - (vx*ay - vy*ax)*(vx*ay - vy*ax); AAA = vx*vx+vy*vy; if ( DELTA < 0. ) { *STATUS=-2; } else if (AAA == 0.){ *STATUS=-3; } else if (DELTA == 0.){ *STATUS=1; POINTS[0] = C0x - vx*(vx*ax + vy*ay)/AAA; POINTS[1] = C0y - vy*(vx*ax + vy*ay)/AAA; POINTS[2] = C0z - vz*(vx*ax + vy*ay)/AAA; }else { *STATUS=0; DELTA = sqrt(DELTA); POINTS[0] = C0x - vx*(vx*ax + vy*ay - DELTA)/AAA; POINTS[1] = C0y - vy*(vx*ax + vy*ay - DELTA)/AAA; POINTS[2] = C0z - vz*(vx*ax + vy*ay - DELTA)/AAA; POINTS[3] = C0x - vx*(vx*ax + vy*ay + DELTA)/AAA; POINTS[4] = C0y - vy*(vx*ax + vy*ay + DELTA)/AAA; POINTS[5] = C0z - vz*(vx*ax + vy*ay + DELTA)/AAA; } return; } //----------end of function PndSttMvdTracking::calculateintersections //--start of function PndSttMvdTracking::WriteMacroSttParallelAssociatedHitsandMvdwithMC void PndSttMvdTracking::WriteMacroSttParallelAssociatedHitsandMvdwithMC( Double_t Ox,Double_t Oy,Double_t R, Double_t primoangolo, Double_t ultimoangolo, UShort_t Nhits, UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Double_t info[][7], UShort_t iTrack, Int_t iNome, Short_t daSttTrackaMCTrack, UShort_t nParalCommon[MAXTRACKSPEREVENT], UShort_t ParalCommonList[MAXMCTRACKS][nmaxSttHits], UShort_t nSpuriParinTrack[MAXTRACKSPEREVENT], UShort_t ParSpuriList[MAXTRACKSPEREVENT][nmaxSttHits], UShort_t nMCParalAlone[MAXTRACKSPEREVENT], UShort_t MCParalAloneList[MAXTRACKSPEREVENT][nmaxSttHits], UShort_t nMvdPixelHitsAssociatedToSttTra, UShort_t ListPixelHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdPixelHitsInTrack], UShort_t nMvdStripHitsAssociatedToSttTra, UShort_t ListStripHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdStripHitsInTrack], UShort_t nMvdPixelCommon, UShort_t *MvdPixelCommonList, UShort_t nMvdPixelSpuriinTrack, UShort_t *MvdPixelSpuriList, UShort_t nMCMvdPixelAlone, UShort_t *MCMvdPixelAloneList, UShort_t nMvdStripCommon, UShort_t *MvdStripCommonList, UShort_t nMvdStripSpuriinTrack, UShort_t *MvdStripSpuriList, UShort_t nMCMvdStripAlone, UShort_t *MCMvdStripAloneList, UShort_t nSkewHitsinTrack, UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], Double_t *SchosenSkew, UShort_t nSkewCommon[MAXTRACKSPEREVENT], UShort_t SkewCommonList[MAXTRACKSPEREVENT][nmaxSttHits], UShort_t nMCSkewAlone[MAXTRACKSPEREVENT], UShort_t MCSkewAloneList[MAXMCTRACKS][nmaxSttHits] ) { Int_t i, j, i1, ii, index, Kincl, nlow, nup, STATUS; Double_t xmin , xmax, ymin, ymax, dx, dy, diff, d1, d2, delta, deltax, deltay, deltaz, deltaS, factor, zmin, zmax, Smin, Smax, S1, S2, z1, z2, y1, y2,x1,x2, vx1, vy1, vz1, C0x1, C0y1, C0z1, aaa, bbb, ccc, angle, minor, major, distance, Rx, Ry, LL, Aellipsis1, Bellipsis1,fi1, fmin, fmax, offset, step, SkewInclWithRespectToS, zpos, zpos1, zpos2, Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; //---------- parallel straws Macro now char nome[300], nome2[300]; sprintf(nome,"MacroAllHitswithMCEvent%dT%d", IVOLTE,iNome); sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); // fprintf(MACRO,"void %s()\n{\n",nome); fprintf(MACRO,"{\n"); xmin=1.e20; xmax=-1.e20; ymin=1.e20; ymax=-1.e20; for( ii=0; ii< Nhits; ii++) { i = ListHitsinTrack[iTrack][ii] ; if (info[i][0]-info[i][3] < xmin) xmin = info[i][0]-info[i][3]; if (info[i][0]+info[i][3] > xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } for( i=0; i< nSkewHitsinTrack; i++) { ii = ListSkewHitsinTrack[iTrack][i] ; aaa=Ox+R*cos(SchosenSkew[ii]); bbb=Oy+R*sin(SchosenSkew[ii]); if (aaa < xmin) xmin = aaa; if (aaa > xmax) xmax = aaa; if (bbb < ymin) ymin = bbb; if (bbb > ymax) ymax = bbb; } // for( ii=0; ii< nMvdStripHit; ii++) { for( i=0; i< nMvdStripHitsAssociatedToSttTra; i++) { ii = ListStripHitsinTrack[iTrack][i]; if (XMvdStrip[ii] < xmin) xmin = XMvdStrip[ii]; if (XMvdStrip[ii] > xmax) xmax = XMvdStrip[ii] ; if (YMvdStrip[ii] < ymin) ymin = YMvdStrip[ii]; if (YMvdStrip[ii] > ymax) ymax = YMvdStrip[ii]; } // for( ii=0; ii< nMvdPixelHit; ii++) { for( i=0; i< nMvdPixelHitsAssociatedToSttTra; i++) { ii = ListPixelHitsinTrack[iTrack][i]; if (XMvdPixel[ii] < xmin) xmin = XMvdPixel[ii]; if (XMvdPixel[ii] > xmax) xmax = XMvdPixel[ii] ; if (YMvdPixel[ii] < ymin) ymin = YMvdPixel[ii]; if (YMvdPixel[ii] > ymax) ymax = YMvdPixel[ii]; } for( i=0; i< nMCMvdPixelAlone; i++) { ii = MCMvdPixelAloneList[i]; if (XMvdPixel[ii] < xmin) xmin = XMvdPixel[ii]; if (XMvdPixel[ii] > xmax) xmax = XMvdPixel[ii] ; if (YMvdPixel[ii] < ymin) ymin = YMvdPixel[ii]; if (YMvdPixel[ii] > ymax) ymax = YMvdPixel[ii]; } for( i=0; i< nMCMvdStripAlone; i++) { ii = MCMvdStripAloneList[i]; if (XMvdStrip[ii] < xmin) xmin = XMvdStrip[ii]; if (XMvdStrip[ii] > xmax) xmax = XMvdStrip[ii] ; if (YMvdStrip[ii] < ymin) ymin = YMvdStrip[ii]; if (YMvdStrip[ii] > ymax) ymax = YMvdStrip[ii]; } for( ii=0; ii< nMCParalAlone[iTrack]; ii++) { i = MCParalAloneList[iTrack][ii] ; if( info[i][0] < xmin) xmin = info[i][0]; if( info[i][0] > xmax) xmax = info[i][0]; if( info[i][1] < ymin) ymin = info[i][1]; if( info[i][1] > ymax) ymax = info[i][1]; } for( ii=0; ii< nMCSkewAlone[iTrack]; ii++) { i = MCSkewAloneList[iTrack][ii] ; if( info[i][0] < xmin) xmin = info[i][0]; if( info[i][0] > xmax) xmax = info[i][0]; if( info[i][1] < ymin) ymin = info[i][1]; if( info[i][1] > ymax) ymax = info[i][1]; } if( xmin > 0. ) xmin = 0.; if( xmax < 0.) xmax = 0.; if( ymin > 0. ) ymin = 0.; if( ymax < 0.) ymax = 0.; deltax = xmax-xmin; deltay = ymax - ymin; if( deltax > deltay) { ymin -= 0.5*(deltax-deltay); ymax = ymin+ deltax; delta = deltax; } else { xmin -= 0.5*(deltay-deltax); xmax = xmin+ deltay; delta= deltay; } xmax = xmax + delta*0.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); // disegna il BiHexagon destro e sinistro delle inner parallel straws. DrawBiHexagonInMacro( VERTICALGAP, MACRO, RStrawDetectorMin, ApotemaMaxInnerParStraw, 4, // color code, 4= blue. "InnerPar" ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. DrawBiHexagonInMacro( VERTICALGAP, MACRO, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, 2, // color code. "Skew" ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. DrawHexagonCircleInMacro( VERTICALGAP, MACRO, ApotemaMinOuterParStraw, RStrawDetectorMax, 4, // color code. "OuterPar" ); //-------------- fprintf(MACRO,"TEllipse* FoundTrack = new TEllipse(%f,%f,%f,%f,%f,%f);\n" ,Ox,Oy,R,R,primoangolo,ultimoangolo); fprintf(MACRO, "FoundTrack->SetLineColor(2);\nFoundTrack->SetFillStyle(0);\nFoundTrack->Draw(\"only\");\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"); //------------- hits paralleli in comune con traccia MC for( ii=0; ii< nParalCommon[iTrack]; ii++) { i = ParalCommonList[iTrack][ii] ; fprintf(MACRO, "TEllipse* CommonParalHit%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nCommonParalHit%d->SetFillStyle(0);\nCommonParalHit%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } //------------- hits paralleli spuri for( ii=0; ii< nSpuriParinTrack[iTrack]; ii++) { i = ParSpuriList[iTrack][ii] ; fprintf(MACRO, "TEllipse* SpurParalHit%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nSpurParalHit%d->SetFillStyle(0);\nSpurParalHit%d->SetLineColor(2);\nSpurParalHit%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i,i); } //------------- hits paralleli MC 'alone' for( ii=0; ii< nMCParalAlone[iTrack]; ii++) { i = MCParalAloneList[iTrack][ii] ; fprintf(MACRO, "TEllipse* AloneParalHit%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nAloneParalHit%d->SetFillStyle(0);\nAloneParalHit%d->SetLineColor(4);\nAloneParalHit%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i,i); } //------------- //------------- hits skew in comune e spuri di traccia MC for( i=0; i< nSkewHitsinTrack; i++) { ii = ListSkewHitsinTrack[iTrack][i]; aaa=Ox+R*cos(SchosenSkew[ii]); bbb=Oy+R*sin(SchosenSkew[ii]); for( int k=0; kSetMarkerColor(1);\nCommonSkewHit%d->Draw();\n" ,ii,ii); goto punco ; } } fprintf(MACRO,"TMarker* SpurSkewHit%d = new TMarker(%f,%f,%d);\n", ii,aaa,bbb,28); fprintf(MACRO,"SpurSkewHit%d->SetMarkerColor(2);\nSpurSkewHit%d->Draw();\n",ii,ii); punco: ; } //------------- hits paralleli MC 'alone' for( ii=0; ii< nMCSkewAlone[iTrack]; ii++) { i = MCSkewAloneList[iTrack][ii] ; fprintf(MACRO, "TMarker* AloneSkewHit%d = new TMarker(%f,%f,%d);\nAloneSkewHit%d->SetMarkerColor(4);\nAloneSkewHit%d->Draw();\n", i,MCSkewAloneX[i],MCSkewAloneY[i],28,i,i); } //------------- now the Strips for( i=0; i< nMvdStripHitsAssociatedToSttTra; i++) { ii = ListStripHitsinTrack[iTrack][i]; x1= XMvdStrip[ii]-sigmaXMvdStrip[ii]; x2= XMvdStrip[ii]+sigmaXMvdStrip[ii]; y1= YMvdStrip[ii]-sigmaYMvdStrip[ii]; y2= YMvdStrip[ii]+sigmaYMvdStrip[ii]; for( int k=0; kSetMarkerColor(1);\nCommonMvdStrip%d->Draw();\n", ii,ii); goto pinco ; } } fprintf(MACRO,"TMarker* SpurMvdStrip%d = new TMarker(%f,%f,%d);\n", ii,XMvdStrip[ii],YMvdStrip[ii],25); fprintf(MACRO,"SpurMvdStrip%d->SetMarkerColor(2);\nSpurMvdStrip%d->Draw();\n",ii,ii); pinco: ; } //-------------- hit Mvd Strips 'Alone' for( ii=0; ii< nMCMvdStripAlone; ii++) { i = MCMvdStripAloneList[ii] ; fprintf(MACRO, "TMarker* AloneMvdStrip%d = new TMarker(%f,%f,%d);\nAloneMvdStrip%d->SetMarkerColor(4);\nAloneMvdStrip%d->Draw();\n", i,XMvdStrip[i],YMvdStrip[i],25,i,i); } //------------- // for( ii=0; ii< nMvdPixelHit; ii++) { for( i=0; i< nMvdPixelHitsAssociatedToSttTra; i++) { ii = ListPixelHitsinTrack[iTrack][i]; x1= XMvdPixel[ii]-sigmaXMvdPixel[ii]; x2= XMvdPixel[ii]+sigmaXMvdPixel[ii]; y1= YMvdPixel[ii]-sigmaYMvdPixel[ii]; y2= YMvdPixel[ii]+sigmaYMvdPixel[ii]; // fprintf(MACRO,"TBox* BP%d = new TBox(%f,%f,%f,%f);\nBP%d->SetFillColor(4);\nBP%d->Draw();\n", // ii,x1,y1,x2,y2,ii,ii); for( int k=0; kSetMarkerColor(1);\nCommonMvdPixel%d->Draw();\n", ii,ii); goto panco ; } } fprintf(MACRO,"TMarker* SpurMvdPixel%d = new TMarker(%f,%f,%d);\n", ii,XMvdPixel[ii],YMvdPixel[ii],26); fprintf(MACRO,"SpurMvdPixel%d->SetMarkerColor(2);\nSpurMvdPixel%d->Draw();\n",ii,ii); panco: ; } //-------------- hit Mvd Pixel 'Alone' for( ii=0; ii< nMCMvdPixelAlone; ii++) { i = MCMvdPixelAloneList[ii] ; fprintf(MACRO, "TMarker* AloneMvdPixel%d = new TMarker(%f,%f,%d);\nAloneMvdPixel%d->SetMarkerColor(4);\nAloneMvdPixel%d->Draw();\n", i,XMvdPixel[i],YMvdPixel[i],26,i,i); } //------------- //----------------- ora la traccia MC corrispondente a questa traccia Stt if( daSttTrackaMCTrack>-1 ) { Int_t icode, im; Double_t Rr, Dd, Fifi, Oxx, Oyy, Cx, Cy, Px, Py, carica ; PndMCTrack* pMC; im=daSttTrackaMCTrack; pMC = (PndMCTrack*) fMCTrackArray->At(im); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if (fabs(carica)<0.1 ) goto fuori ; Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); fprintf(MACRO, //"TEllipse* MC%d = new TEllipse(%f,%f,%f,%f,%f,%f);\nMC%d->SetFillStyle(0);\nMC%d->SetLineColor(3);\nMC%d->Draw(\"only\");\n", // im,Cx,Cy,Rr,Rr,primoangolo,ultimoangolo,im,im,im); "TEllipse* MC%d = new TEllipse(%f,%f,%f,%f,%f,%f);\nMC%d->SetFillStyle(0);\nMC%d->SetLineColor(3);\nMC%d->Draw(\"only\");\n", im,Cx,Cy,Rr,Rr,0.,360.,im,im,im); } }; //----------- fine parte del MC fuori: ; fprintf(MACRO,"}\n"); fclose(MACRO); return ; } //--end of function PndSttMvdTracking::WriteMacroSttParallelAssociatedHitsandMvdwithMC //----------start of function PndSttMvdTracking::WriteMacroParallelHitsGeneral void PndSttMvdTracking::WriteMacroParallelHitsGeneral( Int_t Nhits, Double_t info[][7], UShort_t nTracksFoundSoFar, bool *keepit, Double_t *Ox, Double_t *Oy, Double_t *Radius, Double_t *FI0, Double_t *ultimoangolo, Double_t *primoangolo ) { Int_t i, j, i1, ii, index, Kincl, nlow, nup, STATUS; Double_t xmin , xmax, ymin, ymax, xl, xu, yl, yu, gamma, dx, dy, diff, d1, d2, delta, deltax, deltay, deltaz, deltaS, factor,ff, zmin, zmax, Smin, Smax, S1, S2, z1, z2, y1, y2,x1,x2, vx1, vy1, vz1, C0x1, C0y1, C0z1, aaa, bbb, ccc, rrr, angle, minor, major, distance, Rx, Ry, LL, Aellipsis1, Bellipsis1,fi1, fmin, fmax, offset, step, SkewInclWithRespectToS, zpos, zpos1, zpos2, Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; char nome[300], nome2[300]; //---------- parallel straws Macro now sprintf(nome,"MacroAllHitsEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); xmin=1.e20; xmax=-1.e20; ymin=1.e20; ymax=-1.e20; for( i=0; i< Nhits; i++) { // all straws, anche le skew if (info[i][0]-info[i][3] < xmin) xmin = info[i][0]-info[i][3]; if (info[i][0]+info[i][3] > xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } for( ii=0; ii< nMvdPixelHit; ii++) { if (XMvdPixel[ii] < xmin) xmin = XMvdPixel[ii]; if (XMvdPixel[ii] > xmax) xmax = XMvdPixel[ii] ; if (YMvdPixel[ii] < ymin) ymin = YMvdPixel[ii]; if (YMvdPixel[ii] > ymax) ymax = YMvdPixel[ii]; } for( ii=0; ii< nMvdStripHit; ii++) { if (XMvdStrip[ii] < xmin) xmin = XMvdStrip[ii]; if (XMvdStrip[ii] > xmax) xmax = XMvdStrip[ii] ; if (YMvdStrip[ii] < ymin) ymin = YMvdStrip[ii]; if (YMvdStrip[ii] > ymax) ymax = YMvdStrip[ii]; } 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; ymin=xmin=-1.05*RStrawDetectorMax; ymax=xmax= 1.05*RStrawDetectorMax; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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"); // disegna il BiHexagon destro e sinistro delle inner parallel straws. DrawBiHexagonInMacro( VERTICALGAP, MACRO, RStrawDetectorMin, ApotemaMaxInnerParStraw, 4, // color code, 4= blue. "InnerPar" ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. DrawBiHexagonInMacro( VERTICALGAP, MACRO, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, 2, // color code. "Skew" ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. DrawHexagonCircleInMacro( VERTICALGAP, MACRO, ApotemaMinOuterParStraw, RStrawDetectorMax, 4, // color code. "OuterPar" ); //-------------- for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO,"TEllipse* Paral%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nParal%d->SetFillStyle(0);\nParal%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } else { // skew straws. fprintf(MACRO,"TMarker* Skew%d = new TMarker(%f,%f,%d);\nSkew%d->SetMarkerColor(1);\nSkew%d->Draw();\n", i,info[i][0],info[i][1],28,i,i,i); } } for( ii=0; ii< nMvdStripHit; ii++) { x1= XMvdStrip[ii]-sigmaXMvdStrip[ii]; x2= XMvdStrip[ii]+sigmaXMvdStrip[ii]; y1= YMvdStrip[ii]-sigmaYMvdStrip[ii]; y2= YMvdStrip[ii]+sigmaYMvdStrip[ii]; fprintf(MACRO, "TMarker* Strip%d = new TMarker(%f,%f,%d);\nStrip%d->SetMarkerColor(1);\nStrip%d->Draw();\n", ii,XMvdStrip[ii],YMvdStrip[ii],25,ii,ii,ii); } for( ii=0; ii< nMvdPixelHit; ii++) { x1= XMvdPixel[ii]-sigmaXMvdPixel[ii]; x2= XMvdPixel[ii]+sigmaXMvdPixel[ii]; y1= YMvdPixel[ii]-sigmaYMvdPixel[ii]; y2= YMvdPixel[ii]+sigmaYMvdPixel[ii]; fprintf(MACRO, "TMarker* Pixel%d = new TMarker(%f,%f,%d);\nPixel%d->SetMarkerColor(1);\nPixel%d->Draw();\n", ii,XMvdPixel[ii],YMvdPixel[ii],26,ii,ii); } //------------------------------- plotting all the tracks found for(i=0, ii=0; iSetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", // i,aaa,bbb,rrr,rrr,i,i,i); "TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,%f,%f);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw(\"only\");\n", ii,aaa,bbb,rrr,rrr,primoangolo[i],ultimoangolo[i],ii,ii,ii); } // ----------- fprintf(MACRO,"}\n"); fclose(MACRO); //------------------------------------------------------------------------------------------------------------ //---------- parallel straws Macro now con anche le tracce MC sprintf(nome,"MacroAllHitswithMCEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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"); // disegna il BiHexagon destro e sinistro delle inner parallel straws. DrawBiHexagonInMacro( VERTICALGAP, MACRO, RStrawDetectorMin, ApotemaMaxInnerParStraw, 4, // color code, 4= blue. "InnerPar" ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. DrawBiHexagonInMacro( VERTICALGAP, MACRO, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, 2, // color code. "Skew" ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. DrawHexagonCircleInMacro( VERTICALGAP, MACRO, ApotemaMinOuterParStraw, RStrawDetectorMax, 4, // color code. "OuterPar" ); //-------------- for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO,"TEllipse* Paral%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nParal%d->SetFillStyle(0);\nParal%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } else { // skew straws. fprintf(MACRO,"TMarker* Skew%d = new TMarker(%f,%f,%d);\nSkew%d->SetMarkerColor(1);\nSkew%d->Draw();\n", i,info[i][0],info[i][1],28,i,i,i); } } for( ii=0; ii< nMvdStripHit; ii++) { x1= XMvdStrip[ii]-sigmaXMvdStrip[ii]; x2= XMvdStrip[ii]+sigmaXMvdStrip[ii]; y1= YMvdStrip[ii]-sigmaYMvdStrip[ii]; y2= YMvdStrip[ii]+sigmaYMvdStrip[ii]; fprintf(MACRO, "TMarker* Strip%d = new TMarker(%f,%f,%d);\nStrip%d->SetMarkerColor(1);\nStrip%d->Draw();\n", ii,XMvdStrip[ii],YMvdStrip[ii],25,ii,ii,ii); } for( ii=0; ii< nMvdPixelHit; ii++) { x1= XMvdPixel[ii]-sigmaXMvdPixel[ii]; x2= XMvdPixel[ii]+sigmaXMvdPixel[ii]; y1= YMvdPixel[ii]-sigmaYMvdPixel[ii]; y2= YMvdPixel[ii]+sigmaYMvdPixel[ii]; // fprintf(MACRO,"TBox* BP%d = new TBox(%f,%f,%f,%f);\nBP%d->SetFillColor(4);\nBP%d->Draw();\n", // ii,x1,y1,x2,y2,ii,ii); fprintf(MACRO, "TMarker* Pixel%d = new TMarker(%f,%f,%d);\nPixel%d->SetMarkerColor(1);\nPixel%d->Draw();\n", ii,XMvdPixel[ii],YMvdPixel[ii],26,ii,ii); } //------------------------------- plotting all the tracks found for(i=0, ii=0; iSetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", // i,aaa,bbb,rrr,rrr,i,i,i); "TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,%f,%f);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw(\"only\");\n", ii,aaa,bbb,rrr,rrr,primoangolo[i],ultimoangolo[i],ii,ii,ii); } // ----------- //----------------- ora le traccia MC for(i=0; iAt(i); TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); if( fabs(carica) > 0.1 ){ Rr = aaa*1000./(BFIELD*CVEL); //R(cm) of Helix of track projected in XY plane; // B = 2 Tesla Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); // calcolo per plottare solo la parte rilevante della traccia MC. primo=alfa0 = atan2(Oyy-Cy, Oxx-Cx); for(j=0;j<90;j++){ newalfa = alfa0 - carica*j*PI/45; newx = Cx + Rr*cos(newalfa); newy = Cy + Rr*sin(newalfa); if(newx > xmax || newx < xmin || newy>ymax||newy ultimo ) { primo = ultimo; ultimo = alfa0;}; goto pippo ; } } primo = 0.; ultimo = 2.*PI; pippo: ; fprintf(MACRO, "TEllipse* MC%d = new TEllipse(%f,%f,%f,%f,%f,%f);\nMC%d->SetFillStyle(0);\nMC%d->SetLineColor(3);\nMC%d->Draw(\"only\");\n", i,Cx,Cy,Rr,Rr,primo*180./PI,ultimo*180./PI,i,i,i); } else { // continuation of if( fabs(carica) > 0.1 ) continue; // per ora non plotto le neutre. // this is a neutral particle (green dashed straight line). if( Oxx < xmin || Oxx > xmax || Oyy < ymin || Oyy > ymax) continue; double time, time1, time2; if(fabs(Px) >1.e-10 ){ time2 = (xmax-Oxx)/Px; time1 = (xmin-Oxx)/Px; if(time1<0. && time2<0.) { continue; }else if(time2>0.&&time1<0.){ time = time2; x2 = xmax; }else if (time1>0. && time2<0.) { time = time1; x2 = xmin; } else{ if( time2>time1 ) { time=time2; x2 = xmax; } else { time=time1; x2 = xmin; } } y2 = time*Py + Oyy; fprintf(MACRO,"TLine* MCneut%d = new TLine(%f,%f,%f,%f);\n", i,Oxx,Oyy,x2,y2); fprintf(MACRO,"MCneut%d->SetLineStyle(2);\n",i); fprintf(MACRO,"MCneut%d->SetLineColor(3);\nMCneut%d->Draw(\"only\");\n" ,i,i); } else if( fabs(Px) >1.e-10 ) { if( fabs(Py) >1.e-10 ) continue; time2 = (ymax-Oyy)/Py; time1 = (ymin-Oyy)/Py; if(time1<0. && time2<0.) { continue; }else if(time2>0.&&time1<0.){ y2 = ymax; }else if (time1>0. && time2<0.) { y2 = ymin; } else{ y2= time2>time1 ? ymax : ymin ; } fprintf(MACRO,"TLine* MCneut%d = new TLine(%f,%f,%f,%f);\n", i,Oxx,Oyy,Oxx,y2); fprintf(MACRO,"MCneut%d->SetLineStyle(2);\n",i); fprintf(MACRO,"MCneut%d->SetLineColor(3);\nMCneut%d->Draw(\"only\");\n" ,i,i); } } // end of if( fabs(carica) > 0.1 ) } // end of if ( pMC ) } // end of for(i=0; i=0.){ i = (Int_t) time; sprintf(nome,"MacroAllHitsTime%dEvent%d",i, IVOLTE); } else { i = (Int_t) -time; sprintf(nome,"MacroAllHitsTime-%dEvent%d",i, IVOLTE); } sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); // fprintf(MACRO,"void %s()\n{\n",nome); fprintf(MACRO,"{\n"); xmin=1.e20; xmax=-1.e20; ymin=1.e20; ymax=-1.e20; for( i=0; i< Nhits; i++) { // all straws, anche le skew if (info[i][0]-info[i][3] < xmin) xmin = info[i][0]-info[i][3]; if (info[i][0]+info[i][3] > xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } for( ii=0; ii< nMvdPixelHit; ii++) { if (XMvdPixel[ii] < xmin) xmin = XMvdPixel[ii]; if (XMvdPixel[ii] > xmax) xmax = XMvdPixel[ii] ; if (YMvdPixel[ii] < ymin) ymin = YMvdPixel[ii]; if (YMvdPixel[ii] > ymax) ymax = YMvdPixel[ii]; } for( ii=0; ii< nMvdStripHit; ii++) { if (XMvdStrip[ii] < xmin) xmin = XMvdStrip[ii]; if (XMvdStrip[ii] > xmax) xmax = XMvdStrip[ii] ; if (YMvdStrip[ii] < ymin) ymin = YMvdStrip[ii]; if (YMvdStrip[ii] > ymax) ymax = YMvdStrip[ii]; } 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); if( time>=0.){ i = (Int_t) time; fprintf(MACRO,"TPaveLabel* pt= new TPaveLabel(%f,%f,%f,%f,\"Bkg; Time shift = %d nsec\");\n", xmin+delta*0.05,ymax-delta*0.15,xmin+delta*0.35,ymax- delta*0.05,i); } else { i = (Int_t) -time; fprintf(MACRO,"TPaveLabel* pt= new TPaveLabel(%f,%f,%f,%f,\"Bkg; Time shift = -%d nsec\");\n", xmin+delta*0.05,ymax-delta*0.15,xmin+delta*0.35,ymax- delta*0.05,i); } fprintf(MACRO,"pt->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( i=0; i< Nhits; i++) { if( info[i][3]+dR >STRAWRADIUS || info[i][3]+dR <0.) continue; if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO, "TEllipse* Paral%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nParal%d->SetFillStyle(0);\nParal%d->Draw();\n", i,info[i][0],info[i][1],info[i][3]+dR,info[i][3]+dR,i,i); } } if(fabs(time)<5.){ for( ii=0; ii< nMvdStripHit; ii++) { x1= XMvdStrip[ii]-sigmaXMvdStrip[ii]; x2= XMvdStrip[ii]+sigmaXMvdStrip[ii]; y1= YMvdStrip[ii]-sigmaYMvdStrip[ii]; y2= YMvdStrip[ii]+sigmaYMvdStrip[ii]; fprintf(MACRO, "TMarker* Strip%d = new TMarker(%f,%f,%d);\nStrip%d->SetMarkerColor(1);\nStrip%d->Draw();\n", ii,XMvdStrip[ii],YMvdStrip[ii],25,ii,ii,ii); } for( ii=0; ii< nMvdPixelHit; ii++) { x1= XMvdPixel[ii]-sigmaXMvdPixel[ii]; x2= XMvdPixel[ii]+sigmaXMvdPixel[ii]; y1= YMvdPixel[ii]-sigmaYMvdPixel[ii]; y2= YMvdPixel[ii]+sigmaYMvdPixel[ii]; fprintf(MACRO, "TMarker* Pixel%d = new TMarker(%f,%f,%d);\nPixel%d->SetMarkerColor(1);\nPixel%d->Draw();\n", ii,XMvdPixel[ii],YMvdPixel[ii],26,ii,ii); } } // end of if(fabs(time)>5.) fprintf(MACRO,"}\n"); fclose(MACRO); //------------------------------------------------------------------------------------------------------------ return ; } //----------end of function PndSttMvdTracking::WriteMacroParallelHitsGeneralspecial //----------start of function PndSttMvdTracking::WriteMacroSkewAssociatedHitswithMC void PndSttMvdTracking::WriteMacroSkewAssociatedHitswithMC( Double_t KAPPA, Double_t FI0, Double_t Ox, Double_t Oy, Double_t R, Short_t charge, Double_t info[][7], Double_t WDX[nmaxSttHits], Double_t WDY[nmaxSttHits], Double_t WDZ[nmaxSttHits], UShort_t iTrack, Int_t iNome, // questo e' per il nome delle Macro solamente. UShort_t nSkewHitsinTrack, UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], UShort_t nSkewCommon, UShort_t SkewCommonList[MAXTRACKSPEREVENT][nmaxSttHits], Short_t daTrackFoundaTrackMC, UShort_t nMCSkewAlone, UShort_t MCSkewAloneList[MAXMCTRACKS][nmaxSttHits], UShort_t nPixelHitsinTrack[MAXTRACKSPEREVENT], // output UShort_t ListPixelHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdPixelHitsInTrack], // output UShort_t nStripHitsinTrack[MAXTRACKSPEREVENT], // output UShort_t ListStripHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdStripHitsInTrack], // output UShort_t nMvdPixelCommon, UShort_t *MvdPixelCommonList, UShort_t nMvdPixelSpuriinTrack, UShort_t *MvdPixelSpuriList, UShort_t nMCMvdPixelAlone, UShort_t *MCMvdPixelAloneList, UShort_t nMvdStripCommon, UShort_t *MvdStripCommonList, UShort_t nMvdStripSpuriinTrack, UShort_t *MvdStripSpuriList, UShort_t nMCMvdStripAlone, UShort_t *MCMvdStripAloneList ) { TDatabasePDG *fdbPDG= TDatabasePDG::Instance(); Int_t i, j, i1, ii, iii, index, Kincl, nlow, nup, STATUS, imc, Nmin, Nmax; Double_t xmin , xmax, ymin, ymax, dx, dy, diff, d1, d2, delta, deltax, deltay, deltaz, deltaS, esse,factor, zmin, zmax, zmin2, zmax2, 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, dist[2], Tiltdirection1[2], zl[200],zu[200], POINTS1[6]; //------------------- skew straws hits Macro now char nome2[300],nome[300]; FILE *MACRO; sprintf(nome, "MacroSZwithMvdwithMCEvent%dT%d", IVOLTE,iNome); sprintf(nome2, "%s.C",nome); MACRO = fopen(nome2,"w"); fprintf(MACRO,"{\n"); //KAPPA = 1./166.67 ; FI0 = 1.5*PI; Smin=zmin = 1.e10; Smax=zmax = -zmin; index=0; for( iii=0; iii< nSkewHitsinTrack; iii++) { i = ListSkewHitsinTrack[iTrack][iii] ; aaa = sqrt(WDX[i]*WDX[i]+WDY[i]*WDY[i]+ WDZ[i]*WDZ[i]); vx1 = WDX[i]/aaa; vy1 = WDY[i]/aaa; vz1 = WDZ[i]/aaa; C0x1 = info[i][0]; C0y1 = info[i][1]; C0z1 = info[i][2]; calculateintersections(Ox,Oy,R,C0x1,C0y1,C0z1,info[i][3], vx1,vy1,vz1, &STATUS,POINTS1); if(STATUS < 0 ) continue ; for( ii=0; ii<2; ii++){ j=3*ii; distance = sqrt( (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) ); Rx = POINTS1[j]-Ox ; // x component Radial vector of cylinder of trajectory Ry = POINTS1[1+j]-Oy ; // y direction Radial vector of cylinder of trajectory aaa = sqrt(Rx*Rx+Ry*Ry); SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ; SkewInclWithRespectToS /= R; bbb = sqrt( SkewInclWithRespectToS*SkewInclWithRespectToS + vz1*vz1); // the tilt direction of this ellipse is (1,0) when major axis along Z direction LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; if( distance >= info[i][4] + Aellipsis1) continue; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; } // end of for( ii=0; ii<2; ii++) } // end of for( i=1; i< Nhits; i++) //------ aggiungo in blu eventuali punti della traccia MC che sono non mecciati for( iii=0; iii< nMCSkewAlone; iii++) { i = MCSkewAloneList[iTrack][iii]; aaa = sqrt(WDX[i]*WDX[i]+WDY[i]*WDY[i]+ WDZ[i]*WDZ[i]); vx1 = WDX[i]/aaa; vy1 = WDY[i]/aaa; vz1 = WDZ[i]/aaa; C0x1 = info[i][0]; C0y1 = info[i][1]; C0z1 = info[i][2]; calculateintersections(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; dist[ii] = sqrt( (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) ); } if(dist[0]>dist[1]) ii=1; else ii=0; distance = dist[ii]; // 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]*1.2 ) 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 LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; // } // end of for( ii=0; ii<2; ii++) } // end of for( iii=0; iii< nMCSkewAlone[imaxima]; iii++) //------------------------------- //------ fine aggiunta in blu eventuali punti della traccia MC che sono non mecciati // ora la parte delle Mvd // prima i pixel for(i=0; i ZMvdPixel[ ii ] ) zmin = ZMvdPixel[ ii ]; if( zmax < ZMvdPixel[ii ] ) zmax = ZMvdPixel[ii ]; esse = atan2( YMvdPixel[ ii ]-Oy, XMvdPixel[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; } // poi le strip for(i=0; i ZMvdStrip[ ii ] ) zmin = ZMvdStrip[ ii ]; if( zmax < ZMvdStrip[ii ] ) zmax = ZMvdStrip[ii ]; esse = atan2( YMvdStrip[ ii ]-Oy, XMvdStrip[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; } // ora i pixel 'Alone' for(i=0; i ZMvdPixel[ ii ] ) zmin = ZMvdPixel[ ii ]; if( zmax < ZMvdPixel[ii ] ) zmax = ZMvdPixel[ii ]; esse = atan2( YMvdPixel[ ii ]-Oy, XMvdPixel[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; } // ora le strip 'Alone' for(i=0; i ZMvdStrip[ ii ] ) zmin = ZMvdStrip[ ii ]; if( zmax < ZMvdStrip[ii ] ) zmax = ZMvdStrip[ii ]; esse = atan2( YMvdStrip[ ii ]-Oy, XMvdStrip[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; } if( zmax < zmin ) goto nohits ; if( Smax < Smin ) goto nohits; aaa = Smax-Smin; Smin -= aaa*0.2; Smax += aaa*0.2; aaa = zmax-zmin; zmin -= aaa*0.05; zmax += aaa*0.05; if(Smax > 2.*PI) Smax = 2.*PI; if( Smin < 0.) Smin = 0.; // Smin -= 10.; // Smax += 10.; deltaz = zmax-zmin; deltaS = Smax-Smin; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n", zmin-0.1*deltaz,R*(Smin-.1*deltaS),zmax+0.1*deltaz,R*(Smax+0.1*deltaS)); // fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",zmin,Smin,zmax,Smax); for( ii=0; ii< index; ii++) { fprintf(MACRO,"E%d->Draw();\n",ii); } fprintf(MACRO,"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); zmin-0.01*deltaz,R*(Smin+0.05*deltaS),zmax+0.01*deltaz,R*(Smin+0.05*deltaS), zmin-0.01*deltaz,zmax+0.01*deltaz); fprintf(MACRO,"Assex->SetTitle(\"Z (cm)\");\n"); fprintf(MACRO,"Assex->Draw();\n"); fprintf(MACRO,"TGaxis *Assey = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", zmin+0.05*deltaz,R*(Smin-0.01*deltaS),zmin+0.05*deltaz,R*(Smax+0.01*deltaS), R*(Smin-0.01*deltaS),R*(Smax+0.01*deltaS)); fprintf(MACRO,"Assey->SetTitle(\"Helix crf (cm)\");\n"); fprintf(MACRO,"Assey->Draw();\n"); // -------------------------------- for( iii=0; iii< nSkewHitsinTrack; iii++) { i = ListSkewHitsinTrack[iTrack][iii] ; aaa = sqrt(WDX[i]*WDX[i]+WDY[i]*WDY[i]+ WDZ[i]*WDZ[i]); vx1 = WDX[i]/aaa; vy1 = WDY[i]/aaa; vz1 = WDZ[i]/aaa; C0x1 = info[i][0]; C0y1 = info[i][1]; C0z1 = info[i][2]; calculateintersections(Ox,Oy,R,C0x1,C0y1,C0z1,info[i][3], vx1,vy1,vz1, &STATUS,POINTS1); if(STATUS < 0 ) continue ; for( ii=0; ii<2; ii++){ j=3*ii; distance = sqrt( (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) ); Rx = POINTS1[j]-Ox ; // x component Radial vector of cylinder of trajectory Ry = POINTS1[1+j]-Oy ; // y direction Radial vector of cylinder of trajectory aaa = sqrt(Rx*Rx+Ry*Ry); SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ; SkewInclWithRespectToS /= R; bbb = sqrt( SkewInclWithRespectToS*SkewInclWithRespectToS + vz1*vz1); // the tilt direction of this ellipse is (1,0) when major axis along Z direction if( bbb > 1.e-10){ Tiltdirection1[0] = vz1/bbb; Tiltdirection1[1] = SkewInclWithRespectToS/bbb; } else { Tiltdirection1[0] = 1.; Tiltdirection1[1] = 0.; } LL = fabs(vx1*Rx + vy1*Ry); if( LL < 1.e-10) continue; Aellipsis1 = info[i][3]*aaa/LL; Bellipsis1 = info[i][3]/R; if( distance >= info[i][4] + Aellipsis1) continue; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; Double_t rotation1 = 180.*atan2(Tiltdirection1[1],Tiltdirection1[0])/PI; fprintf(MACRO,"TEllipse* Skew%d_%d = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\nSkew%d_%d->SetFillStyle(0);\n", // i,ii,POINTS1[j+2],fi1,Aellipsis1,Bellipsis1,rotation1,i,ii); i,ii,POINTS1[j+2],R*fi1,Aellipsis1,R*Bellipsis1,rotation1,i,ii); // ------ se lo hit e' spurio marcalo in rosso for( i1=0; i1SetLineColor(2);\n",i,ii); fuori: ; fprintf(MACRO,"Skew%d_%d->Draw();\n",i,ii); index++; } // end of for( ii=0; ii<2; ii++) } // end of for( i=1; i< Nhits; i++) //------ aggiungo in blu eventuali punti della traccia MC che sono non mecciati for( iii=0; iii< nMCSkewAlone; iii++) { i = MCSkewAloneList[iTrack][iii]; aaa = sqrt(WDX[i]*WDX[i]+WDY[i]*WDY[i]+ WDZ[i]*WDZ[i]); vx1 = WDX[i]/aaa; vy1 = WDY[i]/aaa; vz1 = WDZ[i]/aaa; C0x1 = info[i][0]; C0y1 = info[i][1]; C0z1 = info[i][2]; calculateintersections(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; dist[ii] = sqrt( (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) ); } if(dist[0]>dist[1]) ii=1; else ii=0; distance = dist[ii]; j=3*ii; Rx = POINTS1[j]-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; //-------------------------- fi1 = atan2(POINTS1[j+1]-Oy, POINTS1[j]-Ox) ; // atan2 returns radians in (-pi and +pi] if( fi1 < 0.) fi1 += 2.*PI; if( zmin > POINTS1[j+2] - Aellipsis1 ) zmin = POINTS1[j+2] - Aellipsis1; if( zmax < POINTS1[j+2] + Aellipsis1 ) zmax = POINTS1[j+2] + Aellipsis1; if( Smin > fi1 - Bellipsis1 ) Smin = fi1 - Bellipsis1; if( Smax < fi1 + Bellipsis1 ) Smax = fi1 + Bellipsis1; Double_t rotation1 = 180.*atan2(Tiltdirection1[1],Tiltdirection1[0])/PI; fprintf(MACRO,"TEllipse* AloneSkew%d_%d = new TEllipse(%f,%f,%f,%f,0.,360.,%f);\nAloneSkew%d_%d->SetFillStyle(0);\n", // i,ii,POINTS1[j+2],fi1,Aellipsis1,Bellipsis1,rotation1,i,ii); i,ii,POINTS1[j+2],R*fi1,Aellipsis1,R*Bellipsis1,rotation1,i,ii); // ------ marca lo hit in blu fprintf(MACRO,"AloneSkew%d_%d->SetLineColor(4);\n",i,ii); fprintf(MACRO,"AloneSkew%d_%d->Draw();\n",i,ii); index++; } // end of for( iii=0; iii< nMCSkewAlone[imaxima]; iii++) //------------------------------- //------ fine aggiunta in blu eventuali punti della traccia MC che sono non mecciati // ora la parte delle Mvd // prima i pixel for(i=0; i ZMvdPixel[ ii ] ) zmin = ZMvdPixel[ ii ]; if( zmax < ZMvdPixel[ii ] ) zmax = ZMvdPixel[ii ]; esse = atan2( YMvdPixel[ ii ]-Oy, XMvdPixel[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; for( int k=0; kSetMarkerColor(1);\n", // ii,ZMvdPixel[ii],esse,26,ii); ii,ZMvdPixel[ii],R*esse,26,ii); fprintf(MACRO,"CommonPixel%d->Draw();\n",ii); goto punco ; } } fprintf(MACRO,"TMarker* SpuriousPixel%d = new TMarker(%f,%f,%d);\nSpuriousPixel%d->SetMarkerColor(2);\n", // ii,ZMvdPixel[ii],esse,26,ii); ii,ZMvdPixel[ii],R*esse,26,ii); fprintf(MACRO,"SpuriousPixel%d->Draw();\n",ii); punco: ; } // poi le strip for(i=0; i ZMvdStrip[ ii ] ) zmin = ZMvdStrip[ ii ]; if( zmax < ZMvdStrip[ii ] ) zmax = ZMvdStrip[ii ]; esse = atan2( YMvdStrip[ ii ]-Oy, XMvdStrip[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; for( int k=0; kSetMarkerColor(1);\n", // ii,ZMvdStrip[ii],esse,25,ii); ii,ZMvdStrip[ii],R*esse,25,ii); fprintf(MACRO,"CommonStrip%d->Draw();\n",ii); goto ponco ; } } fprintf(MACRO,"TMarker* SpuriousStrip%d = new TMarker(%f,%f,%d);\nSpuriousStrip%d->SetMarkerColor(2);\n", // ii,ZMvdStrip[ii],esse,25,ii); ii,ZMvdStrip[ii],R*esse,25,ii); fprintf(MACRO,"SpuriousStrip%d->Draw();\n",ii); ponco: ; } // ora i pixel 'Alone' for(i=0; i ZMvdPixel[ ii ] ) zmin = ZMvdPixel[ ii ]; if( zmax < ZMvdPixel[ii ] ) zmax = ZMvdPixel[ii ]; esse = atan2( YMvdPixel[ ii ]-Oy, XMvdPixel[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; fprintf(MACRO,"TMarker* AlonePixel%d = new TMarker(%f,%f,%d);\nAlonePixel%d->SetMarkerColor(4);\n", // ii,ZMvdPixel[ii],esse,26,ii); ii,ZMvdPixel[ii],R*esse,26,ii); fprintf(MACRO,"AlonePixel%d->Draw();\n",ii); } // ora le strip 'Alone' for(i=0; i ZMvdStrip[ ii ] ) zmin = ZMvdStrip[ ii ]; if( zmax < ZMvdStrip[ii ] ) zmax = ZMvdStrip[ii ]; esse = atan2( YMvdStrip[ ii ]-Oy, XMvdStrip[ ii ]-Ox); if(esse<0.) esse +=2.*PI; if( Smin > esse ) Smin = esse; if( Smax < esse ) Smax = esse; fprintf(MACRO,"TMarker* AloneStrip%d = new TMarker(%f,%f,%d);\nAloneStrip%d->SetMarkerColor(4);\n", // ii,ZMvdStrip[ii],esse,25,ii); ii,ZMvdStrip[ii],R*esse,25,ii); fprintf(MACRO,"AloneStrip%d->Draw();\n",ii); } // plot della traccia trovata dal finder zmin2=zmin; zmax2=zmax; if( -KAPPA*charge>0.) { // Pz>0. if( zmax <0.) { cout<<"da WriteMacroSkewAssociatedHitswithMC, questa traccia" <<" e' inconsistente col proprio Pz, non plottata!\n"; goto dopp ; } zmin = 0.; } else { // Pz<0. if( zmin >0.) { cout<<"da WriteMacroSkewAssociatedHitswithMC, questa traccia" <<" e' inconsistente col proprio Pz, non plottata!\n"; goto dopp ; } zmax = 0.; } if ( KAPPA >= 0.) { fmin = KAPPA*zmin + FI0; fmax = KAPPA*zmax + FI0; } else { fmax = KAPPA*zmin + FI0; fmin = KAPPA*zmax + FI0; } if( fmax>=0.) { Nmax = (int) (0.5*fmax/ PI); } else { Nmax = ( (int) (0.5*fmax/ PI) ) -1; } if( fmin>=0.) { Nmin = (int) (0.5*fmin/ PI); } else { Nmin = ((int) (0.5*fmin/ PI) )-1; } if(fabs(KAPPA)<1.e-10) { cout<<"da WriteMacroSkewAssociatedHitswithMC, questa traccia Found da PR non plottata" <<" perche' ha fabs(KAPPA)<1.e-10.\n"; goto dopp ; } 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); i-Nmin,z1,0.,z2, R*2.*PI,i-Nmin,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) dopp: ; zmin=zmin2; zmax=zmax2; //----------------- ora la traccia MC corrispondente a questa traccia Stt // for(imc=0; imc-1 ) { Int_t icode ; Double_t Rr, Dd, Fifi, Kakka, Oxx, Oyy, Cx, Cy, Px, Py, carica ; PndMCTrack* pMC; pMC = (PndMCTrack*) fMCTrackArray->At(imc); if ( pMC ) { icode = pMC->GetPdgCode() ; // PDG code of track Oxx = pMC->GetStartVertex().X(); // X of starting point track Oyy = pMC->GetStartVertex().Y(); // Y of starting point track Px = pMC->GetMomentum().X(); Py = pMC->GetMomentum().Y(); aaa = sqrt( Px*Px + Py*Py); Rr = aaa*1000./(BFIELD*CVEL); // R (cm) of Helix of track projected in XY plane; B = 2 Tesla TParticlePDG *fParticle= fdbPDG->GetParticle(icode); if (icode>1000000000) carica = 1.; else carica = fParticle->Charge()/3. ; // charge of track if (fabs(carica)<0.1 ) goto pinco ; Cx = Oxx + Py*1000./(BFIELD*CVEL*carica); Cy = Oyy - Px*1000./(BFIELD*CVEL*carica); Fifi = atan2(Cy, Cx); // MC truth Fifi angle of circle of Helix trajectory if(Fifi<0.) Fifi += 2.*PI; if( fabs( pMC->GetMomentum().Z() )< 1.e-20) Kakka = 99999999.; else Kakka = -carica*0.001*BFIELD*CVEL/pMC->GetMomentum().Z(); KAPPA=Kakka; FI0 = fmod(Fifi+ PI, 2.*PI); if ( KAPPA >= 0.) { fmin = KAPPA*zmin + FI0; fmax = KAPPA*zmax + FI0; } else { fmax = KAPPA*zmin + FI0; fmin = KAPPA*zmax + FI0; } if( fmax>=0.) { Nmax = (int) (0.5*fmax/ PI); } else { Nmax = ( (int) (0.5*fmax/ PI) ) -1; } if( fmin>=0.) { Nmin = (int) (0.5*fmin/ PI); } else { Nmin = ((int) (0.5*fmin/ PI) )-1; } for(i=Nmin; i<= Nmax;i++){ offset = 2.*PI*i; z1 = (i*2.*PI-FI0)/KAPPA; z2 = ((i+1)*2.*PI-FI0)/KAPPA; fprintf(MACRO,"TLine* MC%d_%d = new TLine(%f,%f,%f,%f);\nMC%d_%d->SetLineColor(3);\nMC%d_%d->Draw();\n", // imc,i-Nmin,z1,0.,z2, 2.*PI,imc,i-Nmin,imc,i-Nmin); imc,i-Nmin,z1,0.,z2,R* 2.*PI,imc,i-Nmin,imc,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) } // end of if ( pMC ) } // end of if( imc>-1 ) pinco: ; nohits: ; fprintf(MACRO,"}\n"); fclose(MACRO); } //----------end of function PndSttMvdTracking::WriteMacroSkewAssociatedHitswithMC //----------begin of function PndSttMvdTracking::WriteMacroAllHitsRestanti void PndSttMvdTracking::WriteMacroAllHitsRestanti( UShort_t nSttHit, UShort_t nSttParHit, UShort_t nSttSkewHit, Double_t info[][7], bool *keepit, UShort_t nSttTrackCand, UShort_t nCandHit[MAXTRACKSPEREVENT], UShort_t ListCandHit[MAXTRACKSPEREVENT][nmaxSttHitsInTrack+ nmaxMvdPixelHitsInTrack+ nmaxMvdStripHitsInTrack], Short_t ListCandHitType[MAXTRACKSPEREVENT][nmaxSttHitsInTrack+ nmaxMvdPixelHitsInTrack+ nmaxMvdStripHitsInTrack] ) { // nSttHit = parallel+skew. bool exclusionStt[nmaxSttHits], exclusionPixel[nmaxMvdPixelHits], exclusionStrip[nmaxMvdStripHits]; char nome[300], nome2[300]; int i,j,k; Double_t delta, deltax, deltay, xmin, xmax, ymin, ymax; for(i=0;i xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } for( i=0; i< nMvdPixelHit; i++) { if (XMvdPixel[i] < xmin) xmin = XMvdPixel[i]; if (XMvdPixel[i] > xmax) xmax = XMvdPixel[i]; if (YMvdPixel[i] < ymin) ymin = YMvdPixel[i]; if (YMvdPixel[i] > ymax) ymax = YMvdPixel[i]; } for( i=0; i< nMvdStripHit; i++) { if (XMvdStrip[i] < xmin) xmin = XMvdStrip[i]; if (XMvdStrip[i] > xmax) xmax = XMvdStrip[i]; if (YMvdStrip[i] < ymin) ymin = YMvdStrip[i]; if (YMvdStrip[i] > ymax) ymax = YMvdStrip[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; sprintf(nome,"MacroHitsRestantiEvent%d", IVOLTE); sprintf(nome2,"%s.C",nome); FILE * MACRO = fopen(nome2,"w"); fprintf(MACRO,"void %s()\n{\n",nome); fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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"); // disegna il BiHexagon destro e sinistro delle inner parallel straws. DrawBiHexagonInMacro( VERTICALGAP, MACRO, RStrawDetectorMin, ApotemaMaxInnerParStraw, 4, // color code, 4= blue. "InnerPar" ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. DrawBiHexagonInMacro( VERTICALGAP, MACRO, ApotemaMinSkewStraw, ApotemaMaxSkewStraw, 2, // color code. "Skew" ); //-------------- // disegna il BiHexagon destro e sinistro delle skew straws. DrawHexagonCircleInMacro( VERTICALGAP, MACRO, ApotemaMinOuterParStraw, RStrawDetectorMax, 4, // color code. "OuterPar" ); //-------------- for( i=0; i< nSttHit; i++) { if( !exclusionStt[i]) { // all straws if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nE%d->SetFillStyle(0);\nE%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } else { // skew straws. fprintf(MACRO,"TMarker* SS%d = new TMarker(%f,%f,%d);\nSS%d->SetMarkerColor(1);\nSS%d->Draw();\n", i,info[i][0],info[i][1],28,i,i,i); } } } for( i=0; i< nMvdPixelHit; i++) { if( !exclusionPixel[i]) { // all Pixels fprintf(MACRO, "TMarker* Pixel%d = new TMarker(%f,%f,%d);\nPixel%d->SetMarkerColor(1);\nPixel%d->Draw();\n", i,XMvdPixel[i],YMvdPixel[i],26,i,i); } } for( i=0; i< nMvdStripHit; i++) { if( !exclusionStrip[i]) { // all Pixels fprintf(MACRO, "TMarker* Strip%d = new TMarker(%f,%f,%d);\nStrip%d->SetMarkerColor(1);\nStrip%d->Draw();\n", i,XMvdStrip[i],YMvdStrip[i],25,i,i,i); } } fprintf(MACRO,"}\n"); fclose(MACRO); return; } //----------end of function PndSttMvdTracking::WriteMacroAllHitsRestanti //----------begin of function PndSttMvdTracking::DrawBiHexagonInMacro void PndSttMvdTracking::DrawBiHexagonInMacro( Double_t vgap, FILE * MACRO, Double_t Ami, Double_t Ama, UShort_t color, char *name ) { UShort_t iside; // these are the points defining the sides of the biexhagon on the left(inner). Double_t side_x[] = { -vgap/2., -Ama , -Ama, -vgap/2., -vgap/2., -Ami, -Ami, -vgap/2., -vgap/2.}, side_y[] = {(-0.5*vgap+2.*Ama)/sqrt(3.), Ama/sqrt(3.), -Ama/sqrt(3.), -(-0.5*vgap+2.*Ama)/sqrt(3.), -(-0.5*vgap+2.*Ami)/sqrt(3.), -Ami/sqrt(3.), Ami/sqrt(3.), (-0.5*vgap+2.*Ami)/sqrt(3.), (-0.5*vgap+2.*Ama)/sqrt(3.) }; // Inner straws left for(iside=0;iside<8;iside++){ fprintf(MACRO, "TLine* %sL%d = new TLine(%f,%f,%f,%f);\n",name, iside,side_x[iside],side_y[iside],side_x[iside+1],side_y[iside+1]); fprintf(MACRO,"%sL%d->SetLineColor(%d);\n",name,iside,color); fprintf(MACRO,"%sL%d->SetLineStyle(2);\n",name,iside); fprintf(MACRO,"%sL%d->Draw();\n",name,iside); } // Inner straws right for(iside=0;iside<8;iside++){ fprintf(MACRO, "TLine* %sR%d = new TLine(%f,%f,%f,%f);\n",name, iside,-side_x[iside],side_y[iside],-side_x[iside+1],side_y[iside+1]); fprintf(MACRO,"%sR%d->SetLineColor(%d);\n",name,iside,color); fprintf(MACRO,"%sR%d->SetLineStyle(2);\n",name,iside); fprintf(MACRO,"%sR%d->Draw();\n",name,iside); } } //----------end of function PndSttMvdTracking::DrawBiHexagonInMacro //----------begin of function PndSttMvdTracking::DrawHexagonCircleInMacro void PndSttMvdTracking::DrawHexagonCircleInMacro( Double_t GAP, FILE * MACRO, Double_t ApotemaMin, Double_t Rma, UShort_t color, char *name ) { UShort_t iside; Double_t angle1, angle2; // these are the points defining the sides of the exhagon on the left(outer). Double_t side_x[] = { -GAP/2., -GAP/2. , -ApotemaMin, -ApotemaMin, -GAP/2., -GAP/2. }, side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-0.5*GAP)/sqrt(3.), ApotemaMin/sqrt(3.), -ApotemaMin/sqrt(3.), -(2.*ApotemaMin-0.5*GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)}; // Outer straws left for(iside=0;iside<5;iside++){ fprintf(MACRO, "TLine* %sL%d = new TLine(%f,%f,%f,%f);\n",name, iside,side_x[iside],side_y[iside],side_x[iside+1],side_y[iside+1]); fprintf(MACRO,"%sL%d->SetLineColor(%d);\n",name,iside,color); fprintf(MACRO,"%sL%d->SetLineStyle(2);\n",name,iside); fprintf(MACRO,"%sL%d->Draw();\n",name,iside); } // Outer straws right for(iside=0;iside<5;iside++){ fprintf(MACRO, "TLine* %sR%d = new TLine(%f,%f,%f,%f);\n",name, iside,-side_x[iside],side_y[iside],-side_x[iside+1],side_y[iside+1]); fprintf(MACRO,"%sR%d->SetLineColor(%d);\n",name,iside,color); fprintf(MACRO,"%sR%d->SetLineStyle(2);\n",name,iside); fprintf(MACRO,"%sR%d->Draw();\n",name,iside); } // drawing the left circle. angle1 = atan2 ( side_y[0], side_x[0])*180./PI; angle2 = 360. + atan2 ( side_y[4], side_x[4])*180./PI; fprintf(MACRO,"TEllipse* %sCircleL = new TEllipse(0.,0.,%f,%f,%f,%f);\n", name,Rma,Rma,angle1,angle2); fprintf(MACRO,"%sCircleL->SetFillStyle(0);\n",name); fprintf(MACRO,"%sCircleL->SetLineColor(%d);\n",name,color); fprintf(MACRO,"%sCircleL->Draw(\"only\");\n",name,color); // drawing the right circle. angle2 = atan2 ( side_y[0], -side_x[0])*180./PI; angle1 = atan2 ( side_y[4], -side_x[4])*180./PI; fprintf(MACRO,"TEllipse* %sCircleR = new TEllipse(0.,0.,%f,%f,%f,%f);\n", name,Rma,Rma,angle1,angle2); fprintf(MACRO,"%sCircleR->SetFillStyle(0);\n",name); fprintf(MACRO,"%sCircleR->SetLineColor(%d);\n",name,color); fprintf(MACRO,"%sCircleR->Draw(\"only\");\n",name,color); } //----------end of function PndSttMvdTracking::DrawHexagonCircleInMacro //----------begin of function PndSttMvdTracking::DrawSttDetectorInMacro void PndSttMvdTracking::AssociateFoundTrackstoMC( Double_t info[][7], UShort_t nTracksFoundSoFar, UShort_t nHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool inclusionMC[MAXTRACKSPEREVENT][nmaxSttHits], inclusionExp[MAXTRACKSPEREVENT]; UShort_t ntoMCtrack[MAXTRACKSPEREVENT], toMCtracklist[MAXTRACKSPEREVENT][nmaxSttHits], toMCtrackfrequency[MAXTRACKSPEREVENT][nmaxSttHits]; UShort_t i, j, jtemp,jexp; Short_t itemp, massimo; Int_t enne; for(i=0; i -1){ itemp=-1; massimo = -1; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if( !inclusionExp[jexp]) continue; for(i=0; i< ntoMCtrack[jexp]; i++){ if( !inclusionMC[jexp][i]) continue; if( toMCtrackfrequency[jexp][i]>massimo){ massimo=toMCtrackfrequency[jexp][i]; itemp = toMCtracklist[jexp][i]; jtemp = jexp; } } } if( itemp>-1 ){ daTrackFoundaTrackMC[jtemp]=itemp; inclusionExp[jtemp]=false; for(jexp=0; jexp -1) return; } //----------end of function PndSttMvdTracking::AssociateFoundTrackstoMC //----------begin of function PndSttMvdTracking::AssociateFoundTrackstoMCbis void PndSttMvdTracking::AssociateFoundTrackstoMCbis( Double_t info[][7], UShort_t nTracksFoundSoFar, UShort_t nHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool inclusionMC[MAXTRACKSPEREVENT][nmaxSttHits], inclusionExp[MAXTRACKSPEREVENT]; UShort_t ntoMCtrack[MAXTRACKSPEREVENT], toMCtracklist[MAXTRACKSPEREVENT][nmaxSttHits], toMCtrackfrequency[MAXTRACKSPEREVENT][nmaxSttHits]; UShort_t i, j, enne, jtemp,jexp; Short_t itemp, massimo; for(i=0; i -1){ itemp=-1; massimo = -1; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if( !inclusionExp[jexp]) continue; for(i=0; i< ntoMCtrack[jexp]; i++){ if( !inclusionMC[jexp][i]) continue; if( toMCtrackfrequency[jexp][i]>massimo){ massimo=toMCtrackfrequency[jexp][i]; itemp = toMCtracklist[jexp][i]; jtemp = jexp; } } } if( itemp>-1 ){ daTrackFoundaTrackMC[jtemp]=itemp; inclusionExp[jtemp]=false; for(jexp=0; jexp -1) return; } //----------end of function PndSttMvdTracking::AssociateFoundTrackstoMCbis //----------begin of function PndSttMvdTracking::AssociateFoundTrackstoMCtris void PndSttMvdTracking::AssociateFoundTrackstoMCtris( Double_t info[][7], Double_t Ox[MAXTRACKSPEREVENT], Double_t Oy[MAXTRACKSPEREVENT], Double_t R[MAXTRACKSPEREVENT], UShort_t nTracksFoundSoFar, UShort_t nHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool firstime, inclusionMC[MAXTRACKSPEREVENT][nmaxSttHits], inclusionExp[MAXTRACKSPEREVENT]; UShort_t ntoMCtrack[MAXTRACKSPEREVENT], toMCtrackfrequency[MAXTRACKSPEREVENT][nmaxSttHits]; UShort_t i, j, jtemp,jexp , nmid; Short_t itemp, massimo, toMCtracklist[MAXTRACKSPEREVENT][nmaxSttHits]; Int_t enne; Double_t dx, Cx, Cy, Rr, alfa, beta, gamma, minimo, tanlow[MAXTRACKSPEREVENT], tanmid[MAXTRACKSPEREVENT], tanup[MAXTRACKSPEREVENT], toMCtrackdistance[MAXTRACKSPEREVENT][nmaxSttHits]; int nevento=4; for(i=0; i 2) { dx = info[ ListHitsinTrack[i][0] ][0] - Ox[i]; if(fabs(dx)> 1.e-10){ tanlow[i] = ( info[ ListHitsinTrack[i][0] ][1] - Oy[i] )/dx; } else { tanlow[i] = 999999.; } nmid = (UShort_t) (nHitsinTrack[i]+1)/2; dx = info[ ListHitsinTrack[i][nmid] ][0] - Ox[i]; if(fabs(dx)> 1.e-10){ tanmid[i] = ( info[ ListHitsinTrack[i][nmid] ][1] - Oy[i] )/dx; } else { tanmid[i] = 999999.; } dx = info[ ListHitsinTrack[i][nHitsinTrack[i]-1] ][0] - Ox[i]; if(fabs(dx)> 1.e-10){ tanup[i] = ( info[ ListHitsinTrack[i][nHitsinTrack[i]-1] ][1] - Oy[i] )/dx; } else { tanup[i] = 999999.; } } // end of if( nHitsinTrack[i] > 2) //----------------------------- } // end of for(i=0; i -1){ itemp=-1; massimo = -1; minimo = 999999999999.; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if( !inclusionExp[jexp]) continue; for(i=0; i< ntoMCtrack[jexp]; i++){ if( !inclusionMC[jexp][i]) continue; // if( toMCtrackfrequency[jexp][i]>massimo){ if( toMCtrackdistance[jexp][i]<-0.5) continue; if( toMCtrackdistance[jexp][i]-1 ){ daTrackFoundaTrackMC[jtemp]=itemp; inclusionExp[jtemp]=false; for(jexp=0; jexp -1) return; } //----------end of function PndSttMvdTracking::AssociateFoundTrackstoMCtris //----------begin of function PndSttMvdTracking::AssociateFoundTrackstoMCquater void PndSttMvdTracking::AssociateFoundTrackstoMCquater( bool *keepit, Double_t info[][7], Double_t Ox[MAXTRACKSPEREVENT], Double_t Oy[MAXTRACKSPEREVENT], Double_t R[MAXTRACKSPEREVENT], Double_t X1[MAXTRACKSPEREVENT], Double_t Y1[MAXTRACKSPEREVENT], Double_t X2[MAXTRACKSPEREVENT], Double_t Y2[MAXTRACKSPEREVENT], Double_t X3[MAXTRACKSPEREVENT], Double_t Y3[MAXTRACKSPEREVENT], UShort_t nTracksFoundSoFar, UShort_t nHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], UShort_t nPixelHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListPixelHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdPixelHitsInTrack], Short_t *FromPixeltoMCTrack, UShort_t nStripHitsinTrack[MAXTRACKSPEREVENT], UShort_t ListStripHitsinTrack[MAXTRACKSPEREVENT][nmaxMvdStripHitsInTrack], Short_t *FromStriptoMCTrack, Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { bool firstime, // inclusionMC[nTracksFoundSoFar][nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits], // inclusionExp[nTracksFoundSoFar]; inclusionMC[MAXTRACKSPEREVENT][nmaxSttHits+nmaxMvdPixelHits+nmaxMvdStripHits], inclusionExp[MAXTRACKSPEREVENT]; UShort_t ntoMCtrack[MAXTRACKSPEREVENT], toMCtrackfrequency[MAXTRACKSPEREVENT][nmaxSttHits]; UShort_t i, j, jtemp,jexp , nmid; Short_t itemp, massimo, toMCtracklist[MAXTRACKSPEREVENT][nmaxSttHits]; Int_t enne; Double_t dx, Cx, Cy, Rr, alfa, beta, gamma, minimo, tanlow[MAXTRACKSPEREVENT], tanmid[MAXTRACKSPEREVENT], tanup[MAXTRACKSPEREVENT], toMCtrackdistance[MAXTRACKSPEREVENT][nmaxSttHits]; for(i=0; i 2 || nn == 1) { dx = X1[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanlow[i] = (Y1[i] - Oy[i] )/dx; } else { tanlow[i] = 999999.; } dx = X2[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanmid[i] = ( Y2[i] - Oy[i] )/dx; } else { tanmid[i] = 999999.; } dx = X3[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanup[i] = (Y3[i] - Oy[i] )/dx; } else { tanup[i] = 999999.; } } else if (nn==2) { // continuation of if( nn > 2 || nn == 1) dx = X1[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanlow[i] = (Y1[i] - Oy[i] )/dx; } else { tanlow[i] = 999999.; } dx = X2[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanmid[i] = ( Y2[i] - Oy[i] )/dx; } else { tanmid[i] = 999999.; } tanup[i] = tanmid[i]; } /* else{ // continuation of if( nn > 2), case in which there is only 1 hit in the Track cand. dx = X1[i] - Ox[i]; if(fabs(dx)> 1.e-10){ tanlow[i] = (Y1[i] - Oy[i] )/dx; } else { tanlow[i] = 999999.; } tanmid[i] = tanlow[i]; tanup[i] = tanlow[i]; } // end of if( nn > 2) */ //----------------------------- } // end of for(i=0; i -1){ itemp=-1; massimo = -1; minimo = 999999999999.; for(jexp=0; jexp< nTracksFoundSoFar ;jexp++){ if(!keepit[jexp]) continue; if( !inclusionExp[jexp]) continue; for(i=0; i< ntoMCtrack[jexp]; i++){ if( !inclusionMC[jexp][i]) continue; // if( toMCtrackfrequency[jexp][i]>massimo){ if( toMCtrackdistance[jexp][i]<-0.5) continue; if( toMCtrackdistance[jexp][i]-1 ){ daTrackFoundaTrackMC[jtemp]=itemp; inclusionExp[jtemp]=false; for(jexp=0; jexp -1) return; } //----------end of function PndSttMvdTracking::AssociateFoundTrackstoMCquater //------------------------- begin of function PndSttMvdTracking::SttMatchedSpurious void PndSttMvdTracking::SttMatchedSpurious( bool *keepit, UShort_t ntotalHits, Double_t info[][7], UShort_t nTracksFoundSoFar, // quelle trovate dal PR UShort_t nHitsinTrack[MAXTRACKSPEREVENT], // n. hits PARALLELI+SKEW, dal PR UShort_t ListHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], // dal PR UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], // n. hits skew, dal PR UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHitsInTrack], // dal PR UShort_t nParalCommon[MAXTRACKSPEREVENT], UShort_t ParalCommonList[MAXMCTRACKS][nmaxSttHits], UShort_t nSpuriParinTrack[MAXTRACKSPEREVENT], UShort_t ParSpuriList[MAXTRACKSPEREVENT][nmaxSttHits], UShort_t nSkewCommon[MAXTRACKSPEREVENT], UShort_t SkewCommonList[MAXMCTRACKS][nmaxSttHits], UShort_t nSpuriSkewinTrack[MAXTRACKSPEREVENT], UShort_t SkewSpuriList[MAXTRACKSPEREVENT][nmaxSttHits], UShort_t nHitsInMCTrack[MAXTRACKSPEREVENT], UShort_t nSkewHitsInMCTrack[MAXTRACKSPEREVENT], UShort_t nMCParalAlone[MAXTRACKSPEREVENT], UShort_t MCParalAloneList[MAXTRACKSPEREVENT][nmaxSttHits], UShort_t nMCSkewAlone[MAXTRACKSPEREVENT], UShort_t MCSkewAloneList[MAXTRACKSPEREVENT][nmaxSttHits], Short_t daTrackFoundaTrackMC[MAXTRACKSPEREVENT] ) { UShort_t i, jexp, exphit, iHit, enne[MAXTRACKSPEREVENT][nmaxSttHits]; Short_t emme; for(jexp=0; jexp 2. || (emme != daTrackFoundaTrackMC[jexp]) ) continue; if( !ExclusionListStt[i]) continue; // escludo gli hits con multiple hits for(exphit=0; exphit=3 && IVOLTE<20) cout<<"MatchMvdHitsToSttTracks, i= "<=4){ cout<<"n. punti nel fit "<=3 && IVOLTE <= 20){ sprintf(stringa, "/home/boca/panda/glpk/glpk-4.39/examples/glpsol --min -o soluztrack%dEvent%dstep%d GeneralParallelHitsConformeTraccia%dEvent%d.mcs", 0,IVOLTE,1,0,IVOLTE); } else { sprintf(stringa, "/home/boca/panda/glpk/glpk-4.39/examples/glpsol --min -o soluztrack%dEvent%dstep%d GeneralParallelHitsConformeTraccia%dEvent%d.mcs >& /dev/null", 0, IVOLTE,1,0,IVOLTE,1); } */ m1_result = final_values[0]; m2_result = final_values[1]; q1_result = final_values[2]; q2_result = final_values[3]; if(istampa>2) cout<<"Results : m1 = "<nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack) NpointsInFit=nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack; // NpointsInFit =10; bool mvdhit[nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack]; Double_t ave, avex, avey, cose, sine, M = 50., m_result, q_result, A, alfetta, angle, offsety, rotationangle, Ox[nmaxSttHitsInTrack], Oy[nmaxSttHitsInTrack], Delta[nmaxSttHitsInTrack]; UShort_t i, j, ii, iii, n, nSttHits, nMvdHits; Short_t Status; char nome[300], stringa[300], stringa2[300]; float m1_result,m2_result, q1_result,q2_result, A1_result, A2_result; // -- if(nSkewHitsinTrack==0) { cout<<"from FitSZspace, Evento "< 1.e-10&& DriftRadius[ i ]>0.){ n++; ave += (S[ i ] - FInot)/Z[ i ]; avex += Z[ i ]; avey += (S[ i ] - FInot); } } if( n>0) { ave /=n; avex /=n; avey /=n; // rotationangle = atan(ave); rotationangle = atan2(avey,avex); } else { rotationangle=PI/2.; } // use the trick of increasing the rotation angle by 10 degrees in order to obtain always a positive m // rotationangle -= PI/18.; rotationangle = PI/2.; cose = cos(rotationangle); sine = sin(rotationangle); nSttHits = nMvdHits = 0; for(i=0;i 1.e-10) { *emme=((*emme)*cose+sine)/(cose-(*emme)*sine); return 1; } else { // in this case the equation is 0 = U in the Conformal plane --> x=0 in the XY plane. return -99; } } //----------end of function PndSttMvdTracking::FitSZspace //------------------ begin function PndSttMvdTracking::RefitMvdStt void PndSttMvdTracking::RefitMvdStt( UShort_t nCandHit, UShort_t *ListCandHit, Short_t *ListCandHitType, Double_t info[][7], Double_t rotationangle, // this is between 0. and 2*PI Double_t tv[2], Short_t iexcl, Double_t *ALFA, // output of the fit Double_t *BETA, // output of the fit Double_t *GAMMA,// set at zero always for now bool *status // fit status; true = successful ) { bool TypeConf; UShort_t i, iparallel; UShort_t MAXIMUMHITSINFIT = 20; Short_t exitstatus; Double_t dist2, maxdis2, emme, factor, gamma, qu, ErrorStraw = 0.03, ErrorMvd = 0.01, Xconformal[nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], Yconformal[nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], DriftRadiusconformal[nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack], ErrorDriftRadiusconformal[nmaxSttHitsInTrack+nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack]; *status= false; factor=3.; maxdis2=0.5; // trajectory_vertex[0]=trajectory_vertex[1]=0.; for(i=0, iparallel=0; i2){ cout<<"PndSttMvdTracking::RefitMvdStt, n. Hits (Mvd+Stt || ) = "<2) cout<<"distanza**2 di Strip hit n. "< -1){ if( daTrackFoundaTrackMC[i] == FromPixeltoMCTrack[ ListPixelHitsinTrack[i][j] ] ){ index = i*nMvdPixelHit+nMvdPixelCommon[i]; *(MvdPixelCommonList+index) = ListPixelHitsinTrack[i][j]; nMvdPixelCommon[i]++; } else { index = i*nMvdPixelHit+nMvdPixelSpuriinTrack[i]; *(MvdPixelSpuriList+index) = ListPixelHitsinTrack[i][j]; nMvdPixelSpuriinTrack[i]++; } } } for(j=0;j -1){ if( daTrackFoundaTrackMC[i] == FromStriptoMCTrack[ ListStripHitsinTrack[i][j] ] ){ index = i*nMvdStripHit+nMvdStripCommon[i]; *(MvdStripCommonList+index) = ListStripHitsinTrack[i][j]; nMvdStripCommon[i]++; } else { index = i*nMvdStripHit+nMvdStripSpuriinTrack[i]; *(MvdStripSpuriList+index) = ListStripHitsinTrack[i][j]; nMvdStripSpuriinTrack[i]++; } } } } // end of for(i=0; i -1){ for(i=0;i -1){ for(i=0;i0){ // track must rotate clockwise looking into the beam. anglemax = FI0[i]; anglemin = Fifirst[i]; } else { anglemin = FI0[i]; anglemax = Fifirst[i]; } if(anglemax < anglemin) anglemax += 2.*PI; if(anglemax < anglemin) anglemax==anglemin; // this is just to be super-sure. pippo: ; //-------------------- nPixelHitsinTrack[i]=0; nStripHitsinTrack[i]=0; ngoodmix=0; nn[0]=0; for( imvdcand=0; imvdcandGetBranchId(fMvdPixelBranch)){ ncont++; angle = atan2( YMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]]-Oy[i], XMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]]-Ox[i] ); if(angle<0.) angle += 2.*PI; if( angle>anglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax) { dist=fabs( sqrt( (Ox[i]-XMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]])* (Ox[i]-XMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]]) +(Oy[i]-YMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]])* (Oy[i]-YMvdPixel[ListHitMvdTrackCand[imvdcand][jmvdhit]])) -R[i]); if(distGetBranchId(fMvdPixelBranch); Dist += dist; if( dist anglemin) } else {// at this point this is a Strip hit; already made sure // earlier in the code that there is no third possibility. ncont++; angle = atan2( YMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]]-Oy[i], XMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]]-Ox[i] ); if(angle<0.) angle += 2.*PI; if( angle>anglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]])* (Ox[i]-XMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]]) +(Oy[i]-YMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]])* (Oy[i]-YMvdStrip[ListHitMvdTrackCand[imvdcand][jmvdhit]])) -R[i]); if(distGetBranchId(fMvdStripBranch); Dist += dist; if( dist anglemin) } // end of if(ListHitTypeMvdTrackCand[imvdcand][jmvdhit] } // end of for( jmvdhit=0; jmvdhit0) { DIST[ngoodmix]=Dist/nn[ngoodmix]; ngoodmix++; //--------- stampaggi if(istampa>=3 && IVOLTE<20 ){cout<<"\tquesto Mvd candidato (n. ngoodmix = "<GetBranchId(fMvdPixelBranch)) { cout<<"\tPixel hit n. "<GetBranchId(fMvdStripBranch)){ cout<<"\tStrip hit n. "<anglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdPixel[ListMvdDSPixelHitNotTrackCand[jmvdhit]])* (Ox[i]-XMvdPixel[ListMvdDSPixelHitNotTrackCand[jmvdhit]]) +(Oy[i]-YMvdPixel[ListMvdDSPixelHitNotTrackCand[jmvdhit]])* (Oy[i]-YMvdPixel[ListMvdDSPixelHitNotTrackCand[jmvdhit]]) ) -R[i]); if(distGetBranchId(fMvdPixelBranch); DIST[ngoodmix] += dist; if( dist anglemin ) } // end of for( jmvdhit=0; jmvdhitanglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdStrip[ListMvdDSStripHitNotTrackCand[jmvdhit]])* (Ox[i]-XMvdStrip[ListMvdDSStripHitNotTrackCand[jmvdhit]]) +(Oy[i]-YMvdStrip[ListMvdDSStripHitNotTrackCand[jmvdhit]])* (Oy[i]-YMvdStrip[ListMvdDSStripHitNotTrackCand[jmvdhit]]) ) -R[i]); if(distGetBranchId(fMvdStripBranch); DIST[ngoodmix] += dist; if( dist anglemin ) } // end of for( jmvdhit=0; jmvdhit0) { DIST[ngoodmix] /= nn[ngoodmix]; ngoodmix++; //--------- stampaggi if(istampa>2 && IVOLTE<20 ){cout<<"\tevento n. "<GetBranchId(fMvdPixelBranch)) { cout<<"\tDS Pixel hit n. "<GetBranchId(fMvdStripBranch)){ cout<<"\tDS Strip hit n. "<0) //-------- now the US (upstream) Mvd hits nn[ngoodmix]=0; DIST[ngoodmix] = 0.; nHighQuality[ngoodmix]=0; for( jmvdhit=0; jmvdhitanglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdPixel[ListMvdUSPixelHitNotTrackCand[jmvdhit]])* (Ox[i]-XMvdPixel[ListMvdUSPixelHitNotTrackCand[jmvdhit]]) +(Oy[i]-YMvdPixel[ListMvdUSPixelHitNotTrackCand[jmvdhit]])* (Oy[i]-YMvdPixel[ListMvdUSPixelHitNotTrackCand[jmvdhit]]) ) -R[i]); if(distGetBranchId(fMvdPixelBranch); DIST[ngoodmix] += dist; if( dist anglemin ) } // end of for( jmvdhit=0; jmvdhitanglemax){ angle -= 2.*PI; if( angle>anglemax) angle = anglemax; } else if (angle anglemin && angle < anglemax){ dist=fabs( sqrt( (Ox[i]-XMvdStrip[ListMvdUSStripHitNotTrackCand[jmvdhit]])* (Ox[i]-XMvdStrip[ListMvdUSStripHitNotTrackCand[jmvdhit]]) +(Oy[i]-YMvdStrip[ListMvdUSStripHitNotTrackCand[jmvdhit]])* (Oy[i]-YMvdStrip[ListMvdUSStripHitNotTrackCand[jmvdhit]]) ) -R[i]); if(distGetBranchId(fMvdStripBranch); DIST[ngoodmix] += dist; if( dist anglemin ) } // end of for( jmvdhit=0; jmvdhit0) { DIST[ngoodmix] /= nn[ngoodmix]; ngoodmix++; //--------- stampaggi if(istampa>=3 && IVOLTE<20 ){cout<<"\tevento n. "<0) { for(int icc=0; iccGetBranchId(fMvdPixelBranch)) { cout<<"\tUS Pixel hit n. "<GetBranchId(fMvdStripBranch)){ cout<<"\tUS Strip hit n. "<0) //------- end of using the Mvd which are in no Mvd Track Candidate if(istampa>=3 ){cout<<"da PndSttMvdTracking : appena prima arbitration, IVOLTE = "<< IVOLTE<<", Stt track cand = "< Fi_up_limit) continue; } else if( S[NAssociated] > Fi_up_limit) { if( S[NAssociated]- 2.*PI < Fi_low_limit) continue; } //--------------------------- end check Z[NAssociated] = POINTS1[j+2]; ZDrift[NAssociated] = Aellipsis1*Tiltdirection1[0]; ZErrorafterTilt[NAssociated] = 0.02*aaa*Tiltdirection1[0]/LL; SkewList[NAssociated][0] = i; // n. skew hit in original hit numbering SkewList[NAssociated][1] = ii; // solution 0 or solution 1 were accepted // check if this skew hit doesn't "push out" the most external parallel hit // (see Gianluigi's logbook on page 251) /* Double_t Zh1 = Z[NAssociated] - ZDrift[NAssociated]; Double_t Zh2 = Z[NAssociated] + ZDrift[NAssociated]; Double_t Sh1 = S[NAssociated] - Aellipsis1*Tiltdirection1[1]; Double_t Sh2 = S[NAssociated] + Aellipsis1*Tiltdirection1[1]; Double_t Zlast1 = (Fi_final_helix_referenceframe-Fi_initial_helix_referenceframe)*Zh1 /(Sh1-Fi_initial_helix_referenceframe); Double_t Zlast2 = (Fi_final_helix_referenceframe-Fi_initial_helix_referenceframe)*Zh2 /(Sh2-Fi_initial_helix_referenceframe); if( fabs(Zlast1 - ZCENTER_STRAIGHT) > SEMILENGTH_STRAIGHT && fabs(Zlast2 - ZCENTER_STRAIGHT) > SEMILENGTH_STRAIGHT ) continue; */ NAssociated++; } // end of for( ii=0; ii<2; ii++) } // for( iii=0; iii< NSkewhits; iii++) return NAssociated; } //----------end of function PndSttMvdTracking::AssociateSkewHitsToXYTrack //----------begin of function PndSttMvdTracking::CalculateCircleThru3Points bool PndSttMvdTracking::CalculateCircleThru3Points( Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3, Double_t *Ox, Double_t *Oy, Double_t *R) { // inspiration from http://mathworld.wolfram.com/Circle.html Double_t a, d, e, q1, q2, q3; a = x1*y2 + y1*x3 + x2*y3 - x3*y2 - x2*y1 - x1*y3; if( fabs(a)<1.e-10 ) return false; q1 = x1*x1 + y1*y1; q2 = x2*x2 + y2*y2; q3 = x3*x3 + y3*y3; d = -(q1*y2 + y1*q3 + q2*y3 - q3*y2 - q2*y1 - q1*y3); e = -(x1*q2 + q1*x3 + x2*q3 - x3*q2 - x2*q1 - x1*q3); *Ox = -0.5*d/a; *Oy = -0.5*e/a; *R = sqrt( (x1-*Ox)*(x1-*Ox) + (y1-*Oy)*(y1-*Oy) ); return true; } //----------end of function PndSttMvdTracking::CalculateCircleThru3Points //----------begin of function PndSttMvdTracking::Ordering_Loading_ListTrackCandHit void PndSttMvdTracking::Ordering_Loading_ListTrackCandHit( bool *keepit, UShort_t FirstCandidate, UShort_t LastCandidate, Double_t info[][7] ) { UShort_t i, j, ipar, iskew, ncand; // ordering all the hits belonging to the candidate track, by increasing R; // forming the new track with Mvd+Stt hits for(ncand=FirstCandidate; ncand< LastCandidate; ncand++){ if(!keepit[ncand]) continue; nTrackCandHit[ncand] =nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand]+ nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; UShort_t tempmvdindex[nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] ], tempmvdtype[nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] ], auxIndex[nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] ]; Double_t auxR[nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] ]; // adding the Mvd hits (Pixel and Strips) for(i=0; i< nMvdPixelHitsinTrack[ncand]; i++){ auxR[i] = XMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]* XMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]+ YMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]* YMvdPixel[ ListMvdPixelHitsinTrack[ncand][i] ]; tempmvdindex[i]=ListMvdPixelHitsinTrack[ncand][i]; tempmvdtype[i]=0; auxIndex[i] = i; } for(i=0; i< nMvdStripHitsinTrack[ncand]; i++){ auxR[i+nMvdPixelHitsinTrack[ncand]] = XMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]* XMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]+ YMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]* YMvdStrip[ ListMvdStripHitsinTrack[ncand][i] ]; tempmvdindex[i+nMvdPixelHitsinTrack[ncand]]= ListMvdStripHitsinTrack[ncand][i]; tempmvdtype[i+nMvdPixelHitsinTrack[ncand]]=1; auxIndex[i+nMvdPixelHitsinTrack[ncand]]= i+nMvdPixelHitsinTrack[ncand]; } // ordering the Mvd Hits if( nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand] >0){ Merge_Sort( nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand], auxR, auxIndex); // constructing the first part of the ordered new Track Candidate for(i=0; i< nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; i++){ ListTrackCandHit[ncand][i] = tempmvdindex[ auxIndex[i] ]; ListTrackCandHitType[ncand][i] = tempmvdtype[ auxIndex[i] ]; } } // end of if( nMvdPixelHitsinTrack[ncand]+ // construction of the second part of the ordered new Track Candidate UShort_t tempmvdindex2[nSttParHitsinTrack[ncand]+ nSttSkewHitsinTrack[ncand] ], tempmvdtype2[nSttParHitsinTrack[ncand]+ nSttSkewHitsinTrack[ncand] ], auxIndex2[nSttParHitsinTrack[ncand]+ nSttSkewHitsinTrack[ncand] ]; Double_t auxR2[nSttParHitsinTrack[ncand]+ nSttSkewHitsinTrack[ncand] ]; for(i=0; i0){ Merge_Sort( nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand], auxR2, auxIndex2); for(j=0,ipar=0,iskew=0; j< nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand];j++){ i = j+nMvdPixelHitsinTrack[ncand]+ nMvdStripHitsinTrack[ncand]; ListTrackCandHit[ncand][i] = tempmvdindex2[ auxIndex2[j] ]; ListTrackCandHitType[ncand][i] = tempmvdtype2[ auxIndex2[j] ]; if( ListTrackCandHitType[ncand][i]==2) { ListSttParHitsinTrack[ncand][ipar]=tempmvdindex2[auxIndex2[j]]; ipar++; } else { ListSttSkewHitsinTrack[ncand][iskew]=tempmvdindex2[auxIndex2[j]]; iskew++; } } } // end of if( nSttParHitsinTrack[ncand]+ } // end of for(ncand=FirstCandidate; ncand< LastCandidate; ncand++) return; } //----------end of function PndSttMvdTracking::Ordering_Loading_ListTrackCandHit //----------begin of function PndSttMvdTracking::TrackCleanup bool PndSttMvdTracking::TrackCleanup( Double_t GAP, Double_t Oxx, Double_t Oyy, Double_t Rr, Double_t KAPPA, Double_t FI0, Short_t Charge, Double_t Start[3], UShort_t &nHitsPar, // n. hits parall Stt UShort_t *ListHitsPar, UShort_t &nHitsSkew, // n. hits parall Stt UShort_t *ListHitsSkew, Double_t *auxS, Double_t info[][7], Double_t RStrawDetMin, Double_t ApotemaMaxInnerPar, Double_t ApotemaMinSkew, Double_t ApotemaMaxSkew, Double_t ApotemaMinOuterPar, Double_t RStrawDetMax ) { // this method does 3 things : // // 1) finds the entrance and exit points in the STT parallel and skew volumes of the current track; // 2) eliminates from the track hit list possible spurious hits that are not encompassed // by the entrance and exit point; // 3) eliminates the tracks if the hit sequence is not continuous enough. bool ConsiderLastHit; Short_t flagInnerStt, flagOuterStt; UShort_t ihit, nInnerHits, nOuterHits, nIntersections[2], ListInnerHits[nmaxSttHits], ListOuterHits[nmaxSttHits]; Double_t FiLimitAdmissible, r, Xcross[2], Ycross[2], XintersectionList[12][2], // second index =0 --> inner Hexagon, =1 --> outer. YintersectionList[12][2]; // first index : all the possible intersections // (up to 12 intersections). //------------------------ // calculation of the Maximum FI angle possible (if it is a +ve charge) of the Minimum // for this track, taking into account // that the maximum possible Z of a hit is ZCENTER_STRAIGHT + SEMILENGTH_STRAIGHT; the minimum // Z of a hit is ZCENTER_STRAIGHT - SEMILENGTH_STRAIGHT. if(Charge<0){ if( KAPPA>0.){ FiLimitAdmissible = FI0 + KAPPA*(ZCENTER_STRAIGHT + SEMILENGTH_STRAIGHT) ; } else { FiLimitAdmissible = FI0 + KAPPA*(ZCENTER_STRAIGHT - SEMILENGTH_STRAIGHT) ; } } else { if( KAPPA>0.){ FiLimitAdmissible = FI0 + KAPPA*(ZCENTER_STRAIGHT - SEMILENGTH_STRAIGHT) ; } else { FiLimitAdmissible = FI0 + KAPPA*(ZCENTER_STRAIGHT + SEMILENGTH_STRAIGHT) ; } } // end of if(Charge<0) //----------------------------------------------------------------------------------------------------- // parallel cleanup. //----------------stampe if(istampa>=2&&IVOLTE==0){ cout<<" IVOLTE = "<0 && !SttParalCleanup( GAP, Oxx, Oyy, Rr, Charge, Start, FI0, FiLimitAdmissible, nHitsPar, // input and output. ListHitsPar, // input and output. info, RStrawDetMin, ApotemaMaxInnerPar, ApotemaMinOuterPar, RStrawDetMax ) ){ return false; } //---------------------------------------------------------------------------- // skew cleanup. if ( ! (SttSkewCleanup( GAP, Oxx, Oyy, Rr, Charge, Start, // strarting point of trajectory. FI0, FiLimitAdmissible, nHitsSkew, ListHitsSkew, auxS, ApotemaMinSkew, //distance hexagon side from (0,0) ApotemaMaxSkew, //distance hexagon side from (0,0) 3., 1 ) ) ) { return false; } return true; }; //----------end of function PndSttMvdTracking::TrackCleanup //----------begin of function PndSttMvdTracking::SttParalCleanup bool PndSttMvdTracking::SttParalCleanup( Double_t GAP, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t FI0, Double_t FiLimitAdmissible, UShort_t &nHits, UShort_t *ListHits, Double_t info[][7], Double_t RStrawDetMin, Double_t ApotemaInnerParMax, Double_t ApotemaOuterParMin, Double_t RStrawDetMax ) { // this method does 3 things : // // 1) finds the entrance and exit points in the STT parallel volumes of the current track; // 2) eliminates from the track hit list possible spurious hits that are not encompassed // by the entrance and exit point; // 3) eliminates the tracks if the hit sequence is not continuous enough. // bool ConsiderLastHit; Short_t flagInnerSttR, flagOuterSttR, flagInnerSttL, flagOuterSttL; UShort_t enne, i, ihit, ipurged, nInnerHits, nOuterHits, nIntersections[2], ListInnerHits[nmaxSttHits], ListOuterHits[nmaxSttHits]; Double_t epsilonTheta, fi, r, aux[2], Xcross[2], Ycross[2], XcrossL[2], YcrossL[2], XcrossR[2], YcrossR[2], XintersectionList[5], // there is also the last boundary FiLimitAdmissible YintersectionList[5]; // take into account. //------------------------ // elimination of hits outside the physical FI range (FiLimitAdmissible) due to finite length of // Straws. epsilonTheta = STRAWRADIUS/Rr; // some extra slac for being conservative. for(i=0, ipurged=0; i< nHits; i++){ fi = atan2( info[ListHits[i]][1]-Oyy,info[ListHits[i]][0]-Oxx); if (fi<0.) fi+=2.*PI; if(Charge <0) { if( fi > FI0){ if( fi>FiLimitAdmissible+epsilonTheta) continue; } else { fi += 2.*PI; if( fi >FiLimitAdmissible+epsilonTheta ) continue; } // end of if( fi > FI0) } else { // continuation of if(Charge <0) if( fi > FI0){ fi -= 2.*PI; } // end of if( fi > FI0) if (fi < FiLimitAdmissible-epsilonTheta) continue; } // end of if(Charge <0) ListHits[ipurged]=ListHits[i]; ipurged++; } // end of for(i=0, ipurged=0; i< nHits; i++) nHits = ipurged; if(nHits==0){ nInnerHits=0; nOuterHits=0; goto jumpa ; // don't discard track yet, see if it should have parallel hits. } //------------------ // separation of inner Parallel Stt hits from outer Parallel Stt hits. SeparateInnerOuterParallel( // input nHits, ListHits, info, 2.*ApotemaInnerParMax/sqrt(3.), // output &nInnerHits, ListInnerHits, &nOuterHits, ListOuterHits ); //------------------------------------------ // find the entrance and exit of the track in the Inner Left Parallel Straw region. // This region is bounded by two Hexagons, and it has the target gap in the middle. // flag meaning : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; jumpa: ; flagInnerSttL=FindTrackEntranceExitbiHexagonLeft( GAP, Oxx, Oyy, Rr, Charge, Start, RStrawDetMin, ApotemaInnerParMax, XcrossL, YcrossL ); // find the entrance and exit of the track in the Inner Right Parallel Straw region. // This region is bounded by two Hexagons, and it has the target gap in the middle. flagInnerSttR=FindTrackEntranceExitbiHexagonRight( GAP, Oxx, Oyy, Rr, Charge, Start, RStrawDetMin, ApotemaInnerParMax, XcrossR, YcrossR ); //----------------- // working in the hypothesis that his is a track coming from Vertex at (0,0). if( flagInnerSttL != 0 && flagInnerSttR != 0 ){ nHits=0; // eliminate all the hits from hit list. return true; } if( flagInnerSttL == 0 && flagInnerSttR != 0 ){ XintersectionList[0]=XcrossL[0]; YintersectionList[0]=YcrossL[0]; XintersectionList[1]=XcrossL[1]; YintersectionList[1]=YcrossL[1]; if(fabs(FiLimitAdmissible-FI0) < 2.*PI){ // in this case the point // corresponding to FiLimitAdmissible can play a role in the // determination of the limiting points of the hits. aux[0]=Oxx+Rr*cos(FiLimitAdmissible); aux[1]=Oyy+Rr*sin(FiLimitAdmissible); XintersectionList[2]=aux[0]; YintersectionList[2]=aux[1]; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 3,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); if( fabs(aux[0]-Xcross[0])<1.e-5 && fabs(aux[1]-Ycross[0])<1.e-5 ){ nHits=0; // eliminate all the hits from hit list. return true; } if (nInnerHits == 0) return false; } else { // continuation of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) if (nInnerHits == 0) return false; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 2,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); } // end of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) } else if ( flagInnerSttL != 0 && flagInnerSttR == 0 ){ XintersectionList[0]=XcrossR[0]; YintersectionList[0]=YcrossR[0]; XintersectionList[1]=XcrossR[1]; YintersectionList[1]=YcrossR[1]; if(fabs(FiLimitAdmissible-FI0) < 2.*PI){ // in this case the point // corresponding to FiLimitAdmissible can play a role in the // determination of the limiting points of the hits. aux[0]=Oxx+Rr*cos(FiLimitAdmissible); aux[1]=Oyy+Rr*sin(FiLimitAdmissible); XintersectionList[2]=aux[0]; YintersectionList[2]=aux[1]; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 3,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); if( fabs(aux[0]-Xcross[0])<1.e-5 && fabs(aux[1]-Ycross[0])<1.e-5 ){ nHits=0; // eliminate all the hits from hit list. return true; } if (nInnerHits == 0) return false; } else { // continuation of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) if (nInnerHits == 0) return false; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 2,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); } // end of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) } else { // continuation of if (flagInnerSttL == 0 && flagInnerSttR != 0) // case in which flagInnerSttL==0 && flagInnerSttR==0. for(i=0;i<2;i++){ XintersectionList[i]=XcrossR[i]; XintersectionList[i+2]=XcrossL[i]; YintersectionList[i]=YcrossR[i]; YintersectionList[i+2]=YcrossL[i]; } if(fabs(FiLimitAdmissible-FI0) < 2.*PI){ // in this case the point // corresponding to FiLimitAdmissible can play a role in the // determination of the limiting points of the hits. aux[0]=Oxx+Rr*cos(FiLimitAdmissible); aux[1]=Oyy+Rr*sin(FiLimitAdmissible); XintersectionList[4]=aux[0]; YintersectionList[4]=aux[1]; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 5,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); if( fabs(aux[0]-Xcross[0])<1.e-5 && fabs(aux[1]-Ycross[0])<1.e-5 ){ nHits=0; // eliminate all the hits from hit list. return true; } if (nInnerHits == 0) return false; } else { // continuation of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) if (nInnerHits == 0) return false; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 4,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); } // end of if(fabs(FiLimitAdmissible-FI0) < 2.*PI). } // end of if( flagInnerSttL == 0 && flagInnerSttR != 0 ){ //------------- cleanup of the spurious tracks first using the inner parallel straws. // at this point the n. of inner hits cannot be 0 for a true track. if ( BadTrack_ParStt( Oxx, Oyy, Rr, Charge, Xcross, // Xcross[0]=point of entrance; Xcross[1]=point of exit. Ycross, nInnerHits, ListInnerHits, info, 1.5*DiameterStrawTube, // cut of proximity between hits. 1 // maximum allowed # consecutive hits with distance > cut. ) ){ return false; } //----------------------------------------------------- //------------ // find the entrance and exit of the track in the Outer Parallel Straw region, Left side. // This region is bounded by a Hexagon (inner), a Circle (outer) and it has // the target gap in the middle. // Returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections. flagOuterSttL=FindTrackEntranceExitHexagonCircleLeft( Oxx, Oyy, Rr, Charge, Start, ApotemaMinOuterParStraw, RStrawDetMax, GAP, XcrossL, YcrossL ); //------------ // find the entrance and exit of the track in the Outer Parallel Straw region, Right side. // This region is bounded by a Hexagon (inner), a Circle (outer) and it has // the target gap in the middle. flagOuterSttR=FindTrackEntranceExitHexagonCircleRight( Oxx, Oyy, Rr, Charge, Start, ApotemaMinOuterParStraw, RStrawDetMax, GAP, XcrossR, YcrossR ); //-------------------------------- if( flagOuterSttL != 0 && flagOuterSttR != 0 ){ if(nOuterHits==0){ return true; } else { nHits -= nOuterHits; // eliminate the Outer hits. return true; } } if( flagOuterSttL == 0 && flagOuterSttR != 0 ){ XintersectionList[0]=XcrossL[0]; YintersectionList[0]=YcrossL[0]; XintersectionList[1]=XcrossL[1]; YintersectionList[1]=YcrossL[1]; if(fabs(FiLimitAdmissible-FI0) < 2.*PI){ // in this case the point // corresponding to FiLimitAdmissible can play a role in the // determination of the limiting points of the hits. aux[0]=Oxx+Rr*cos(FiLimitAdmissible); aux[1]=Oyy+Rr*sin(FiLimitAdmissible); XintersectionList[2]=aux[0]; YintersectionList[2]=aux[1]; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 3,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); if( fabs(aux[0]-Xcross[0])<1.e-5 && fabs(aux[1]-Ycross[0])<1.e-5 ){ nHits -= nOuterHits; // eliminate the Outer hits. return true; } if (nOuterHits == 0) return false; } else { // continuation of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) if (nOuterHits == 0) return false; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 2,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); } // end of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) } else if ( flagOuterSttL != 0 && flagOuterSttR == 0 ){ XintersectionList[0]=XcrossR[0]; YintersectionList[0]=YcrossR[0]; XintersectionList[1]=XcrossR[1]; YintersectionList[1]=YcrossR[1]; if(fabs(FiLimitAdmissible-FI0) < 2.*PI){ // in this case the point // corresponding to FiLimitAdmissible can play a role in the // determination of the limiting points of the hits. aux[0]=Oxx+Rr*cos(FiLimitAdmissible); aux[1]=Oyy+Rr*sin(FiLimitAdmissible); XintersectionList[2]=aux[0]; YintersectionList[2]=aux[1]; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 3,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); if( fabs(aux[0]-Xcross[0])<1.e-5 && fabs(aux[1]-Ycross[0])<1.e-5 ){ nHits -= nOuterHits; // eliminate the Outer hits. return true; } if (nOuterHits == 0) return false; } else { // continuation of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) if (nOuterHits == 0) return false; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 2,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); } // end of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) } else { // this is when flagOuterSttL == 0 && flagOuterSttR == 0. for(i=0;i<2;i++){ XintersectionList[i]=XcrossR[i]; XintersectionList[i+2]=XcrossL[i]; YintersectionList[i]=YcrossR[i]; YintersectionList[i+2]=YcrossL[i]; } if(fabs(FiLimitAdmissible-FI0) < 2.*PI){ // in this case the point // corresponding to FiLimitAdmissible can play a role in the // determination of the limiting points of the hits. aux[0]=Oxx+Rr*cos(FiLimitAdmissible); aux[1]=Oyy+Rr*sin(FiLimitAdmissible); XintersectionList[4]=aux[0]; YintersectionList[4]=aux[1]; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 5,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); if( fabs(aux[0]-Xcross[0])<1.e-5 && fabs(aux[1]-Ycross[0])<1.e-5 ){ nHits=0; return true; } if (nOuterHits == 0) return false; } else { // continuation of if(fabs(FiLimitAdmissible-FI0) < 2.*PI) if (nOuterHits == 0) return false; ChooseEntranceExitbis( Oxx, Oyy, Charge, FI0, 4,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); } // end of if(fabs(FiLimitAdmissible-FI0) < 2.*PI). //-------------------- stampe if(istampa>=2&&IVOLTE<20){ cout<<"OUTER, caso R&L true, IVOLTE = "< cut. ) ) return false; // } // end of if(nOuterHits==0) //---------------------------------------------------------------------------- // if the code comes here it means that the track is acceptable. return true; }; //----------end of function PndSttMvdTracking::SttParalCleanup //----------begin of function PndSttMvdTracking::SttSkewCleanup bool PndSttMvdTracking::SttSkewCleanup( Double_t GAP, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t FI0, Double_t FiLimitAdmissible, UShort_t &nHits, UShort_t *ListHits, Double_t *S, Double_t RminStrawSkew, Double_t RmaxStrawSkew, Double_t cut, UShort_t maxnum ) { bool ConsiderLastHit; Short_t flagSttL, flagSttR; UShort_t i, ipurged, ibad, nIntersections[2]; Double_t cut2, epsilonTheta, fi, FiStart, r, Sprevious, aux[2], Distance[nmaxSttHits+1], Xcross[2], Ycross[2], XcrossL[2], YcrossL[2], XcrossR[2], YcrossR[2], XintersectionList[5], // second index =0 --> inner Hexagon, =1 --> outer. YintersectionList[5]; // first index : all the possible intersections // (up to 12 intersections). cut2=cut*cut; //------------------------ // elimination of hits outside the physical FI range (FiLimitAdmissible) due to finite length of // Straws. epsilonTheta = STRAWRADIUS/Rr; // some extra slac for being conservative. if(istampa>=2&&IVOLTE<20) cout<<"\n\nevt "<=2&&IVOLTE<20)cout<<"\thit // n. "< FI0){ if( fi>FiLimitAdmissible+epsilonTheta) continue; } else { fi += 2.*PI; if( fi > FiLimitAdmissible+epsilonTheta ) continue; } // end of if( fi > FI0) } else { // continuation of if(Charge <0) if( fi > FI0){ fi -= 2.*PI; } // end of if( fi > FI0) if (fi < FiLimitAdmissible-epsilonTheta) continue; } // end of if(Charge <0) if(istampa>=2&&IVOLTE<20)cout<<"\t\thit preso!"<cut2){ if(Distance[i]>16.*cut2){ return false; } ibad++; } } // end of do (i=1; icut2 ){ if( Distance[nHits]>16.*cut2){ return false; } ibad++; } } // end of if(IsInsideArc( if( ibad > maxnum){ return false; } return true; }; //----------end of function PndSttMvdTracking::SttSkewCleanup //----------begin of function PndSttMvdTracking::BadTrack_ParStt bool PndSttMvdTracking::BadTrack_ParStt( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Xcross[2], // Xcross[0]=point of entrance; // Xcross[1]=point of exit. Double_t Ycross[2], // bool ConsiderLastHit, UShort_t nHits, UShort_t* ListHits, Double_t info[][7], Double_t cut, UShort_t maxnum ) { UShort_t ibad, ihit, ninarc; Double_t cut2, Xprevious, Yprevious, S, Distance[nmaxSttHits+1]; cut2=cut*cut; ibad=0; Xprevious=Xcross[0]; Yprevious=Ycross[0]; for(ihit=0,ninarc=0; ihitcut2){ if(Distance[ihit]>16.*cut2){ return true; } ibad++; } } // end of for(ihit=0,NinDet=0 ;ihit=2)cout<<"from BadTrack_ParStt, Stt || hit n. (original notation) "<< ListHits[nHits-1]<<", Distance to boundary = "<cut2 ){ if( Distance[nHits]>16.*cut2){ return true; } ibad++; } } // end of if( IsInsideArc if(istampa>=2)cout<<"from BadTrack_ParStt, ibad "< length_segmentq || (x-P2x)*(x-P2x)+(y-P2y)*(y-P2y) > length_segmentq ) continue; status = true; XintersectionList[*Nintersections] = x; YintersectionList[*Nintersections] = y; (*Nintersections)++; } // end of for(ipossibility=-1; ... return status; } //----------end of function PndSttMvdTracking::IntersectionCircle_Segment //----------begin of function PndSttMvdTracking::IntersectionsWithGapSemicircle UShort_t PndSttMvdTracking::IntersectionsWithGapSemicircle( Double_t Oxx, Double_t Oyy, Double_t Rr, Double_t GAP, bool left, // if true --> Left semicircle; false --> Right semicircle. Double_t Rma, Double_t *XintersectionList, Double_t *YintersectionList ) { UShort_t nIntersectionsCircle; Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; nIntersectionsCircle=0; aaa = sqrt(Oxx*Oxx+Oyy*Oyy); // preliminary condition for having intersections between trajectory and Circle. if( !( aaa >= Rr + Rma || Rr >= aaa + Rma) && !( aaa + Rr <= Rma || aaa - Rr >= Rma ) ) { // now the calculation FI0 = atan2(-Oyy,-Oxx); cosFi = (aaa*aaa + Rr*Rr - Rma*Rma)/(2.*Rr*aaa); if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.; Fi = acos(cosFi); // (x1, y1) and (x2,y2) are the intersections between the trajectory // and the circle, in the laboratory reference frame. x1 = Oxx+Rr*cos(FI0 - Fi); y1 = Oyy+Rr*sin(FI0 - Fi); theta1 = atan2(y1,x1); // in this way theta1 is between -PI and PI. x2 = Oxx+Rr*cos(FI0 + Fi); y2 = Oyy+Rr*sin(FI0 + Fi); theta2 = atan2(y2,x2); // in this way theta2 is between -PI and PI. // Theta1, Theta2 = angle of the edges of the outer circle + Gap in the laboratory frame. // Theta1, Theta2 are angles between -PI/2 and +PI/2. if(!left){ // Right (looking into the beam) Semicircle. Theta2 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.); Theta1 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.); if( Theta1<= theta1 && theta1 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x1; YintersectionList[nIntersectionsCircle]=y1; nIntersectionsCircle++; } if( Theta1<= theta2 && theta2 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x2; YintersectionList[nIntersectionsCircle]=y2; nIntersectionsCircle++; } } else { // Left (looking into the beam) Semicircle. Theta2 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.); Theta1 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.); if( Theta1<= theta1 || theta1 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x1; YintersectionList[nIntersectionsCircle]=y1; nIntersectionsCircle++; } if( Theta1<= theta2 || theta2 <= Theta2 ){ XintersectionList[nIntersectionsCircle]=x2; YintersectionList[nIntersectionsCircle]=y2; nIntersectionsCircle++; } } } // end of if (!( a >= Rr + Rma || ..... //---------------------------- end of calculation intersection with outer circle. return nIntersectionsCircle; } //----------end of function PndSttMvdTracking::IntersectionsWithGapSemicircle //----------star of function PndSttMvdTracking::IsInternal bool PndSttMvdTracking::IsInternal( Double_t Px, // point Double_t Py, Double_t Xtraslation, Double_t Ytraslation, Double_t Theta ) { // for explanations see Gianluigi's logbook on page 278-280. if( (Xtraslation-Px)*sin(Theta) + (Py-Ytraslation)*cos(Theta) >= 0. ) return true; else return false; } //----------end of function PndSttMvdTracking::IsInternal //----------star of function PndSttMvdTracking::ChooseEntranceExitbis void PndSttMvdTracking::ChooseEntranceExitbis( Double_t Oxx, Double_t Oyy, Short_t Charge, Double_t FiStart, UShort_t nIntersections, Double_t *XintersectionList, // second index =1 -->inner polygon; Double_t *YintersectionList, // second index =2 -->outer polygon. Double_t Xcross[2], // output Double_t Ycross[2] // output ) { UShort_t i, j; // this method works under the hypothesis that there are at least 2 intersections. if (nIntersections<2) return; if(Charge > 0) { UShort_t auxIndex[100]; Double_t fi[100]; for( i=0;i FiStart) fi[i] -= 2.*PI; if( fi[i] > FiStart) fi[i] = FiStart; auxIndex[i]=i; } // end of for( i=0;i2.*ApotemaMaxInnerParStraw/sqrt(3.) ){ // outer Parallel hit. ListOuterHits[ *nOuterHits ] = ListHits[ihit]; (*nOuterHits)++; }else{ ListInnerHits[ *nInnerHits ] = ListHits[ihit]; (*nInnerHits)++; } } } //----------end of function PndSttMvdTracking::SeparateInnerOuterParallel //----------begin of function PndSttMvdTracking::FindTrackEntranceExitbiHexagonLeft Short_t PndSttMvdTracking::FindTrackEntranceExitbiHexagonLeft( Double_t vgap, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t ApotemaMax, Double_t Xcross[2], Double_t Ycross[2] ) { Short_t flag; UShort_t nIntersections; Double_t FiStart, XintersectionList[16], // all the possible intersections YintersectionList[16]; // (up to 16 intersections). // The following is the form of the Left (looking from downstream into the beam) biHexagon // geometrical shape considered in this method : // /* /| / | / | / / / / / / / / | | | | | | | | \ \ \ \ \ \ \ \ \ | \ | \| */ // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; flag=IntersectionsWithClosedbiHexagonLeft( vgap, Oxx, Oyy, Rr, ApotemaMin, // Apotema of the inner Hexagon, ApotemaMax,// Apotema of the outer Hexagon. &nIntersections, XintersectionList, // XintersectionList[..][0] --> inner polygon, // XintersectionList[..][1] --> outer polygon. YintersectionList ); // IMPORTANT : // this is true because here it is assumed that the track comes from (0,0,0) // otherwise the code must be changed! if (!(flag == 0)) return flag; if( nIntersections<2) return -1; //------- among all possible intersection find the entrance point (Xcross[0], Ycross[0]) // and the exit point (Xcross[1], Ycross[1]) of this track. //-------- the starting point of the track. FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); if(FiStart<0.) FiStart+= 2.*PI; if(FiStart<0.) FiStart =0.; // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. // so at this point, the intersections are at least 2. ChooseEntranceExitbis( Oxx, Oyy, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); //---------------------- return flag; } //----------end of function PndSttMvdTracking::FindTrackEntranceExitbiHexagonLeft //----------begin of function PndSttMvdTracking::FindTrackEntranceExitbiHexagonRight Short_t PndSttMvdTracking::FindTrackEntranceExitbiHexagonRight( Double_t vgap, Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t ApotemaMax, Double_t Xcross[2], Double_t Ycross[2] ) { Short_t flag; UShort_t nIntersections; Double_t FiStart, XintersectionList[16], // all the possible intersections YintersectionList[16]; // (up to 16 intersections). // The following is the form of the Left (looking from downstream into the beam) biHexagon // geometrical shape considered in this method : // /* |\ | \ | \ \ \ \ \ | | | | / / / / / / | / | / |/ */ // finding all possible intersections with inner parallel straw region. // The inner parallel straw region is delimited by two hexagons. // flag meaning : // -1 --> track outside outer perimeter; // 0 --> at least 1 intersection with polygon, therefore a possible entry and an exit; // 1 --> track contained completely between the two polygons; flag=IntersectionsWithClosedbiHexagonRight( vgap, Oxx, Oyy, Rr, ApotemaMin, // Apotema of the inner Hexagon, ApotemaMax,// Apotema of the outer Hexagon. &nIntersections, XintersectionList, YintersectionList ); // IMPORTANT : // this is true because here it is assumed that the track comes from (0,0,0) // otherwise the code must be changed! if (!(flag == 0)) return flag; if( nIntersections<2) return -1; //------- among all possible intersection find the entrance point (Xcross[0], Ycross[0]) // and the exit point (Xcross[1], Ycross[1]) of this track. //-------- the starting point of the track. FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); if(FiStart<0.) FiStart+= 2.*PI; if(FiStart<0.) FiStart =0.; // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. // so at this point, the intersections are at least 2. ChooseEntranceExitbis( Oxx, Oyy, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); //---------------------- return flag; } //----------end of function PndSttMvdTracking::FindTrackEntranceExitbiHexagonRight //----------begin of function PndSttMvdTracking::FindTrackEntranceExitHexagonCircleLeft Short_t PndSttMvdTracking::FindTrackEntranceExitHexagonCircleLeft( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t Rma, // outer radius of the Circle. Double_t GAP, Double_t Xcross[2], Double_t Ycross[2] ) { // This methods finds the intersections between a trajectory coming from (0,0) and // parameters Oxx, Oyy, Rr with the closed geometrical figure (in XY) formed by the STT // external Left semicircle and the Outer STT parallel Left straw semi-Hexagon + Gap for // the pellet target target. // It returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections. Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; //------------------ UShort_t nIntersectionsCircle, nIntersections; Double_t FiStart, XintersectionList[12], YintersectionList[12]; // all the possible intersections (up to 12 intersections). // finding all possible intersections with inner parallel straw region. Double_t Side_x[] = { -GAP/2., -GAP/2. , -ApotemaMin, -ApotemaMin, -GAP/2., -GAP/2. }, Side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-0.5*GAP)/sqrt(3.), ApotemaMin/sqrt(3.), -ApotemaMin/sqrt(3.), -(2.*ApotemaMin-0.5*GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)}, a[] = {1., -1./sqrt(3.), 1., 1./sqrt(3.), 1.}, b[] = {0., 1., 0., 1., 0.}, c[] = {GAP/2., -2.*ApotemaMin/sqrt(3.),ApotemaMin, 2.*ApotemaMin/sqrt(3.), GAP/2.}; nIntersections=IntersectionsWithOpenPolygon( Oxx, Oyy, Rr, 5, // n. Sides of open Polygon. a, // coefficient of formula : aX + bY + c = 0 defining the Polygon sides. b, c, Side_x, // X,Y coordinate of the Sides vertices (in sequence, following Side_y, // the Polygon along. XintersectionList, // XintersectionList YintersectionList // YintersectionList. ); //------------------------------------------------------------------------- // finding intersections of trajectory [assumed to originate from (0,0) ] // with outer semicircle, the Left part. nIntersectionsCircle=IntersectionsWithGapSemicircle( Oxx, // input from trajectory Oyy, // input from trajectory Rr, // input from trajectory GAP, // input, vertical gap in XY plane of STT detector. true, // true --> Left semicircle, false --> Right semicircle. Rma, // radius of external Stt containment. &XintersectionList[nIntersections], // output, X list of intersections (maximal 2). &YintersectionList[nIntersections] ); nIntersections += nIntersectionsCircle; //-------- the starting point of the track. if(nIntersections<2) return -1; FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); if(FiStart<0.) FiStart+= 2.*PI; if(FiStart<0.) FiStart =0.; // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. ChooseEntranceExitbis( Oxx, Oyy, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); return 0; } //----------end of function PndSttMvdTracking::FindTrackEntranceExitHexagonCircleLeft //----------begin of function PndSttMvdTracking::FindTrackEntranceExitHexagonCircleRight Short_t PndSttMvdTracking::FindTrackEntranceExitHexagonCircleRight( Double_t Oxx, Double_t Oyy, Double_t Rr, Short_t Charge, Double_t Start[3], Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0). Double_t Rma, // outer radius of the Circle. Double_t GAP, Double_t Xcross[2], Double_t Ycross[2] ) { // This methods finds the intersections between a trajectory coming from (0,0) and // parameters Oxx, Oyy, Rr with the closed geometrical figure (in XY) formed by the STT // external Left semicircle and the Outer STT parallel Left straw semi-Hexagon + Gap for // the pellet target target. // It returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections. Double_t cosFi, theta1, theta2, Theta1, Theta2, aaa, Fi, FI0, x1, x2, y1, y2; //------------------ UShort_t nIntersectionsCircle, nIntersections; Double_t FiStart, XintersectionList[12], YintersectionList[12]; // all the possible intersections (up to 10 intersections). // finding all possible intersections with inner parallel straw region. Double_t Side_x[] = { GAP/2., GAP/2. , ApotemaMin, ApotemaMin, GAP/2., GAP/2. }, Side_y[] = { sqrt(Rma*Rma-GAP*GAP/4.), (2.*ApotemaMin-GAP)/sqrt(3.), ApotemaMin/sqrt(3.), -ApotemaMin/sqrt(3.), -(2.*ApotemaMin-GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)}, a[] = {1., 1./sqrt(3.), 1., -1./sqrt(3.), 1.}, b[] = {0., 1., 0., 1., 0.}, c[] = {-GAP/2.,-2.*ApotemaMin/sqrt(3.),-ApotemaMin, 2.*ApotemaMin/sqrt(3.), -GAP/2.}; nIntersections=IntersectionsWithOpenPolygon( Oxx, Oyy, Rr, 5, // n. Sides of open Polygon. a, // coefficient of formula : aX + bY + c = 0 defining the Polygon sides. b, c, Side_x, // X,Y coordinate of the Sides vertices (in sequence, following Side_y, // the Polygon along. XintersectionList, // XintersectionList YintersectionList // YintersectionList. ); //------------------------------------------------------------------------- // finding intersections of trajectory [assumed to originate from (0,0) ] // with outer semicircle, the Left part. nIntersectionsCircle=IntersectionsWithGapSemicircle( Oxx, // input from trajectory Oyy, // input from trajectory Rr, // input from trajectory GAP, // input, vertical gap in XY plane of STT detector. false, // true --> Left semicircle, false --> Right semicircle. Rma, // radius of external Stt containment. &XintersectionList[nIntersections], // output, X list of intersections (maximal 2). &YintersectionList[nIntersections] ); nIntersections += nIntersectionsCircle; //-------- the starting point of the track. if(nIntersections<2) return -1; FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx); if(FiStart<0.) FiStart+= 2.*PI; if(FiStart<0.) FiStart =0.; // this method selects the entrance and exit points of the trajectory among all // geometrical intersections of the circular trajectory with the straw particular // volume. ChooseEntranceExitbis( Oxx, Oyy, Charge, FiStart, nIntersections, XintersectionList, YintersectionList, Xcross, // output Ycross // output ); return 0; } //----------end of function PndSttMvdTracking::FindTrackEntranceExitHexagonCircleRight //----------begin of function PndSttMvdTracking::IsInsideCircle bool PndSttMvdTracking::IsInsideArc( Double_t Oxx, Double_t Oyy, Short_t Charge, Double_t Xcross[2], Double_t Ycross[2], Double_t f // f should be between 0 and 2PI. ) { Double_t f1, f2; // Xcross[0],Ycross[0] is the point of entrance. f1 = atan2(Ycross[0]-Oyy, Xcross[0]-Oxx); if(f1<0.) f1+= 2.*PI; if(f1<0.) f1= 0.; f2 = atan2(Ycross[1]-Oyy, Xcross[1]-Oxx); if(f2<0.) f2+= 2.*PI; if(f2<0.) f2= 0.; if(Charge<0){ if(f1 > f2 ) f2 +=2.*PI; if(f1 > f2 ) f2 = f1; if( f>f1){ if(f 0. if(f1 < f2 ) f1 +=2.*PI; if(f1 < f2 ) f1 = f2; if( f>f2){ if(f0.) { XintersectionList[0] = gap; YintersectionList[0] = sqrt(delta) + Oyy; XintersectionList[1] = gap; YintersectionList[1] = -sqrt(delta) + Oyy; nintersections+=2; } // intersections of trajectory with plane X = -gap; delta = Rr*Rr - (-gap-Oxx)*(-gap-Oxx); if(delta>0.) { XintersectionList[nintersections] = -gap; YintersectionList[nintersections] = sqrt(delta) + Oyy; XintersectionList[nintersections+1] = -gap; YintersectionList[nintersections+1] = -sqrt(delta) + Oyy; nintersections+=2; } // intersections of trajectory with plane Z = gap; XintersectionList[nintersections] = Oxx + Rr*cos(fi0+kappa*gap); YintersectionList[nintersections] = Oyy + Rr*sin(fi0+kappa*gap); nintersections++; // intersections of trajectory with plane Z = -gap; XintersectionList[nintersections] = Oxx + Rr*cos(fi0-kappa*gap); YintersectionList[nintersections] = Oyy + Rr*sin(fi0-kappa*gap); nintersections++; ChooseEntranceExitbis( Oxx, Oyy, Charge, fi0, nintersections,// n. intersection in input. XintersectionList, YintersectionList, Xcross, // output Ycross // output ); // here I suppose that the Mvd system CONSERVATIVELY ends at Rmax=5 cm. if( fabs(Ycross[0]) > 5. ) return true; else return false; } //----------end of function PndSttMvdTracking::IsTargetInPipe ClassImp(PndSttMvdTracking)