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