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