From 1391e63398cd74db041d3a824ee184c4d6c8c485 Mon Sep 17 00:00:00 2001 From: martinez Date: Thu, 9 Oct 2003 13:43:57 +0000 Subject: [PATCH] New Stepmanager function. (Khalil Boudjemline Gines Martinez) New coordinate system in Slat segmentation (Eric Dumonteil) Config file option param operational --- MUON/AliMUON.cxx | 5 +- MUON/AliMUONData.cxx | 14 ++ MUON/AliMUONData.h | 7 +- MUON/AliMUONDigitizerv1.cxx | 18 ++ MUON/AliMUONDisplay.cxx | 19 +- MUON/AliMUONEventReconstructor.cxx | 4 +- MUON/AliMUONHit.cxx | 43 +++- MUON/AliMUONHit.h | 16 +- MUON/AliMUONSegmentationSlat.cxx | 14 +- MUON/AliMUONv1.cxx | 326 +++++++++++++---------------- MUON/AliMUONv1.h | 34 +-- MUON/Config_MUON_test.C | 20 +- MUON/MUONrawclusters.C | 10 + 13 files changed, 301 insertions(+), 229 deletions(-) diff --git a/MUON/AliMUON.cxx b/MUON/AliMUON.cxx index d34c95898ab..addbdb55ee8 100644 --- a/MUON/AliMUON.cxx +++ b/MUON/AliMUON.cxx @@ -147,8 +147,9 @@ AliMUON::AliMUON(const char *name, const char *title) } // Chamber stCH (0, 1) in } // Station st (0...) - fMaxStepGas=0.01; - fMaxStepAlu=0.1; + // Negatives values are ignored by geant3 CONS200 in the calculation of the tracking parameters + fMaxStepGas=0.1; + fMaxStepAlu=0.1; fMaxDestepGas=-1; fMaxDestepAlu=-1; diff --git a/MUON/AliMUONData.cxx b/MUON/AliMUONData.cxx index e3b70b483a4..fa6a74ebd3a 100644 --- a/MUON/AliMUONData.cxx +++ b/MUON/AliMUONData.cxx @@ -150,6 +150,20 @@ void AliMUONData::AddHit(Int_t fIshunt, Int_t track, Int_t iChamber, phi, length, destep); } //____________________________________________________________________________ +void AliMUONData::AddHit(Int_t fIshunt, Int_t track, Int_t iChamber, + Int_t idpart, Float_t X, Float_t Y, Float_t Z, + Float_t tof, Float_t momentum, Float_t theta, + Float_t phi, Float_t length, Float_t destep, + Float_t Xref,Float_t Yref,Float_t Zref) +{ + TClonesArray &lhits = *fHits; + new(lhits[fNhits++]) AliMUONHit(fIshunt, track, iChamber, + idpart, X, Y, Z, + tof, momentum, theta, + phi, length, destep, + Xref,Yref,Zref); +} +//____________________________________________________________________________ void AliMUONData::AddLocalTrigger(Int_t *localtr) { // add a MUON Local Trigger to the list diff --git a/MUON/AliMUONData.h b/MUON/AliMUONData.h index f54c0cdf00f..abbb0a94fc7 100644 --- a/MUON/AliMUONData.h +++ b/MUON/AliMUONData.h @@ -1,7 +1,7 @@ #ifndef ALIMUONDATA_H #define ALIMUONDATA_H // -/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ /* $Id$ */ @@ -42,6 +42,11 @@ class AliMUONData : public TNamed { Int_t idpart, Float_t X, Float_t Y, Float_t Z, Float_t tof, Float_t momentum, Float_t theta, Float_t phi, Float_t length, Float_t destep); + virtual void AddHit(Int_t fIshunt, Int_t track, Int_t iChamber, + Int_t idpart, Float_t X, Float_t Y, Float_t Z, + Float_t tof, Float_t momentum, Float_t theta, + Float_t phi, Float_t length, Float_t destep, + Float_t Xref,Float_t Yref,Float_t Zref); virtual void AddGlobalTrigger(Int_t *singlePlus, Int_t *singleMinus, Int_t *singleUndef, Int_t *pairUnlike, Int_t *pairLike); diff --git a/MUON/AliMUONDigitizerv1.cxx b/MUON/AliMUONDigitizerv1.cxx index ceafc77e6a7..882c55b06ac 100644 --- a/MUON/AliMUONDigitizerv1.cxx +++ b/MUON/AliMUONDigitizerv1.cxx @@ -77,7 +77,25 @@ Bool_t AliMUONDigitizerv1::ExistTransientDigit(AliMUONTransientDigit * mTD) //------------------------------------------------------------------------ Bool_t AliMUONDigitizerv1::Init() { + // Initialization + if (GetDebug()>2) Info("Init","AliMUONDigitizerv1::Init() starts"); + + //Loaders (We assume input0 to be the output too) + AliRunLoader * runloader; // Input loader + AliLoader * gime; + + // Getting runloader + runloader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0)); + if (runloader == 0x0) { + Error("Init","RunLoader is not in input file 0"); + return kFALSE; // RunDigitizer is not working. + } + // Getting MUONloader + gime = runloader->GetLoader("MUONLoader"); + gime->LoadHits("READ"); + gime->LoadDigits("RECREATE"); + return kTRUE; } diff --git a/MUON/AliMUONDisplay.cxx b/MUON/AliMUONDisplay.cxx index 5c8a018d7e4..1dc4250c045 100644 --- a/MUON/AliMUONDisplay.cxx +++ b/MUON/AliMUONDisplay.cxx @@ -622,6 +622,7 @@ void AliMUONDisplay::DrawTitle(Option_t *option) Float_t dy = ymax-ymin; + AliRunLoader * RunLoader = AliRunLoader::GetRunLoader("MUONFolder"); if (strlen(option) == 0) { TPaveText *title = new TPaveText(xmin +0.01*dx, ymax-0.09*dy, xmin +0.5*dx, ymax-0.01*dy); @@ -632,7 +633,7 @@ void AliMUONDisplay::DrawTitle(Option_t *option) title->Draw(); char ptitle[100]; sprintf(ptitle, "Alice event:%d Run:%d Chamber:%d Cathode:%d", - gAlice->GetHeader()->GetEvent(), + RunLoader->GetEventNumber(), gAlice->GetHeader()->GetRun(), fChamber, fCathode); @@ -1119,18 +1120,20 @@ void AliMUONDisplay::SetView(Float_t theta, Float_t phi, Float_t psi) //_____________________________________________________________________________ void AliMUONDisplay::ShowNextEvent(Int_t delta) { + + AliRunLoader * RunLoader = AliRunLoader::GetRunLoader("MUONFolder"); + // Display (current event_number + delta) // delta = 1 shown next event // delta = -1 show previous event if (delta) { - gAlice->Clear(); - Int_t currentEvent = gAlice->GetHeader()->GetEvent(); - Int_t newEvent = currentEvent + delta; - gAlice->GetEvent(newEvent); - fEvent=newEvent; - if (!gAlice->TreeD()) return; + //RunLoader->CleanDetectors(); + //RunLoader->CleanKinematics(); + Int_t currentEvent = RunLoader->GetEventNumber(); + Int_t newEvent = currentEvent + delta; + RunLoader->GetEvent(newEvent); + fEvent=newEvent; } - LoadDigits(fChamber, fCathode); fPad->cd(); Draw(); diff --git a/MUON/AliMUONEventReconstructor.cxx b/MUON/AliMUONEventReconstructor.cxx index af219e38f23..c79f995e21d 100644 --- a/MUON/AliMUONEventReconstructor.cxx +++ b/MUON/AliMUONEventReconstructor.cxx @@ -509,13 +509,15 @@ void AliMUONEventReconstructor::AddHitsForRecFromGEANT(TTree *TH) Int_t ihit, nhits=0; nhits = (Int_t) muondata->Hits()->GetEntriesFast(); AliMUONHit* mHit=0x0; + for(ihit=0; ihit(muondata->Hits()->At(ihit)); Int_t ipart = TMath::Abs ((Int_t) mHit->Particle()); //AZ if (NewHitForRecFromGEANT(mHit,track, hit, 1) && ipart == 13 //if (NewHitForRecFromGEANT(mHit,itrack-1, hit, 1) && ipart == 13 - && itrack <= 2) chamBits |= BIT(mHit->Chamber()-1); //AZ - set bit + && itrack <= 2 && !BIT(mHit->Chamber()-1) ) chamBits |= BIT(mHit->Chamber()-1); //AZ - set bit } + if (chamBits&3 && chamBits>>2&3 && chamBits>>4&3 && chamBits>>6&3 && chamBits>>8&3 && ((chamBits>>6&3)==3 || (chamBits>>8&3)==3)) fMuons += 1; //AZ diff --git a/MUON/AliMUONHit.cxx b/MUON/AliMUONHit.cxx index d59ea595806..00d6d90fde3 100644 --- a/MUON/AliMUONHit.cxx +++ b/MUON/AliMUONHit.cxx @@ -40,9 +40,14 @@ AliMUONHit::AliMUONHit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits): fPy = hits[12]; fPz = hits[13]; fAge = hits[14]; + fXref = 0.; + fYref = 0.; + fZref = 0.; } //___________________________________________ -AliMUONHit::AliMUONHit(Int_t shunt, Int_t track, Int_t iChamber, Int_t idpart, Float_t X, Float_t Y, Float_t Z, Float_t tof, Float_t momentum, Float_t theta, Float_t phi, Float_t length, Float_t destep): +AliMUONHit::AliMUONHit(Int_t shunt, Int_t track, Int_t iChamber, Int_t idpart, + Float_t X, Float_t Y, Float_t Z, Float_t tof, Float_t momentum, + Float_t theta, Float_t phi, Float_t length, Float_t destep): AliHit(shunt, track) { // Constructor @@ -59,12 +64,34 @@ AliMUONHit::AliMUONHit(Int_t shunt, Int_t track, Int_t iChamber, Int_t idpart, F fPHlast = 0; fPTot = momentum; fAge = tof; + fXref = 0.; + fYref = 0.; + fZref = 0.; } - - - - - - - +//----------------------------------------------------------------------------------------------- +AliMUONHit::AliMUONHit(Int_t shunt, Int_t track, Int_t iChamber, Int_t idpart, + Float_t X, Float_t Y, Float_t Z, Float_t tof, Float_t momentum, + Float_t theta, Float_t phi, Float_t length, Float_t destep, + Float_t Xref,Float_t Yref,Float_t Zref): + AliHit(shunt, track) +{ +// Constructor + fChamber = iChamber; + fParticle = idpart; + fX = X; + fY = Y; + fZ = Z; + fTheta = theta; + fPhi = phi; + fTlength = length; + fEloss = destep; + fPHfirst = 0; + fPHlast = 0; + fPTot = momentum; + fAge = tof; + fXref = Xref; + fYref = Yref; + fZref = Zref; +} +//----------------------------------------------------------------------------------------------- diff --git a/MUON/AliMUONHit.h b/MUON/AliMUONHit.h index edd1e4d70db..44c2cf6bdf0 100644 --- a/MUON/AliMUONHit.h +++ b/MUON/AliMUONHit.h @@ -14,6 +14,10 @@ class AliMUONHit : public AliHit { AliMUONHit() {} AliMUONHit(Int_t fIshunt, Int_t track, Int_t *vol, Float_t *hits); AliMUONHit(Int_t fIshunt, Int_t track, Int_t iChamber, Int_t idpart, Float_t X, Float_t Y, Float_t Z, Float_t tof, Float_t momentum, Float_t theta, Float_t phi, Float_t length, Float_t destep); + AliMUONHit(Int_t fIshunt, Int_t track, Int_t iChamber, Int_t idpart, + Float_t X, Float_t Y, Float_t Z, Float_t tof, Float_t momentum, + Float_t theta, Float_t phi, Float_t length, Float_t destep, + Float_t Xref, Float_t Yref, Float_t Zref); virtual ~AliMUONHit() {} Int_t Chamber() {return fChamber;} Float_t Particle() {return fParticle;} @@ -32,6 +36,11 @@ class AliMUONHit : public AliHit { Float_t Cy() {return fPy/fPTot;} Float_t Cz() {return fPz/fPTot;} + Float_t Xref() {return fXref;} + Float_t Yref() {return fYref;} + Float_t Zref() {return fZref;} + + private: Int_t fChamber; // Chamber number Float_t fParticle; // Geant3 particle type @@ -43,11 +52,16 @@ class AliMUONHit : public AliHit { Int_t fPHfirst; // first padhit Int_t fPHlast; // last padhit - Float_t fPTot; // hit local momentum P + Float_t fPTot; // Local momentum P of the track when entering in the chamber Float_t fPx; // Px Float_t fPy; // Py Float_t fPz; // Pz + Float_t fXref; // X position of hit in the center of the chamber (without angle effect) + Float_t fYref; // Y position of hit in the center of the chamber (without angle effect) + Float_t fZref; // Z position of hit in the center of the chamber (without angle effect) + + ClassDef(AliMUONHit,1) //Hit object for MUON }; #endif diff --git a/MUON/AliMUONSegmentationSlat.cxx b/MUON/AliMUONSegmentationSlat.cxx index 9add3aeb9ac..9e819681592 100644 --- a/MUON/AliMUONSegmentationSlat.cxx +++ b/MUON/AliMUONSegmentationSlat.cxx @@ -137,15 +137,17 @@ void AliMUONSegmentationSlat::GlobalToLocal( // positive side is shifted up // by half the overlap zlocal = z-fChamber->Z(); - zlocal = (x>0) ? zlocal-2.*fDz : zlocal+2.*fDz; + +// zlocal = (x>0) ? zlocal-2.*fDz : zlocal+2.*fDz; + zlocal = (x>0) ? zlocal+2.*fDz : zlocal-2.*fDz; //Change? + // Set the signs for the symmetry transformation and transform to first quadrant SetSymmetry(x); Float_t xabs=TMath::Abs(x); - Int_t ifirst = (zlocal < Float_t(0))? 0:1; -// + // Find slat number - for (i=ifirst; i= fYPosition[i]-eps) && (y <= fYPosition[i]+fSlatY+eps)) break; } @@ -203,7 +205,7 @@ LocalToGlobal(Int_t islat, Float_t xlocal, Float_t ylocal, Float_t &x, Float_ x = (xlocal+fXPosition[islat])*fSym; y=(ylocal+fYPosition[islat]); - z = (TMath::Even(islat)) ? -fDz : fDz ; + z = (TMath::Even(islat)) ? fDz : -fDz ; //Change for new referential z = (x>0) ? z+2.*fDz : z-2.*fDz ; z+=fChamber->Z(); @@ -277,7 +279,7 @@ GetPadC(Int_t ix, Int_t iy, Float_t &x, Float_t &y, Float_t &z) x=x*TMath::Sign(1,ix); // z-position - z = (TMath::Even(islat)) ? -fDz : fDz ; + z = (TMath::Even(islat)) ? fDz : -fDz ; //Change for new referential z = (x>0) ? z+2.*fDz : z-2.*fDz ; z += fChamber->Z(); } diff --git a/MUON/AliMUONv1.cxx b/MUON/AliMUONv1.cxx index 2cd88bec14f..ea11291fe59 100644 --- a/MUON/AliMUONv1.cxx +++ b/MUON/AliMUONv1.cxx @@ -27,6 +27,7 @@ #include #include #include +#include #include "AliCallf77.h" #include "AliConst.h" @@ -44,21 +45,22 @@ ClassImp(AliMUONv1) //___________________________________________ AliMUONv1::AliMUONv1() : AliMUON() + ,fTrackMomentum(), fTrackPosition() { // Constructor - fChambers = 0; - fStations = 0; - fStepManagerVersionOld = kFALSE; - fStepManagerVersionNew = kFALSE; - fStepManagerVersionTest = kFALSE; - - fStepMaxInActiveGas = 2.0; -} - - + fChambers = 0; + fStations = 0; + fStepManagerVersionOld = kFALSE; + fStepMaxInActiveGas = 0.6; + fStepSum = 0x0; + fDestepSum = 0x0; + fElossRatio = 0x0; + fAngleEffect10 = 0x0; + fAngleEffectNorma= 0x0; +} //___________________________________________ AliMUONv1::AliMUONv1(const char *name, const char *title) - : AliMUON(name,title) + : AliMUON(name,title), fTrackMomentum(), fTrackPosition() { // Constructor // By default include all stations @@ -69,10 +71,35 @@ AliMUONv1::AliMUONv1(const char *name, const char *title) factory.Build(this, title); fStepManagerVersionOld = kFALSE; - fStepManagerVersionNew = kFALSE; - fStepManagerVersionTest = kFALSE; - fStepMaxInActiveGas = 2.0; + fStepMaxInActiveGas = 0.6; + + fStepSum = new Float_t [AliMUONConstants::NCh()]; + fDestepSum = new Float_t [AliMUONConstants::NCh()]; + for (Int_t i=0; iSetParameter(0,1.02138); + fElossRatio->SetParameter(1,-9.54149e-02); + fElossRatio->SetParameter(2,+7.83433e-02); + fElossRatio->SetParameter(3,-9.98208e-03); + fElossRatio->SetParameter(4,+3.83279e-04); + + // Angle effect in tracking chambers at theta =10 degres as a function of ElossRatio (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis) (in micrometers) + fAngleEffect10 = new TF1("AngleEffect10","[0]+[1]*x+[2]*x*x",0.5,3.0); + fAngleEffect10->SetParameter(0, 1.90691e+02); + fAngleEffect10->SetParameter(1,-6.62258e+01); + fAngleEffect10->SetParameter(2,+1.28247e+01); + // Angle effect: Normalisation form theta=10 degres to theta between 0 and 10 (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis) + // Angle with respect to the wires assuming that chambers are perpendicular to the z axis. + fAngleEffectNorma = new TF1("AngleEffectNorma","[0]+[1]*x+[2]*x*x+[3]*x*x*x",0.0,10.0); + fAngleEffectNorma->SetParameter(0,4.148); + fAngleEffectNorma->SetParameter(1,-6.809e-01); + fAngleEffectNorma->SetParameter(2,5.151e-02); + fAngleEffectNorma->SetParameter(3,-1.490e-03); } //___________________________________________ @@ -1596,190 +1623,139 @@ void AliMUONv1::Init() //cp } -//___________________________________________ + +//_______________________________________________________________________________ +Int_t AliMUONv1::GetChamberId(Int_t volId) const +{ +// Check if the volume with specified volId is a sensitive volume (gas) +// of some chamber and returns the chamber number; +// if not sensitive volume - return 0. +// --- + + for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++) + if (volId==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) return i; + + return 0; +} +//_______________________________________________________________________________ void AliMUONv1::StepManager() { - if (fStepManagerVersionOld) { + if (fStepManagerVersionOld) { StepManagerOld(); return; } - if (fStepManagerVersionNew) { - StepManagerNew(); - return; - } - - if (fStepManagerVersionTest) { - StepManagerTest(); - return; - } - - - // Volume id - Int_t copy, id; - Int_t idvol; - Int_t iChamber=0; - // Particule id, pos and mom vectors, - // theta, phi angles with respect the normal of the chamber, - // spatial step, delta_energy and time of flight - Int_t ipart; - TLorentzVector pos, mom; - Float_t theta, phi, tof; - Float_t destep, step; - const Float_t kBig = 1.e10; // Only charged tracks if( !(gMC->TrackCharge()) ) return; - + // Only charged tracks + // Only gas gap inside chamber // Tag chambers and record hits when track enters - idvol=-1; + Int_t idvol=-1; + Int_t iChamber=0; + Int_t id=0; + Int_t copy; + const Float_t kBig = 1.e10; + id=gMC->CurrentVolID(copy); + // printf("id == %d \n",id); for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++) { if(id==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) { iChamber = i; idvol = i-1; } } - if (idvol == -1) return; - - // printf(">>>> This Chamber %d\n",iChamber); - - // record hits when track enters ... - if( gMC->IsTrackEntering()) gMC->SetMaxStep(fStepMaxInActiveGas); - - if (gMC->TrackStep() > 0.) { - // Get current particle id (ipart), track position (pos) and momentum (mom) - gMC->TrackPosition(pos); - gMC->TrackMomentum(mom); - ipart = gMC->TrackPid(); - theta = mom.Theta()*kRaddeg; // theta of track - phi = mom.Phi() *kRaddeg; // phi of the track - tof = gMC->TrackTime(); // Time of flight - // - // momentum loss and steplength in last step - destep = gMC->Edep(); - step = gMC->TrackStep(); - - //new hit - GetMUONData()->AddHit(fIshunt, gAlice->GetCurrentTrackNumber(), iChamber, ipart, - pos.X(), pos.Y(), pos.Z(), tof, mom.P(), - theta, phi, step, destep); - } - // Track left chamber ... - if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){ - gMC->SetMaxStep(kBig); + if (idvol == -1) { + return; } -} - - -Int_t AliMUONv1::GetChamberId(Int_t volId) const -{ -// Check if the volume with specified volId is a sensitive volume (gas) -// of some chamber and returns the chamber number; -// if not sensitive volume - return 0. -// --- - for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++) - if (volId==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) return i; - - return 0; -} -//__ - - - -void AliMUONv1::StepManagerTest() -{ - return; -} -//________________________________________ -void AliMUONv1::StepManagerNew() -{ - -// // Volume id -// Int_t copy, id; -// Int_t idvol; -// Int_t iChamber=0; -// // Particule id, pos and mom vectors, -// // theta, phi angles with respect the normal of the chamber, -// // spatial step, delta_energy and time of flight -// Int_t ipart; -// TLorentzVector pos, mom; -// Float_t theta, phi, tof; -// Float_t destep, step; -// const Float_t kBig = 1.e10; - -// // Only charged tracks -// if( !(gMC->TrackCharge()) ) return; - -// // Only gas gap inside chamber -// // Tag chambers and record hits when track enters -// idvol=-1; -// id=gMC->CurrentVolID(copy); -// for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++) { -// if(id==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) { -// iChamber = i; -// idvol = i-1; -// } -// } -// static Float_t Sstep[20]; // Sum of steps per chamber -// // static Float_t Sdestep[20]; // Sum of eloss per chamber -// Float_t GAP; -// Float_t TEST; - -// if (idvol == -1) return; - -// // printf(">>>> This Chamber %d\n",iChamber); - -// // record hits when track enters ... -// //if( gMC->IsTrackEntering()) gMC->SetMaxStep(fStepMaxInActiveGas); - -// if (gMC->TrackStep() > 0.) { -// // Get current particle id (ipart), track position (pos) and momentum (mom) -// gMC->TrackPosition(pos); -// gMC->TrackMomentum(mom); -// ipart = gMC->TrackPid(); // Particle -// theta = mom.Theta()*kRaddeg; // theta of track -// phi = mom.Phi() *kRaddeg; // phi of the track -// tof = gMC->TrackTime(); // Time of flight -// // -// // momentum loss and steplength in last step -// destep = gMC->Edep(); -// step = gMC->TrackStep(); - -// Sstep[iChamber]+=step; -// // Sdestep[iChamber]+=destep; + if( gMC->IsTrackEntering() ) { + Float_t theta = fTrackMomentum.Theta(); + Float_t phi = fTrackMomentum.Phi(); + Float_t theta_wires = TMath::Abs( TMath::ASin( TMath::Sin(theta) * TMath::Sin(phi) ) ); + if ( (theta_wires*kRaddeg>=10) ) gMC->SetMaxStep(fStepMaxInActiveGas); + } -// } +// if (GetDebug()) { +// Float_t z = ( (AliMUONChamber*)(*fChambers)[idvol])->Z() ; +// Info("StepManager Step","Active volume found %d chamber %d Z chamber is %f ",idvol,iChamber, z); +// } + // Particule id and mass, + Int_t ipart = gMC->TrackPid(); + Float_t mass = gMC->TrackMass(); + + fDestepSum[idvol]+=gMC->Edep(); + // Get current particle id (ipart), track position (pos) and momentum (mom) + if ( fStepSum[idvol]==0.0 ) gMC->TrackMomentum(fTrackMomentum); + fStepSum[idvol]+=gMC->TrackStep(); -// step = Sstep[iChamber]; // Total step >= gap -// // destep = Sdestep[iChamber]; // Total eloss - - -// // Track left chamber ... -// if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){ -// gMC->SetMaxStep(kBig); - -// Sstep[iChamber]=0; // Reset for the next event -// //Sdestep[iChamber]=0; // Reset for the next event +// if (GetDebug()) { +// Info("StepManager Step","iChamber %d, Particle %d, theta %f phi %f mass %f StepSum %f eloss %g", +// iChamber,ipart, fTrackMomentum.Theta()*kRaddeg, fTrackMomentum.Phi()*kRaddeg, mass, fStepSum[idvol], gMC->Edep()); +// Info("StepManager Step","Track Momentum %f %f %f", fTrackMomentum.X(), fTrackMomentum.Y(), fTrackMomentum.Z()) ; +// gMC->TrackPosition(fTrackPosition); +// Info("StepManager Step","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ; +// } -// if (iChamber>=1 && iChamber<=2) GAP=0.4; -// if (iChamber>=11 && iChamber<=14) GAP=0.2; -// if (iChamber>=3 && iChamber<=10) GAP=0.5; - -// TF1 *ELOSS1 = new TF1("Gauss1","exp(-((x-4.13727e+01)**2)/(2*1.42223e+01**2))",0,75); -// TF1 *ELOSS2 = new TF1("Gauss2","exp(-((x+6.83795e+02)**2)/(2*4.48415e+02**2))",75,350); -// TEST=gRandom->Rndm(); -// if (TEST <=0.89) destep=ELOSS1->GetRandom(); -// else destep=ELOSS2->GetRandom(); -// destep*=pow(10,-6)*0.0274; -// destep*=GAP/0.5; - -// // One hit per chamber -// GetMUONData()->AddHit(fIshunt, gAlice->GetCurrentTrackNumber(), iChamber, ipart, -// pos.X()-(step/2*sin(theta*kDegrad)*cos(phi*kDegrad)), pos.Y()-(step/2*sin(theta*kDegrad)*sin(phi*kDegrad)), pos.Z()-GAP/2, tof, mom.P(),theta, phi, step, destep); + // Track left chamber or StepSum larger than fStepMaxInActiveGas + if ( gMC->IsTrackExiting() || + gMC->IsTrackStop() || + gMC->IsTrackDisappeared()|| + (fStepSum[idvol]>fStepMaxInActiveGas) ) { + + if ( gMC->IsTrackExiting() || + gMC->IsTrackStop() || + gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig); -// } + gMC->TrackPosition(fTrackPosition); + Float_t theta = fTrackMomentum.Theta(); + Float_t phi = fTrackMomentum.Phi(); + + TLorentzVector BackToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi), + fStepSum[idvol]/2.*sin(theta)*sin(phi), + fStepSum[idvol]/2.*cos(theta),0.0 ); + // if (GetDebug()) + // Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ; + // if (GetDebug()) + // Info("StepManager Exit ","Track BackToWire %f %f %f",BackToWire.X(),BackToWire.Y(),BackToWire.Z()) ; + fTrackPosition-=BackToWire; + + //-------------- Angle effect + // Ratio between energy loss of particle and Mip as a function of BetaGamma of particle (Energy/Mass) + + Float_t Beta_x_Gamma = fTrackMomentum.P()/mass;// pc/mc2 + Float_t SigmaEffect_10degrees; + Float_t SigmaEffect_thetadegrees; + Float_t ELossParticle_ELossMip; + Float_t YAngleEffect=0.; + Float_t theta_wires = TMath::Abs( TMath::ASin( TMath::Sin(theta) * TMath::Sin(phi) ) ); + + if ( (Beta_x_Gamma >3.2) && (theta_wires*kRaddeg<=10) ) { + Beta_x_Gamma=TMath::Log(Beta_x_Gamma); + ELossParticle_ELossMip = fElossRatio->Eval(Beta_x_Gamma); + // 10 degrees is a reference for a model (arbitrary) + SigmaEffect_10degrees=fAngleEffect10->Eval(ELossParticle_ELossMip);// in micrometers + // Angle with respect to the wires assuming that chambers are perpendicular to the z axis. + SigmaEffect_thetadegrees = SigmaEffect_10degrees/fAngleEffectNorma->Eval(theta_wires*kRaddeg); // For 5mm gap + if ( (iChamber==1) || (iChamber==2) ) + SigmaEffect_thetadegrees/=(1.09833e+00+1.70000e-02*theta_wires*kRaddeg); // The gap is different (4mm) + YAngleEffect=1.e-04*gRandom->Gaus(0,SigmaEffect_thetadegrees); // Error due to the angle effect in cm + } + + + // One hit per chamber + GetMUONData()->AddHit(fIshunt, gAlice->GetCurrentTrackNumber(), iChamber, ipart, + fTrackPosition.X(), fTrackPosition.Y()+YAngleEffect, fTrackPosition.Z(), 0.0, + fTrackMomentum.P(),theta, phi, fStepSum[idvol], fDestepSum[idvol], + fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()); +// if (GetDebug()){ +// Info("StepManager Exit","Particle exiting from chamber %d",iChamber); +// Info("StepManager Exit","StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]); +// Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ; +// } + fStepSum[idvol] =0; // Reset for the next event + fDestepSum[idvol]=0; // Reset for the next event + } } //___________________________________________ diff --git a/MUON/AliMUONv1.h b/MUON/AliMUONv1.h index a59b1a6842f..02ca246ae7a 100644 --- a/MUON/AliMUONv1.h +++ b/MUON/AliMUONv1.h @@ -10,8 +10,11 @@ // Manager and hits classes for set:MUON version 0 // ///////////////////////////////////////////////////////// +#include "TLorentzVector.h" + #include "AliMUON.h" -#include + +class TF1; class AliMUONv1 : public AliMUON { public: @@ -24,28 +27,29 @@ public: virtual Int_t IsVersion() const {return 1;} virtual void StepManager(); void StepManagerOld(); - void StepManagerNew(); - void StepManagerTest(); - - - void SetStepManagerVersionOld(Bool_t Opt) { fStepManagerVersionOld = Opt; } - void SetStepManagerVersionNew(Bool_t Opt) - { fStepManagerVersionNew = Opt; } - void SetStepManagerVersionTest(Bool_t Opt) - { fStepManagerVersionTest = Opt; } void SetStepMaxInActiveGas(Float_t StepMax) {fStepMaxInActiveGas = StepMax; } protected: - Int_t* fStations; //! allow to externally set which station to create + Int_t* fStations; //! allow to externally set which station to create Bool_t fStepManagerVersionOld; // Version of StepManager, Default is false - Bool_t fStepManagerVersionNew; // Version of StepManager, Default is false - Bool_t fStepManagerVersionTest; // Version of StepManager, Default is false - Float_t fStepMaxInActiveGas; // Step mas in active gas default 0.6cm + Float_t fStepMaxInActiveGas; // Step max in active gas default 0.6cm virtual Int_t GetChamberId(Int_t volId) const; - + // StepManager + Float_t * fStepSum; //! + Float_t * fDestepSum; //! + // Momentum of the particle entering in the active gas of chamber + TLorentzVector fTrackMomentum; + // Position of the particle exiting the active gas of chamber + TLorentzVector fTrackPosition; + // Ratio of particle mean eloss with respect MIP's + TF1 * fElossRatio; + // Angle effect in tracking chambers at theta =10 degres as a function of ElossRatio (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis) (in micrometers) + TF1 * fAngleEffect10; + // Angle effect: Normalisation form theta=10 degres to theta between 0 and 10 (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis) + TF1 * fAngleEffectNorma; private: ClassDef(AliMUONv1,1) // MUON Detector class Version 1 diff --git a/MUON/Config_MUON_test.C b/MUON/Config_MUON_test.C index 1a8642b5c8e..81a5d7852fb 100644 --- a/MUON/Config_MUON_test.C +++ b/MUON/Config_MUON_test.C @@ -6,7 +6,7 @@ void Config(char directory[100]="", char option[6]="box") { // // Config file for MUON test - // Gines MARITNEZ, Subatech, mai 2003, august 2003 + // Gines MARTINEZ, Subatech, mai 2003, august 2003 // //===================================================================== @@ -87,8 +87,8 @@ void Config(char directory[100]="", char option[6]="box") if (!strcmp(option,"box")) { AliGenBox * gener = new AliGenBox(1); gener->SetMomentumRange(20.,20.1); - gener->SetPhiRange(45., 45.01); - gener->SetThetaRange(4.000,4.001); + gener->SetPhiRange(-180, 180.01); + gener->SetThetaRange(171.000,178.001); gener->SetPart(13); // Muons gener->SetOrigin(0.,0., 0.); //vertex position gener->SetSigma(0.0, 0.0, 0.0); //Sigma in (X,Y,Z) (cm) on IP position @@ -112,9 +112,8 @@ void Config(char directory[100]="", char option[6]="box") gener->SetMomentumRange(0,999); gener->SetPtRange(0,100.); gener->SetPhiRange(-180, 180); - gener->SetYRange(2.5,4); gener->SetCutOnChild(1); - gener->SetChildThetaRange(2.0,9); + gener->SetChildThetaRange(171.0,178.0); gener->SetOrigin(0,0,0); //vertex position gener->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position gener->SetForceDecay(kDiMuon); gener->SetTrackingFlag(1); @@ -123,21 +122,18 @@ void Config(char directory[100]="", char option[6]="box") //============================================================= // Field (L3 0.4 T) - AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k4kG); - gAlice->SetField(field); - //gAlice->SetField(2,1) ; //(-999,2); - + AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k4kG); + gAlice->SetField(field); //=================== Alice BODY parameters ============================= AliBODY *BODY = new AliBODY("BODY","Alice envelop"); //=================== MUON Subsystem =========================== - AliMUONv1 *MUON = new AliMUONv1("MUON","default"); - + cout << ">>> Config_MUON_test.C: Creating AliMUONv1 ..."<LoadDigits("READ"); MUONLoader->LoadRecPoints("UPDATE"); + + // Testing if RawClusterisation has already been done + RunLoader->GetEvent(0); + if (MUONLoader->TreeR()) { + if (muondata->IsRawClusterBranchesInTree()) { + MUONLoader->UnloadRecPoints(); + MUONLoader->LoadRecPoints("RECREATE"); + printf("Recreating RecPoints files\n"); + } + } // // Loop over events // -- 2.39.3