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 "TH2F.h"
  46 #include "AliGenEventHeader.h"
  47 #include "AliHeader.h"
  48 #include "TFile.h"
  49 #include "TAxis.h"
  50 #include "AliCDBId.h"
  51 #include "AliCDBMetaData.h"
  52 #include "TROOT.h"
  53 #include "AliFMDParameters.h"
  54 #include "AliLog.h"
  55 #include "TList.h"
  56 #include "AliFMDAnaParameters.h"
  57 #include "AliFMDAnaCalibBackgroundCorrection.h"
  58 #include "AliTrackReference.h"
  59 #include "AliFMDStripIndex.h"
  60
  61 ClassImp(AliFMDBackgroundCorrection)
  62 //_____________________________________________________________________
  63 AliFMDBackgroundCorrection::AliFMDBackgroundCorrection() :
  64   TNamed(),
  65   fCorrectionArray(),
  66   fPrimaryList()
  67 {}
  68
  69 //_____________________________________________________________________
  70 AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG(Bool_t hits_not_trackref) :
  71   AliFMDInput(),
  72   fPrimaryArray(),
  73   fHitArray(),
  74   fPrim(0),
  75   fHits(0),
  76   fZvtxCut(0),
  77   fNvtxBins(0),
  78   fPrevTrack(-1),
  79   fPrevDetector(-1),
  80   fPrevRing('Q'),
  81   fNbinsEta(100)
  82 {
  83   AddLoad(kTracks);
  84   if(hits_not_trackref)
  85     AddLoad(kHits);
  86   else
  87     AddLoad(kTrackRefs);
  88   AddLoad(kKinematics);
  89   AddLoad(kHeader);
  90 }
  91
  92 //_____________________________________________________________________
  93
  94 void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Bool_t from_hits,
  95                                                               const Int_t nvtxbins,
  96                                                               Float_t zvtxcut,
  97                                                               const Int_t nBinsEta,
  98                                                               Bool_t storeInAlien,
  99                                                               Int_t runNo,
 100                                                               Int_t endRunNo,
 101                                                               const Char_t* filename,
 102                                                               Bool_t simulate,
 103                                                               Int_t nEvents,
 104                                                               Bool_t inFile,
 105                                                               const Char_t* infilename) {
 106
 107   //TGrid::Connect("alien:",0,0,"t");
 108   if(simulate)
 109     Simulate(nEvents);
 110   else {
 111     //AliCDBManager::Instance()->SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
 112     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
 113     AliCDBManager::Instance()->SetRun(runNo);
 114
 115 #if defined(__CINT__)
 116     gSystem->Load("liblhapdf");
 117     gSystem->Load("libEGPythia6");
 118     gSystem->Load("libpythia6");
 119     gSystem->Load("libAliPythia6");
 120     gSystem->Load("libgeant321");
 121 #endif
 122       //
 123   }
 124
 125   //Setting up the geometry
 126   //-----------------------------------------------
 127   if (AliGeomManager::GetGeometry() == NULL)
 128     AliGeomManager::LoadGeometry();
 129
 130   AliFMDGeometry* geo = AliFMDGeometry::Instance();
 131   geo->Init();
 132   geo->InitTransformations();
 133
 134   AliInfo("Processing hits and primaries ");
 135
 136   AliFMDInputBG input(from_hits);
 137
 138   if(!inFile) {
 139
 140     input.SetVtxCutZ(zvtxcut);
 141     input.SetNvtxBins(nvtxbins);
 142     input.SetNbinsEta(nBinsEta);
 143     input.Run();
 144   }
 145
 146   AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
 147   TObjArray* hitArray ;
 148   TObjArray* primaryArray;
 149   if(inFile) {
 150     TFile* infile = TFile::Open(infilename);
 151     hitArray     = new TObjArray();
 152     primaryArray = new TObjArray();
 153
 154     for(Int_t det =1; det<=3;det++)
 155       {
 156         TObjArray* detArrayHits = new TObjArray();
 157         detArrayHits->SetName(Form("FMD%d",det));
 158         hitArray->AddAtAndExpand(detArrayHits,det);
 159         Int_t nRings = (det==1 ? 1 : 2);
 160         for(Int_t ring = 0;ring<nRings;ring++)
 161           {
 162
 163             Char_t ringChar = (ring == 0 ? 'I' : 'O');
 164             TObjArray* vtxArrayHits = new TObjArray();
 165             vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
 166             detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
 167             for(Int_t v=0; v<nvtxbins;v++)
 168               {
 169
 170                 TH2F* hHits          = (TH2F*)infile->Get(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v));
 171
 172
 173               vtxArrayHits->AddAtAndExpand(hHits,v);
 174
 175             }
 176         }
 177
 178     }
 179
 180     for(Int_t iring = 0; iring<2;iring++) {
 181       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 182       TObjArray* ringArray = new TObjArray();
 183       ringArray->SetName(Form("FMD_%c",ringChar));
 184       primaryArray->AddAtAndExpand(ringArray,iring);
 185       for(Int_t v=0; v<nvtxbins;v++) {
 186
 187         TH2F* hPrimary       = (TH2F*)infile->Get(Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
 188         ringArray->AddAtAndExpand(hPrimary,v);
 189       }
 190     }
 191
 192
 193   }
 194   else {
 195     hitArray     = input.GetHits();
 196     primaryArray = input.GetPrimaries();
 197   }
 198   fCorrectionArray.SetName("FMD_bg_correction");
 199   fCorrectionArray.SetOwner();
 200
 201   TList* primaryList     = new TList();
 202   primaryList->SetName("primaries");
 203
 204   TList* hitList     = new TList();
 205   hitList->SetName("hits");
 206   TList* corrList    = new TList();
 207   corrList->SetName("corrections");
 208
 209   AliFMDAnaCalibBackgroundCorrection* background = new AliFMDAnaCalibBackgroundCorrection();
 210
 211   for(Int_t det= 1; det <=3; det++) {
 212     Int_t nRings = (det==1 ? 1 : 2);
 213
 214     TObjArray* detArrayCorrection = new TObjArray();
 215     detArrayCorrection->SetName(Form("FMD%d",det));
 216     fCorrectionArray.AddAtAndExpand(detArrayCorrection,det);
 217
 218
 219     for(Int_t iring = 0; iring<nRings; iring++) {
 220       TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
 221       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 222       TObjArray* vtxArrayCorrection = new TObjArray();
 223       vtxArrayCorrection->SetName(Form("FMD%d%c",det,ringChar));
 224       detArrayCorrection->AddAtAndExpand(vtxArrayCorrection,iring);
 225
 226       for(Int_t vertexBin=0;vertexBin<nvtxbins;vertexBin++) {
 227         TObjArray* detArray  = (TObjArray*)hitArray->At(det);
 228         TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
 229         TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
 230         hitList->Add(hHits);
 231         TH2F* hPrimary  = (TH2F*)primRingArray->At(vertexBin);
 232         primaryList->Add(hPrimary);
 233         TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
 234         hCorrection->Divide(hPrimary);
 235         corrList->Add(hCorrection);
 236         hCorrection->SetTitle(hCorrection->GetName());
 237         vtxArrayCorrection->AddAtAndExpand(hCorrection,vertexBin);
 238         background->SetBgCorrection(det,ringChar,vertexBin,hCorrection);
 239       }
 240
 241     }
 242   }
 243
 244   TAxis refAxis(nvtxbins,-1*zvtxcut,zvtxcut);
 245   if(inFile) {
 246     TFile* infile = TFile::Open(infilename);
 247     TAxis* refaxis = (TAxis*)infile->Get("vertexbins");
 248     background->SetRefAxis(refaxis);
 249
 250   }
 251   else background->SetRefAxis(&refAxis);
 252
 253   TFile*  fout = new TFile(filename,"RECREATE");
 254   refAxis.Write("vertexbins");
 255
 256   hitList->Write();
 257   primaryList->Write();
 258   corrList->Write();
 259
 260   TObjArray* container = new TObjArray();
 261   container->SetOwner();
 262   container->AddAtAndExpand(&refAxis,0);
 263   container->AddAtAndExpand(&fCorrectionArray,1);
 264   container->AddAtAndExpand(hitArray,2);
 265   container->AddAtAndExpand(primaryArray,3);
 266
 267
 268   if(storeInAlien) {
 269     AliCDBManager* cdb = AliCDBManager::Instance();
 270     cdb->SetDefaultStorage("local://$ALICE_ROOT");
 271     AliCDBId      id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
 272
 273     AliCDBMetaData* meta = new AliCDBMetaData;
 274     meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
 275     meta->SetAliRootVersion(gROOT->GetVersion());
 276     meta->SetBeamPeriod(1);
 277     meta->SetComment("Background Correction for FMD");
 278     meta->SetProperty("key1", background );
 279     cdb->Put(background, id, meta);
 280
 281   }
 282
 283   fout->Close();
 284
 285
 286   }
 287
 288 //_____________________________________________________________________
 289 void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
 290
 291   AliSimulation sim ;
 292   sim.SetRunNumber(0);
 293   TGrid::Connect("alien:",0,0,"t");
 294   sim.SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
 295   sim.SetConfigFile("Config.C");
 296   sim.SetRunQA("FMD:");
 297   TStopwatch timer;
 298   timer.Start();
 299   sim.RunSimulation(nEvents);
 300   timer.Stop();
 301   timer.Print();
 302
 303 }
 304
 305 //_____________________________________________________________________
 306 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit( AliFMDHit* h, TParticle* p) {
 307
 308   if(!h)
 309     return kTRUE;
 310   Bool_t retval = ProcessEvent(h->Detector(),h->Ring(),h->Sector(),h->Strip(),h->Track(),h->Q());
 311
 312   return retval;
 313 }
 314 //_____________________________________________________________________
 315 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef( AliTrackReference* tr, TParticle* p) {
 316
 317   if(!tr)
 318     return kTRUE;
 319   UShort_t det,sec,strip;
 320   Char_t   ring;
 321   AliFMDStripIndex::Unpack(tr->UserId(),det,ring,sec,strip);
 322   Int_t    nTrack = tr->GetTrack();
 323   TDatabasePDG* pdgDB = TDatabasePDG::Instance();
 324   TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
 325   Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
 326   Bool_t retval  = ProcessEvent(det,ring,sec,strip,nTrack,charge);
 327   return retval;
 328
 329 }
 330 //_____________________________________________________________________
 331 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
 332                                                                Char_t   ring,
 333                                                                UShort_t sec,
 334                                                                UShort_t strip,
 335                                                                Int_t    nTrack,
 336                                                                Float_t  charge)
 337 {
 338
 339   AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
 340   TArrayF vertex;
 341   genHeader->PrimaryVertex(vertex);
 342
 343   if(TMath::Abs(vertex.At(2)) > fZvtxCut)
 344     return kTRUE;
 345
 346   Double_t delta = 2*fZvtxCut/fNvtxBins;
 347   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 348   Int_t vertexBin = (Int_t)vertexBinDouble;
 349
 350   if(charge !=  0 && (nTrack != fPrevTrack || det != fPrevDetector || ring != fPrevRing)) {
 351
 352     Double_t x,y,z;
 353     AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
 354     fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
 355
 356     Int_t iring = (ring == 'I' ? 0 : 1);
 357
 358     TObjArray* detArray  = (TObjArray*)fHitArray.At(det);
 359     TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
 360     TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
 361
 362     Float_t phi   = TMath::ATan2(y,x);
 363     if(phi<0)
 364       phi = phi+2*TMath::Pi();
 365     Float_t theta = TMath::ATan2(TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2)),z+vertex.At(2));
 366     Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 367     hHits->Fill(eta,phi);
 368     fHits++;
 369     fPrevDetector = det;
 370     fPrevRing     = ring;
 371   }
 372
 373   fPrevTrack = nTrack;
 374   return kTRUE;
 375
 376 }
 377
 378 //_____________________________________________________________________
 379 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init() {
 380
 381   fPrimaryArray.SetOwner();
 382   fPrimaryArray.SetName("FMD_primary");
 383
 384   for(Int_t iring = 0; iring<2;iring++) {
 385     Char_t ringChar = (iring == 0 ? 'I' : 'O');
 386     TObjArray* ringArray = new TObjArray();
 387     ringArray->SetName(Form("FMD_%c",ringChar));
 388     fPrimaryArray.AddAtAndExpand(ringArray,iring);
 389     Int_t nSec = (iring == 1 ? 40 : 20);
 390     for(Int_t v=0; v<fNvtxBins;v++) {
 391
 392       TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
 393                                       Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
 394                                       fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
 395       hPrimary->Sumw2();
 396       ringArray->AddAtAndExpand(hPrimary,v);
 397     }
 398   }
 399
 400
 401   fHitArray.SetOwner();
 402   fHitArray.SetName("FMD_hits");
 403
 404   for(Int_t det =1; det<=3;det++)
 405     {
 406       TObjArray* detArrayHits = new TObjArray();
 407       detArrayHits->SetName(Form("FMD%d",det));
 408       fHitArray.AddAtAndExpand(detArrayHits,det);
 409       Int_t nRings = (det==1 ? 1 : 2);
 410       for(Int_t ring = 0;ring<nRings;ring++)
 411         {
 412           Int_t nSec = (ring == 1 ? 40 : 20);
 413           Char_t ringChar = (ring == 0 ? 'I' : 'O');
 414           TObjArray* vtxArrayHits = new TObjArray();
 415           vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
 416           detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
 417           for(Int_t v=0; v<fNvtxBins;v++)
 418             {
 419
 420               TH2F* hHits          = new TH2F(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v),
 421                                               Form("hHits_FMD%d%c_vtx%d",det,ringChar,v),
 422                                               fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
 423               hHits->Sumw2();
 424               vtxArrayHits->AddAtAndExpand(hHits,v);
 425
 426             }
 427         }
 428
 429     }
 430
 431
 432
 433   AliFMDInput::Init();
 434
 435
 436   return kTRUE;
 437 }
 438 //_____________________________________________________________________
 439
 440 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event ) {
 441
 442   Bool_t retVal = AliFMDInput::Begin(event);
 443
 444   AliStack* partStack=fLoader->Stack();
 445
 446   Int_t nTracks=partStack->GetNtrack();
 447   AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
 448   TArrayF vertex;
 449   genHeader->PrimaryVertex(vertex);
 450
 451   if(TMath::Abs(vertex.At(2)) > fZvtxCut)
 452     return kTRUE;
 453
 454   Double_t delta = 2*fZvtxCut/fNvtxBins;
 455   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 456   Int_t vertexBin = (Int_t)vertexBinDouble;
 457
 458   for(Int_t j=0;j<nTracks;j++)
 459         {
 460           TParticle* p  = partStack->Particle(j);
 461           TDatabasePDG* pdgDB = TDatabasePDG::Instance();
 462           TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
 463           Float_t charge = (pdgPart ? pdgPart->Charge() : 0);
 464           Float_t phi = TMath::ATan2(p->Py(),p->Px());
 465
 466           if(phi<0)
 467             phi = phi+2*TMath::Pi();
 468           if(p->Theta() == 0) continue;
 469               Float_t eta   = -1*TMath::Log(TMath::Tan(0.5*p->Theta()));
 470
 471               Bool_t primary = partStack->IsPhysicalPrimary(j);
 472               //(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);
 473               if(primary && charge !=0) {
 474
 475                 fPrim++;
 476                 TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
 477                 TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
 478
 479                 TH2F* hPrimaryInner  = (TH2F*)innerArray->At(vertexBin);
 480                 TH2F* hPrimaryOuter  = (TH2F*)outerArray->At(vertexBin);
 481                 hPrimaryInner->Fill(eta,phi);
 482                 hPrimaryOuter->Fill(eta,phi);
 483               }
 484         }
 485
 486   return retVal;
 487 }
 488 //_____________________________________________________________________
 489 //
 490 // EOF
 491 //
 492
 493