changes in the MagF constructor
[u/mrichter/AliRoot.git] / HMPID / HESDfromKin.C
1 AliRunLoader *gAL=0; 
2 Int_t gEvt=0; Int_t gMaxEvt=0;
3 TObjArray *pNmean,*pQthre;
4 TTree *gEsdTr;
5 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6 void HESDfromKin(const char *name="default")
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();
18     AliESDEvent *pEsd = new AliESDEvent();   
19     TFile *pEsdFl=TFile::Open("AliESDs.root","recreate"); 
20     gEsdTr=new TTree("esdTree","Sim ESD from kinematics"); 
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        
25   }  else return;  
26
27   if(!OpenCalib()) {Printf("Problems in OpenCalib!Bye.");return;}
28     
29   TString ttl=name;
30   Bool_t htaCheck=ttl.Contains("HTA");
31 //  if(!htaCheck) SimEsd(pHL,pEsd); else SimEsdHidden(pHL,pEsd);
32   SimEsd(pHL,pEsd,htaCheck);
33   
34   pEsdFl->cd();
35   pEsdFl->Write();pEsdFl->Close();        
36 }
37 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
38 void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck)
39 {
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 }
57
58   InitGRP();
59 //  AliMagF *magFieldMap = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
60
61 //  AliHMPIDTracker::SetFieldMap(gAL->GetAliRun()->Field(),kTRUE);
62   AliHMPIDTracker pTracker;
63   AliHMPID *pH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID");
64   Int_t iNevt=gAL->GetNumberOfEvents();
65   pEsd->SetMagneticField(AliHMPIDTracker::GetBz());
66   for(Int_t iEvt=0;iEvt<iNevt;iEvt++){//events loop
67     gAL->GetEvent(iEvt);    
68     pHL->TreeR()->GetEntry(0);
69     AliStack *pStack=gAL->Stack();
70     Int_t nTrkHMPID=0;
71
72     for(Int_t i=0;i<pStack->GetNtrack();i++){
73       if(!pStack->IsPhysicalPrimary(i)) continue;
74       TParticle *pTrack=pStack->Particle(i); 
75       if(pTrack->GetPDG()->Charge()==0) continue;
76       Printf("track n. %i",i);
77       AliESDtrack trk(pTrack); 
78       Float_t xPc,yPc,xRa,yRa,thRa,phRa;
79       Int_t iCh=pTracker.IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa);         //get chamber intersected by this track 
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       }
86       nTrkHMPID++;
87       trk.SetHMPIDcluIdx   (iCh,99999);                                                          //chamber not found, mip not yet considered
88       
89       if(phRa<0) phRa += TMath::TwoPi(); // to be verified
90       
91       trk.SetHMPIDtrk(xPc,yPc,thRa,phRa);                                                        //store initial infos
92       pEsd->AddTrack(&trk);
93     
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
98 //      Printf("status %i",status);
99       if(status !=0) continue;
100
101     
102     }// track loop
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     
107     gEsdTr->Fill();
108     pEsd->Reset();
109   }// event loop
110   Printf("Events processed %i",iEvt);
111   if(htaCheck) {
112     fout->Write();
113     fout->Close();
114     delete fout;
115   }
116   gAL->UnloadHeader();  gAL->UnloadKinematics();
117 }//SimEsd()
118 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
119 Bool_t OpenCalib()
120 {
121   AliCDBManager* pCDB = AliCDBManager::Instance();
122   pCDB->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
123   pCDB->SetRun(0);
124   AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0);
125   AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0);
126   
127   if(!pQthreEnt || !pNmeanEnt) return kFALSE;
128   
129   pNmean=(TObjArray*)pNmeanEnt->GetObject(); 
130   pQthre=(TObjArray*)pQthreEnt->GetObject(); 
131
132   if(!pQthre || !pNmean) return kFALSE;  
133   return kTRUE;
134 }
135 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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   }
326   Printf("------------------------------");
327   Printf(" Summary for B: %s",s.Data());
328   Printf("------------------------------");
329   AliMagF* fld = new AliMagF("MagneticFieldMap", s.Data(), fcL3, fcDip, map, btype,beamenergy,path);
330   TGeoGlobalMagField::Instance()->SetField( fld );
331   TGeoGlobalMagField::Instance()->Lock();
332   //
333   return kTRUE;
334 }