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