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       AliESDtrack trk(pTrack);
  77       Float_t xPc,yPc,xRa,yRa,thRa,phRa;
  78       Int_t iCh=pTracker.IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa);         //get chamber intersected by this track
  79       if(iCh<0) {
  80         trk.SetHMPIDtrk(0,0,0,0);                                                                //no intersection found
  81         trk.SetHMPIDcluIdx   (99,99999);                                                         //chamber not found, mip not yet considered
  82         trk.SetHMPIDsignal(AliHMPIDRecon::kNotPerformed);                                        //ring reconstruction not yet performed
  83         continue;                                                           //no intersection at all, go after next track
  84       }
  85       nTrkHMPID++;
  86       trk.SetHMPIDcluIdx   (iCh,99999);                                                          //chamber not found, mip not yet considered
  87
  88       if(phRa<0) phRa += TMath::TwoPi(); // to be verified
  89
  90       trk.SetHMPIDtrk(xRa,yRa,thRa,phRa);                                                        //store initial infos
  91       pEsd->AddTrack(&trk);
  92
  93       Int_t status;
  94       if(!htaCheck) status = pTracker.Recon         (pEsd,pH->CluLst(),pNmean,pQthre);
  95       else          status = pTracker.ReconHiddenTrk(pEsd,pH->CluLst(),pNmean,pQthre);
  96
  97 //      Printf("status %i",status);
  98       if(status !=0) continue;
  99
 100
 101     }// track loop
 102
 103     if(!(iEvt%50)) Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID);
 104 //      Printf("Number of events processed: %i with tracks %i in HMPID",iEvt,nTrkHMPID);
 105
 106     gEsdTr->Fill();
 107     pEsd->Reset();
 108   }// event loop
 109   Printf("Events processed %i",iEvt);
 110   if(htaCheck) {
 111     fout->Write();
 112     fout->Close();
 113     delete fout;
 114   }
 115   gAL->UnloadHeader();  gAL->UnloadKinematics();
 116 }//SimEsd()
 117 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 118 Bool_t OpenCalib()
 119 {
 120   AliCDBManager* pCDB = AliCDBManager::Instance();
 121   pCDB->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 122   pCDB->SetRun(0);
 123   AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0);
 124   AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0);
 125
 126   if(!pQthreEnt || !pNmeanEnt) return kFALSE;
 127
 128   pNmean=(TObjArray*)pNmeanEnt->GetObject();
 129   pQthre=(TObjArray*)pQthreEnt->GetObject();
 130
 131   if(!pQthre || !pNmean) return kFALSE;
 132   return kTRUE;
 133 }
 134 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 135 Bool_t InitGRP() {
 136   //------------------------------------
 137   // Initialization of the GRP entry
 138   //------------------------------------
 139
 140   AliGRPObject *fGRPData;
 141
 142   AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
 143
 144   if (entry) {
 145
 146     TMap* m = dynamic_cast<TMap*>(entry->GetObject());  // old GRP entry
 147
 148     if (!m) {
 149        Printf("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
 150        m->Print();
 151        fGRPData = new AliGRPObject();
 152        fGRPData->ReadValuesFromMap(m);
 153     }
 154
 155     else {
 156        Printf("Found an AliGRPObject in GRP/GRP/Data, reading it");
 157        fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject());  // new GRP entry
 158        entry->SetOwner(0);
 159     }
 160
 161     AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
 162   }
 163
 164   if (!fGRPData) {
 165      Printf("No GRP entry found in OCDB!");
 166      return kFALSE;
 167   }
 168
 169   TString lhcState = fGRPData->GetLHCState();
 170   if (lhcState==AliGRPObject::GetInvalidString()) {
 171     Printf("GRP/GRP/Data entry:  missing value for the LHC state ! Using UNKNOWN");
 172     lhcState = "UNKNOWN";
 173   }
 174
 175   TString beamType = fGRPData->GetBeamType();
 176   if (beamType==AliGRPObject::GetInvalidString()) {
 177     Printf("GRP/GRP/Data entry:  missing value for the beam type ! Using UNKNOWN");
 178     beamType = "UNKNOWN";
 179   }
 180
 181   Float_t beamEnergy = fGRPData->GetBeamEnergy();
 182   if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
 183     Printf("GRP/GRP/Data entry:  missing value for the beam energy ! Using 0");
 184     beamEnergy = 0;
 185   }
 186   // energy is provided in MeV*120
 187   beamEnergy /= 120E3;
 188
 189   TString runType = fGRPData->GetRunType();
 190   if (runType==AliGRPObject::GetInvalidString()) {
 191     Printf("GRP/GRP/Data entry:  missing value for the run type ! Using UNKNOWN");
 192     runType = "UNKNOWN";
 193   }
 194
 195   Int_t activeDetectors = fGRPData->GetDetectorMask();
 196   if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
 197     Printf("GRP/GRP/Data entry:  missing value for the detector mask ! Using 1074790399");
 198     activeDetectors = 1074790399;
 199   }
 200
 201   fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
 202   printf("qqqqqqqqqqqqqqqqqqqqqqq %s %s %f %s %d\n", lhcState.Data(), beamType.Data(), beamEnergy, runType.Data(), activeDetectors);
 203   fRunInfo->Dump();
 204
 205   //*** Dealing with the magnetic field map
 206   if ( TGeoGlobalMagField::Instance()->IsLocked() ) {Printf("Running with the externally locked B field !");}
 207   else {
 208     // Construct the field map out of the information retrieved from GRP.
 209     Bool_t ok = kTRUE;
 210     // L3
 211     Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
 212     if (l3Current == AliGRPObject::GetInvalidFloat()) {
 213       Prtinf("GRP/GRP/Data entry:  missing value for the L3 current !");
 214       ok = kFALSE;
 215     }
 216
 217     Char_t l3Polarity = fGRPData->GetL3Polarity();
 218     if (l3Polarity == AliGRPObject::GetInvalidChar()) {
 219       Printf("GRP/GRP/Data entry:  missing value for the L3 polarity !");
 220       ok = kFALSE;
 221     }
 222
 223     // Dipole
 224     Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
 225     if (diCurrent == AliGRPObject::GetInvalidFloat()) {
 226       Printf("GRP/GRP/Data entry:  missing value for the dipole current !");
 227       ok = kFALSE;
 228     }
 229
 230     Char_t diPolarity = fGRPData->GetDipolePolarity();
 231     if (diPolarity == AliGRPObject::GetInvalidChar()) {
 232       Printf("GRP/GRP/Data entry:  missing value for the dipole polarity !");
 233       ok = kFALSE;
 234     }
 235
 236     if (ok) {
 237       if ( !SetFieldMap(l3Current, diCurrent, l3Polarity ? -1:1, diPolarity ? -1:1) )
 238         AliFatal("Failed to creat a B field map ! Exiting...");
 239       Printf("Running with the B field constructed out of GRP !");
 240     }
 241     else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
 242
 243   }
 244
 245   return kTRUE;
 246 }
 247
 248 //_____________________________________________________________________________
 249 //_____________________________________________________________________________
 250   Bool_t SetFieldMap(Float_t l3Cur=30000., Float_t diCur=6000.,
 251                      Float_t l3Pol=1., Float_t diPol=1., Float_t beamenergy=7000.,
 252                      const Char_t* beamtype="pp",
 253                      const Char_t* path="$(ALICE_ROOT)/data/maps/mfchebKGI_sym.root")
 254 {
 255   //------------------------------------------------
 256   // The magnetic field map, defined externally...
 257   // L3 current 30000 A  -> 0.5 T
 258   // L3 current 12000 A  -> 0.2 T
 259   // dipole current 6000 A
 260   // The polarities must be the same
 261   //------------------------------------------------
 262   const Float_t l3NominalCurrent1=30000.; // (A)
 263   const Float_t l3NominalCurrent2=12000.; // (A)
 264   const Float_t diNominalCurrent =6000. ; // (A)
 265
 266   const Float_t tolerance=0.03; // relative current tolerance
 267   const Float_t zero=77.;       // "zero" current (A)
 268   //
 269   TString s=(l3Pol < 0) ? "L3: -" : "L3: +";
 270   //
 271   AliMagF::BMap_t map = AliMagF::k5kG;
 272   //
 273   double fcL3,fcDip;
 274   //
 275   l3Cur = TMath::Abs(l3Cur);
 276   if (TMath::Abs(l3Cur-l3NominalCurrent1)/l3NominalCurrent1 < tolerance) {
 277     fcL3 = l3Cur/l3NominalCurrent1;
 278     map  = AliMagF::k5kG;
 279     s   += "0.5 T;  ";
 280   } else if (TMath::Abs(l3Cur-l3NominalCurrent2)/l3NominalCurrent2 < tolerance) {
 281     fcL3 = l3Cur/l3NominalCurrent2;
 282     map  = AliMagF::k2kG;
 283     s   += "0.2 T;  ";
 284   } else if (l3Cur <= zero) {
 285     fcL3 = 0;
 286     map  = AliMagF::k5kGUniform;
 287     s   += "0.0 T;  ";
 288     fUniformField=kTRUE;        // track with the uniform (zero) B field
 289   } else {
 290     AliError(Form("Wrong L3 current (%f A)!",l3Cur));
 291     return kFALSE;
 292   }
 293   //
 294   diCur = TMath::Abs(diCur);
 295   if (TMath::Abs(diCur-diNominalCurrent)/diNominalCurrent < tolerance) {
 296     // 3% current tolerance...
 297     fcDip = diCur/diNominalCurrent;
 298     s    += "Dipole ON";
 299   } else if (diCur <= zero) { // some small current..
 300     fcDip = 0.;
 301     s    += "Dipole OFF";
 302   } else {
 303     AliError(Form("Wrong dipole current (%f A)!",diCur));
 304     return kFALSE;
 305   }
 306   //
 307   if (l3Pol!=diPol && (map==AliMagF::k5kG || map==AliMagF::k2kG) && fcDip!=0) {
 308     AliError("L3 and Dipole polarities must be the same");
 309     return kFALSE;
 310   }
 311   //
 312   if (l3Pol<0) fcL3  = -fcL3;
 313   if (diPol<0) fcDip = -fcDip;
 314   //
 315   AliMagF::BeamType_t btype = AliMagF::kNoBeamField;
 316   TString btypestr = beamtype;
 317   btypestr.ToLower();
 318   TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
 319   TPRegexp ionBeam("(lead|pb|ion|a)\\s*-?\\s*\\1");
 320   if (btypestr.Contains(ionBeam)) btype = AliMagF::kBeamTypeAA;
 321   else if (btypestr.Contains(protonBeam)) btype = AliMagF::kBeamTypepp;
 322   else {
 323     Printf(Form("Cannot determine the beam type from %s, assume no LHC magnet field",beamtype));
 324   }
 325
 326   AliMagF* fld = new AliMagF("MagneticFieldMap", s.Data(), 2, fcL3, fcDip, 10., map, path,
 327                              btype,beamenergy);
 328   TGeoGlobalMagField::Instance()->SetField( fld );
 329   TGeoGlobalMagField::Instance()->Lock();
 330   //
 331   return kTRUE;
 332 }