FMD/AliFMDBackgroundCorrection.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2004, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // Thil class computes background corrections for the FMD. The correction is computed
  17 // in eta,phi cells and the objects stored can be put into alien to use with analysis.
  18 // It is based on the AliFMDInput class that is used to loop over hits and primaries.
  19 //
  20 // Author: Hans Hjersing Dalsgaard, NBI, hans.dalsgaard@cern.ch
  21 //
  22 //
  23
  24 #include "AliSimulation.h"
  25 #include "TStopwatch.h"
  26 #include "iostream"
  27 #include "TGrid.h"
  28 #include "AliRunLoader.h"
  29 #include "AliGeomManager.h"
  30 #include "AliFMDGeometry.h"
  31 #include "AliStack.h"
  32 #include "TParticle.h"
  33 #include "TDatabasePDG.h"
  34 #include "TParticlePDG.h"
  35 #include "AliRun.h"
  36 #include "AliFMDBackgroundCorrection.h"
  37 #include "TSystem.h"
  38 #include "AliCDBManager.h"
  39 #include "TTree.h"
  40 #include "TClonesArray.h"
  41 #include "TBranch.h"
  42 #include "AliFMDHit.h"
  43 #include "AliLoader.h"
  44 #include "AliFMD.h"
  45 #include "AliGenEventHeader.h"
  46 #include "AliHeader.h"
  47 #include "TFile.h"
  48 #include "TAxis.h"
  49 #include "AliCDBId.h"
  50 #include "AliCDBMetaData.h"
  51 #include "TROOT.h"
  52 #include "AliFMDParameters.h"
  53 #include "AliLog.h"
  54 #include "TList.h"
  55 #include "AliFMDAnaParameters.h"
  56 #include "AliFMDAnaCalibBackgroundCorrection.h"
  57 #include "AliTrackReference.h"
  58 #include "AliFMDStripIndex.h"
  59
  60 ClassImp(AliFMDBackgroundCorrection)
  61 //_____________________________________________________________________
  62 AliFMDBackgroundCorrection::AliFMDBackgroundCorrection() :
  63   TNamed(),
  64   fCorrectionArray(),
  65   fPrimaryList()
  66 {}
  67
  68 //_____________________________________________________________________
  69 AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG(Bool_t hits_not_trackref) :
  70   AliFMDInput("galice.root"),
  71   fPrimaryArray(),
  72   fHitArray(),
  73   fHitMap(),
  74   fPrim(0),
  75   fHits(0),
  76   fZvtxCut(0),
  77   fNvtxBins(0),
  78   fPrevTrack(-1),
  79   fPrevDetector(-1),
  80   fPrevRing('Q'),
  81   fPrevSec(-1),
  82   fNbinsEta(100)
  83 {
  84   AddLoad(kTracks);
  85   if(hits_not_trackref)
  86     AddLoad(kHits);
  87   else
  88     AddLoad(kTrackRefs);
  89   AddLoad(kKinematics);
  90   AddLoad(kHeader);
  91
  92
  93 }
  94
  95 //_____________________________________________________________________
  96
  97 void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Bool_t from_hits,
  98                                                               const Int_t nvtxbins,
  99                                                               Float_t zvtxcut,
 100                                                               const Int_t nBinsEta,
 101                                                               Bool_t storeInAlien,
 102                                                               Int_t runNo,
 103                                                               Int_t endRunNo,
 104                                                               const Char_t* filename,
 105                                                               Bool_t simulate,
 106                                                               Int_t nEvents,
 107                                                               Bool_t inFile,
 108                                                               const Char_t* infilename) {
 109
 110   //TGrid::Connect("alien:",0,0,"t");
 111   if(simulate)
 112     Simulate(nEvents);
 113   else {
 114     //AliCDBManager::Instance()->SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
 115     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
 116     AliCDBManager::Instance()->SetRun(runNo);
 117
 118 #if defined(__CINT__)
 119     gSystem->Load("liblhapdf");
 120     gSystem->Load("libEGPythia6");
 121     gSystem->Load("libpythia6");
 122     gSystem->Load("libAliPythia6");
 123     gSystem->Load("libgeant321");
 124 #endif
 125       //
 126   }
 127
 128   //Setting up the geometry
 129   //-----------------------------------------------
 130   if (AliGeomManager::GetGeometry() == NULL)
 131     AliGeomManager::LoadGeometry();
 132
 133   AliFMDGeometry* geo = AliFMDGeometry::Instance();
 134   geo->Init();
 135   geo->InitTransformations();
 136
 137   AliInfo("Processing hits and primaries ");
 138
 139   AliFMDInputBG input(from_hits);
 140
 141   if(!inFile) {
 142
 143     input.SetVtxCutZ(zvtxcut);
 144     input.SetNvtxBins(nvtxbins);
 145     input.SetNbinsEta(nBinsEta);
 146     input.Run();
 147   }
 148
 149   AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
 150   TObjArray* hitArray ;
 151   TObjArray* primaryArray;
 152   if(inFile) {
 153     TFile* infile = TFile::Open(infilename);
 154     hitArray     = new TObjArray();
 155     primaryArray = new TObjArray();
 156
 157     for(Int_t det =1; det<=3;det++)
 158       {
 159         TObjArray* detArrayHits = new TObjArray();
 160         detArrayHits->SetName(Form("FMD%d",det));
 161         hitArray->AddAtAndExpand(detArrayHits,det);
 162         Int_t nRings = (det==1 ? 1 : 2);
 163         for(Int_t ring = 0;ring<nRings;ring++)
 164           {
 165
 166             Char_t ringChar = (ring == 0 ? 'I' : 'O');
 167             TObjArray* vtxArrayHits = new TObjArray();
 168             vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
 169             detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
 170             for(Int_t v=0; v<nvtxbins;v++)
 171               {
 172
 173                 TH2F* hHits          = (TH2F*)infile->Get(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v));
 174
 175
 176               vtxArrayHits->AddAtAndExpand(hHits,v);
 177
 178             }
 179         }
 180
 181     }
 182
 183     for(Int_t iring = 0; iring<2;iring++) {
 184       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 185       TObjArray* ringArray = new TObjArray();
 186       ringArray->SetName(Form("FMD_%c",ringChar));
 187       primaryArray->AddAtAndExpand(ringArray,iring);
 188       for(Int_t v=0; v<nvtxbins;v++) {
 189
 190         TH2F* hPrimary       = (TH2F*)infile->Get(Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
 191         ringArray->AddAtAndExpand(hPrimary,v);
 192       }
 193     }
 194
 195
 196   }
 197   else {
 198     hitArray     = input.GetHits();
 199     primaryArray = input.GetPrimaries();
 200   }
 201   fCorrectionArray.SetName("FMD_bg_correction");
 202   fCorrectionArray.SetOwner();
 203
 204   TList* primaryList     = new TList();
 205   primaryList->SetName("primaries");
 206
 207   TList* hitList     = new TList();
 208   hitList->SetName("hits");
 209   TList* corrList    = new TList();
 210   corrList->SetName("corrections");
 211
 212   AliFMDAnaCalibBackgroundCorrection* background = new AliFMDAnaCalibBackgroundCorrection();
 213
 214   for(Int_t det= 1; det <=3; det++) {
 215     Int_t nRings = (det==1 ? 1 : 2);
 216
 217     TObjArray* detArrayCorrection = new TObjArray();
 218     detArrayCorrection->SetName(Form("FMD%d",det));
 219     fCorrectionArray.AddAtAndExpand(detArrayCorrection,det);
 220
 221
 222     for(Int_t iring = 0; iring<nRings; iring++) {
 223       TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
 224       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 225       TObjArray* vtxArrayCorrection = new TObjArray();
 226       vtxArrayCorrection->SetName(Form("FMD%d%c",det,ringChar));
 227       detArrayCorrection->AddAtAndExpand(vtxArrayCorrection,iring);
 228
 229       for(Int_t vertexBin=0;vertexBin<nvtxbins;vertexBin++) {
 230         TObjArray* detArray  = (TObjArray*)hitArray->At(det);
 231         TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
 232         TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
 233         hitList->Add(hHits);
 234         TH2F* hPrimary  = (TH2F*)primRingArray->At(vertexBin);
 235         primaryList->Add(hPrimary);
 236         TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
 237         hCorrection->Divide(hPrimary);
 238         corrList->Add(hCorrection);
 239         hCorrection->SetTitle(hCorrection->GetName());
 240         vtxArrayCorrection->AddAtAndExpand(hCorrection,vertexBin);
 241         background->SetBgCorrection(det,ringChar,vertexBin,hCorrection);
 242       }
 243
 244     }
 245   }
 246
 247   TAxis refAxis(nvtxbins,-1*zvtxcut,zvtxcut);
 248   if(inFile) {
 249     TFile* infile = TFile::Open(infilename);
 250     TAxis* refaxis = (TAxis*)infile->Get("vertexbins");
 251     background->SetRefAxis(refaxis);
 252
 253   }
 254   else background->SetRefAxis(&refAxis);
 255
 256   TFile*  fout = new TFile(filename,"RECREATE");
 257   refAxis.Write("vertexbins");
 258
 259   hitList->Write();
 260   primaryList->Write();
 261   corrList->Write();
 262
 263   TObjArray* container = new TObjArray();
 264   container->SetOwner();
 265   container->AddAtAndExpand(&refAxis,0);
 266   container->AddAtAndExpand(&fCorrectionArray,1);
 267   container->AddAtAndExpand(hitArray,2);
 268   container->AddAtAndExpand(primaryArray,3);
 269
 270
 271   if(storeInAlien) {
 272     AliCDBManager* cdb = AliCDBManager::Instance();
 273     cdb->SetDefaultStorage("local://$ALICE_ROOT");
 274     AliCDBId      id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
 275
 276     AliCDBMetaData* meta = new AliCDBMetaData;
 277     meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
 278     meta->SetAliRootVersion(gROOT->GetVersion());
 279     meta->SetBeamPeriod(1);
 280     meta->SetComment("Background Correction for FMD");
 281     meta->SetProperty("key1", background );
 282     cdb->Put(background, id, meta);
 283
 284   }
 285
 286   fout->Close();
 287
 288
 289   }
 290
 291 //_____________________________________________________________________
 292 void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
 293
 294   AliSimulation sim ;
 295   sim.SetRunNumber(0);
 296   TGrid::Connect("alien:",0,0,"t");
 297   sim.SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
 298   sim.SetConfigFile("Config.C");
 299   sim.SetRunQA("FMD:");
 300   TStopwatch timer;
 301   timer.Start();
 302   sim.RunSimulation(nEvents);
 303   timer.Stop();
 304   timer.Print();
 305
 306 }
 307
 308 //_____________________________________________________________________
 309 Bool_t
 310 AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h,
 311                                                       TParticle* /*p*/)
 312 {
 313
 314   if(!h)
 315     return kTRUE;
 316   Bool_t retval = ProcessEvent(h->Detector(),
 317                                h->Ring(),
 318                                h->Sector(),
 319                                h->Strip(),
 320                                h->Track(),
 321                                h->Q());
 322
 323   return retval;
 324 }
 325 //_____________________________________________________________________
 326 Bool_t
 327 AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr,
 328                                                            TParticle* p)
 329 {
 330   if(!tr)
 331     return kTRUE;
 332   UShort_t det,sec,strip;
 333   Char_t   ring;
 334   AliFMDStripIndex::Unpack(tr->UserId(),det,ring,sec,strip);
 335   Int_t         nTrack  = tr->GetTrack();
 336   TDatabasePDG* pdgDB   = TDatabasePDG::Instance();
 337   TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
 338   Float_t       charge  = (pdgPart ? pdgPart->Charge() : 0);
 339   Bool_t        retval  = ProcessEvent(det,ring,sec,strip,nTrack,charge);
 340   return retval;
 341
 342 }
 343 //_____________________________________________________________________
 344 Bool_t
 345 AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
 346                                                         Char_t   ring,
 347                                                         UShort_t sec,
 348                                                         UShort_t strip,
 349                                                         Int_t    nTrack,
 350                                                         Float_t  charge)
 351 {
 352
 353
 354
 355   if(charge !=  0 &&
 356      ((nTrack != fPrevTrack) ||
 357       (det != fPrevDetector) ||
 358       (ring != fPrevRing)    ||
 359       (sec != fPrevSec))) {
 360     fHitMap.operator()(det,ring,sec,strip) = 1;
 361     fHits++;
 362
 363   }
 364
 365   fPrevDetector = det;
 366   fPrevRing     = ring;
 367   fPrevSec      = sec;
 368   fPrevTrack    = nTrack;
 369
 370   return kTRUE;
 371
 372 }
 373
 374 //_____________________________________________________________________
 375 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init()
 376 {
 377   fPrimaryArray.SetOwner();
 378   fPrimaryArray.SetName("FMD_primary");
 379
 380   fPrimaryMapInner.SetBins(fNbinsEta, -6,6, 20, 0, 2*TMath::Pi());
 381   fPrimaryMapOuter.SetBins(fNbinsEta, -6,6, 40, 0, 2*TMath::Pi());
 382   fPrimaryMapInner.SetName("fPrimaryMapInner");
 383   fPrimaryMapInner.SetName("fPrimaryMapOuter");
 384
 385   fPrimaryMapInner.Sumw2();
 386   fPrimaryMapOuter.Sumw2();
 387   for(Int_t iring = 0; iring<2;iring++) {
 388     Char_t ringChar = (iring == 0 ? 'I' : 'O');
 389     TObjArray* ringArray = new TObjArray();
 390     ringArray->SetName(Form("FMD_%c",ringChar));
 391     fPrimaryArray.AddAtAndExpand(ringArray,iring);
 392     Int_t nSec = (iring == 1 ? 40 : 20);
 393     for(Int_t v=0; v<fNvtxBins;v++) {
 394
 395       TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
 396                                       Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
 397                                       fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
 398       hPrimary->Sumw2();
 399       ringArray->AddAtAndExpand(hPrimary,v);
 400     }
 401   }
 402
 403
 404   fHitArray.SetOwner();
 405   fHitArray.SetName("FMD_hits");
 406
 407   for(Int_t det =1; det<=3;det++) {
 408     TObjArray* detArrayHits = new TObjArray();
 409     detArrayHits->SetName(Form("FMD%d",det));
 410     fHitArray.AddAtAndExpand(detArrayHits,det);
 411     Int_t nRings = (det==1 ? 1 : 2);
 412     for(Int_t ring = 0;ring<nRings;ring++) {
 413       Int_t nSec = (ring == 1 ? 40 : 20);
 414       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 415       TObjArray* vtxArrayHits = new TObjArray();
 416       vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
 417       detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
 418       for(Int_t v=0; v<fNvtxBins;v++) {
 419         TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 420                                Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 421                                fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
 422         hHits->Sumw2();
 423         vtxArrayHits->AddAtAndExpand(hHits,v);
 424
 425       }
 426     }
 427   }
 428
 429   AliFMDInput::Init();
 430
 431   return kTRUE;
 432 }
 433 //_____________________________________________________________________
 434
 435 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event )
 436 {
 437
 438   Bool_t             retVal    = AliFMDInput::Begin(event);
 439   AliStack*          partStack = fLoader->Stack();
 440   Int_t              nTracks   = partStack->GetNtrack();
 441
 442
 443   for(Int_t j=0;j<nTracks;j++) {
 444     TParticle* p           = partStack->Particle(j);
 445     TDatabasePDG* pdgDB    = TDatabasePDG::Instance();
 446     TParticlePDG* pdgPart  = pdgDB->GetParticle(p->GetPdgCode());
 447     Float_t       charge   = (pdgPart ? pdgPart->Charge() : 0);
 448     Float_t       phi      = TMath::ATan2(p->Py(),p->Px());
 449
 450     if(phi<0) phi = phi+2*TMath::Pi();
 451     // if(p->Theta() == 0) continue;
 452     Float_t eta   = p->Eta();
 453
 454     // std::cout<<-1*TMath::Log(TMath::Tan(0.5*p->Theta()))<<std::endl;
 455
 456     Bool_t primary = partStack->IsPhysicalPrimary(j);
 457     //(charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01);
 458     if(primary && charge !=0) {
 459
 460
 461
 462       fPrim++;
 463
 464       fPrimaryMapInner.Fill(eta,phi);
 465       fPrimaryMapOuter.Fill(eta,phi);
 466     }
 467   }
 468
 469   return retVal;
 470 }
 471 //_____________________________________________________________________
 472 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::End()  {
 473
 474   Bool_t retval = AliFMDInput::End();
 475
 476   AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
 477   TArrayF vertex;
 478   genHeader->PrimaryVertex(vertex);
 479
 480   if(TMath::Abs(vertex.At(2)) > fZvtxCut)
 481     return kTRUE;
 482
 483   Double_t delta           = 2*fZvtxCut/fNvtxBins;
 484   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 485   Int_t    vertexBin       = (Int_t)vertexBinDouble;
 486   //Primaries
 487   TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
 488   TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
 489
 490   TH2F* hPrimaryInner  = (TH2F*)innerArray->At(vertexBin);
 491   TH2F* hPrimaryOuter  = (TH2F*)outerArray->At(vertexBin);
 492
 493   hPrimaryInner->Add(&fPrimaryMapInner);
 494   hPrimaryOuter->Add(&fPrimaryMapOuter);
 495
 496   //Hits
 497   for(UShort_t det=1;det<=3;det++) {
 498     Int_t nRings = (det==1 ? 1 : 2);
 499     for (UShort_t ir = 0; ir < nRings; ir++) {
 500       Char_t   ring = (ir == 0 ? 'I' : 'O');
 501       UShort_t nsec = (ir == 0 ? 20  : 40);
 502       UShort_t nstr = (ir == 0 ? 512 : 256);
 503
 504       for(UShort_t sec =0; sec < nsec;  sec++) {
 505
 506         for(UShort_t strip = 0; strip < nstr; strip++) {
 507
 508           if(fHitMap.operator()(det,ring,sec,strip) > 0) {
 509
 510             Double_t x,y,z;
 511             AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
 512             fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
 513
 514             Int_t iring = (ring == 'I' ? 0 : 1);
 515
 516             TObjArray* detArray  = (TObjArray*)fHitArray.At(det);
 517             TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
 518             TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
 519
 520             Float_t   phi   = TMath::ATan2(y,x);
 521             if(phi<0) phi   = phi+2*TMath::Pi();
 522             Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
 523             Float_t   theta = TMath::ATan2(r,z-vertex.At(2));
 524             Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 525             hHits->Fill(eta,phi);
 526           }
 527
 528         }
 529       }
 530     }
 531   }
 532
 533   fPrimaryMapInner.Reset();
 534   fPrimaryMapOuter.Reset();
 535   fHitMap.Reset(0);
 536
 537   return retval;
 538 }
 539
 540 //
 541 // EOF
 542 //
 543
 544