PWG2/FORWARD/analysis/AliFMDAnalysisTaskCollector.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id:$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Here some comments on what this task does                                //
  21
  22 #include <TChain.h>
  23 #include <TF1.h>
  24 #include <TFile.h>
  25 #include <TInterpreter.h>
  26 #include <TList.h>
  27 #include <TMath.h>
  28 #include <TROOT.h>
  29 #include <TSystem.h>
  30 #include <iostream>
  31
  32 #include "AliFMDAnalysisTaskCollector.h"
  33 #include "AliAnalysisManager.h"
  34 #include "AliESDFMD.h"
  35 #include "AliESDEvent.h"
  36 #include "AliAODEvent.h"
  37 #include "AliAODHandler.h"
  38 #include "AliMCEventHandler.h"
  39 #include "AliStack.h"
  40 #include "AliMultiplicity.h"
  41 #include "AliESDVertex.h"
  42 #include "AliFMDAnaParameters.h"
  43 #include "AliFMDAnaCalibBackgroundCorrection.h"
  44 #include "AliFMDAnaCalibEnergyDistribution.h"
  45 #include "AliFMDAnaCalibEventSelectionEfficiency.h"
  46 #include "AliFMDAnaCalibSharingEfficiency.h"
  47
  48 //#include "AliFMDGeometry.h"
  49
  50
  51 ClassImp(AliFMDAnalysisTaskCollector)
  52
  53 //____________________________________________________________________
  54 Double_t  AliFMDAnalysisTaskCollector::TripleLandau(Double_t *x, Double_t *par) {
  55
  56   Double_t energy        = x[0];
  57   Double_t constant = par[0];
  58   Double_t mpv      = par[1];
  59   Double_t sigma    = par[2];
  60   Double_t alpha    = par[3];
  61   Double_t beta     = par[4];
  62
  63   Double_t f = constant*(TMath::Landau(energy,mpv,sigma,kTRUE)+
  64                          alpha*TMath::Landau(energy,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
  65                          beta*TMath::Landau(energy,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE) );
  66
  67   return f;
  68 }
  69 //____________________________________________________________________
  70
  71 AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
  72 : fDebug(0),
  73   fOutputList(0),
  74   fArray(0),
  75   fZvtxDist(0),
  76   fEvents(0),
  77   fEmptyEvents(0),
  78   fClusters(0),
  79   fClustersEmpty(0),
  80   fFirstEvent(kTRUE),
  81   fParam(0)
  82 {
  83   // Default constructor
  84
  85
  86 }
  87 //____________________________________________________________________
  88 AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
  89     AliAnalysisTaskSE(name),
  90     fDebug(0),
  91     fOutputList(0),
  92     fArray(0),
  93     fZvtxDist(0),
  94     fEvents(0),
  95     fEmptyEvents(0),
  96     fClusters(0),
  97     fClustersEmpty(0),
  98     fFirstEvent(kTRUE),
  99     fParam(0)
 100 {
 101   // Default constructor
 102
 103   DefineOutput(1, TList::Class());
 104   fParam = AliFMDAnaParameters::Instance();
 105 }
 106 //____________________________________________________________________
 107 void AliFMDAnalysisTaskCollector::UserCreateOutputObjects()
 108 {
 109   // Create the output container
 110   printf("AnalysisTaskFMD::CreateOutPutData() \n");
 111   fFirstEvent = kTRUE;
 112   fOutputList = new TList();//(TList*)GetOutputData(0);
 113   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 114
 115   fArray     = new TObjArray();
 116   fArray->SetName("FMD");
 117   fArray->SetOwner();
 118   TH1F* hEdist = 0;
 119   TH1F* hEmptyEdist = 0;
 120   TH1F* hRingEdist = 0;
 121   for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
 122     TObjArray* etaArray = new TObjArray();
 123     fArray->AddAtAndExpand(etaArray,nEta);
 124     for(Int_t det =1; det<=3;det++)
 125       {
 126         TObjArray* detArray = new TObjArray();
 127         detArray->SetName(Form("FMD%d",det));
 128         etaArray->AddAtAndExpand(detArray,det);
 129         Int_t nRings = (det==1 ? 1 : 2);
 130         for(Int_t ring = 0;ring<nRings;ring++)
 131           {
 132             Char_t ringChar = (ring == 0 ? 'I' : 'O');
 133             hEdist = new TH1F(Form("FMD%d%c_etabin%d",det,ringChar,nEta),Form("FMD%d%c_etabin%d",det,ringChar,nEta),200,0,6);
 134             hEdist->SetXTitle("#Delta E / E_{MIP}");
 135             fOutputList->Add(hEdist);
 136             detArray->AddAtAndExpand(hEdist,ring);
 137
 138           }
 139       }
 140
 141   }
 142
 143   for(Int_t det =1; det<=3;det++)
 144     {
 145       Int_t nRings = (det==1 ? 1 : 2);
 146       for(Int_t ring = 0;ring<nRings;ring++)
 147         {
 148           Char_t ringChar = (ring == 0 ? 'I' : 'O');
 149           hRingEdist = new TH1F(Form("ringFMD%d%c",det,ringChar),Form("ringFMD%d%c",det,ringChar),200,0,6);
 150           hRingEdist->SetXTitle("#Delta E / E_{MIP}");
 151           fOutputList->Add(hRingEdist);
 152           hEmptyEdist = new TH1F(Form("emptyFMD%d%c",det,ringChar),Form("emptyFMD%d%c",det,ringChar),200,0,6);
 153           hEmptyEdist->SetXTitle("#Delta E / E_{MIP}");
 154           fOutputList->Add(hEmptyEdist);
 155
 156         }
 157     }
 158   fZvtxDist  = new TH1F("ZvtxDist","Vertex distribution",100,-30,30);
 159   fZvtxDist->SetXTitle("z vertex");
 160   fOutputList->Add(fZvtxDist);
 161 }
 162
 163 //____________________________________________________________________
 164 void AliFMDAnalysisTaskCollector::UserExec(Option_t */*option*/)
 165 {
 166
 167   AliESDEvent* esd = (AliESDEvent*)InputEvent();
 168   AliESD* old = esd->GetAliESDOld();
 169   if (old) {
 170     esd->CopyFromOldESD();
 171   }
 172   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 173   pars->SetTriggerStatus(esd);
 174   if (fFirstEvent) {
 175     pars->SetParametersFromESD(esd);
 176     pars->PrintStatus();
 177     fFirstEvent = kFALSE;
 178   }
 179
 180   TString triggers = esd->GetFiredTriggerClasses();
 181   //if(!triggers.IsNull()) return;
 182   //Bool_t trigger = pars->IsEventTriggered(esd);
 183
 184   Bool_t physics = pars->IsEventTriggered(AliFMDAnaParameters::kMB1);
 185   Bool_t empty   = pars->IsEventTriggered(AliFMDAnaParameters::kEMPTY);
 186   //std::cout<<physics<<"   "<<empty<<std::endl;
 187   //if(!trigger)
 188   //  physics = kFALSE;
 189   Double_t vertex[3] ={0,0,0};
 190
 191   Bool_t vtxStatus =   pars->GetVertex(esd,vertex);
 192
 193   if(!vtxStatus)
 194     physics = kFALSE;
 195
 196   if(physics) {
 197     fZvtxDist->Fill(vertex[2]);
 198     if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
 199       physics = kFALSE;
 200   }
 201   //std::cout<<"Bananer "<<vtxStatus<<"    "<<physics<<std::endl;
 202   AliESDFMD* fmd = esd->GetFMDData();
 203   if (!fmd) return;
 204   if(physics)
 205     fEvents++;
 206   else if(empty)
 207     fEmptyEvents++;
 208
 209   if(!physics && !empty)
 210     return;
 211
 212   const AliMultiplicity* spdMult = esd->GetMultiplicity();
 213   if(physics)
 214     fClusters+= (spdMult->GetNumberOfSingleClusters() + spdMult->GetNumberOfTracklets());
 215
 216   if(empty)
 217     fClustersEmpty+= (spdMult->GetNumberOfSingleClusters() + spdMult->GetNumberOfTracklets());
 218
 219   if(empty)
 220     std::cout<<spdMult->GetNumberOfSingleClusters()<<"  "<<spdMult->GetNumberOfTracklets()<<std::endl;
 221
 222   TH1F* Edist = 0;
 223   TH1F* emptyEdist = 0;
 224   TH1F* ringEdist = 0;
 225   for(UShort_t det=1;det<=3;det++) {
 226
 227     Int_t nRings = (det==1 ? 1 : 2);
 228     for (UShort_t ir = 0; ir < nRings; ir++) {
 229       Char_t   ring = (ir == 0 ? 'I' : 'O');
 230       emptyEdist = (TH1F*)fOutputList->FindObject(Form("emptyFMD%d%c",det,ring));
 231       ringEdist = (TH1F*)fOutputList->FindObject(Form("ringFMD%d%c",det,ring));
 232       UShort_t nsec = (ir == 0 ? 20  : 40);
 233       UShort_t nstr = (ir == 0 ? 512 : 256);
 234
 235       for(UShort_t sec =0; sec < nsec;  sec++)  {
 236         for(UShort_t strip = 0; strip < nstr; strip++) {
 237
 238
 239           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
 240           if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
 241
 242           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
 243
 244           Int_t nEta  = pars->GetEtaBin(eta);
 245
 246           //      std::cout<<det<<"  "<<ring<<"   "<<sec<<"    "<<strip<<"   "<<vertex[2]<<"   "<<eta<<"   "<<nEta<<std::endl;
 247           if(physics) {
 248             Edist = (TH1F*)fOutputList->FindObject(Form("FMD%d%c_etabin%d",det,ring,nEta));
 249             Edist->Fill(mult);
 250             ringEdist->Fill(mult);
 251           }
 252           else if(empty) {
 253             emptyEdist->Fill(mult);
 254
 255           }
 256           else {
 257             AliWarning("Something is wrong - wrong trigger");
 258             continue;
 259           }
 260
 261
 262         }
 263       }
 264     }
 265
 266     PostData(1, fOutputList);
 267
 268   }
 269 }
 270 //____________________________________________________________________
 271
 272 void AliFMDAnalysisTaskCollector::Terminate(Option_t */*option*/) {
 273   std::cout<<"Analysed "<<fEvents<<" events and "<<fEmptyEvents<<" empty"<<std::endl;
 274   std::cout<<fClusters / fEvents << " clusters per physics event and "<<fClustersEmpty / fEmptyEvents<< " clusters per empty event"<<std::endl;
 275
 276   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 277
 278
 279     for(Int_t det = 1; det<=3; det++) {
 280       Int_t nRings  =  (det == 1 ? 1 : 2);
 281       for(Int_t ring = 0;ring<nRings; ring++) {
 282
 283         Char_t ringChar = (ring == 0 ? 'I' : 'O');
 284         TH1F* hRingEdist = (TH1F*)fOutputList->FindObject(Form("ringFMD%d%c",det,ringChar));
 285         TH1F* hEmptyEdist = (TH1F*)fOutputList->FindObject(Form("emptyFMD%d%c",det,ringChar));
 286         if(fEmptyEvents)
 287           hEmptyEdist->Scale(1./(Float_t)fEmptyEvents);
 288
 289         if(fEvents)
 290           hRingEdist->Scale(1./(Float_t)fEvents);
 291         for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
 292           TH1F* hEdist = (TH1F*)fOutputList->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
 293           if(fEvents)
 294             hEdist->Scale(1./(Float_t)fEvents);
 295
 296
 297       }
 298     }
 299   }
 300
 301 }
 302 //____________________________________________________________________
 303 void AliFMDAnalysisTaskCollector::ReadFromFile(const Char_t* filename, Bool_t store, Int_t speciesOption)  {
 304
 305   //speciesOption:
 306   //0: p+p Landau fit
 307   //1: Pb+Pb triple landau convolution fit
 308
 309   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 310   //pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
 311
 312   TFile fin(filename,"UPDATE");
 313
 314   TList* list = (TList*)fin.Get("energyDist");
 315
 316   AliFMDAnaCalibEnergyDistribution* EnergyDist = new AliFMDAnaCalibEnergyDistribution();
 317
 318   EnergyDist->SetNetaBins(pars->GetNetaBins());
 319   EnergyDist->SetEtaLimits(pars->GetEtaMin(),pars->GetEtaMax());
 320
 321   TF1* fitFunc = 0;
 322
 323   for(Int_t det = 1; det<=3; det++) {
 324     Int_t nRings  =  (det == 1 ? 1 : 2);
 325     for(Int_t ring = 0;ring<nRings; ring++) {
 326       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 327       TH1F* hEmptyEdist = (TH1F*)list->FindObject(Form("emptyFMD%d%c",det,ringChar));
 328       TH1F* hRingEdist = (TH1F*)list->FindObject(Form("ringFMD%d%c",det,ringChar));
 329       fitFunc = FitEnergyDistribution(hEmptyEdist, speciesOption) ;
 330       if(fitFunc)
 331         fitFunc->Write(Form("emptyFMD%d%c_fitfunc",det,ringChar),TObject::kWriteDelete);
 332       fitFunc = FitEnergyDistribution(hRingEdist, speciesOption) ;
 333       if(fitFunc)
 334         fitFunc->Write(Form("FMD%d%c_fitfunc",det,ringChar),TObject::kWriteDelete);
 335
 336
 337       EnergyDist->SetEmptyEnergyDistribution(det,ringChar,hEmptyEdist);
 338       EnergyDist->SetRingEnergyDistribution(det,ringChar,hRingEdist);
 339       for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
 340         TH1F* hEdist = (TH1F*)list->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
 341
 342         fitFunc = FitEnergyDistribution(hEdist, speciesOption) ;
 343         if(fitFunc)
 344           fitFunc->Write(Form("FMD%d%c_etabin%d_fitfunc",det,ringChar,nEta),TObject::kWriteDelete);
 345         EnergyDist->SetEnergyDistribution(det,ringChar,nEta,hEdist);
 346
 347       }
 348
 349     }
 350
 351   }
 352
 353   fin.Write();
 354   fin.Close();
 355
 356   if(store) {
 357     TFile fcalib(pars->GetPath(pars->GetEdistID() ),"RECREATE");
 358     EnergyDist->Write(AliFMDAnaParameters::GetEdistID());
 359     fcalib.Close();
 360   }
 361
 362
 363
 364 }
 365 //____________________________________________________________________
 366 TF1* AliFMDAnalysisTaskCollector::FitEnergyDistribution(TH1F* hEnergy, Int_t /*speciesOption*/) {
 367   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 368   TF1* fitFunc = 0;
 369   if(hEnergy->GetEntries() != 0) {
 370
 371     hEnergy->GetXaxis()->SetRangeUser(0.2,hEnergy->GetXaxis()->GetXmax());
 372
 373     if(pars->GetCollisionSystem() == 0)
 374       fitFunc =  new TF1("FMDfitFunc","landau",hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,3);
 375     if(pars->GetCollisionSystem() == 1) {
 376       fitFunc = new TF1("FMDfitFunc",TripleLandau,hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,5,5);
 377       fitFunc->SetParNames("constant","MPV","sigma","2-Mip weight","3-Mip weight");
 378       fitFunc->SetParameters(10,0.8,0.1,0.05,0.01);
 379       fitFunc->SetParLimits(1,0.6,1.2);
 380       fitFunc->SetParLimits(3,0,1);
 381       fitFunc->SetParLimits(4,0,1);
 382
 383     }
 384     hEnergy->Fit(fitFunc,"","",hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,3);
 385
 386   }
 387   return fitFunc;
 388
 389 }
 390 //____________________________________________________________________
 391 //
 392 // EOF
 393 //
 394 // EOF