FMD/analysis/AliFMDAnalysisTaskGenerateBackground.cxx

   1 #include "AliFMDAnalysisTaskGenerateBackground.h"
   2 #include "AliESDEvent.h"
   3 #include "iostream"
   4 #include "AliESDFMD.h"
   5 #include "TH2F.h"
   6 #include "AliTrackReference.h"
   7 #include "AliStack.h"
   8 #include "AliFMDAnaParameters.h"
   9 #include "AliFMDStripIndex.h"
  10 #include "AliStack.h"
  11 #include "AliMCParticle.h"
  12 #include "AliMCEvent.h"
  13 //#include "AliFMDGeometry.h"
  14 #include "TArray.h"
  15 #include "AliGenEventHeader.h"
  16 #include "AliHeader.h"
  17 #include "AliFMDAnaCalibBackgroundCorrection.h"
  18 //#include "AliCDBManager.h"
  19 //#include "AliCDBId.h"
  20 //#include "AliCDBMetaData.h"
  21 #include "TSystem.h"
  22 #include "TROOT.h"
  23 #include "TAxis.h"
  24 ClassImp(AliFMDAnalysisTaskGenerateBackground)
  25
  26 //_____________________________________________________________________
  27 AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground():
  28 AliAnalysisTaskSE()
  29 {
  30   // Default constructor
  31 }
  32 //_____________________________________________________________________
  33 AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground(const char* name):
  34   AliAnalysisTaskSE(name),
  35   fZvtxCut(10),
  36   fNvtxBins(10),
  37   fNbinsEta(200)
  38 {
  39
  40   DefineOutput(1, TList::Class());
  41   DefineOutput(2, TList::Class());
  42   DefineOutput(3, TH1F::Class());
  43   DefineOutput(4, TList::Class());
  44 }
  45 //_____________________________________________________________________
  46 void AliFMDAnalysisTaskGenerateBackground::UserCreateOutputObjects()
  47 {
  48 // Create the output containers
  49 //
  50
  51   std::cout<<"Creating output objects"<<std::endl;
  52   for(Int_t iring = 0; iring<2;iring++) {
  53     Char_t ringChar = (iring == 0 ? 'I' : 'O');
  54     Int_t nSec = (iring == 1 ? 40 : 20);
  55     for(Int_t v=0; v<fNvtxBins;v++) {
  56
  57       TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
  58                                       Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
  59                                       fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
  60       hPrimary->Sumw2();
  61       fListOfPrimaries.Add(hPrimary);
  62     }
  63   }
  64
  65
  66   for(Int_t det =1; det<=3;det++) {
  67     Int_t nRings = (det==1 ? 1 : 2);
  68     for(Int_t ring = 0;ring<nRings;ring++) {
  69       Int_t nSec = (ring == 1 ? 40 : 20);
  70       Char_t ringChar = (ring == 0 ? 'I' : 'O');
  71       TH1F* doubleHits = new TH1F(Form("DoubleHits_FMD%d%c",det,ringChar),
  72                                   Form("DoubleHits_FMD%d%c",det,ringChar),fNbinsEta, -6,6);
  73       TH1F* allHits = new TH1F(Form("allHits_FMD%d%c",det,ringChar),
  74                                Form("allHits_FMD%d%c",det,ringChar), fNbinsEta, -6,6);
  75
  76       doubleHits->Sumw2();
  77       allHits->Sumw2();
  78       fListOfHits.Add(allHits);
  79       fListOfHits.Add(doubleHits);
  80
  81       for(Int_t v=0; v<fNvtxBins;v++) {
  82         TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
  83                                Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
  84                                fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
  85         hHits->Sumw2();
  86         fListOfHits.Add(hHits);
  87
  88       }
  89     }
  90   }
  91
  92   fVertexBins.SetName("VertexBins");
  93   fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
  94
  95 }
  96 //_____________________________________________________________________
  97 void AliFMDAnalysisTaskGenerateBackground::Init()
  98 {
  99   fLastTrackByStrip.Reset(-1);
 100
 101
 102 }
 103 //_____________________________________________________________________
 104 void AliFMDAnalysisTaskGenerateBackground::UserExec(Option_t */*option*/)
 105 {
 106
 107   fLastTrackByStrip.Reset(-1);
 108   fHitsByStrip.Reset(0);
 109   AliMCEvent* mcevent = MCEvent();
 110   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 111
 112   AliMCParticle* particle = 0;
 113   AliStack* stack = mcevent->Stack();
 114
 115   UShort_t det,sec,strip;
 116   Char_t   ring;
 117
 118   Int_t nTracks = mcevent->GetNumberOfTracks();
 119   AliHeader* header            = mcevent->Header();
 120   AliGenEventHeader* genHeader = header->GenEventHeader();
 121
 122   TArrayF vertex;
 123   genHeader->PrimaryVertex(vertex);
 124
 125   if(TMath::Abs(vertex.At(2)) > fZvtxCut)
 126     return;
 127   for(Int_t i = 0 ;i<nTracks;i++) {
 128     particle = mcevent->GetTrack(i);
 129
 130     if(!particle)
 131       continue;
 132
 133     Double_t delta           = 2*fZvtxCut/fNvtxBins;
 134     Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 135     Int_t    vertexBin       = (Int_t)vertexBinDouble;
 136
 137     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
 138
 139
 140       TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
 141       TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin));
 142       hPrimaryInner->Fill(particle->Eta(),particle->Phi());
 143       hPrimaryOuter->Fill(particle->Eta(),particle->Phi());
 144     }
 145
 146     for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
 147
 148       AliTrackReference* ref = particle->GetTrackReference(j);
 149
 150       if(ref->DetectorId() != AliTrackReference::kFMD)
 151         continue;
 152       AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
 153       Float_t thisStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip);
 154       if(particle->Charge() != 0 && i != thisStripTrack ) {
 155         //Double_t x,y,z;
 156         //AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
 157         //fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
 158         Float_t phi = pars->GetPhiFromSector(det,ring,sec);
 159         Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));
 160         //Float_t   phi   = TMath::ATan2(y,x);
 161         //if(phi<0) phi   = phi+2*TMath::Pi();
 162         //      Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
 163         //Float_t   theta = TMath::ATan2(r,z-vertex.At(2));
 164         //Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 165         TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
 166         hHits->Fill(eta,phi);
 167         Float_t nstrips = (ring =='O' ? 256 : 512);
 168         fHitsByStrip.operator()(det,ring,sec,strip) +=1;
 169         TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
 170         TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
 171
 172         if(fHitsByStrip.operator()(det,ring,sec,strip) == 1)
 173           allHits->Fill(eta);
 174
 175         doubleHits->Fill(eta);
 176         /*if(fHitsByStrip.operator()(det,ring,sec,strip) == 2){
 177           TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
 178           doubleHits->Fill(eta,2);
 179           }*/
 180         //if(fHitsByStrip.operator()(det,ring,sec,strip) > 1){
 181         //  doubleHits->Fill(eta);
 182         //      }
 183
 184
 185         fLastTrackByStrip.operator()(det,ring,sec,strip) = (Float_t)i;
 186         if(strip >0)
 187           fLastTrackByStrip.operator()(det,ring,sec,strip-1) = (Float_t)i;
 188         if(strip < (nstrips - 1))
 189           fLastTrackByStrip.operator()(det,ring,sec,strip+1) = (Float_t)i;
 190       }
 191     }
 192
 193   }
 194
 195
 196   PostData(1, &fListOfHits);
 197   PostData(2, &fListOfPrimaries);
 198   PostData(3, &fVertexBins);
 199 }
 200 //_____________________________________________________________________
 201 void AliFMDAnalysisTaskGenerateBackground::Terminate(Option_t */*option*/)
 202 {
 203   /*  TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
 204   TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
 205
 206   doubleHits->Divide(allHits);
 207   GenerateCorrection();
 208   PostData(1, &fListOfHits);
 209   PostData(4, &fListOfCorrection);*/
 210
 211 }
 212 //_____________________________________________________________________
 213 void AliFMDAnalysisTaskGenerateBackground::GenerateCorrection() {
 214
 215   fBackground = new AliFMDAnaCalibBackgroundCorrection();
 216
 217   for(Int_t det= 1; det <=3; det++) {
 218     Int_t nRings = (det==1 ? 1 : 2);
 219
 220     for(Int_t iring = 0; iring<nRings; iring++) {
 221       Char_t ring = (iring == 0 ? 'I' : 'O');
 222       TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
 223       TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
 224       fBackground->SetDoubleHitCorrection(det,ring,doubleHits);
 225       doubleHits->Divide(allHits);
 226       for(Int_t vertexBin=0;vertexBin<fNvtxBins  ;vertexBin++) {
 227         TH2F* hHits          = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
 228         TH2F* hPrimary  = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",ring,vertexBin));
 229         TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
 230         hCorrection->Divide(hPrimary);
 231         hCorrection->SetTitle(hCorrection->GetName());
 232         fListOfCorrection.Add(hCorrection);
 233         fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
 234
 235
 236       }
 237
 238     }
 239   }
 240   TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 241   fBackground->SetRefAxis(&refAxis);
 242
 243 }
 244 //_____________________________________________________________________
 245 void AliFMDAnalysisTaskGenerateBackground::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t runNo) {
 246
 247   TFile infile(filename);
 248   TH1F* hVertex = (TH1F*)infile.Get("VertexBins");
 249   fZvtxCut = hVertex->GetXaxis()->GetXmax();
 250   fNvtxBins = hVertex->GetXaxis()->GetNbins();
 251   fVertexBins.SetName("VertexBins");
 252   fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 253
 254   TList* listOfHits = (TList*)infile.Get("Hits");
 255   TList* listOfPrim = (TList*)infile.Get("Primaries");
 256
 257   for(Int_t det =1; det<=3;det++)
 258       {
 259         Int_t nRings = (det==1 ? 1 : 2);
 260         for(Int_t ring = 0;ring<nRings;ring++)
 261           {
 262             Char_t ringChar = (ring == 0 ? 'I' : 'O');
 263             TH1F* allHits = (TH1F*)listOfHits->FindObject(Form("allHits_FMD%d%c",det,ringChar));
 264             TH1F* doubleHits = (TH1F*)listOfHits->FindObject(Form("DoubleHits_FMD%d%c",det,ringChar));
 265             fListOfHits.Add(allHits);
 266             fListOfHits.Add(doubleHits);
 267             for(Int_t v=0; v<fNvtxBins;v++)
 268               {
 269
 270                 TH2F* hHits          = (TH2F*)listOfHits->FindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v));
 271                 fListOfHits.Add(hHits);
 272               }
 273           }
 274       }
 275   for(Int_t iring = 0; iring<2;iring++) {
 276     Char_t ringChar = (iring == 0 ? 'I' : 'O');
 277     for(Int_t v=0; v<fNvtxBins;v++) {
 278
 279       TH2F* hPrimary       = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
 280       fListOfPrimaries.Add(hPrimary);
 281
 282     }
 283   }
 284   GenerateCorrection();
 285
 286   TFile fout("backgroundFromFile.root","recreate");
 287   fListOfHits.Write();
 288   fListOfPrimaries.Write();
 289   fListOfCorrection.Write();
 290   fVertexBins.Write();
 291
 292   fout.Close();
 293
 294   if(storeInOCDB) {
 295     TFile fcalib("$ALICE_ROOT/FMD/Correction/Background/background.root","RECREATE");
 296     fBackground->Write(AliFMDAnaParameters::fkBackgroundID);
 297     fcalib.Close();
 298     /*  AliCDBManager* cdb = AliCDBManager::Instance();
 299     cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 300     AliCDBId      id(AliFMDAnaParameters::GetBackgroundPath(),runNo,999999999);
 301
 302     AliCDBMetaData* meta = new AliCDBMetaData;
 303     meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
 304     meta->SetAliRootVersion(gROOT->GetVersion());
 305     meta->SetBeamPeriod(1);
 306     meta->SetComment("Background Correction for FMD");
 307     meta->SetProperty("key1", fBackground );
 308     cdb->Put(fBackground, id, meta);
 309     */
 310   }
 311
 312 }
 313 //_____________________________________________________________________
 314 //
 315 // EOF
 316 //