ANALYSIS/AliAnalysisTaskAODCentralityMaker.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, 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 /////////////////////////////////////////////////////////////
  17 //
  18 // AliAnalysisTaskSE to make AOD centrality
  19 // Author: Alberica Toia, CERN, Alberica.Toia@cern.ch
  20 //
  21 /////////////////////////////////////////////////////////////
  22
  23 #include <TROOT.h>
  24 #include <TSystem.h>
  25
  26 #include "AliVEvent.h"
  27 #include "AliESDEvent.h"
  28 #include "AliAODEvent.h"
  29 #include "AliAODCentrality.h"
  30 #include "AliAnalysisManager.h"
  31 #include "AliESDInputHandler.h"
  32 #include "AliESDZDC.h"
  33 #include "AliESDFMD.h"
  34 #include "AliESDVZERO.h"
  35 #include "AliMultiplicity.h"
  36 #include "AliAODHandler.h"
  37 #include "AliAODHeader.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODVertex.h"
  40 #include "AliMCEvent.h"
  41 #include "AliMCEventHandler.h"
  42 #include "AliMCParticle.h"
  43 #include "AliStack.h"
  44 #include "AliHeader.h"
  45 #include "AliGenEventHeader.h"
  46 #include "AliGenHijingEventHeader.h"
  47 #include "AliAnalysisTaskAODCentralityMaker.h"
  48
  49 ClassImp(AliAnalysisTaskAODCentralityMaker)
  50
  51
  52 //________________________________________________________________________
  53 AliAnalysisTaskAODCentralityMaker::AliAnalysisTaskAODCentralityMaker():
  54 AliAnalysisTaskSE(),
  55 fAODCentrality(0),
  56 fDeltaAODFileName("AliAOD.Centrality.root"),
  57 fAODHeader        (0),
  58 fIsMCInput        (0)
  59 {
  60   // Default constructor
  61 }
  62
  63 //________________________________________________________________________
  64 AliAnalysisTaskAODCentralityMaker::AliAnalysisTaskAODCentralityMaker(const char *name):
  65 AliAnalysisTaskSE(name),
  66 fAODCentrality(0),
  67 fDeltaAODFileName("AliAOD.Centrality.root"),
  68 fAODHeader        (0),
  69 fIsMCInput        (0)
  70 {
  71   // Standard constructor
  72 }
  73
  74
  75 //________________________________________________________________________
  76 AliAnalysisTaskAODCentralityMaker::~AliAnalysisTaskAODCentralityMaker()
  77 {
  78   // Destructor
  79 }
  80
  81 //________________________________________________________________________
  82 void AliAnalysisTaskAODCentralityMaker::Init()
  83 {
  84   // Initialization
  85   if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker::Init() \n");
  86   AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFileName.Data());
  87
  88   return;
  89 }
  90
  91 //________________________________________________________________________
  92 void AliAnalysisTaskAODCentralityMaker::UserCreateOutputObjects()
  93 {
  94
  95   // Create the output container
  96   //
  97   if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker::UserCreateOutPutData() \n");
  98   // Support both the case when the AOD + deltaAOD are produced in an ESD
  99   // analysis or if the deltaAOD is produced on an analysis on AOD's. (A.G. 27/04/09)
 100   if(!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) {
 101     Fatal("UserCreateOutputObjects", "This task needs an AOD handler");
 102     return;
 103   }
 104   TString filename = fDeltaAODFileName;
 105   // When running on standard AOD to produce deltas, IsStandardAOD is never set,
 106   // If AODEvent is NULL, new branches have to be added to the new file(s) (A.G. 15/01/10)
 107   if(!IsStandardAOD() && AODEvent()) filename = "";
 108
 109   fAODCentrality = new AliAODCentrality();
 110   fAODCentrality->SetName("AODCentrality");
 111   AddAODBranch("AliAODCentrality", &fAODCentrality, filename);
 112
 113
 114   fAODHeader = new AliAODHeader();
 115   AddAODBranch("AliAODHeader", &fAODHeader, filename);
 116   return;
 117 }
 118
 119
 120 //________________________________________________________________________
 121 void AliAnalysisTaskAODCentralityMaker::UserExec(Option_t */*option*/)
 122 {
 123   AliVEvent*   event = InputEvent();
 124   AliESDEvent* esd   = dynamic_cast<AliESDEvent*>(event);
 125
 126   Float_t beamEnergy = esd->GetBeamEnergy();
 127   Int_t   nTracks    = event->GetNumberOfTracks();
 128   Int_t   nPmdTracks = esd->GetNumberOfPmdTracks();
 129
 130   // ***** V0 info
 131   AliESDVZERO* esdV0 = esd->GetVZEROData();
 132   Double_t multV0A = esdV0->GetMTotV0A();
 133   Double_t multV0C = esdV0->GetMTotV0C();
 134
 135   // ***** Trigger selection
 136   char     trigClass[100];      //  fired trigger classes
 137   TString triggerClass = esd->GetFiredTriggerClasses();
 138   sprintf(trigClass,"%s",triggerClass.Data());
 139
 140   // ***** vertex info
 141   const AliESDVertex *vertex = esd->GetPrimaryVertexSPD();
 142   Double_t xVertex = vertex->GetX();
 143   Double_t yVertex = vertex->GetY();
 144   Double_t zVertex = vertex->GetZ();
 145   Bool_t vertexer3d;
 146
 147   if(vertex->IsFromVertexer3D()) vertexer3d = kTRUE;
 148   else vertexer3d = kFALSE;
 149   Double_t vertex3[3];
 150   vertex->GetXYZ(vertex3);
 151
 152   // ***** CB info (tracklets, clusters, chips)
 153   const AliMultiplicity *mult = esd->GetMultiplicity();
 154   Int_t nTracklets = mult->GetNumberOfTracklets();
 155   Int_t nSingleClusters;
 156   Int_t nClusters[6];
 157
 158   for(Int_t ilay = 0; ilay < 6; ilay++){
 159     nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 160   }
 161   nSingleClusters = mult->GetNumberOfSingleClusters();
 162
 163   Int_t nChips[2];
 164   for(Int_t ilay = 0; ilay < 2; ilay++){
 165     nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 166   }
 167
 168   // ***** FMD info
 169   AliESDFMD *fmd = esd->GetFMDData();
 170   Float_t totalMultA = 0;
 171   Float_t totalMultC = 0;
 172   const Float_t fmdLowCut = 0.4;
 173
 174   for(UShort_t det = 1;det <= 3; det++) {
 175       Int_t nRings = (det==1 ? 1 : 2);
 176       for (UShort_t ir = 0; ir < nRings; ir++) {
 177           Char_t   ring = (ir == 0 ? 'I' : 'O');
 178           UShort_t nsec = (ir == 0 ? 20  : 40);
 179           UShort_t nstr = (ir == 0 ? 512 : 256);
 180           for(UShort_t sec =0; sec < nsec;  sec++)  {
 181               for(UShort_t strip = 0; strip < nstr; strip++) {
 182                   Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
 183                   if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
 184                   Float_t nParticles=0;
 185                   if(fmdMult > fmdLowCut) {
 186                       nParticles = 1.;
 187                   }
 188
 189                   if (det<3) totalMultA = totalMultA + nParticles;
 190                   else totalMultC = totalMultC + nParticles;
 191
 192               }
 193           }
 194       }
 195   }
 196   Float_t multFMDA = totalMultA;
 197   Float_t multFMDC = totalMultC;
 198
 199   // ***** ZDC info
 200   AliESDZDC *esdZDC = esd->GetESDZDC();
 201   UInt_t esdFlag =  esdZDC->GetESDQuality();
 202
 203   Float_t znCEnergy  = (Float_t) (esdZDC->GetZDCN1Energy());
 204   Float_t zpCEnergy  = (Float_t) (esdZDC->GetZDCP1Energy());
 205   Float_t znAEnergy  = (Float_t) (esdZDC->GetZDCN2Energy());
 206   Float_t zpAEnergy  = (Float_t) (esdZDC->GetZDCP2Energy());
 207   Float_t zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
 208   Float_t zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
 209
 210   Double_t bZDC      = esdZDC->GetImpactParameter();
 211   Int_t    nPartZDC  = esdZDC->GetZDCParticipants();
 212   Double_t bZDCA     = esdZDC->GetImpactParamSideA();
 213   Int_t    nPartZDCA = esdZDC->GetZDCPartSideA();
 214   Double_t bZDCC     = esdZDC->GetImpactParamSideC();
 215   Int_t    nPartZDCC = esdZDC->GetZDCPartSideC();
 216
 217   const Double_t * towZNC = esdZDC->GetZN1TowerEnergy();
 218   const Double_t * towZPC = esdZDC->GetZP1TowerEnergy();
 219   const Double_t * towZNA = esdZDC->GetZN2TowerEnergy();
 220   const Double_t * towZPA = esdZDC->GetZP2TowerEnergy();
 221   //
 222   Float_t  znCtower[5]; //  ZNC 5 tower signals
 223   Float_t  zpCtower[5]; //  ZPC 5 tower signals
 224   Float_t  znAtower[5]; //  ZNA 5 tower signals
 225   Float_t  zpAtower[5]; //  ZPA 5 tower signals
 226   Float_t  centrZNC[2]; //  centroid over ZNC
 227   Float_t  centrZNA[2]; //  centroid over ZNA
 228
 229   for(Int_t it = 0; it < 5; it++){
 230     znCtower[it] = (Float_t) (towZNC[it]);
 231     zpCtower[it] = (Float_t) (towZPC[it]);
 232     znAtower[it] = (Float_t) (towZNA[it]);
 233     zpAtower[it] = (Float_t) (towZPA[it]);
 234   }
 235
 236   Double_t xyZNC[2] = {-99.,-99.};
 237   Double_t xyZNA[2] = {-99.,-99.};
 238
 239   esdZDC->GetZNCentroidInPbPb(beamEnergy, xyZNC, xyZNA);
 240   for(Int_t it = 0; it < 2; it++){
 241       centrZNC[it] = xyZNC[it];
 242       centrZNA[it] = xyZNA[it];
 243   }
 244
 245   // ***** MC info
 246   Double_t bMC          = 0.;
 247   Int_t specNeutronProj = 0;
 248   Int_t specProtonProj  = 0;
 249   Int_t specNeutronTarg = 0;
 250   Int_t specProtonTarg  = 0;
 251   Int_t nPartTargMC     = 0;
 252   Int_t nPartProjMC     = 0;
 253   Int_t nnColl          = 0;
 254   Int_t nnwColl         = 0;
 255   Int_t nwNColl         = 0;
 256   Int_t nwNwColl        = 0;
 257
 258   if(fIsMCInput){
 259
 260     AliMCEvent* mcEvent = MCEvent();
 261     if (!mcEvent) {
 262       printf("   Could not retrieve MC event!!!\n");
 263       return;
 264     }
 265
 266     Int_t nMyTracks_gen = 0;
 267     AliStack *stack = 0x0; // needed for MC studies
 268     stack = MCEvent()->Stack();
 269     for (Int_t iTrack = 0; iTrack < MCEvent()->GetNumberOfTracks(); iTrack++) {
 270       //get properties of mc particle
 271       AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(iTrack);
 272       // Primaries only
 273       if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
 274       //charged tracks only
 275       if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
 276       //same cuts as on ESDtracks
 277       //          if(TMath::Abs(mcP->Eta())>0.9)continue;
 278       //          if(mcP->Pt()<0.2)continue;
 279       //          if(mcP->Pt()>200)continue;
 280
 281       nMyTracks_gen ++;
 282     }
 283
 284     AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 285     if(!genHeader){
 286       printf("  Event generator header not available!!!\n");
 287       return;
 288     }
 289
 290
 291     if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
 292         bMC = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
 293         specNeutronProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsn();
 294         specProtonProj  = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsp();
 295         specNeutronTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsn();
 296         specProtonTarg  = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsp();
 297         nPartTargMC = Int_t (208.-(specNeutronTarg+specProtonTarg));
 298         nPartProjMC = Int_t (208.-(specNeutronProj+specProtonProj));
 299         nnColl   = ((AliGenHijingEventHeader*) genHeader)->NN();
 300         nnwColl  = ((AliGenHijingEventHeader*) genHeader)->NNw();
 301         nwNColl  = ((AliGenHijingEventHeader*) genHeader)->NwN();
 302         nwNwColl = ((AliGenHijingEventHeader*) genHeader)->NwNw();
 303     }
 304   }
 305
 306   fAODCentrality->SetTrigClass         (trigClass[100]);
 307   fAODCentrality->SetxVertex           (xVertex       );
 308   fAODCentrality->SetyVertex           (yVertex       );
 309   fAODCentrality->SetzVertex           (zVertex       );
 310   fAODCentrality->SetVertexer3d        (vertexer3d    );
 311   fAODCentrality->SetbMC               (bMC);
 312   fAODCentrality->SetNpartTargMC       (nPartTargMC);
 313   fAODCentrality->SetNpartProjMC       (nPartProjMC);
 314   fAODCentrality->SetNNColl            (nnColl);
 315   fAODCentrality->SetNNwColl           (nnwColl);
 316   fAODCentrality->SetNwNColl           (nwNColl);
 317   fAODCentrality->SetNwNwColl          (nwNwColl);
 318   fAODCentrality->SetNTracklets        (nTracklets);
 319   fAODCentrality->SetNSingleClusters   (nSingleClusters);
 320   fAODCentrality->SetNClusters         (
 321                                            nClusters[0],
 322                                            nClusters[1],
 323                                            nClusters[2],
 324                                            nClusters[3],
 325                                            nClusters[4],
 326                                            nClusters[5]);
 327   fAODCentrality->SetNChips            (
 328                                            nChips[0],
 329                                            nChips[1]);
 330   fAODCentrality->SetbZDC              (bZDC);
 331   fAODCentrality->SetNpartZDC          (nPartZDC);
 332   fAODCentrality->SetbZDCA             (bZDCA);
 333   fAODCentrality->SetNpartZDCA         (nPartZDCA);
 334   fAODCentrality->SetbZDCC             (bZDCC);
 335   fAODCentrality->SetNpartZDCC         (nPartZDCC);
 336   fAODCentrality->SetESDFlag      (esdFlag);
 337   fAODCentrality->SetZNCEnergy    (znCEnergy);
 338   fAODCentrality->SetZPCEnergy    (zpCEnergy);
 339   fAODCentrality->SetZNAEnergy    (znAEnergy);
 340   fAODCentrality->SetZPAEnergy    (zpAEnergy);
 341   fAODCentrality->SetZEM1Energy   (zem1Energy);
 342   fAODCentrality->SetZEM2Energy   (zem2Energy);
 343   fAODCentrality->SetZNCtower          (
 344                                            znCtower[0],
 345                                            znCtower[1],
 346                                            znCtower[2],
 347                                            znCtower[3],
 348                                            znCtower[4]);
 349   fAODCentrality->SetZPCtower          (
 350                                          zpCtower[0],
 351                                          zpCtower[1],
 352                                          zpCtower[2],
 353                                          zpCtower[3],
 354                                          zpCtower[4]);
 355
 356   fAODCentrality-> SetZNAtower          (
 357                              znAtower[0],
 358                              znAtower[1],
 359                              znAtower[2],
 360                              znAtower[3],
 361                              znAtower[4]);
 362   fAODCentrality-> SetZPAtower          (
 363                                             zpAtower[0],
 364                                             zpAtower[1],
 365                                             zpAtower[2],
 366                                             zpAtower[3],
 367                                             zpAtower[4]);
 368   fAODCentrality-> SetCentrZNC          (
 369                                             centrZNC[0],
 370                                             centrZNC[1]);
 371   fAODCentrality-> SetCentrZNA          (
 372                                             centrZNA[0],
 373                                             centrZNA[1]);
 374   fAODCentrality-> SetNTracks           (nTracks);
 375   fAODCentrality-> SetNPmdTracks        (nPmdTracks);
 376   fAODCentrality-> SetMultV0A           (multV0A);
 377   fAODCentrality-> SetMultV0C           (multV0C);
 378   fAODCentrality-> SetMultFMDA          (multFMDA);
 379   fAODCentrality-> SetMultFMDC          (multFMDC);
 380
 381 //
 382 // Header Replication
 383 //
 384   AliAODHeader* hdr = AODEvent()->GetHeader();
 385   *fAODHeader =  *hdr;
 386 //
 387   return;
 388 }
 389
 390 //________________________________________________________________________
 391 void AliAnalysisTaskAODCentralityMaker::Terminate(Option_t */*option*/)
 392 {
 393   // Terminate analysis
 394   //
 395   if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker: Terminate() \n");
 396 }