PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateCorrection.cxx

   1 #include "AliFMDAnalysisTaskGenerateCorrection.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 "AliMultiplicity.h"
  17 #include "AliHeader.h"
  18 #include "AliFMDAnaCalibBackgroundCorrection.h"
  19 #include "AliFMDAnaCalibEventSelectionEfficiency.h"
  20 #include "AliGenPythiaEventHeader.h"
  21 //#include "AliCDBManager.h"
  22 //#include "AliCDBId.h"
  23 //#include "AliCDBMetaData.h"
  24 #include "TSystem.h"
  25 #include "TROOT.h"
  26 #include "TAxis.h"
  27 ClassImp(AliFMDAnalysisTaskGenerateCorrection)
  28
  29 //_____________________________________________________________________
  30 AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection():
  31 AliAnalysisTaskSE(),
  32   fListOfHits(),
  33   fListOfPrimaries(),
  34   fListOfCorrection(),
  35   fVertexBins(),
  36   fLastTrackByStrip(0),
  37   fHitsByStrip(0),
  38   fZvtxCut(10),
  39   fNvtxBins(10),
  40   fNbinsEta(200),
  41   fBackground(0),
  42   fEventSelectionEff(0)
  43 {
  44   // Default constructor
  45 }
  46 //_____________________________________________________________________
  47 AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection(const char* name):
  48   AliAnalysisTaskSE(name),
  49   fListOfHits(),
  50   fListOfPrimaries(),
  51   fListOfCorrection(),
  52   fVertexBins(),
  53   fLastTrackByStrip(0),
  54   fHitsByStrip(0),
  55   fZvtxCut(10),
  56   fNvtxBins(10),
  57   fNbinsEta(200),
  58   fBackground(0),
  59   fEventSelectionEff(0)
  60 {
  61
  62   DefineOutput(1, TList::Class());
  63   DefineOutput(2, TList::Class());
  64   DefineOutput(3, TH1F::Class());
  65   DefineOutput(4, TList::Class());
  66 }
  67 //_____________________________________________________________________
  68 void AliFMDAnalysisTaskGenerateCorrection::UserCreateOutputObjects()
  69 {
  70 // Create the output containers
  71 //
  72
  73   std::cout<<"Creating output objects"<<std::endl;
  74   for(Int_t v=0; v<fNvtxBins;v++) {
  75
  76     TH2F* hSPDhits       = new TH2F(Form("hSPDhits_vtx%d",v),
  77                                     Form("hSPDhits_vtx%d",v),
  78                                     fNbinsEta, -4,6, 20, 0,2*TMath::Pi());
  79     hSPDhits->Sumw2();
  80     fListOfHits.Add(hSPDhits);
  81
  82     for(Int_t iring = 0; iring<2;iring++) {
  83       Char_t ringChar = (iring == 0 ? 'I' : 'O');
  84       Int_t nSec = (iring == 1 ? 40 : 20);
  85
  86       TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
  87                                       Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
  88                                       fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
  89       hPrimary->Sumw2();
  90       fListOfPrimaries.Add(hPrimary);
  91
  92
  93     }
  94   }
  95
  96
  97
  98
  99   for(Int_t det =1; det<=3;det++) {
 100     Int_t nRings = (det==1 ? 1 : 2);
 101     for(Int_t ring = 0;ring<nRings;ring++) {
 102       Int_t nSec = (ring == 1 ? 40 : 20);
 103       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 104       TH1F* doubleHits = new TH1F(Form("DoubleHits_FMD%d%c",det,ringChar),
 105                                   Form("DoubleHits_FMD%d%c",det,ringChar),fNbinsEta, -4,6);
 106       TH1F* allHits = new TH1F(Form("allHits_FMD%d%c",det,ringChar),
 107                                Form("allHits_FMD%d%c",det,ringChar), fNbinsEta, -4,6);
 108
 109       doubleHits->Sumw2();
 110       allHits->Sumw2();
 111       fListOfHits.Add(allHits);
 112       fListOfHits.Add(doubleHits);
 113
 114       for(Int_t v=0; v<fNvtxBins;v++) {
 115         TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 116                                Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 117                                fNbinsEta, -4,6, nSec, 0, 2*TMath::Pi());
 118         hHits->Sumw2();
 119         fListOfHits.Add(hHits);
 120
 121       }
 122     }
 123   }
 124
 125   TH1F* hEventsSelected  = new TH1F("EventsSelected","EventsSelected",fNvtxBins,0,fNvtxBins);
 126   TH1F* hEventsAll    = new TH1F("EventsAll","EventsAll",fNvtxBins,0,fNvtxBins);
 127   TH1F* hEventsAllNSD    = new TH1F("EventsAllNSD","EventsAllNSD",fNvtxBins,0,fNvtxBins);
 128   TH1F* hEventsSelectedVtx  = new TH1F("EventsSelectedVtx","EventsSelectedVtx",fNvtxBins,0,fNvtxBins);
 129   TH1F* hEventsSelectedTrigger  = new TH1F("EventsSelectedTrigger","EventsSelectedTrigger",fNvtxBins,0,fNvtxBins);
 130
 131   TH1F* hEventsSelectedNSDVtx  = new TH1F("EventsSelectedNSDVtx","EventsSelectedNSDVtx",fNvtxBins,0,fNvtxBins);
 132   TH1F* hEventsSelectedNSD     = new TH1F("EventsSelectedNSD","EventsSelectedNSD",fNvtxBins,0,fNvtxBins);
 133
 134   TH1F* hXvtx = new TH1F("hXvtx","x vertex distribution",100,-2,2);
 135   TH1F* hYvtx = new TH1F("hYvtx","y vertex distribution",100,-2,2);
 136   TH1F* hZvtx = new TH1F("hZvtx","z vertex distribution",4*fNvtxBins,-4*fZvtxCut,4*fZvtxCut);
 137
 138   fListOfPrimaries.Add(hXvtx);
 139   fListOfPrimaries.Add(hYvtx);
 140   fListOfPrimaries.Add(hZvtx);
 141
 142   hEventsSelected->Sumw2();
 143   hEventsAll->Sumw2();
 144   hEventsAllNSD->Sumw2();
 145
 146   fListOfHits.Add(hEventsSelected);
 147   fListOfHits.Add(hEventsSelectedVtx);
 148   fListOfHits.Add(hEventsSelectedTrigger);
 149   fListOfHits.Add(hEventsSelectedNSDVtx);
 150   fListOfHits.Add(hEventsSelectedNSD);
 151
 152
 153   fListOfPrimaries.Add(hEventsAll);
 154   fListOfPrimaries.Add(hEventsAllNSD);
 155
 156   for(Int_t v=0; v<fNvtxBins;v++) {
 157
 158     for(Int_t iring = 0; iring<2;iring++) {
 159       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 160       Int_t nSec = (iring == 1 ? 40 : 20);
 161       TH2F* hAnalysed       = new TH2F(Form("Analysed_FMD%c_vtx%d",ringChar,v),
 162                                        Form("Analysed_FMD%c_vtx%d",ringChar,v),
 163                                        fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
 164
 165       hAnalysed->Sumw2();
 166       fListOfPrimaries.Add(hAnalysed);
 167
 168       TH2F* hAnalysedNSD       = new TH2F(Form("AnalysedNSD_FMD%c_vtx%d",ringChar,v),
 169                                           Form("AnalysedNSD_FMD%c_vtx%d",ringChar,v),
 170                                           fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
 171
 172       hAnalysedNSD->Sumw2();
 173       fListOfPrimaries.Add(hAnalysedNSD);
 174
 175       TH2F* hInel           = new TH2F(Form("Inel_FMD%c_vtx%d",ringChar,v),
 176                                        Form("Inel_FMD%c_vtx%d",ringChar,v),
 177                                        fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
 178
 179       hInel->Sumw2();
 180       fListOfPrimaries.Add(hInel);
 181
 182       TH2F* hNSD           = new TH2F(Form("NSD_FMD%c_vtx%d",ringChar,v),
 183                                       Form("NSD_FMD%c_vtx%d",ringChar,v),
 184                                       fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
 185
 186       hNSD->Sumw2();
 187       fListOfPrimaries.Add(hNSD);
 188
 189     }
 190   }
 191
 192   fVertexBins.SetName("VertexBins");
 193   fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 194
 195 }
 196 //_____________________________________________________________________
 197 void AliFMDAnalysisTaskGenerateCorrection::Init()
 198 {
 199   fLastTrackByStrip.Reset(-1);
 200
 201
 202 }
 203 //_____________________________________________________________________
 204 void AliFMDAnalysisTaskGenerateCorrection::UserExec(Option_t */*option*/)
 205 {
 206
 207   fLastTrackByStrip.Reset(-1);
 208   fHitsByStrip.Reset(0);
 209   AliMCEvent* mcevent = MCEvent();
 210
 211   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 212
 213
 214   AliESDEvent* esdevent = (AliESDEvent*)InputEvent();
 215
 216   pars->SetTriggerStatus(esdevent);
 217
 218   Double_t esdvertex[3];
 219   Bool_t vtxStatus =  pars->GetVertex(esdevent,esdvertex);
 220
 221   AliMCParticle* particle = 0;
 222   AliStack* stack = mcevent->Stack();
 223
 224   UShort_t det,sec,strip;
 225   Char_t   ring;
 226
 227   Int_t nTracks                = mcevent->GetNumberOfTracks();
 228   AliHeader* header            = mcevent->Header();
 229   AliGenEventHeader* genHeader = header->GenEventHeader();
 230
 231   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 232   if (!pythiaGenHeader) {
 233     std::cout<<" no pythia header! - NSD selection unusable"<<std::endl;
 234     //return;
 235   }
 236   Bool_t nsd = kTRUE;
 237   if(pythiaGenHeader) {
 238     Int_t pythiaType = pythiaGenHeader->ProcessType();
 239
 240     if(pythiaType==92 || pythiaType==93)
 241       nsd = kFALSE;
 242   }
 243   //if(pythiaType==94){
 244   //  std::cout<<"double diffractive"<<std::endl;
 245   //  return;
 246   //}
 247
 248   TArrayF vertex;
 249   genHeader->PrimaryVertex(vertex);
 250
 251   TH1F* hXvtx = (TH1F*)fListOfPrimaries.FindObject("hXvtx");
 252   hXvtx->Fill(vertex.At(0));
 253   TH1F* hYvtx = (TH1F*)fListOfPrimaries.FindObject("hYvtx");
 254   hYvtx->Fill(vertex.At(1));
 255   TH1F* hZvtx = (TH1F*)fListOfPrimaries.FindObject("hZvtx");
 256   hZvtx->Fill(vertex.At(2));
 257
 258
 259
 260
 261   Double_t delta           = 2*fZvtxCut/fNvtxBins;
 262   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 263   Int_t    vertexBin       = (Int_t)vertexBinDouble;
 264
 265   // Vertex determination correction
 266   TH1F* hEventsSelected           = (TH1F*)fListOfHits.FindObject("EventsSelected");
 267   TH1F* hEventsSelectedVtx        = (TH1F*)fListOfHits.FindObject("EventsSelectedVtx");
 268   TH1F* hEventsSelectedTrigger    = (TH1F*)fListOfHits.FindObject("EventsSelectedTrigger");
 269   TH1F* hEventsSelectedNSDVtx     = (TH1F*)fListOfHits.FindObject("EventsSelectedNSDVtx");
 270   TH1F* hEventsSelectedNSD        = (TH1F*)fListOfHits.FindObject("EventsSelectedNSD");
 271   TH1F* hEventsAll                = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
 272   TH1F* hEventsAllNSD             = (TH1F*)fListOfPrimaries.FindObject("EventsAllNSD");
 273
 274   // TH1F* hTriggered      = (TH1F*)fListOfHits.FindObject("Triggered");
 275   //  TH1F* hTriggeredAll   = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
 276
 277   Bool_t vtxFound = kTRUE;
 278   if(!vtxStatus)
 279     vtxFound = kFALSE;
 280
 281   //if(TMath::Abs(vertex.At(2)) > fZvtxCut) {
 282   //  vtxFound = kFALSE;
 283
 284   //}
 285   Bool_t isTriggered    = pars->IsEventTriggered(AliFMDAnaParameters::kMB1);
 286   Bool_t isTriggeredNSD = pars->IsEventTriggered(AliFMDAnaParameters::kNSD);
 287   if(vtxFound && isTriggered) hEventsSelected->Fill(vertexBin);
 288
 289   if(vtxFound) hEventsSelectedVtx->Fill(vertexBin);
 290   if(isTriggered) hEventsSelectedTrigger->Fill(vertexBin);
 291
 292   if(vtxFound && isTriggeredNSD) hEventsSelectedNSDVtx->Fill(vertexBin);
 293   if(isTriggeredNSD) hEventsSelectedNSD->Fill(vertexBin);
 294
 295
 296   hEventsAll->Fill(vertexBin);
 297   if(nsd) hEventsAllNSD->Fill(vertexBin);
 298
 299   //  if(!vtxFound || !isTriggered) return;
 300
 301   if(TMath::Abs(vertex.At(2)) > fZvtxCut) {
 302     return;
 303   }
 304
 305   for(Int_t i = 0 ;i<nTracks;i++) {
 306     particle = (AliMCParticle*) mcevent->GetTrack(i);
 307
 308     if(!particle)
 309       continue;
 310
 311     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
 312
 313
 314       TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
 315       TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin));
 316       hPrimaryInner->Fill(particle->Eta(),particle->Phi());
 317       hPrimaryOuter->Fill(particle->Eta(),particle->Phi());
 318       TH2F* hAnalysedInner = (TH2F*)fListOfPrimaries.FindObject( Form("Analysed_FMD%c_vtx%d",'I',vertexBin));
 319       TH2F* hAnalysedOuter = (TH2F*)fListOfPrimaries.FindObject( Form("Analysed_FMD%c_vtx%d",'O',vertexBin));
 320       TH2F* hAnalysedNSDInner = (TH2F*)fListOfPrimaries.FindObject( Form("AnalysedNSD_FMD%c_vtx%d",'I',vertexBin));
 321       TH2F* hAnalysedNSDOuter = (TH2F*)fListOfPrimaries.FindObject( Form("AnalysedNSD_FMD%c_vtx%d",'O',vertexBin));
 322       TH2F* hInelInner = (TH2F*)fListOfPrimaries.FindObject( Form("Inel_FMD%c_vtx%d",'I',vertexBin));
 323       TH2F* hInelOuter = (TH2F*)fListOfPrimaries.FindObject( Form("Inel_FMD%c_vtx%d",'O',vertexBin));
 324       TH2F* hNSDInner = (TH2F*)fListOfPrimaries.FindObject( Form("NSD_FMD%c_vtx%d",'I',vertexBin));
 325       TH2F* hNSDOuter = (TH2F*)fListOfPrimaries.FindObject( Form("NSD_FMD%c_vtx%d",'O',vertexBin));
 326
 327       //if(vtxFound && isTriggered) {
 328       if(isTriggeredNSD) {
 329         hAnalysedNSDInner->Fill(particle->Eta(),particle->Phi());
 330         hAnalysedNSDOuter->Fill(particle->Eta(),particle->Phi());
 331       }
 332       if(isTriggered) {
 333         hAnalysedInner->Fill(particle->Eta(),particle->Phi());
 334         hAnalysedOuter->Fill(particle->Eta(),particle->Phi());
 335       }
 336       hInelInner->Fill(particle->Eta(),particle->Phi());
 337       hInelOuter->Fill(particle->Eta(),particle->Phi());
 338
 339       if(nsd) {
 340         hNSDInner->Fill(particle->Eta(),particle->Phi());
 341         hNSDOuter->Fill(particle->Eta(),particle->Phi());
 342       }
 343     }
 344
 345     for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
 346
 347       AliTrackReference* ref = particle->GetTrackReference(j);
 348
 349       if(ref->DetectorId() != AliTrackReference::kFMD)
 350         continue;
 351
 352       AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
 353       Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
 354       if(particle->Charge() != 0 && i != thisStripTrack ) {
 355
 356         Float_t phi = pars->GetPhiFromSector(det,ring,sec);
 357         Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));
 358
 359         TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
 360         hHits->Fill(eta,phi);
 361         Float_t nstrips = (ring =='O' ? 256 : 512);
 362         fHitsByStrip(det,ring,sec,strip) +=1;
 363         TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
 364         TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
 365
 366         if(fHitsByStrip(det,ring,sec,strip) == 1)
 367           allHits->Fill(eta);
 368
 369         doubleHits->Fill(eta);
 370
 371         fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
 372         if(strip >0)
 373           fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
 374         if(strip < (nstrips - 1))
 375           fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
 376       }
 377     }
 378
 379   }
 380
 381   //SPD part HHD
 382   TH2F* hSPDMult = (TH2F*)fListOfHits.FindObject(Form("hSPDhits_vtx%d", vertexBin));
 383
 384   const AliMultiplicity* spdmult = esdevent->GetMultiplicity();
 385   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
 386     hSPDMult->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
 387
 388   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
 389     hSPDMult->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),spdmult->GetPhiSingle(j));
 390
 391
 392   PostData(1, &fListOfHits);
 393   PostData(2, &fListOfPrimaries);
 394   PostData(3, &fVertexBins);
 395 }
 396 //_____________________________________________________________________
 397 void AliFMDAnalysisTaskGenerateCorrection::Terminate(Option_t */*option*/)
 398 {
 399   /*  TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
 400   TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
 401
 402   doubleHits->Divide(allHits);
 403   GenerateCorrection();
 404   PostData(1, &fListOfHits);
 405   PostData(4, &fListOfCorrection);*/
 406
 407 }
 408 //_____________________________________________________________________
 409 void AliFMDAnalysisTaskGenerateCorrection::GenerateCorrection() {
 410
 411   fBackground         = new AliFMDAnaCalibBackgroundCorrection();
 412   fEventSelectionEff  = new AliFMDAnaCalibEventSelectionEfficiency();
 413
 414   //Event selection
 415   TH1F* hEventsSelected           = (TH1F*)fListOfHits.FindObject("EventsSelected");
 416   TH1F* hEventsSelectedVtx        = (TH1F*)fListOfHits.FindObject("EventsSelectedVtx");
 417   TH1F* hEventsSelectedTrigger    = (TH1F*)fListOfHits.FindObject("EventsSelectedTrigger");
 418   TH1F* hEventsSelectedNSDVtx     = (TH1F*)fListOfHits.FindObject("EventsSelectedNSDVtx");
 419   TH1F* hEventsSelectedNSD        = (TH1F*)fListOfHits.FindObject("EventsSelectedNSD");
 420
 421   TH1F* hEventsAll                = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
 422   TH1F* hEventsAllNSD             = (TH1F*)fListOfPrimaries.FindObject("EventsAllNSD");
 423   TH1F* hEventsSelectedVtxDivByTr = (TH1F*)hEventsSelectedVtx->Clone("hEventsSelectedVtxDivByTr");
 424   TH1F* hEventsSelectedNSDVtxDivByNSD = (TH1F*)hEventsSelectedNSDVtx->Clone("hEventsSelectedNSDVtxDivByNSD");
 425
 426   fListOfHits.Add(hEventsSelectedVtxDivByTr);
 427   fListOfHits.Add(hEventsSelectedNSDVtxDivByNSD);
 428   hEventsSelectedNSDVtxDivByNSD->Divide(hEventsSelectedNSD);
 429   hEventsSelectedVtxDivByTr->Divide(hEventsSelectedTrigger);
 430
 431   for(Int_t v=0;v<fNvtxBins ; v++) {
 432     TH2F* hAnalysedInner = (TH2F*)fListOfPrimaries.FindObject(Form("Analysed_FMD%c_vtx%d",'I',v));
 433     TH2F* hAnalysedOuter = (TH2F*)fListOfPrimaries.FindObject(Form("Analysed_FMD%c_vtx%d",'O',v));
 434     TH2F* hAnalysedNSDInner = (TH2F*)fListOfPrimaries.FindObject( Form("AnalysedNSD_FMD%c_vtx%d",'I',v));
 435     TH2F* hAnalysedNSDOuter = (TH2F*)fListOfPrimaries.FindObject( Form("AnalysedNSD_FMD%c_vtx%d",'O',v));
 436
 437     TH2F* hInelInner = (TH2F*)fListOfPrimaries.FindObject(Form("Inel_FMD%c_vtx%d",'I',v));
 438     TH2F* hInelOuter = (TH2F*)fListOfPrimaries.FindObject(Form("Inel_FMD%c_vtx%d",'O',v));
 439     TH2F* hNSDInner  = (TH2F*)fListOfPrimaries.FindObject( Form("NSD_FMD%c_vtx%d",'I',v));
 440     TH2F* hNSDOuter  = (TH2F*)fListOfPrimaries.FindObject( Form("NSD_FMD%c_vtx%d",'O',v));
 441     // hAnalysedInner->Scale((hEventsSelected->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelected->GetBinContent(v+1)));
 442     //hAnalysedOuter->Scale((hEventsSelected->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelected->GetBinContent(v+1)));
 443
 444     hAnalysedInner->Scale((hEventsSelectedTrigger->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelectedTrigger->GetBinContent(v+1)));
 445
 446     hAnalysedOuter->Scale((hEventsSelectedTrigger->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelectedTrigger->GetBinContent(v+1)));
 447     hAnalysedNSDInner->Scale((hEventsSelectedNSD->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelectedNSD->GetBinContent(v+1)));
 448
 449     hAnalysedNSDOuter->Scale((hEventsSelectedNSD->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelectedNSD->GetBinContent(v+1)));
 450
 451
 452
 453     hInelInner->Scale((hEventsAll->GetBinContent(v+1) == 0 ? 0 : 1/hEventsAll->GetBinContent(v+1)));
 454     hInelOuter->Scale((hEventsAll->GetBinContent(v+1) == 0 ? 0 : 1/hEventsAll->GetBinContent(v+1)));
 455
 456     hNSDInner->Scale((hEventsAllNSD->GetBinContent(v+1) == 0 ? 0 : 1/hEventsAllNSD->GetBinContent(v+1)));
 457     hNSDOuter->Scale((hEventsAllNSD->GetBinContent(v+1) == 0 ? 0 : 1/hEventsAllNSD->GetBinContent(v+1)));
 458
 459
 460
 461     //hAnalysedInner->Scale(1./0.84);  //numbers for real data, run104892
 462     //hAnalysedOuter->Scale(1./0.84);
 463     hAnalysedInner->Divide(hInelInner);
 464     hAnalysedOuter->Divide(hInelOuter);
 465
 466     hAnalysedNSDInner->Divide(hNSDInner);
 467     hAnalysedNSDOuter->Divide(hNSDOuter);
 468
 469     fEventSelectionEff->SetCorrection("INEL",v,'I',hAnalysedInner);
 470     fEventSelectionEff->SetCorrection("INEL",v,'O',hAnalysedOuter);
 471     //NSD
 472     fEventSelectionEff->SetCorrection("NSD",v,'I',hAnalysedNSDInner);
 473     fEventSelectionEff->SetCorrection("NSD",v,'O',hAnalysedNSDOuter);
 474
 475
 476 }
 477
 478   Float_t vtxEff = 1;
 479   if(hEventsSelectedTrigger->GetEntries())
 480     vtxEff = hEventsSelectedVtx->GetEntries() / hEventsSelectedTrigger->GetEntries();
 481   fEventSelectionEff->SetVtxToTriggerRatio(vtxEff);
 482
 483   //  hEventsAll->Divide(hEventsAll,hEventsSelected,1,1,"B");
 484   hEventsSelectedVtx->Divide(hEventsAll);
 485   hEventsSelectedTrigger->Divide(hEventsAll);
 486
 487   hEventsSelectedNSDVtx->Divide(hEventsAllNSD);
 488   hEventsSelectedNSD->Divide(hEventsAllNSD);
 489
 490   for(Int_t i = 1; i<=hEventsSelected->GetNbinsX(); i++) {
 491     if(hEventsSelected->GetBinContent(i) == 0 )
 492       continue;
 493     Float_t b    = hEventsSelected->GetBinContent(i);
 494     // Float_t b    = hEventsSelected->GetBinContent(i)*hEventsSelectedTrigger->GetBinContent(i);
 495     Float_t db   = hEventsSelected->GetBinError(i);
 496     Float_t sum  = hEventsAll->GetBinContent(i);
 497     Float_t dsum = hEventsAll->GetBinError(i);
 498     Float_t a    = sum-b;
 499     Float_t da   = TMath::Sqrt(TMath::Power(db,2) + TMath::Power(dsum,2));
 500
 501     Float_t cor  = sum / b;
 502     Float_t ecor = TMath::Sqrt(TMath::Power(da,2) + TMath::Power(a/(b*db),2)) / b;
 503
 504     hEventsAll->SetBinContent(i,cor);
 505     hEventsAll->SetBinError(i,ecor);
 506
 507   }
 508
 509   //TH1F* hEventTest = (TH1F*)hEventsAll->Clone("hEventTest");
 510
 511
 512
 513   fEventSelectionEff->SetCorrection(hEventsAll);
 514
 515   for(Int_t det= 1; det <=3; det++) {
 516     Int_t nRings = (det==1 ? 1 : 2);
 517
 518     for(Int_t iring = 0; iring<nRings; iring++) {
 519       Char_t ring = (iring == 0 ? 'I' : 'O');
 520       TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
 521       TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
 522       allHits->Divide(doubleHits);
 523
 524       fBackground->SetDoubleHitCorrection(det,ring,allHits);
 525
 526       for(Int_t vertexBin=0;vertexBin<fNvtxBins  ;vertexBin++) {
 527         TH2F* hHits          = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
 528         TH2F* hPrimary  = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",ring,vertexBin));
 529         TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
 530         hCorrection->Divide(hPrimary);
 531
 532
 533         hCorrection->SetTitle(hCorrection->GetName());
 534         fListOfCorrection.Add(hCorrection);
 535         fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
 536
 537
 538       }
 539
 540     }
 541   }
 542   for(Int_t vertexBin=0;vertexBin<fNvtxBins  ;vertexBin++) {
 543     TH2F* hPrimary  = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
 544     TH2F* hSPDMult = (TH2F*)fListOfHits.FindObject(Form("hSPDhits_vtx%d", vertexBin));
 545     if(!hSPDMult) continue;
 546
 547     TH2F* hCorrection = (TH2F*)hSPDMult->Clone(Form("SPD_vtxbin_%d_correction",vertexBin));
 548     hCorrection->SetTitle(hCorrection->GetName());
 549     fListOfCorrection.Add(hCorrection);
 550     hCorrection->Divide(hPrimary);
 551     fBackground->SetBgCorrection(0,'Q',vertexBin,hCorrection);
 552
 553     TH1F* hAlive = new TH1F(Form("hAliveSPD_vtxbin%d",vertexBin),Form("hAliveSPD_vtxbin%d",vertexBin),hSPDMult->GetNbinsX(),hSPDMult->GetXaxis()->GetXmin(), hSPDMult->GetXaxis()->GetXmax());
 554     TH1F* hPresent = new TH1F(Form("hPresentSPD_vtxbin%d",vertexBin),Form("hPresentSPD_vtxbin%d",vertexBin),hSPDMult->GetNbinsX(),hSPDMult->GetXaxis()->GetXmin(), hSPDMult->GetXaxis()->GetXmax());
 555     for(Int_t xx = 1; xx <=hSPDMult->GetNbinsX(); xx++) {
 556
 557       if(TMath::Abs(hCorrection->GetXaxis()->GetBinCenter(xx)) > 2)
 558         continue;
 559       for(Int_t yy = 1; yy <=hSPDMult->GetNbinsY(); yy++) {
 560         if(hCorrection->GetBinContent(xx,yy) > 0.1)
 561           hAlive->Fill(hCorrection->GetXaxis()->GetBinCenter(xx));
 562         hPresent->Fill(hCorrection->GetXaxis()->GetBinCenter(xx));
 563
 564       }
 565     }
 566     TH1F* hDeadCorrection = (TH1F*)hAlive->Clone(Form("hSPDDeadCorrection_vtxbin%d",vertexBin));
 567     hDeadCorrection->Divide(hPresent);
 568     fBackground->SetSPDDeadCorrection(vertexBin,hDeadCorrection);
 569     fListOfCorrection.Add(hDeadCorrection);
 570   }
 571
 572
 573
 574
 575
 576
 577   TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 578   fBackground->SetRefAxis(&refAxis);
 579
 580 }
 581 //_____________________________________________________________________
 582 void AliFMDAnalysisTaskGenerateCorrection::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /*runNo*/) {
 583
 584   TFile infile(filename);
 585   TH1F* hVertex = (TH1F*)infile.Get("VertexBins");
 586   fZvtxCut = hVertex->GetXaxis()->GetXmax();
 587   fNvtxBins = hVertex->GetXaxis()->GetNbins();
 588   fVertexBins.SetName("VertexBins");
 589   fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 590
 591   TList* listOfHits = (TList*)infile.Get("Hits");
 592   TList* listOfPrim = (TList*)infile.Get("Primaries");
 593
 594   TH1F* hEventsSelected           = (TH1F*)listOfHits->FindObject("EventsSelected");
 595   TH1F* hEventsSelectedVtx        = (TH1F*)listOfHits->FindObject("EventsSelectedVtx");
 596   TH1F* hEventsSelectedTrigger    = (TH1F*)listOfHits->FindObject("EventsSelectedTrigger");
 597   TH1F* hEventsSelectedNSDVtx     = (TH1F*)listOfHits->FindObject("EventsSelectedNSDVtx");
 598   TH1F* hEventsSelectedNSD        = (TH1F*)listOfHits->FindObject("EventsSelectedNSD");
 599   TH1F* hEventsAll                = (TH1F*)listOfPrim->FindObject("EventsAll");
 600   TH1F* hEventsAllNSD             = (TH1F*)listOfPrim->FindObject("EventsAllNSD");
 601
 602   fListOfHits.Add(hEventsSelected);
 603   fListOfHits.Add(hEventsSelectedVtx);
 604   fListOfHits.Add(hEventsSelectedNSD);
 605   fListOfHits.Add(hEventsSelectedNSDVtx);
 606   fListOfHits.Add(hEventsSelectedTrigger);
 607   fListOfPrimaries.Add(hEventsAll);
 608   fListOfPrimaries.Add(hEventsAllNSD);
 609
 610   TH1F* hXvtx = (TH1F*)listOfPrim->FindObject("hXvtx");
 611   TH1F* hYvtx = (TH1F*)listOfPrim->FindObject("hYvtx");
 612   TH1F* hZvtx = (TH1F*)listOfPrim->FindObject("hZvtx");
 613   fListOfPrimaries.Add(hXvtx);
 614   fListOfPrimaries.Add(hYvtx);
 615   fListOfPrimaries.Add(hZvtx);
 616
 617   for(Int_t det =1; det<=3;det++)
 618     {
 619       Int_t nRings = (det==1 ? 1 : 2);
 620       for(Int_t ring = 0;ring<nRings;ring++)
 621         {
 622           Char_t ringChar = (ring == 0 ? 'I' : 'O');
 623           TH1F* allHits = (TH1F*)listOfHits->FindObject(Form("allHits_FMD%d%c",det,ringChar));
 624           TH1F* doubleHits = (TH1F*)listOfHits->FindObject(Form("DoubleHits_FMD%d%c",det,ringChar));
 625           fListOfHits.Add(allHits);
 626           fListOfHits.Add(doubleHits);
 627           for(Int_t v=0; v<fNvtxBins;v++)
 628             {
 629
 630               TH2F* hHits          = (TH2F*)listOfHits->FindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v));
 631               fListOfHits.Add(hHits);
 632             }
 633         }
 634     }
 635   for(Int_t v=0; v<fNvtxBins;v++) {
 636     TH2F* hSPDHits          = (TH2F*)listOfHits->FindObject(Form("hSPDhits_vtx%d", v));
 637     fListOfHits.Add(hSPDHits);
 638
 639     for(Int_t iring = 0; iring<2;iring++) {
 640       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 641
 642
 643       TH2F* hPrimary       = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
 644       TH2F* hAnalysed      = (TH2F*)listOfPrim->FindObject(Form("Analysed_FMD%c_vtx%d",ringChar,v));
 645       TH2F* hAnalysedNSD   = (TH2F*)listOfPrim->FindObject(Form("AnalysedNSD_FMD%c_vtx%d",ringChar,v));
 646       TH2F* hInel          = (TH2F*)listOfPrim->FindObject(Form("Inel_FMD%c_vtx%d",ringChar,v));
 647       TH2F* hNSD           = (TH2F*)listOfPrim->FindObject(Form("NSD_FMD%c_vtx%d",ringChar,v));
 648
 649       fListOfPrimaries.Add(hPrimary);
 650       fListOfPrimaries.Add(hAnalysed);
 651       fListOfPrimaries.Add(hInel);
 652       fListOfPrimaries.Add(hAnalysedNSD);
 653       fListOfPrimaries.Add(hNSD);
 654     }
 655   }
 656   GenerateCorrection();
 657
 658   TFile fout("backgroundFromFile.root","recreate");
 659   fListOfHits.Write();
 660   fListOfPrimaries.Write();
 661   fListOfCorrection.Write();
 662   fVertexBins.Write();
 663
 664   fout.Close();
 665   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 666   if(storeInOCDB) {
 667     TFile fbg(pars->GetPath(pars->GetBackgroundID()),"RECREATE");
 668     fBackground->Write(AliFMDAnaParameters::GetBackgroundID());
 669     fbg.Close();
 670     TFile feselect(pars->GetPath(pars->GetEventSelectionEffID()),"RECREATE");
 671     fEventSelectionEff->Write(AliFMDAnaParameters::GetEventSelectionEffID());
 672     feselect.Close();
 673
 674   }
 675
 676
 677
 678
 679 }
 680 //_____________________________________________________________________
 681 //
 682 // EOF
 683 //