PWGLF/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(const Double_t *x, Double_t *par) {
  55   //A convolution of three landaus to fit three MIP peaks
  56
  57   Double_t energy        = x[0];
  58   Double_t constant = par[0];
  59   Double_t mpv      = par[1];
  60   Double_t sigma    = par[2];
  61   Double_t alpha    = par[3];
  62   Double_t beta     = par[4];
  63
  64   Double_t f = constant*(TMath::Landau(energy,mpv,sigma,kTRUE)+
  65                          alpha*TMath::Landau(energy,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
  66                          beta*TMath::Landau(energy,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE) );
  67
  68   return f;
  69 }
  70 //____________________________________________________________________
  71
  72 AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
  73   : //fDebug(0),
  74   fOutputList(0),
  75   fArray(0),
  76   fZvtxDist(0),
  77   fEvents(0),
  78   fEmptyEvents(0),
  79   fClusters(0),
  80   fClustersEmpty(0),
  81   fFirstEvent(kTRUE),
  82   fParam(0)
  83 {
  84   // Default constructor
  85
  86
  87 }
  88 //____________________________________________________________________
  89 AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
  90     AliAnalysisTaskSE(name),
  91     //fDebug(0),
  92     fOutputList(0),
  93     fArray(0),
  94     fZvtxDist(0),
  95     fEvents(0),
  96     fEmptyEvents(0),
  97     fClusters(0),
  98     fClustersEmpty(0),
  99     fFirstEvent(kTRUE),
 100     fParam(0)
 101 {
 102   // Default constructor
 103
 104   DefineOutput(1, TList::Class());
 105   fParam = AliFMDAnaParameters::Instance();
 106 }
 107 //____________________________________________________________________
 108 void AliFMDAnalysisTaskCollector::UserCreateOutputObjects()
 109 {
 110   // Create the output container
 111   printf("AnalysisTaskFMD::CreateOutPutData() \n");
 112   fFirstEvent = kTRUE;
 113   fOutputList = new TList();//(TList*)GetOutputData(0);
 114   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 115
 116   fArray     = new TObjArray();
 117   fArray->SetName("FMD");
 118   fArray->SetOwner();
 119   TH1F* hEdist = 0;
 120   TH1F* hEmptyEdist = 0;
 121   TH1F* hRingEdist = 0;
 122   for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
 123     TObjArray* etaArray = new TObjArray();
 124     fArray->AddAtAndExpand(etaArray,nEta);
 125     for(Int_t det =1; det<=3;det++)
 126       {
 127         TObjArray* detArray = new TObjArray();
 128         detArray->SetName(Form("FMD%d",det));
 129         etaArray->AddAtAndExpand(detArray,det);
 130         Int_t nRings = (det==1 ? 1 : 2);
 131         for(Int_t ring = 0;ring<nRings;ring++)
 132           {
 133             Char_t ringChar = (ring == 0 ? 'I' : 'O');
 134             hEdist = new TH1F(Form("FMD%d%c_etabin%d",det,ringChar,nEta),Form("FMD%d%c_etabin%d",det,ringChar,nEta),200,0,6);
 135             hEdist->SetXTitle("#Delta E / E_{MIP}");
 136             fOutputList->Add(hEdist);
 137             detArray->AddAtAndExpand(hEdist,ring);
 138
 139           }
 140       }
 141
 142   }
 143
 144   for(Int_t det =1; det<=3;det++)
 145     {
 146       Int_t nRings = (det==1 ? 1 : 2);
 147       for(Int_t ring = 0;ring<nRings;ring++)
 148         {
 149           Char_t ringChar = (ring == 0 ? 'I' : 'O');
 150           hRingEdist = new TH1F(Form("ringFMD%d%c",det,ringChar),Form("ringFMD%d%c",det,ringChar),200,0,6);
 151           hRingEdist->SetXTitle("#Delta E / E_{MIP}");
 152           fOutputList->Add(hRingEdist);
 153           hEmptyEdist = new TH1F(Form("emptyFMD%d%c",det,ringChar),Form("emptyFMD%d%c",det,ringChar),200,0,6);
 154           hEmptyEdist->SetXTitle("#Delta E / E_{MIP}");
 155           fOutputList->Add(hEmptyEdist);
 156
 157         }
 158     }
 159   fZvtxDist  = new TH1F("ZvtxDist","Vertex distribution",100,-30,30);
 160   fZvtxDist->SetXTitle("z vertex");
 161   fOutputList->Add(fZvtxDist);
 162 }
 163
 164 //____________________________________________________________________
 165 void AliFMDAnalysisTaskCollector::UserExec(Option_t */*option*/)
 166 {
 167   //Collect data for fitting
 168   AliESDEvent* esd = (AliESDEvent*)InputEvent();
 169   AliESD* old = esd->GetAliESDOld();
 170   if (old) {
 171     esd->CopyFromOldESD();
 172   }
 173   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 174   pars->SetTriggerStatus(esd);
 175   if (fFirstEvent) {
 176     pars->SetParametersFromESD(esd);
 177     pars->PrintStatus();
 178     fFirstEvent = kFALSE;
 179   }
 180
 181   TString triggers = esd->GetFiredTriggerClasses();
 182   //if(!triggers.IsNull()) return;
 183   //Bool_t trigger = pars->IsEventTriggered(esd);
 184
 185   Bool_t physics = pars->IsEventTriggered(AliFMDAnaParameters::kMB1);
 186   Bool_t empty   = pars->IsEventTriggered(AliFMDAnaParameters::kEMPTY);
 187   //std::cout<<physics<<"   "<<empty<<std::endl;
 188   //if(!trigger)
 189   //  physics = kFALSE;
 190   Double_t vertex[3] ={0,0,0};
 191
 192   Bool_t vtxStatus =   pars->GetVertex(esd,vertex);
 193
 194   if(!vtxStatus)
 195     physics = kFALSE;
 196
 197   if(physics) {
 198     fZvtxDist->Fill(vertex[2]);
 199     if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
 200       physics = kFALSE;
 201   }
 202   //std::cout<<"Bananer "<<vtxStatus<<"    "<<physics<<std::endl;
 203   AliESDFMD* fmd = esd->GetFMDData();
 204   if (!fmd) return;
 205   if(physics)
 206     fEvents++;
 207   else if(empty)
 208     fEmptyEvents++;
 209
 210   if(!physics && !empty)
 211     return;
 212
 213   const AliMultiplicity* spdMult = esd->GetMultiplicity();
 214   if(physics)
 215     fClusters+= (spdMult->GetNumberOfSingleClusters() + spdMult->GetNumberOfTracklets());
 216
 217   if(empty)
 218     fClustersEmpty+= (spdMult->GetNumberOfSingleClusters() + spdMult->GetNumberOfTracklets());
 219
 220   if(empty)
 221     std::cout<<spdMult->GetNumberOfSingleClusters()<<"  "<<spdMult->GetNumberOfTracklets()<<std::endl;
 222
 223   TH1F* edist = 0;
 224   TH1F* emptyEdist = 0;
 225   TH1F* ringEdist = 0;
 226   for(UShort_t det=1;det<=3;det++) {
 227
 228     Int_t nRings = (det==1 ? 1 : 2);
 229     for (UShort_t ir = 0; ir < nRings; ir++) {
 230       Char_t   ring = (ir == 0 ? 'I' : 'O');
 231       emptyEdist = (TH1F*)fOutputList->FindObject(Form("emptyFMD%d%c",det,ring));
 232       ringEdist = (TH1F*)fOutputList->FindObject(Form("ringFMD%d%c",det,ring));
 233       UShort_t nsec = (ir == 0 ? 20  : 40);
 234       UShort_t nstr = (ir == 0 ? 512 : 256);
 235
 236       for(UShort_t sec =0; sec < nsec;  sec++)  {
 237         for(UShort_t strip = 0; strip < nstr; strip++) {
 238
 239
 240           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
 241           if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
 242
 243           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
 244
 245           Int_t nEta  = pars->GetEtaBin(eta);
 246
 247           //      std::cout<<det<<"  "<<ring<<"   "<<sec<<"    "<<strip<<"   "<<vertex[2]<<"   "<<eta<<"   "<<nEta<<std::endl;
 248           if(physics) {
 249             edist = (TH1F*)fOutputList->FindObject(Form("FMD%d%c_etabin%d",det,ring,nEta));
 250             edist->Fill(mult);
 251             ringEdist->Fill(mult);
 252           }
 253           else if(empty) {
 254             emptyEdist->Fill(mult);
 255
 256           }
 257           else {
 258             AliWarning("Something is wrong - wrong trigger");
 259             continue;
 260           }
 261
 262
 263         }
 264       }
 265     }
 266
 267     PostData(1, fOutputList);
 268
 269   }
 270 }
 271 //____________________________________________________________________
 272
 273 void AliFMDAnalysisTaskCollector::Terminate(Option_t */*option*/) {
 274   //Terminate
 275
 276   std::cout<<"Analysed "<<fEvents<<" events and "<<fEmptyEvents<<" empty"<<std::endl;
 277   std::cout<<fClusters / fEvents << " clusters per physics event and "<<fClustersEmpty / fEmptyEvents<< " clusters per empty event"<<std::endl;
 278
 279   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 280
 281
 282     for(Int_t det = 1; det<=3; det++) {
 283       Int_t nRings  =  (det == 1 ? 1 : 2);
 284       for(Int_t ring = 0;ring<nRings; ring++) {
 285
 286         Char_t ringChar = (ring == 0 ? 'I' : 'O');
 287         TH1F* hRingEdist = (TH1F*)fOutputList->FindObject(Form("ringFMD%d%c",det,ringChar));
 288         TH1F* hEmptyEdist = (TH1F*)fOutputList->FindObject(Form("emptyFMD%d%c",det,ringChar));
 289         if(fEmptyEvents)
 290           hEmptyEdist->Scale(1./(Float_t)fEmptyEvents);
 291
 292         if(fEvents)
 293           hRingEdist->Scale(1./(Float_t)fEvents);
 294         for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
 295           TH1F* hEdist = (TH1F*)fOutputList->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
 296           if(fEvents)
 297             hEdist->Scale(1./(Float_t)fEvents);
 298
 299
 300       }
 301     }
 302   }
 303
 304 }
 305 //____________________________________________________________________
 306 void AliFMDAnalysisTaskCollector::ReadFromFile(const Char_t* filename, Bool_t store, Int_t speciesOption)  {
 307
 308   //speciesOption:
 309   //0: p+p Landau fit
 310   //1: Pb+Pb triple landau convolution fit
 311
 312   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 313   //pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
 314
 315   TFile fin(filename,"UPDATE");
 316
 317   TList* list = (TList*)fin.Get("energyDist");
 318
 319   AliFMDAnaCalibEnergyDistribution* energyDist = new AliFMDAnaCalibEnergyDistribution();
 320
 321   energyDist->SetNetaBins(pars->GetNetaBins());
 322   energyDist->SetEtaLimits(pars->GetEtaMin(),pars->GetEtaMax());
 323
 324   TF1* fitFunc = 0;
 325
 326   for(Int_t det = 1; det<=3; det++) {
 327     Int_t nRings  =  (det == 1 ? 1 : 2);
 328     for(Int_t ring = 0;ring<nRings; ring++) {
 329       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 330       TH1F* hEmptyEdist = (TH1F*)list->FindObject(Form("emptyFMD%d%c",det,ringChar));
 331       TH1F* hRingEdist = (TH1F*)list->FindObject(Form("ringFMD%d%c",det,ringChar));
 332       fitFunc = FitEnergyDistribution(hEmptyEdist, speciesOption) ;
 333       if(fitFunc)
 334         fitFunc->Write(Form("emptyFMD%d%c_fitfunc",det,ringChar),TObject::kWriteDelete);
 335       fitFunc = FitEnergyDistribution(hRingEdist, speciesOption) ;
 336       if(fitFunc)
 337         fitFunc->Write(Form("FMD%d%c_fitfunc",det,ringChar),TObject::kWriteDelete);
 338
 339
 340       energyDist->SetEmptyEnergyDistribution(det,ringChar,hEmptyEdist);
 341       energyDist->SetRingEnergyDistribution(det,ringChar,hRingEdist);
 342       for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
 343         TH1F* hEdist = (TH1F*)list->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
 344
 345         fitFunc = FitEnergyDistribution(hEdist, speciesOption) ;
 346         if(fitFunc)
 347           fitFunc->Write(Form("FMD%d%c_etabin%d_fitfunc",det,ringChar,nEta),TObject::kWriteDelete);
 348         energyDist->SetEnergyDistribution(det,ringChar,nEta,hEdist);
 349
 350       }
 351
 352     }
 353
 354   }
 355
 356   fin.Write();
 357   fin.Close();
 358
 359   if(store) {
 360     TFile fcalib(pars->GetPath(pars->GetEdistID() ),"RECREATE");
 361     energyDist->Write(AliFMDAnaParameters::GetEdistID());
 362     fcalib.Close();
 363   }
 364
 365
 366
 367 }
 368 //____________________________________________________________________
 369 TF1* AliFMDAnalysisTaskCollector::FitEnergyDistribution(TH1F* hEnergy, Int_t /*speciesOption*/) {
 370   //Fit energy distribution
 371   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 372   TF1* fitFunc = 0;
 373   if(hEnergy->GetEntries() != 0) {
 374
 375     hEnergy->GetXaxis()->SetRangeUser(0.2,hEnergy->GetXaxis()->GetXmax());
 376
 377     if(pars->GetCollisionSystem() == 0)
 378       fitFunc =  new TF1("FMDfitFunc","landau",hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,3);
 379     if(pars->GetCollisionSystem() == 1) {
 380       fitFunc = new TF1("FMDfitFunc",TripleLandau,hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,5,5);
 381       fitFunc->SetParNames("constant","MPV","sigma","2-Mip weight","3-Mip weight");
 382       fitFunc->SetParameters(10,0.8,0.1,0.05,0.01);
 383       fitFunc->SetParLimits(1,0.6,1.2);
 384       fitFunc->SetParLimits(3,0,1);
 385       fitFunc->SetParLimits(4,0,1);
 386
 387     }
 388     hEnergy->Fit(fitFunc,"","",hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,3);
 389
 390   }
 391   return fitFunc;
 392
 393 }
 394 //____________________________________________________________________
 395 //
 396 // EOF
 397 //
 398 // EOF