]>
Commit | Line | Data |
---|---|---|
d8770a0f | 1 | AliRunLoader *gAL=0; |
2 | Int_t gEvt=0; Int_t gMaxEvt=0; | |
3df18479 | 3 | TObjArray *pNmean,*pQthre; |
d8770a0f | 4 | TTree *gEsdTr; |
5 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
e4290ede | 6 | void HESDfromKin(const char *name="default") |
d8770a0f | 7 | {//simulate ESD from kinematics |
8 | ||
9 | if(gSystem->IsFileInIncludePath("galice.root")){// tries to open session | |
10 | if(gAlice) delete gAlice; //in case we execute this in aliroot delete default AliRun object | |
11 | gAL=AliRunLoader::Open(); //try to open galice.root from current dir | |
12 | gAL->LoadgAlice(); //take new AliRun object from galice.root | |
13 | if(gAL->LoadHeader()) return; | |
14 | if(gAL->LoadKinematics()) return; | |
15 | ||
16 | AliLoader *pHL=gAL->GetDetectorLoader("HMPID"); | |
17 | pHL->LoadRecPoints(); | |
6e06db1d | 18 | AliESDEvent *pEsd = new AliESDEvent(); |
d8770a0f | 19 | TFile *pEsdFl=TFile::Open("AliESDs.root","recreate"); |
20 | gEsdTr=new TTree("esdTree","Sim ESD from kinematics"); | |
e4290ede | 21 | pEsd->CreateStdContent(); pEsd->WriteToTree(gEsdTr); //clm: new ESD write schema: see Task Force meeting 20th June, 2007 |
22 | gEsdTr->GetUserInfo()->Add(pEsd); //clm: TList has to be created for ReadFromTree method -- this was not needed by the old ESD | |
23 | ||
24 | ||
d8770a0f | 25 | } else return; |
26 | ||
3df18479 | 27 | if(!OpenCalib()) {Printf("Problems in OpenCalib!Bye.");return;} |
d8770a0f | 28 | |
c0d7adf8 | 29 | TString ttl=name; |
30 | Bool_t htaCheck=ttl.Contains("HTA"); | |
9785d5fb | 31 | // if(!htaCheck) SimEsd(pHL,pEsd); else SimEsdHidden(pHL,pEsd); |
32 | SimEsd(pHL,pEsd,htaCheck); | |
e4290ede | 33 | |
34 | pEsdFl->cd(); | |
d8770a0f | 35 | pEsdFl->Write();pEsdFl->Close(); |
36 | } | |
37 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
9785d5fb | 38 | void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck) |
d8770a0f | 39 | { |
9785d5fb | 40 | if(htaCheck) { |
41 | TFile *fout = new TFile("HTA.root","recreate"); | |
42 | TH1F *hdC = new TH1F("dC" ,";delta Cerenkov (rad)",100,-0.2,0.2); | |
43 | TH1F *hCer = new TH1F("Cer" ,"Theta Cerenkov (rad)",250,0.,0.75); | |
44 | TH2F *htvsp = new TH2F("tvsp",";momentum (GeV/c);theta Cerenkov (rad)",100,0.,5.,1000,0.,0.75); | |
45 | TH1F *hdth = new TH1F("dth" ,";Delta theta Trk (mrad)",100,-250,250); | |
46 | TH1F *hdph = new TH1F("dph" ,";Delta phi Trk (mrad)",100,-500,500); | |
47 | Double_t rd=TMath::RadToDeg(); | |
48 | Printf("----------------------------------------------"); | |
49 | Printf("| SimHTA:Utility to embed ESD from kinematics|"); | |
50 | Printf("| with Hidden Track Algorithm (HTA) |"); | |
51 | Printf("----------------------------------------------"); | |
52 | } else { | |
53 | Printf("-----------------------------------------------"); | |
54 | Printf("| SimESD: Utility to embed ESD from kinematics|"); | |
55 | Printf("-----------------------------------------------"); | |
56 | } | |
efc29373 | 57 | |
58 | InitGRP(); | |
59 | // AliMagF *magFieldMap = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); | |
60 | ||
61 | // AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE); | |
e56b695f | 62 | AliHMPIDTracker pTracker; |
d8770a0f | 63 | AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID"); |
d8770a0f | 64 | Int_t iNevt=gAL->GetNumberOfEvents(); |
1b46ef29 | 65 | pEsd->SetMagneticField(AliHMPIDTracker::GetBz()); |
d8770a0f | 66 | for(Int_t iEvt=0;iEvt<iNevt;iEvt++){//events loop |
d8770a0f | 67 | gAL->GetEvent(iEvt); |
68 | pHL->TreeR()->GetEntry(0); | |
69 | AliStack *pStack=gAL->Stack(); | |
9785d5fb | 70 | Int_t nTrkHMPID=0; |
2b919fd5 | 71 | |
d8770a0f | 72 | for(Int_t i=0;i<pStack->GetNtrack();i++){ |
9785d5fb | 73 | if(!pStack->IsPhysicalPrimary(i)) continue; |
d8770a0f | 74 | TParticle *pTrack=pStack->Particle(i); |
1b46ef29 | 75 | if(pTrack->GetPDG()->Charge()==0) continue; |
97a4d538 | 76 | Printf("track n. %i",i); |
e4290ede | 77 | AliESDtrack trk(pTrack); |
39cd22e6 | 78 | Float_t xPc,yPc,xRa,yRa,thRa,phRa; |
e56b695f | 79 | Int_t iCh=pTracker.IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa); //get chamber intersected by this track |
39cd22e6 | 80 | if(iCh<0) { |
81 | trk.SetHMPIDtrk(0,0,0,0); //no intersection found | |
82 | trk.SetHMPIDcluIdx (99,99999); //chamber not found, mip not yet considered | |
83 | trk.SetHMPIDsignal(AliHMPIDRecon::kNotPerformed); //ring reconstruction not yet performed | |
84 | continue; //no intersection at all, go after next track | |
85 | } | |
9785d5fb | 86 | nTrkHMPID++; |
39cd22e6 | 87 | trk.SetHMPIDcluIdx (iCh,99999); //chamber not found, mip not yet considered |
e56b695f | 88 | |
89 | if(phRa<0) phRa += TMath::TwoPi(); // to be verified | |
90 | ||
97a4d538 | 91 | trk.SetHMPIDtrk(xPc,yPc,thRa,phRa); //store initial infos |
d8770a0f | 92 | pEsd->AddTrack(&trk); |
9785d5fb | 93 | |
6574baa2 | 94 | Int_t status; |
95 | if(!htaCheck) status = pTracker.Recon (pEsd,pH->CluLst(),pNmean,pQthre); | |
96 | else status = pTracker.ReconHiddenTrk(pEsd,pH->CluLst(),pNmean,pQthre); | |
97 | ||
2b919fd5 | 98 | // Printf("status %i",status); |
6574baa2 | 99 | if(status !=0) continue; |
9785d5fb | 100 | |
ac66a50d | 101 | |
6574baa2 | 102 | }// track loop |
2b919fd5 | 103 | |
104 | if(!(iEvt%50)) Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID); | |
105 | // Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID); | |
106 | ||
d8770a0f | 107 | gEsdTr->Fill(); |
108 | pEsd->Reset(); | |
109 | }// event loop | |
110 | Printf("Events processed %i",iEvt); | |
204299b8 | 111 | if(htaCheck) { |
112 | fout->Write(); | |
113 | fout->Close(); | |
114 | delete fout; | |
115 | } | |
d8770a0f | 116 | gAL->UnloadHeader(); gAL->UnloadKinematics(); |
e56b695f | 117 | }//SimEsd() |
d8770a0f | 118 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
3df18479 | 119 | Bool_t OpenCalib() |
d8770a0f | 120 | { |
121 | AliCDBManager* pCDB = AliCDBManager::Instance(); | |
efc29373 | 122 | pCDB->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); |
123 | pCDB->SetRun(0); | |
d8770a0f | 124 | AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0); |
125 | AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0); | |
126 | ||
3df18479 | 127 | if(!pQthreEnt || !pNmeanEnt) return kFALSE; |
d8770a0f | 128 | |
129 | pNmean=(TObjArray*)pNmeanEnt->GetObject(); | |
3df18479 | 130 | pQthre=(TObjArray*)pQthreEnt->GetObject(); |
131 | ||
132 | if(!pQthre || !pNmean) return kFALSE; | |
133 | return kTRUE; | |
d8770a0f | 134 | } |
135 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
efc29373 | 136 | Bool_t InitGRP() { |
137 | //------------------------------------ | |
138 | // Initialization of the GRP entry | |
139 | //------------------------------------ | |
140 | ||
141 | AliGRPObject *fGRPData; | |
142 | ||
143 | AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data"); | |
144 | ||
145 | if (entry) { | |
146 | ||
147 | TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry | |
148 | ||
149 | if (!m) { | |
150 | Printf("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject"); | |
151 | m->Print(); | |
152 | fGRPData = new AliGRPObject(); | |
153 | fGRPData->ReadValuesFromMap(m); | |
154 | } | |
155 | ||
156 | else { | |
157 | Printf("Found an AliGRPObject in GRP/GRP/Data, reading it"); | |
158 | fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry | |
159 | entry->SetOwner(0); | |
160 | } | |
161 | ||
162 | AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data"); | |
163 | } | |
164 | ||
165 | if (!fGRPData) { | |
166 | Printf("No GRP entry found in OCDB!"); | |
167 | return kFALSE; | |
168 | } | |
169 | ||
170 | TString lhcState = fGRPData->GetLHCState(); | |
171 | if (lhcState==AliGRPObject::GetInvalidString()) { | |
172 | Printf("GRP/GRP/Data entry: missing value for the LHC state ! Using UNKNOWN"); | |
173 | lhcState = "UNKNOWN"; | |
174 | } | |
175 | ||
176 | TString beamType = fGRPData->GetBeamType(); | |
177 | if (beamType==AliGRPObject::GetInvalidString()) { | |
178 | Printf("GRP/GRP/Data entry: missing value for the beam type ! Using UNKNOWN"); | |
179 | beamType = "UNKNOWN"; | |
180 | } | |
181 | ||
182 | Float_t beamEnergy = fGRPData->GetBeamEnergy(); | |
183 | if (beamEnergy==AliGRPObject::GetInvalidFloat()) { | |
184 | Printf("GRP/GRP/Data entry: missing value for the beam energy ! Using 0"); | |
185 | beamEnergy = 0; | |
186 | } | |
187 | // energy is provided in MeV*120 | |
188 | beamEnergy /= 120E3; | |
189 | ||
190 | TString runType = fGRPData->GetRunType(); | |
191 | if (runType==AliGRPObject::GetInvalidString()) { | |
192 | Printf("GRP/GRP/Data entry: missing value for the run type ! Using UNKNOWN"); | |
193 | runType = "UNKNOWN"; | |
194 | } | |
195 | ||
196 | Int_t activeDetectors = fGRPData->GetDetectorMask(); | |
197 | if (activeDetectors==AliGRPObject::GetInvalidUInt()) { | |
198 | Printf("GRP/GRP/Data entry: missing value for the detector mask ! Using 1074790399"); | |
199 | activeDetectors = 1074790399; | |
200 | } | |
201 | ||
202 | fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors); | |
203 | printf("qqqqqqqqqqqqqqqqqqqqqqq %s %s %f %s %d\n", lhcState.Data(), beamType.Data(), beamEnergy, runType.Data(), activeDetectors); | |
204 | fRunInfo->Dump(); | |
205 | ||
206 | //*** Dealing with the magnetic field map | |
207 | if ( TGeoGlobalMagField::Instance()->IsLocked() ) {Printf("Running with the externally locked B field !");} | |
208 | else { | |
209 | // Construct the field map out of the information retrieved from GRP. | |
210 | Bool_t ok = kTRUE; | |
211 | // L3 | |
212 | Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0); | |
213 | if (l3Current == AliGRPObject::GetInvalidFloat()) { | |
214 | Prtinf("GRP/GRP/Data entry: missing value for the L3 current !"); | |
215 | ok = kFALSE; | |
216 | } | |
217 | ||
218 | Char_t l3Polarity = fGRPData->GetL3Polarity(); | |
219 | if (l3Polarity == AliGRPObject::GetInvalidChar()) { | |
220 | Printf("GRP/GRP/Data entry: missing value for the L3 polarity !"); | |
221 | ok = kFALSE; | |
222 | } | |
223 | ||
224 | // Dipole | |
225 | Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0); | |
226 | if (diCurrent == AliGRPObject::GetInvalidFloat()) { | |
227 | Printf("GRP/GRP/Data entry: missing value for the dipole current !"); | |
228 | ok = kFALSE; | |
229 | } | |
230 | ||
231 | Char_t diPolarity = fGRPData->GetDipolePolarity(); | |
232 | if (diPolarity == AliGRPObject::GetInvalidChar()) { | |
233 | Printf("GRP/GRP/Data entry: missing value for the dipole polarity !"); | |
234 | ok = kFALSE; | |
235 | } | |
236 | ||
237 | if (ok) { | |
238 | if ( !SetFieldMap(l3Current, diCurrent, l3Polarity ? -1:1, diPolarity ? -1:1) ) | |
239 | AliFatal("Failed to creat a B field map ! Exiting..."); | |
240 | Printf("Running with the B field constructed out of GRP !"); | |
241 | } | |
242 | else AliFatal("B field is neither set nor constructed from GRP ! Exitig..."); | |
243 | ||
244 | } | |
245 | ||
246 | return kTRUE; | |
247 | } | |
248 | ||
249 | //_____________________________________________________________________________ | |
250 | //_____________________________________________________________________________ | |
251 | Bool_t SetFieldMap(Float_t l3Cur=30000., Float_t diCur=6000., | |
252 | Float_t l3Pol=1., Float_t diPol=1., Float_t beamenergy=7000., | |
253 | const Char_t* beamtype="pp", | |
254 | const Char_t* path="$(ALICE_ROOT)/data/maps/mfchebKGI_sym.root") | |
255 | { | |
256 | //------------------------------------------------ | |
257 | // The magnetic field map, defined externally... | |
258 | // L3 current 30000 A -> 0.5 T | |
259 | // L3 current 12000 A -> 0.2 T | |
260 | // dipole current 6000 A | |
261 | // The polarities must be the same | |
262 | //------------------------------------------------ | |
263 | const Float_t l3NominalCurrent1=30000.; // (A) | |
264 | const Float_t l3NominalCurrent2=12000.; // (A) | |
265 | const Float_t diNominalCurrent =6000. ; // (A) | |
266 | ||
267 | const Float_t tolerance=0.03; // relative current tolerance | |
268 | const Float_t zero=77.; // "zero" current (A) | |
269 | // | |
270 | TString s=(l3Pol < 0) ? "L3: -" : "L3: +"; | |
271 | // | |
272 | AliMagF::BMap_t map = AliMagF::k5kG; | |
273 | // | |
274 | double fcL3,fcDip; | |
275 | // | |
276 | l3Cur = TMath::Abs(l3Cur); | |
277 | if (TMath::Abs(l3Cur-l3NominalCurrent1)/l3NominalCurrent1 < tolerance) { | |
278 | fcL3 = l3Cur/l3NominalCurrent1; | |
279 | map = AliMagF::k5kG; | |
280 | s += "0.5 T; "; | |
281 | } else if (TMath::Abs(l3Cur-l3NominalCurrent2)/l3NominalCurrent2 < tolerance) { | |
282 | fcL3 = l3Cur/l3NominalCurrent2; | |
283 | map = AliMagF::k2kG; | |
284 | s += "0.2 T; "; | |
285 | } else if (l3Cur <= zero) { | |
286 | fcL3 = 0; | |
287 | map = AliMagF::k5kGUniform; | |
288 | s += "0.0 T; "; | |
289 | fUniformField=kTRUE; // track with the uniform (zero) B field | |
290 | } else { | |
291 | AliError(Form("Wrong L3 current (%f A)!",l3Cur)); | |
292 | return kFALSE; | |
293 | } | |
294 | // | |
295 | diCur = TMath::Abs(diCur); | |
296 | if (TMath::Abs(diCur-diNominalCurrent)/diNominalCurrent < tolerance) { | |
297 | // 3% current tolerance... | |
298 | fcDip = diCur/diNominalCurrent; | |
299 | s += "Dipole ON"; | |
300 | } else if (diCur <= zero) { // some small current.. | |
301 | fcDip = 0.; | |
302 | s += "Dipole OFF"; | |
303 | } else { | |
304 | AliError(Form("Wrong dipole current (%f A)!",diCur)); | |
305 | return kFALSE; | |
306 | } | |
307 | // | |
308 | if (l3Pol!=diPol && (map==AliMagF::k5kG || map==AliMagF::k2kG) && fcDip!=0) { | |
309 | AliError("L3 and Dipole polarities must be the same"); | |
310 | return kFALSE; | |
311 | } | |
312 | // | |
313 | if (l3Pol<0) fcL3 = -fcL3; | |
314 | if (diPol<0) fcDip = -fcDip; | |
315 | // | |
316 | AliMagF::BeamType_t btype = AliMagF::kNoBeamField; | |
317 | TString btypestr = beamtype; | |
318 | btypestr.ToLower(); | |
319 | TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1"); | |
320 | TPRegexp ionBeam("(lead|pb|ion|a)\\s*-?\\s*\\1"); | |
321 | if (btypestr.Contains(ionBeam)) btype = AliMagF::kBeamTypeAA; | |
322 | else if (btypestr.Contains(protonBeam)) btype = AliMagF::kBeamTypepp; | |
323 | else { | |
324 | Printf(Form("Cannot determine the beam type from %s, assume no LHC magnet field",beamtype)); | |
325 | } | |
97a4d538 | 326 | Printf("------------------------------"); |
327 | Printf(" Summary for B: %s",s.Data()); | |
328 | Printf("------------------------------"); | |
efc29373 | 329 | AliMagF* fld = new AliMagF("MagneticFieldMap", s.Data(), 2, fcL3, fcDip, 10., map, path, |
330 | btype,beamenergy); | |
331 | TGeoGlobalMagField::Instance()->SetField( fld ); | |
332 | TGeoGlobalMagField::Instance()->Lock(); | |
333 | // | |
334 | return kTRUE; | |
335 | } |