#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; } // ------------------------------------------------------------------------- PndSttMvdTracking::PndSttMvdTracking(Int_t verbose) : FairTask("STT Stt-Mvd Tracking") { fPersistence = kTRUE; fVerbose = verbose; istampa = verbose; } // ------------------------------------------------------------------------- // ----- 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.); */ 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 SttTrack array // tracce che vengono dal fit di Pavia dell'elica // fSttTrackArray = (TClonesArray*) ioman->GetObject("STTTrack"); // if ( ! fSttTrackArray) // { // cout << "-E- PndSttMvdTracking::Init: No SttTrack 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 dopo digi 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); // Create and register output array /* fHelixHitArray = new TClonesArray("PndSttHelixHit"); ioman->Register("SttHelixHit","STT",fHelixHitArray, fPersistence); */ // CHECK added /* PndSttMapCreator *mapper = new PndSttMapCreator(fSttParameters); fSttTubeArray = mapper->FillTubeArray(); */ // ------------------------- get the Mvd hits 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("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[MAXTRACKSPEREVENT], SttSZfit[MAXTRACKSPEREVENT], Mvdhits[MAXTRACKSPEREVENT]; Short_t Candidato, Charge, 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, nTrackCandHit[MAXTRACKSPEREVENT], ListTrackCandHit[MAXTRACKSPEREVENT][nmaxSttHits+ nmaxMvdPixelHitsInTrack+ nmaxMvdStripHitsInTrack], ListTrackCandHitType[MAXTRACKSPEREVENT][nmaxSttHits+ // type = 0 --> Mvd Pixel nmaxMvdPixelHitsInTrack+ // type = 1 --> Mvd Strip nmaxMvdStripHitsInTrack], // type = 2 --> Stt Parallel // type = 3 --> Stt Straw // nHitMvdTrackCand, // nMvdPixelHit, // nMvdStripHit, nMvdMCPoint, // nMvdTrackCand, nSttHit, nSttParHit, nSttSkewHit, nSttHelixTrack, nSttMCPoint, nSttTrackCand, TemporarynSkewHitsinTrack, ListAllParHits[nmaxSttHits], ListAllSkewHits[nmaxSttHits], nSttHitsinTrack[MAXTRACKSPEREVENT], nSttParHitsinTrack[MAXTRACKSPEREVENT], nSttSkewHitsinTrack[MAXTRACKSPEREVENT], ListSttHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], ListSttHitsinTrackType[MAXTRACKSPEREVENT][nmaxSttHits], ListSttParHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], ListSttSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], nMvdPixelHitsAssociatedToSttTrack[MAXTRACKSPEREVENT], ListMvdPixelHitsAssociatedToSttTrack[MAXTRACKSPEREVENT][nmaxMvdPixelHits], nMvdStripHitsAssociatedToSttTrack[MAXTRACKSPEREVENT], ListMvdStripHitsAssociatedToSttTrack[MAXTRACKSPEREVENT][nmaxMvdStripHits] ; UShort_t ipinco, ncand, n; 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, highqualitycut, Pxini, px, Pyini, py, Pzini, qop, x, y, s[2], versor[2], z[2], zeta0, zeta1, zdrift[2], zerror[2], primoangolo[MAXTRACKSPEREVENT], ultimoangolo[MAXTRACKSPEREVENT]; 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++; //---------- fetching the MC truth tracks nMCTracks = fMCTrackArray->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"<=1&& IVOLTE<20) cout<<"da PndSttMvdTracking : n. punti MC = "<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 N. totale Hits in STT : "<GetEntriesFast(); nMvdStripHit = fMvdStripHitArray->GetEntriesFast(); if (nMvdPixelHit+nMvdStripHit ==0){ cout<<"da PndSttMvdTracking : N. of MvdPixelHit=nMvdStripHit=0, return!\n"<= nmaxMvdPixelHits) { cout<<"da PndSttMvdTracking : N. of MvdPixelHit = "<= nmaxMvdPixelHits ("<= nmaxMvdStripHits) { cout<<"da PndSttMvdTracking : N. of MvdStripHit = "<= nmaxMvdStripHits ("<GetEntriesFast(); // if (nMvdMCPoint ==0){ // cout<<"da PndSttMvdTracking : N. di MvdMCPoint =0, return!\n"<2&& IVOLTE<20) cout<<"N. MC Points delle Mvd = "<GetEntriesFast(); /* if (nMvdTrackCand ==0){ cout<<"da PndSttMvdTracking : N. of MvdTrackCand =0, return!\n"<= MAXMVDTRACKSPEREVENT) { cout<<"da PndSttMvdTracking : N. of MvdTrackCand = "<< nMvdTrackCand<<" and it is >= MAXMVDTRACKSPEREVENT (="<=1 && 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 (="<=3 && IVOLTE<20){ cout<<"N. totale di PndTrackCand del PR solo = " <At(i); ipunto= pSttHit->GetRefIndex();// modo giusto di estrarre il punto MC corrispondente. 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(); info[i][6]= pSttMCPoint[ipunto]->GetTrackID(); 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.; ListAllSkewHits[nSttSkewHit]=i; nSttSkewHit++; } //--------------- inizio stampaggi, stampe di controllo if (istampa >= 2) { cout <<"da PndSttMvdTracking, Stt iHit "<< i << " e n. punto MC ottenuto con RefIndex = " <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(); // cout<<"\tPixel n. "<=2 && IVOLTE<20) //------------------------------------------ stampaggi // --------------------------------------------- get MC Points of MVD Short_t MCPointtoMCTrackID[nMvdMCPoint]; Double_t XMvdMCPoint[nMvdMCPoint], YMvdMCPoint[nMvdMCPoint], ZMvdMCPoint[nMvdMCPoint]; for( i= 0; i< nMvdMCPoint; i++){ pMvdMCPoint = (PndSdsMCPoint*) fMvdMCPointArray->At(i); if(istampa>2&& IVOLTE<20) cout<<"Il punto n. "<GetTrackID()<Position(position); XMvdMCPoint[i]=position.X(); YMvdMCPoint[i]=position.Y(); ZMvdMCPoint[i]=position.Z(); MCPointtoMCTrackID[i]= pMvdMCPoint->GetTrackID(); } // ---------------------------------------------------------------------- //--------------- recupero le Helix-Lia PndSttTrack trovate dopo il PR // delle STT + fit di Lia /* nSttHelixTrack = fSttTrackArray->GetEntriesFast(); if(istampa>2 && IVOLTE<20){ cout<<"N. totale di PndSttTrack dopo PR+fit Lia = " <At(i); // ------ estraggo il PndTrackCand [che era uscito dal PR] // sul quale HelixFit e' stato applicato Candidato = pSttHelixTrack->GetTrackCandIndex(); // parameters of the helix: d0, phi0, Rad, tanlambda, z0 Dist =pSttHelixTrack->GetDist(); Phi = pSttHelixTrack->GetPhi(); Rad=pSttHelixTrack->GetRad(); TanL=pSttHelixTrack->GetTanL(); Z=pSttHelixTrack->GetZ(); Charge= pSttHelixTrack->GetCharge(); iflag = pSttHelixTrack->GetFlag(); } */ //--------------- fine del recupero le Helix-Lia PndSttTrack trovate dopo il PR delle STT + fit di Lia // ------------------------------------------ estraggo le altre info della trackcand del MVD bool inMvdTrackCandPixel[nMvdPixelHit], inMvdTrackCandStrip[nMvdStripHit]; for( i= 0; i< nMvdTrackCand ; i++){ 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 if(nHitMvdTrackCand[i]>= nmaxMvdPixelHitsInTrack+nmaxMvdStripHitsInTrack){ cout<<"from PndSttMvdTracking : # Mvd hits of this Mvd Track Cand is >= than the maximum allowed (" <2 && IVOLTE<20) cout<<"Evento n. "<GetSortedHit(j); ListHitMvdTrackCand[i][j] = pndtrackcandhit.GetHitId(); // questo e' il n. Hit nativo che posso usare // per estrarre tutte le info che voglio. ListHitTypeMvdTrackCand[i][j] = pndtrackcandhit.GetDetId(); // questo e' il n. Hit nativo che posso usare // per estrarre tutte le info che voglio. if( ListHitTypeMvdTrackCand[i][j]==FairRootManager::Instance()->GetBranchId("MVDHitsPixel")){ inMvdTrackCandPixel[ ListHitMvdTrackCand[i][j] ]= true; } else if( ListHitTypeMvdTrackCand[i][j]==FairRootManager::Instance()->GetBranchId("MVDHitsStrip")){ inMvdTrackCandStrip[ ListHitMvdTrackCand[i][j] ]= true; } if(istampa>2 && IVOLTE<20){ cout<<" hit n. "<=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 // for the calculation of the approximate Fi of the first hit in thos track pndtrackcandhit = pSttTrackCand->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.; } if(istampa > 2&& IVOLTE<20){ cout<<" da PndSttMvdtracking, Evento n. "<=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 nTotalCandidates = nSttTrackCand; // this may change in the future for(ncand=0; ncand< nTotalCandidates; ncand++){ nTrackCandHit[ncand] =nSttHitsinTrack[ncand]+ nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand]; UShort_t tempmvdindex[nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] ], tempmvdtype[nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] ], auxIndex[nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] ]; Double_t auxR[nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] ]; // adding the Mvd hits (Pixel and Strips) for(i=0; i< nMvdPixelHitsAssociatedToSttTrack[ncand]; i++){ auxR[i] = sqrt( XMvdPixel[ ListMvdPixelHitsAssociatedToSttTrack[ncand][i] ]* XMvdPixel[ ListMvdPixelHitsAssociatedToSttTrack[ncand][i] ]+ YMvdPixel[ ListMvdPixelHitsAssociatedToSttTrack[ncand][i] ]* YMvdPixel[ ListMvdPixelHitsAssociatedToSttTrack[ncand][i] ] ); tempmvdindex[i]=ListMvdPixelHitsAssociatedToSttTrack[ncand][i]; tempmvdtype[i]=0; auxIndex[i] = i; } for(i=0; i< nMvdStripHitsAssociatedToSttTrack[ncand]; i++){ auxR[i+nMvdPixelHitsAssociatedToSttTrack[ncand]] = sqrt( XMvdStrip[ ListMvdStripHitsAssociatedToSttTrack[ncand][i] ]* XMvdStrip[ ListMvdStripHitsAssociatedToSttTrack[ncand][i] ]+ YMvdStrip[ ListMvdStripHitsAssociatedToSttTrack[ncand][i] ]* YMvdStrip[ ListMvdStripHitsAssociatedToSttTrack[ncand][i] ] ); tempmvdindex[i+nMvdPixelHitsAssociatedToSttTrack[ncand]]= ListMvdStripHitsAssociatedToSttTrack[ncand][i]; tempmvdtype[i+nMvdPixelHitsAssociatedToSttTrack[ncand]]=1; auxIndex[i+nMvdPixelHitsAssociatedToSttTrack[ncand]]= i+nMvdPixelHitsAssociatedToSttTrack[ncand]; } // ordering the Mvd Hits if( nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] >0){ Merge_Sort( nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand], auxR, auxIndex); // constructing the ordered new Track Candidate now for(i=0; i< nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand]; i++){ ListTrackCandHit[ncand][i] = tempmvdindex[ auxIndex[i] ]; ListTrackCandHitType[ncand][i] = tempmvdtype[ auxIndex[i] ]; } } // end of if( nMvdPixelHitsAssociatedToSttTrack[ncand]+ for(i=0; i=2) cout<<"da PndSttMvdTracking, IVOLTE = "<< IVOLTE<<", ncand = " < 1) 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]]; } } } if( diff > 0.5 ) { trajectory_vertex[0]=trajectory_vertex[1]=0.; iexcl=-1; } RefitMvdStt( &nTrackCandHit[ncand], //this is both input and output after the new fit &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(istampa>=2) { cout<<"PndSttMvdTracking, prima di MatchMvdHitsToSttTracksagain, IVOLTE = "<< IVOLTE<<", ncand "< Fi_low_limit and // possibly > 2PI. ); } // 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( Mvdhits, delta, highqualitycut, nSttTrackCand, Ox, Oy, R, FI0, Fi_low_limit, // because here we deal with hits in Mvd region CHARGE, nMvdPixelHitsAssociatedToSttTrack, // input and output ListMvdPixelHitsAssociatedToSttTrack, // input and output nMvdStripHitsAssociatedToSttTrack, // input and output ListMvdStripHitsAssociatedToSttTrack // 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++) { //-----stampaggi if(istampa>=3) { cout<<"PndSttMvdTracking, dopo di MatchMvdHitsToSttTracksagain, IVOLTE = "<< IVOLTE<<", ncand "<=2) { cout<<"PndSttMvdTracking, prima di AssociateSkewHitsToXYTrack, IVOLTE = "<< IVOLTE<<", ncand "<=2) { 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 = nMvdPixelHitsAssociatedToSttTrack[ncand]+ j=nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[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[0])*(s[1]-FI0[0]); double modulo0 = ZED[0]*ZED[0] + (s[0]-FI0[0])*(s[0]-FI0[0]); double modulo1 = ZED[1]*ZED[1] + (s[1]-FI0[0])*(s[1]-FI0[0]); 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, S, // S can be modified by +-2*PI if necessary. &FI0[ncand], // this remains unchanged. CHARGE[ncand] ); //--------------------- 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; //------------------------------------------- // 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) EliminateSpuriousSZ( &nMvdPixelHitsAssociatedToSttTrack[ncand], // input and output &ListMvdPixelHitsAssociatedToSttTrack[ncand][0],// input and output &nMvdStripHitsAssociatedToSttTrack[ncand], // input and output &ListMvdStripHitsAssociatedToSttTrack[ncand][0],// input and output &nSttSkewHitsinTrack[ncand], // input and output &ListSttSkewHitsinTrack[ncand][0],// input and output Sbis, // input ZEDbis, // input DriftRadiusbis, // input ErrorDriftRadiusbis, // input &SchosenPixel[ncand][0], &SchosenStrip[ncand][0], &SchosenSkew[ncand][0], &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]) } // end of for(ncand=0; ncand< nTotalCandidates; ncand++) //--- redo association of parallel Stt straw hits to this track, after better refit. CollectParSttHitsagain( Mvdhits, delta, highqualitycut, info, nSttParHit, ListAllParHits, nSttTrackCand, Ox, Oy, R, KAPPA, FI0, Fi_low_limit, Fi_up_limit, nSttParHitsinTrack, // input and output ListSttParHitsinTrack // input and output ); //---------------------------------------------------------- // Questo va razionalizzato maggiormente dopo Groningen!! Riordino degli hits, per la seconda volta. // ordering all the hits belonging to the candidate track, by increasing R; // forming the new track with Mvd+Stt hits nTotalCandidates = nSttTrackCand; // this may change in the future for(ncand=0; ncand< nTotalCandidates; ncand++){ nTrackCandHit[ncand] =nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand]+ nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand]; UShort_t tempmvdindex[nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] ], tempmvdtype[nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] ], auxIndex[nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] ]; Double_t auxR[nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] ]; // adding the Mvd hits (Pixel and Strips) for(i=0; i< nMvdPixelHitsAssociatedToSttTrack[ncand]; i++){ auxR[i] = XMvdPixel[ ListMvdPixelHitsAssociatedToSttTrack[ncand][i] ]* XMvdPixel[ ListMvdPixelHitsAssociatedToSttTrack[ncand][i] ]+ YMvdPixel[ ListMvdPixelHitsAssociatedToSttTrack[ncand][i] ]* YMvdPixel[ ListMvdPixelHitsAssociatedToSttTrack[ncand][i] ]; tempmvdindex[i]=ListMvdPixelHitsAssociatedToSttTrack[ncand][i]; tempmvdtype[i]=0; auxIndex[i] = i; } for(i=0; i< nMvdStripHitsAssociatedToSttTrack[ncand]; i++){ auxR[i+nMvdPixelHitsAssociatedToSttTrack[ncand]] = XMvdStrip[ ListMvdStripHitsAssociatedToSttTrack[ncand][i] ]* XMvdStrip[ ListMvdStripHitsAssociatedToSttTrack[ncand][i] ]+ YMvdStrip[ ListMvdStripHitsAssociatedToSttTrack[ncand][i] ]* YMvdStrip[ ListMvdStripHitsAssociatedToSttTrack[ncand][i] ]; tempmvdindex[i+nMvdPixelHitsAssociatedToSttTrack[ncand]]= ListMvdStripHitsAssociatedToSttTrack[ncand][i]; tempmvdtype[i+nMvdPixelHitsAssociatedToSttTrack[ncand]]=1; auxIndex[i+nMvdPixelHitsAssociatedToSttTrack[ncand]]= i+nMvdPixelHitsAssociatedToSttTrack[ncand]; } // ordering the Mvd Hits if( nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand] >0){ Merge_Sort( nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand], auxR, auxIndex); // constructing the first part of the ordered new Track Candidate for(i=0; i< nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand]; i++){ ListTrackCandHit[ncand][i] = tempmvdindex[ auxIndex[i] ]; ListTrackCandHitType[ncand][i] = tempmvdtype[ auxIndex[i] ]; } } // end of if( nMvdPixelHitsAssociatedToSttTrack[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; j< nSttParHitsinTrack[ncand]+nSttSkewHitsinTrack[ncand];j++){ i = j+nMvdPixelHitsAssociatedToSttTrack[ncand]+ nMvdStripHitsAssociatedToSttTrack[ncand]; ListTrackCandHit[ncand][i] = tempmvdindex2[ auxIndex2[j] ]; ListTrackCandHitType[ncand][i] = tempmvdtype2[ auxIndex2[j] ]; } } // end of if( nSttParHitsinTrack[ncand]+ } // end of for(ncand=0; ncand< nTotalCandidates; ncand++) //-------------------- end of ordering //----------stampaggi if(istampa>2){ for(int it=0; it0 && nSttTrackCand > 0 ){ AssociateFoundTrackstoMCtris( info, Ox, Oy, R, nSttTrackCand, nSttParHitsinTrack, ListSttParHitsinTrack, nSttSkewHitsinTrack, ListSttSkewHitsinTrack, daTrackFoundaTrackMC ); } // endo of if( nMCTracks >0 && nSttTrackCand > 0 ) if(istampa>2){ for(i=0;i2){ for(ncand=0; ncand< nSttTrackCand; ncand++){ //----------------- ora la traccia MC corrispondente a questa traccia Stt cout<<"da 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 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 Cx = Oxx + Pyy*1000./(BFIELD*CVEL*carica); Cy = Oyy - Pxx*1000./(BFIELD*CVEL*carica); cout<<"da PndSttMvdTracking, evento (cominciando da 0) n. "<=3 ) { fprintf(HANDLE, "\n Evento %d NTotaleTracceMC %d ------\n",IVOLTE, nMCTracksaccettabili); } int ii; Double_t HoughFiii; for (ii=0; ii=3 ;ii++){ fprintf(HANDLE,"----------------------------------------------------------\n"); i=daTrackFoundaTrackMC[ii]; if( i <0 ) { fprintf(HANDLE, " No TracciaMC associated to found track n. %d in pattern recognition, with %d Hits ||, %d skew hits, %f Radius \n " ,ii,nSttParHitsinTrack[ii], nSttSkewHitsinTrack[ii], R[ii] ); continue; } if( ( !SttSZfit[ii] )&&(resultFitSZagain[ii] !=1 )) { fprintf(HANDLE, " TracciaMC %d; sua FoundTrack associata (n. %d) NONsoddisfaRequisitiMinimi perche' in Z-S e' fallito il fit \n",i,ii); continue; } if(fabs(KAPPA[ii])<1.e-20 ){ fprintf(HANDLE, " TracciaMC %d; sua FoundTrack associata (n. %d) NONsoddisfaRequisitiMinimi perche' KAPPA troppo piccolo; KAPPA = %g\n", i,ii,KAPPA[ii]); continue; } dista=sqrt( Ox[ii]*Ox[ii]+Oy[ii]*Oy[ii] ); if(fabs(dista)<1.e-20 ){ fprintf(HANDLE, " TracciaMC %d; sua FoundTrack associata (n. %d) NONsoddisfaRequisitiMinimi perche' centro Helix Cilinder trovato dista solo %g da (0,0)\n", i,ii,dista); continue; } pMCtr = (PndMCTrack*) fMCTrackArray->At(i); if ( ! pMCtr ){ fprintf(HANDLE, " MC track n. %d doesn't have pointer to MC Track TClones Array\n", i); continue; } fprintf(HANDLE, " TracciaMC %d ParHitsMC %d ParMecc %d ParMeccSpuri %d SkewHitsMC %d SkewMecc %d SkewMeccSpuri %d\n", i, nHitsInMCTrack[ii], nParalCommon[ii], nSpuriParinTrack[ii], nSkewHitsInMCTrack[ii], nSkewCommon[ii], nSpuriSkewinTrack[ii] ) ; fprintf(HANDLE, " e corrisponde a track found n. %d\n", ii ); fprintf(HANDLE, " AVENDO %d hits paralleli e %d hits skew non mecciati dalla corrisponde track found\n" , nMCParalAlone[ii],nMCSkewAlone[ii] ); 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 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 ;ii++) //--------------ghosts if( istampa>=3){ int NParghost=0, NParhitsghost=0,icc; for(icc=0; icc=1) //------------------ 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; for(ncand=0, ipinco = 0; ncand< nTotalCandidates; ncand++){ if(!status[ncand]) continue; Oxx = Ox[ncand]; Oyy = Oy[ncand]; dis=sqrt( Oxx*Oxx+Oyy*Oyy ); if( dis < 1.e-20) continue; Ptras = R[ncand]*0.003*BFIELD; Pxini = -CHARGE[ncand]*Ptras*Oyy/dis; Pyini = CHARGE[ncand]*Ptras*Oxx/dis; // starting point not necessarily at x=0., y=0. x= Ox[ncand] + R[ncand]*cos(FI0[ncand]); y= Oy[ncand] + R[ncand]*sin(FI0[ncand]); TVector3 posSeed(x,y,0.); // the starting point if(fabs(KAPPA[ncand])>1.e-20 ){ Pzini = -CHARGE[ncand]*0.003*BFIELD/KAPPA[ncand]; } else { Pzini = 999999.; } // 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] ); for(j=0; j< nTrackCandHit[ncand]; j++){ switch (ListTrackCandHitType[ncand][j]){ case 0: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId("MVDHitsPixel"),(Int_t)ListTrackCandHit[ncand][j],j); break; case 1: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId("MVDHitsStrip"),(Int_t)ListTrackCandHit[ncand][j],j); break; case 2: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId("STTHit"),(Int_t)ListTrackCandHit[ncand][j],j); break; case 3: pTrckCand->AddHit(FairRootManager::Instance()->GetBranchId("STTHit"),(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 ); 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>=1) cout<<" evento = "<=2) if(iplotta && IVOLTE<10){ for( i= 0; i< nSttTrackCand; i++){ int npunti=-1+nSttParHitsinTrack[i]+nSttSkewHitsinTrack[i]+ nMvdPixelHitsAssociatedToSttTrack[i]+nMvdStripHitsAssociatedToSttTrack[i]; if (ListTrackCandHitType[i][npunti] == 0){ // Mvd Pixel ultimoangolo[i] = atan2( YMvdPixel[ ListTrackCandHit[i][npunti] ]-Oy[i], XMvdPixel[ ListTrackCandHit[i][npunti] ]-Ox[i]); } else if (ListTrackCandHitType[i][npunti] == 1){ // Mvd Strip ultimoangolo[i] = atan2( YMvdStrip[ ListTrackCandHit[i][npunti] ]-Oy[i], XMvdStrip[ ListTrackCandHit[i][npunti] ]-Ox[i]); } else if( ListTrackCandHitType[i][npunti] == 2 ){ // it is a parallel straw hit PndSttInfoXYZParal ( info, ListTrackCandHit[i][npunti], Ox[i], Oy[i], R[i], KAPPA[i], FI0[i], CHARGE[i], Posiz1 ); ultimoangolo[i] = atan2( Posiz1[1]-Oy[i],Posiz1[0]-Ox[i]); } else if ( ListTrackCandHitType[i][npunti] == 3 ){ // it is a skew straw hit ultimoangolo[i] = SchosenSkew[i][ ListTrackCandHit[i][npunti] ]; } if( ultimoangolo[i]<0.) ultimoangolo[i]+= 2.*PI; double primo; primoangolo[i] = fmod(FI0[i],2.*PI); if(CHARGE[i]>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) { 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, daTrackFoundaTrackMC[i], nParalCommon, ParalCommonList, nSpuriParinTrack, ParSpuriList, nMCParalAlone, MCParalAloneList, nMvdPixelHitsAssociatedToSttTrack[i], ListMvdPixelHitsAssociatedToSttTrack, nMvdStripHitsAssociatedToSttTrack[i], ListMvdStripHitsAssociatedToSttTrack, 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]+ ... if( nSttSkewHitsinTrack[i]+nMvdPixelHitsAssociatedToSttTrack[i]+ nMvdStripHitsAssociatedToSttTrack[i]>0){ WriteMacroSkewAssociatedHitswithMC( KAPPA[i],FI0[i], Ox[i], Oy[i], R[i], info, WDX,WDY,WDZ, i,0, nSttSkewHitsinTrack[i], ListSttSkewHitsinTrack, nSkewCommon[i], SkewCommonList, daTrackFoundaTrackMC[i], nMCSkewAlone[i], MCSkewAloneList, nMvdPixelHitsAssociatedToSttTrack, ListMvdPixelHitsAssociatedToSttTrack, nMvdStripHitsAssociatedToSttTrack, ListMvdStripHitsAssociatedToSttTrack, 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 for( i= 0; i< nSttTrackCand; i++) i=0; WriteMacroParallelHitsGeneral( nSttHit, info, nSttTrackCand,Ox,Oy,R, FI0, ultimoangolo, primoangolo ); WriteMacroAllHitsRestanti( nSttHit, nSttParHit, nSttSkewHit, info, nSttTrackCand, 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 PndSttTrackFinderReal::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][nmaxSttHits], Double_t info[][7], UShort_t iTrack, 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 ListMvdPixelHitsAssociatedToSttTrack[MAXTRACKSPEREVENT][nmaxMvdPixelHits], UShort_t nMvdStripHitsAssociatedToSttTra, UShort_t ListMvdStripHitsAssociatedToSttTrack[MAXTRACKSPEREVENT][nmaxMvdStripHits], 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][nmaxSttHits], Double_t SchosenSkew[nmaxSttHits], 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]; // Ox = (D+R)*cos(Fi); // Oy = (D+R)*sin(Fi); //cout<<"da MacroTrackparalleletc. Ox, Oy "< xmax) xmax = info[i][0]+info[i][3]; if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; } for( ii=0; ii< nSkewHitsinTrack; ii++) { i = ListSkewHitsinTrack[iTrack][ii] ; aaa=Ox+R*cos(SchosenSkew[i]); bbb=Oy+R*sin(SchosenSkew[i]); if (aaa < xmin) xmin = aaa; if (aaa > xmax) xmax = aaa; if (bbb < ymin) ymin = bbb; if (bbb > ymax) ymax = bbb; // if (info[i][0]-info[i][3] < xmin) xmin = info[i][0]-info[i][3]; // if (info[i][0]+info[i][3] > xmax) xmax = info[i][0]+info[i][3]; // if (info[i][1]-info[i][3] < ymin) ymin = info[i][1]-info[i][3]; // if (info[i][1]+info[i][3] > ymax) ymax = info[i][1]+info[i][3]; // if(iTrack==9&&IVOLTE==2){ cout<<"skew hit n. "<Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); 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 = ListMvdStripHitsAssociatedToSttTrack[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 = ListMvdPixelHitsAssociatedToSttTrack[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 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 fprintf(MACRO,"}\n"); fclose(MACRO); return ; } //--end of function PndSttTrackFinderReal::WriteMacroSttParallelAssociatedHitsandMvdwithMC //----------start of function PndSttMvdTracking::WriteMacroParallelHitsGeneral void PndSttMvdTracking::WriteMacroParallelHitsGeneral( Int_t Nhits, Double_t info[][7], UShort_t nTracksFoundSoFar, Double_t *Ox, Double_t *Oy, Double_t *Radius, Double_t *FI0, Double_t *ultimoangolo, Double_t *primoangolo /* UShort_t nMvdPixelHit Double_t* XMvdPixel, Double_t* sigmaXMvdPixel, Double_t* YMvdPixel, Double_t* sigmaYMvdPixel, Double_t* ZMvdPixel, Double_t* sigmaZMvdPixel, UShort_t nMvdStripHit, Double_t* XMvdStrip, Double_t* sigmaXMvdStrip, Double_t* YMvdStrip, Double_t* sigmaYMvdStrip, Double_t* ZMvdStrip, Double_t* sigmaZMvdStrip */ ) { 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; fprintf(MACRO,"TCanvas* my= new TCanvas();\nmy->Range(%f,%f,%f,%f);\n",xmin,ymin,xmax,ymax); fprintf(MACRO,"TGaxis *Assex = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n",xmin,0.,xmax,0.,xmin,xmax); fprintf(MACRO,"Assex->Draw();\n"); fprintf(MACRO,"TGaxis *Assey = new TGaxis(%f,%f,%f,%f,%f,%f,510);\n", 0.,ymin,0.,ymax,ymin,ymax); fprintf(MACRO,"Assey->Draw();\n"); for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nE%d->SetFillStyle(0);\nE%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } 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( 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,"TBox* BS%d = new TBox(%f,%f,%f,%f);\nBS%d->SetFillColor(2);\nBS%d->Draw();\n", // ii,x1,y1,x2,y2,ii,ii); fprintf(MACRO,"TMarker* BS%d = new TMarker(%f,%f,%d);\nBS%d->SetMarkerColor(1);\nBS%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* BP%d = new TMarker(%f,%f,%d);\nBP%d->SetMarkerColor(1);\nBP%d->Draw();\n", ii,XMvdPixel[ii],YMvdPixel[ii],26,ii,ii); } //------------------------------- plotting all the tracks found for(i=0; iSetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", // i,aaa,bbb,rrr,rrr,i,i,i); "TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,%f,%f);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw(\"only\");\n", i,aaa,bbb,rrr,rrr,primoangolo[i],ultimoangolo[i],i,i,i); } // ----------- 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"); for( i=0; i< Nhits; i++) { if( info[i][5] == 1 ) { // parallel straws fprintf(MACRO,"TEllipse* E%d = new TEllipse(%f,%f,%f,%f,0.,360.);\nE%d->SetFillStyle(0);\nE%d->Draw();\n", i,info[i][0],info[i][1],info[i][3],info[i][3],i,i); } 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( 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,"TBox* BS%d = new TBox(%f,%f,%f,%f);\nBS%d->SetFillColor(2);\nBS%d->Draw();\n", // ii,x1,y1,x2,y2,ii,ii); fprintf(MACRO,"TMarker* BS%d = new TMarker(%f,%f,%d);\nBS%d->SetMarkerColor(1);\nBS%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* BP%d = new TMarker(%f,%f,%d);\nBP%d->SetMarkerColor(1);\nBP%d->Draw();\n", ii,XMvdPixel[ii],YMvdPixel[ii],26,ii,ii); } //------------------------------- plotting all the tracks found for(i=0; iSetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw();\n", // i,aaa,bbb,rrr,rrr,i,i,i); "TEllipse* ris%d=new TEllipse(%f,%f,%f,%f,%f,%f);\nris%d->SetFillStyle(0);\nris%d->SetLineColor(2);\nris%d->Draw(\"only\");\n", i,aaa,bbb,rrr,rrr,primoangolo[i],ultimoangolo[i],i,i,i); } // ----------- //----------------- 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 if( fabs(carica) > 0.1 ){ 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 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,primoangolo[i],ultimoangolo[i],i,i,i); "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); } } } // end of for(i=0; i= 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 if ( KAPPA >= 0.) { fmin = KAPPA*zmin + FI0; fmax = KAPPA*zmax + FI0; } else { fmax = KAPPA*zmin + FI0; fmin = KAPPA*zmax + FI0; } if( fmax>=0.) { Nmax = (int) (0.5*fmax/ PI); } else { Nmax = ( (int) (0.5*fmax/ PI) ) -1; } if( fmin>=0.) { Nmin = (int) (0.5*fmin/ PI); } else { Nmin = ((int) (0.5*fmin/ PI) )-1; } for(i=Nmin; i<= Nmax;i++){ offset = 2.*PI*i; z1 = (i*2.*PI-FI0)/KAPPA; z2 = ((i+1)*2.*PI-FI0)/KAPPA; fprintf(MACRO,"TLine* FOUND%d = new TLine(%f,%f,%f,%f);\nFOUND%d->SetLineColor(2);\nFOUND%d->Draw();\n", // i-Nmin,z1,0.,z2, 2.*PI,i-Nmin,i-Nmin); i-Nmin,z1,0.,z2, R*2.*PI,i-Nmin,i-Nmin); } // end of for(i=Nmin; i<= Nmax;++) //----------------- 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 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 ) 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], UShort_t nSttTrackCand, UShort_t nTrackCandHit[MAXTRACKSPEREVENT], UShort_t ListTrackCandHit[MAXTRACKSPEREVENT][nmaxSttHits+ nmaxMvdPixelHitsInTrack+ nmaxMvdStripHitsInTrack], UShort_t ListTrackCandHitType[MAXTRACKSPEREVENT][nmaxSttHits+ nmaxMvdPixelHitsInTrack+ nmaxMvdStripHitsInTrack] ) { // nSttHit = parallel+skew. bool exclusionStt[nSttHit], exclusionPixel[nMvdPixelHit], exclusionStrip[nMvdStripHit]; 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"); 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* BP%d = new TMarker(%f,%f,%d);\nBP%d->SetMarkerColor(1);\nBP%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* BS%d = new TMarker(%f,%f,%d);\nBS%d->SetMarkerColor(1);\nBS%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::AssociateFoundTrackstoMC 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[nTracksFoundSoFar][nmaxSttHits], inclusionExp[nTracksFoundSoFar]; UShort_t ntoMCtrack[nTracksFoundSoFar], toMCtracklist[nTracksFoundSoFar][nmaxSttHits], toMCtrackfrequency[nTracksFoundSoFar][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::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[nTracksFoundSoFar][nmaxSttHits], inclusionExp[nTracksFoundSoFar]; UShort_t ntoMCtrack[nTracksFoundSoFar], toMCtracklist[nTracksFoundSoFar][nmaxSttHits], toMCtrackfrequency[nTracksFoundSoFar][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[nTracksFoundSoFar][nmaxSttHits], inclusionExp[nTracksFoundSoFar]; UShort_t ntoMCtrack[nTracksFoundSoFar], toMCtrackfrequency[nTracksFoundSoFar][nmaxSttHits]; UShort_t i, j, jtemp,jexp , nmid; Short_t enne, itemp, massimo, toMCtracklist[nTracksFoundSoFar][nmaxSttHits]; Double_t dx, Cx, Cy, Rr, alfa, beta, gamma, minimo, tanlow[MAXTRACKSPEREVENT], tanmid[MAXTRACKSPEREVENT], tanup[MAXTRACKSPEREVENT], toMCtrackdistance[nTracksFoundSoFar][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::SttMatchedSpurious void PndSttMvdTracking::SttMatchedSpurious( 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][nmaxSttHits], // dal PR UShort_t nSkewHitsinTrack[MAXTRACKSPEREVENT], // n. hits skew, dal PR UShort_t ListSkewHitsinTrack[MAXTRACKSPEREVENT][nmaxSttHits], // 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, emme, enne[MAXTRACKSPEREVENT][nmaxSttHits]; 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){ printf("from FitHelixCylinder printout dopo glp_main -------------------------------\n"); printf(" number of structural variables %d\n",NStructVar); int ica; for(ica=0;ica=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<<"Risults : m1 = "< 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 = x+*qu . return -99; } } //----------end of function PndSttMvdTracking::FitSZspace //------------------ begin function PndSttMvdTracking::RefitMvdStt void PndSttMvdTracking::RefitMvdStt( UShort_t *nTrackCandHit, UShort_t *ListTrackCandHit, UShort_t *ListTrackCandHitType, 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[*nTrackCandHit], Yconformal[*nTrackCandHit], DriftRadiusconformal[*nTrackCandHit], ErrorDriftRadiusconformal[*nTrackCandHit]; *status= false; factor=3.; maxdis2=0.5; // trajectory_vertex[0]=trajectory_vertex[1]=0.; for(i=0, iparallel=0; i<*nTrackCandHit && iparallel< MAXIMUMHITSINFIT; i++){ if(i==iexcl) continue; if( ListTrackCandHitType[i] == 0 ){ // mvd pixels //----- translate the little circumference in XY representing // approximately the sensitive area into the conformal space // circumference dist2 = (XMvdPixel[ListTrackCandHit[i]]-tv[0])* (XMvdPixel[ListTrackCandHit[i]]-tv[0])+ (YMvdPixel[ListTrackCandHit[i]]-tv[1])* (YMvdPixel[ListTrackCandHit[i]]-tv[1]); if(dist22){ cout<<"PndSttMvdTracking::RefitMvdStt, n. Hits (Mvd+Stt || ) = "<2) cout<<"distanza**2 di Strip hit n. "< -1){ if( daTrackFoundaTrackMC[i] == FromPixeltoMCTrack[ ListMvdPixelHitsAssociatedToSttTrack[i][j] ] ){ index = i*nMvdPixelHit+nMvdPixelCommon[i]; *(MvdPixelCommonList+index) = ListMvdPixelHitsAssociatedToSttTrack[i][j]; nMvdPixelCommon[i]++; } else { index = i*nMvdPixelHit+nMvdPixelSpuriinTrack[i]; *(MvdPixelSpuriList+index) = ListMvdPixelHitsAssociatedToSttTrack[i][j]; nMvdPixelSpuriinTrack[i]++; } } } for(j=0;j -1){ if( daTrackFoundaTrackMC[i] == FromStriptoMCTrack[ ListMvdStripHitsAssociatedToSttTrack[i][j] ] ){ index = i*nMvdStripHit+nMvdStripCommon[i]; *(MvdStripCommonList+index) = ListMvdStripHitsAssociatedToSttTrack[i][j]; nMvdStripCommon[i]++; } else { index = i*nMvdStripHit+nMvdStripSpuriinTrack[i]; *(MvdStripSpuriList+index) = ListMvdStripHitsAssociatedToSttTrack[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. //-------------------- nMvdPixelHitsAssociatedToSttTrack[i]=0; nMvdStripHitsAssociatedToSttTrack[i]=0; ngoodmix=0; nn[0]=0; for( imvdcand=0; imvdcandGetBranchId("MVDHitsPixel")){ 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("MVDHitsPixel"); Dist += dist; if( dist anglemin) } else if (ListHitTypeMvdTrackCand[imvdcand][jmvdhit]== FairRootManager::Instance()->GetBranchId("MVDHitsStrip")){ 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("MVDHitsStrip"); 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("MVDHitsPixel")) { cout<<"\tPixel hit n. "<GetBranchId("MVDHitsStrip")){ 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("MVDHitsPixel"); 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("MVDHitsStrip"); 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("MVDHitsPixel")) { cout<<"\tDS Pixel hit n. "<GetBranchId("MVDHitsStrip")){ 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("MVDHitsPixel"); 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("MVDHitsStrip"); 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("MVDHitsPixel")) { cout<<"\tUS Pixel hit n. "<GetBranchId("MVDHitsStrip")){ 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 = "<0) } // end of for(i=0; i0){ // track must rotate clockwise looking into the beam. anglemax = FI0[itrack]; anglemin = Fifirst[itrack]; } else { anglemin = FI0[itrack]; anglemax = Fifirst[itrack]; } if(anglemax < anglemin) anglemax += 2.*PI; if(anglemax < anglemin) anglemax==anglemin; // this is just to be super-sure. // find if this track goes downstream or upstream ndownstream=0; for(j=0; j0.){ ndownstream++ ; } } for(j=0; j0.){ ndownstream++ ; } } if(ndownstream>ntot-ndownstream) downstream=true; else downstream=false; // loop over the Mvd Pixel and try to attach new Pixels to each candidate track naddpix=0; for(ipix=0; ipix=0.) || // ((!downstream) && ZMvdPixel[ipix]<0.) // ){ angle = atan2(YMvdPixel[ipix]-Oy[itrack], XMvdPixel[ipix]-Ox[itrack]); 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[itrack]-XMvdPixel[ipix])*(Ox[itrack]-XMvdPixel[ipix]) +(Oy[itrack]-YMvdPixel[ipix])*(Oy[itrack]-YMvdPixel[ipix]) ) -R[itrack]); if(dist anglemin) // } // end of if ( (downstream && out: ; } // end of for(ipix=0; ipix0){ for(j=0;j0) naddstr=0; for(istr=0; istr=0.) || // ( (!downstream) && ZMvdStrip[istr]<0.) // ){ angle = atan2(YMvdStrip[istr]-Oy[itrack], XMvdStrip[istr]-Ox[itrack]); 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[itrack]-XMvdStrip[istr])*(Ox[itrack]-XMvdStrip[istr]) +(Oy[itrack]-YMvdStrip[istr])*(Oy[itrack]-YMvdStrip[istr]) ) -R[itrack]); if(dist anglemin) // } // end of if ( (downstream && out2: ; } // end of for(istr=0; istr0){ for(j=0;j0) } // end of for(itrack=0; itrack 1.e-20){ deltaZ = 2.*PI/KAPPA[itrack]; Zpos = (angle - FI0[itrack])/KAPPA[itrack]; if( fabs(Zpos-info[ihit][2])>1.5*info[ihit][4] && fabs(Zpos+deltaZ-info[ihit][2])>1.5*info[ihit][4] && fabs(Zpos-deltaZ-info[ihit][2])>1.5*info[ihit][4] ) continue; } if( angle>Fi_up_limit[itrack]){ angle -= 2.*PI; if(angle < Fi_low_limit[itrack]) continue; } else if (angle < Fi_low_limit[itrack]){ angle += 2.*PI; if(angle > Fi_up_limit[itrack]) continue; } dist1 = fabs( sqrt( (Ox[itrack]-info[ihit][0])*(Ox[itrack]-info[ihit][0])+ (Oy[itrack]-info[ihit][1])*(Oy[itrack]-info[ihit][1]) ) - R[itrack] ); dist = fabs( dist1 - info[ihit][3] ); if(distAt(MCTrack); 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 *Cx = Oxx + Pyy*1000./(BFIELD*CVEL*carica); *Cy = Oyy - Pxx*1000./(BFIELD*CVEL*carica); } else { *Rr = -1.; } return; } //------------------------- end of function PndSttMvdTracking::getMCInfo //------------------------- begin of function PndSttMvdTracking::FindDistance Double_t PndSttMvdTracking::FindDistance( Double_t Ox, // center from wich distance is calculated Double_t Oy, // center from wich distance is calculated Double_t R, Double_t tanlow, Double_t tanmid, Double_t tanup, Double_t alfa, // intersection circumference parameter Double_t beta, // intersection circumference parameter Double_t gamma // intersection circumference parameter ) { UShort_t i, n; Double_t Delta, m[3], q, dist, dist1, dist2, distlow, distmid, distup, totaldist, x1, x2, y1, y2; int nevento=1; //if(IVOLTE==1) { cout<<"controprova, alfa = "< 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 ClassImp(PndSttMvdTracking)