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