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