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