Printf("| SimESD: Utility to embed ESD from kinematics|");
Printf("-----------------------------------------------");
}
- AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE);
+
+ InitGRP();
+// AliMagF *magFieldMap = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+
+// AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE);
+ AliHMPIDTracker pTracker;
AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID");
Int_t iNevt=gAL->GetNumberOfEvents();
pEsd->SetMagneticField(AliHMPIDTracker::GetBz());
pHL->TreeR()->GetEntry(0);
AliStack *pStack=gAL->Stack();
Int_t nTrkHMPID=0;
+
for(Int_t i=0;i<pStack->GetNtrack();i++){
if(!pStack->IsPhysicalPrimary(i)) continue;
TParticle *pTrack=pStack->Particle(i);
if(pTrack->GetPDG()->Charge()==0) continue;
+ Printf("track n. %i",i);
AliESDtrack trk(pTrack);
Float_t xPc,yPc,xRa,yRa,thRa,phRa;
- Int_t iCh=AliHMPIDTracker::IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa); //get chamber intersected by this track
+ Int_t iCh=pTracker.IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa); //get chamber intersected by this track
if(iCh<0) {
trk.SetHMPIDtrk(0,0,0,0); //no intersection found
trk.SetHMPIDcluIdx (99,99999); //chamber not found, mip not yet considered
}
nTrkHMPID++;
trk.SetHMPIDcluIdx (iCh,99999); //chamber not found, mip not yet considered
- trk.SetHMPIDtrk(xRa,yRa,thRa,phRa); //store initial infos
- pEsd->AddTrack(&trk);
- }// track loop
-
- if(!(iEvt%50)) Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID);
- if(!htaCheck) AliHMPIDTracker::Recon(pEsd,pH->CluLst(),pNmean,pQthre);
-// Perform PID
-
-
- for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){ //ESD tracks loop
- AliESDtrack *pTrk = pEsd->GetTrack(iTrk); //get reconstructed track
+ if(phRa<0) phRa += TMath::TwoPi(); // to be verified
- Float_t radX,radY,thetaTrk,phiTrk;
- pTrk->GetHMPIDtrk(radX,radY,thetaTrk,phiTrk);
-
- if(htaCheck) {
- Int_t iCh = pTrk->GetHMPIDcluIdx(); //chamber not found, mip not yet considered
- iCh/=1000000;
- // Printf("simulated track theta %f phi %f",thetaTrk*rd,phiTrk*rd);
- TObjArray *pClus = pH->CluLst();
-
- if(AliHMPIDTracker::ReconHiddenTrk(pEsd,pClus,pNmean,pQthre)!=0) continue;
-
- Double_t thetaCerSim = TMath::ACos(pTrack->Energy()/(1.292*pTrack->P()));
- Printf("-----------------------------------------------------------");
- Printf(" theta Cerenkov simulated %f",thetaCerSim);
- Printf(" theta Cerenkov reconstructed %f",pTrk->GetHMPIDsignal());
- Float_t thetaFit,phiFit;
- pTrk->GetHMPIDtrk(radX,radY,thetaFit,phiFit);
-// Printf("reconstr. track theta %f phi %f",thetaFit*rd,phiFit*rd);
-
- // fill useful histos
- hdC->Fill(pTrk->GetHMPIDsignal()-thetaCerSim);
- hCer->Fill(pTrk->GetHMPIDsignal());
- htvsp->Fill(pTrk->P(),pTrk->GetHMPIDsignal());
- hdth->Fill((thetaFit-thetaTrk)*1000);
- hdph->Fill((phiFit-phiTrk)*1000);
- }
-
- AliHMPIDPid pID;
- Double_t prob[5];
- pID.FindPid(pTrk,5,prob);
- pTrk->SetHMPIDpid(prob);
-// Printf(" Prob e- %6.2f mu %6.2f pi %6.2f k %6.2f p %6.2f",prob[0]*100,prob[1]*100,prob[2]*100,prob[3]*100,prob[4]*100);
- }
-
- gEsdTr->Fill();
- pEsd->Reset();
+ trk.SetHMPIDtrk(xPc,yPc,thRa,phRa); //store initial infos
+ pEsd->AddTrack(&trk);
- }// event loop
- Printf("Events processed %i",iEvt);
- gAL->UnloadHeader(); gAL->UnloadKinematics();
-}//Esd()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void SimEsdHidden(AliLoader *pHL,AliESDEvent *pEsd)
-{
- TFile *fout = new TFile("HTA.root","recreate");
- TH1F *hdC = new TH1F("dC" ,";delta Cerenkov (rad)",100,-0.2,0.2);
- TH1F *hCer = new TH1F("Cer" ,"Theta Cerenkov (rad)",250,0.,0.75);
- TH2F *htvsp = new TH2F("tvsp",";momentum (GeV/c);theta Cerenkov (rad)",100,0.,5.,1000,0.,0.75);
- TH1F *hdth = new TH1F("dth" ,";Delta theta Trk (mrad)",100,-250,250);
- TH1F *hdph = new TH1F("dph" ,";Delta phi Trk (mrad)",100,-500,500);
- Double_t rd=TMath::RadToDeg();
- Printf("----------------------------------------------");
- Printf("| SimHTA:Utility to embed ESD from kinematics|");
- Printf("| with Hidden Track Algorithm (HTA) |");
- Printf("----------------------------------------------");
- AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE);
- AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID");
- Int_t iNevt=gAL->GetNumberOfEvents();
- pEsd->SetMagneticField(AliHMPIDTracker::GetBz());
-
- Printf("Number of events to process: %i",iNevt);
-
- for(Int_t iEvt=0;iEvt<iNevt;iEvt++){//events loop
- if(!(iEvt%1)) Printf("Events processed %i",iEvt);
- gAL->GetEvent(iEvt);
- pHL->TreeR()->GetEntry(0);
- AliStack *pStack=gAL->Stack();
+ Int_t status;
+ if(!htaCheck) status = pTracker.Recon (pEsd,pH->CluLst(),pNmean,pQthre);
+ else status = pTracker.ReconHiddenTrk(pEsd,pH->CluLst(),pNmean,pQthre);
+
+// Printf("status %i",status);
+ if(status !=0) continue;
+
- for(Int_t i=0;i<pStack->GetNtrack();i++){
-
- (!pStack->IsPhysicalPrimary(i)) continue;
- TParticle *pTrack=pStack->Particle(i);
-
- if(pTrack->GetPDG()->Charge()==0) continue;
- //find the chamber that intersects HMPID
- pTrack->Print();
- AliESDtrack trk(pTrack);
- Float_t xPc,yPc,xRa,yRa,thRa,phRa;
- Int_t iCh=AliHMPIDTracker::IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa); //get chamber intersected by this track
- if(iCh<0) continue; //no intersection at all, go after next track
- trk.SetHMPIDtrk(xRa,yRa,thRa,phRa); //store initial infos
- Float_t radX,radY,thetaTrk,phiTrk;
- trk.GetHMPIDtrk(radX,radY,thetaTrk,phiTrk);
-// Printf("simulated track theta %f phi %f",thetaTrk*rd,phiTrk*rd);
- TObjArray *pClus = pH->CluLst();
-
- if(AliHMPIDTracker::ReconHiddenTrk(iCh,0,&trk,(TClonesArray *)pClus->At(iCh),pNmean,pQthre)!=0) continue;
-
- pEsd->AddTrack(&trk);
- Double_t thetaCerSim = TMath::ACos(pTrack->Energy()/(1.292*pTrack->P()));
-// Printf(" theta Cerenkov simulated %f",thetaCerSim);
-// Printf(" theta Cerenkov reconstructed %f",trk.GetHMPIDsignal());
- Float_t thetaFit,phiFit;
- trk.GetHMPIDtrk(radX,radY,thetaFit,phiFit);
-// Printf("reconstr. track theta %f phi %f",thetaFit*rd,phiFit*rd);
-
- // fill useful histos
- hdC->Fill(trk.GetHMPIDsignal()-thetaCerSim);
- hCer->Fill(trk.GetHMPIDsignal());
- htvsp->Fill(pTrack->P(),trk.GetHMPIDsignal());
- hdth->Fill((thetaFit-thetaTrk)*1000);
- hdph->Fill((phiFit-phiTrk)*1000);
}// track loop
- pEsd->SetMagneticField(AliHMPIDTracker::GetBz());
+
+ if(!(iEvt%50)) Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID);
+// Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID);
+
gEsdTr->Fill();
pEsd->Reset();
}// event loop
Printf("Events processed %i",iEvt);
- fout->Write();
- fout->Close();
- delete fout;
+ if(htaCheck) {
+ fout->Write();
+ fout->Close();
+ delete fout;
+ }
gAL->UnloadHeader(); gAL->UnloadKinematics();
-}//EsdHidden()
+}//SimEsd()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Bool_t OpenCalib()
{
AliCDBManager* pCDB = AliCDBManager::Instance();
- pCDB->SetDefaultStorage("local://$HOME");
+ pCDB->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ pCDB->SetRun(0);
AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0);
AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0);
return kTRUE;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Bool_t InitGRP() {
+ //------------------------------------
+ // Initialization of the GRP entry
+ //------------------------------------
+
+ AliGRPObject *fGRPData;
+
+ AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
+
+ if (entry) {
+
+ TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry
+
+ if (!m) {
+ Printf("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
+ m->Print();
+ fGRPData = new AliGRPObject();
+ fGRPData->ReadValuesFromMap(m);
+ }
+
+ else {
+ Printf("Found an AliGRPObject in GRP/GRP/Data, reading it");
+ fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry
+ entry->SetOwner(0);
+ }
+
+ AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
+ }
+
+ if (!fGRPData) {
+ Printf("No GRP entry found in OCDB!");
+ return kFALSE;
+ }
+
+ TString lhcState = fGRPData->GetLHCState();
+ if (lhcState==AliGRPObject::GetInvalidString()) {
+ Printf("GRP/GRP/Data entry: missing value for the LHC state ! Using UNKNOWN");
+ lhcState = "UNKNOWN";
+ }
+
+ TString beamType = fGRPData->GetBeamType();
+ if (beamType==AliGRPObject::GetInvalidString()) {
+ Printf("GRP/GRP/Data entry: missing value for the beam type ! Using UNKNOWN");
+ beamType = "UNKNOWN";
+ }
+
+ Float_t beamEnergy = fGRPData->GetBeamEnergy();
+ if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
+ Printf("GRP/GRP/Data entry: missing value for the beam energy ! Using 0");
+ beamEnergy = 0;
+ }
+ // energy is provided in MeV*120
+ beamEnergy /= 120E3;
+
+ TString runType = fGRPData->GetRunType();
+ if (runType==AliGRPObject::GetInvalidString()) {
+ Printf("GRP/GRP/Data entry: missing value for the run type ! Using UNKNOWN");
+ runType = "UNKNOWN";
+ }
+
+ Int_t activeDetectors = fGRPData->GetDetectorMask();
+ if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
+ Printf("GRP/GRP/Data entry: missing value for the detector mask ! Using 1074790399");
+ activeDetectors = 1074790399;
+ }
+
+ fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
+ printf("qqqqqqqqqqqqqqqqqqqqqqq %s %s %f %s %d\n", lhcState.Data(), beamType.Data(), beamEnergy, runType.Data(), activeDetectors);
+ fRunInfo->Dump();
+
+ //*** Dealing with the magnetic field map
+ if ( TGeoGlobalMagField::Instance()->IsLocked() ) {Printf("Running with the externally locked B field !");}
+ else {
+ // Construct the field map out of the information retrieved from GRP.
+ Bool_t ok = kTRUE;
+ // L3
+ Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
+ if (l3Current == AliGRPObject::GetInvalidFloat()) {
+ Prtinf("GRP/GRP/Data entry: missing value for the L3 current !");
+ ok = kFALSE;
+ }
+
+ Char_t l3Polarity = fGRPData->GetL3Polarity();
+ if (l3Polarity == AliGRPObject::GetInvalidChar()) {
+ Printf("GRP/GRP/Data entry: missing value for the L3 polarity !");
+ ok = kFALSE;
+ }
+
+ // Dipole
+ Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
+ if (diCurrent == AliGRPObject::GetInvalidFloat()) {
+ Printf("GRP/GRP/Data entry: missing value for the dipole current !");
+ ok = kFALSE;
+ }
+
+ Char_t diPolarity = fGRPData->GetDipolePolarity();
+ if (diPolarity == AliGRPObject::GetInvalidChar()) {
+ Printf("GRP/GRP/Data entry: missing value for the dipole polarity !");
+ ok = kFALSE;
+ }
+
+ if (ok) {
+ if ( !SetFieldMap(l3Current, diCurrent, l3Polarity ? -1:1, diPolarity ? -1:1) )
+ AliFatal("Failed to creat a B field map ! Exiting...");
+ Printf("Running with the B field constructed out of GRP !");
+ }
+ else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
+
+ }
+
+ return kTRUE;
+}
+
+//_____________________________________________________________________________
+//_____________________________________________________________________________
+ Bool_t SetFieldMap(Float_t l3Cur=30000., Float_t diCur=6000.,
+ Float_t l3Pol=1., Float_t diPol=1., Float_t beamenergy=7000.,
+ const Char_t* beamtype="pp",
+ const Char_t* path="$(ALICE_ROOT)/data/maps/mfchebKGI_sym.root")
+{
+ //------------------------------------------------
+ // The magnetic field map, defined externally...
+ // L3 current 30000 A -> 0.5 T
+ // L3 current 12000 A -> 0.2 T
+ // dipole current 6000 A
+ // The polarities must be the same
+ //------------------------------------------------
+ const Float_t l3NominalCurrent1=30000.; // (A)
+ const Float_t l3NominalCurrent2=12000.; // (A)
+ const Float_t diNominalCurrent =6000. ; // (A)
+
+ const Float_t tolerance=0.03; // relative current tolerance
+ const Float_t zero=77.; // "zero" current (A)
+ //
+ TString s=(l3Pol < 0) ? "L3: -" : "L3: +";
+ //
+ AliMagF::BMap_t map = AliMagF::k5kG;
+ //
+ double fcL3,fcDip;
+ //
+ l3Cur = TMath::Abs(l3Cur);
+ if (TMath::Abs(l3Cur-l3NominalCurrent1)/l3NominalCurrent1 < tolerance) {
+ fcL3 = l3Cur/l3NominalCurrent1;
+ map = AliMagF::k5kG;
+ s += "0.5 T; ";
+ } else if (TMath::Abs(l3Cur-l3NominalCurrent2)/l3NominalCurrent2 < tolerance) {
+ fcL3 = l3Cur/l3NominalCurrent2;
+ map = AliMagF::k2kG;
+ s += "0.2 T; ";
+ } else if (l3Cur <= zero) {
+ fcL3 = 0;
+ map = AliMagF::k5kGUniform;
+ s += "0.0 T; ";
+ fUniformField=kTRUE; // track with the uniform (zero) B field
+ } else {
+ AliError(Form("Wrong L3 current (%f A)!",l3Cur));
+ return kFALSE;
+ }
+ //
+ diCur = TMath::Abs(diCur);
+ if (TMath::Abs(diCur-diNominalCurrent)/diNominalCurrent < tolerance) {
+ // 3% current tolerance...
+ fcDip = diCur/diNominalCurrent;
+ s += "Dipole ON";
+ } else if (diCur <= zero) { // some small current..
+ fcDip = 0.;
+ s += "Dipole OFF";
+ } else {
+ AliError(Form("Wrong dipole current (%f A)!",diCur));
+ return kFALSE;
+ }
+ //
+ if (l3Pol!=diPol && (map==AliMagF::k5kG || map==AliMagF::k2kG) && fcDip!=0) {
+ AliError("L3 and Dipole polarities must be the same");
+ return kFALSE;
+ }
+ //
+ if (l3Pol<0) fcL3 = -fcL3;
+ if (diPol<0) fcDip = -fcDip;
+ //
+ AliMagF::BeamType_t btype = AliMagF::kNoBeamField;
+ TString btypestr = beamtype;
+ btypestr.ToLower();
+ TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
+ TPRegexp ionBeam("(lead|pb|ion|a)\\s*-?\\s*\\1");
+ if (btypestr.Contains(ionBeam)) btype = AliMagF::kBeamTypeAA;
+ else if (btypestr.Contains(protonBeam)) btype = AliMagF::kBeamTypepp;
+ else {
+ Printf(Form("Cannot determine the beam type from %s, assume no LHC magnet field",beamtype));
+ }
+ Printf("------------------------------");
+ Printf(" Summary for B: %s",s.Data());
+ Printf("------------------------------");
+ AliMagF* fld = new AliMagF("MagneticFieldMap", s.Data(), 2, fcL3, fcDip, 10., map, path,
+ btype,beamenergy);
+ TGeoGlobalMagField::Instance()->SetField( fld );
+ TGeoGlobalMagField::Instance()->Lock();
+ //
+ return kTRUE;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff