Bug in magnetic field solved.
[u/mrichter/AliRoot.git] / HMPID / HESDfromKin.C
CommitLineData
d8770a0f 1AliRunLoader *gAL=0;
2Int_t gEvt=0; Int_t gMaxEvt=0;
3df18479 3TObjArray *pNmean,*pQthre;
d8770a0f 4TTree *gEsdTr;
5//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
e4290ede 6void 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 38void 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 119Bool_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 136Bool_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}