PWGUD/DIFFRACTIVE/example/AliAnalysisTaskCDex.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 // Example analysis for diffractive studies
  17 //
  18 // Author:
  19 //  Felix Reidt <Felix.Reidt@cern.ch>
  20
  21
  22 #include <TH1.h>
  23 #include <TH2.h>
  24 #include <TList.h>
  25
  26
  27 #include "AliAODInputHandler.h"
  28 #include "AliAODHandler.h"
  29 #include "AliESDInputHandler.h"
  30 #include "AliPIDResponse.h"
  31 #include "AliPhysicsSelection.h"
  32
  33 #include "AliCDMesonBaseStripped.h"
  34 #include "AliCDMesonUtilsStripped.h"
  35 #include "AliCDMesonTracks.h"
  36 #include "AliAnalysisTaskCDex.h"
  37
  38
  39 //------------------------------------------------------------------------------
  40 AliAnalysisTaskCDex::AliAnalysisTaskCDex(const char* name):
  41         AliAnalysisTaskSE(name)
  42         , fDoAOD(kFALSE)
  43         , fMaxVtxDst(0.5) // value to be checked with the vertex study histograms
  44
  45         , fESDEvent(0x0)
  46         , fAODEvent(0x0)
  47         , fPIDResponse(0x0)
  48         , fTracks(new AliCDMesonTracks())
  49         , fVtxDst(-1.)
  50         , fVtxZ(-20)
  51         , fResidualTracks(0)
  52         , fResidualTracklets(0)
  53         , fMCprocessType(0)
  54         , fMCprocess(-1)
  55
  56         , fRun(-999)
  57         , fCurrentGapCondition(0)
  58
  59         , fHist(0x0)
  60
  61         , fv0ntrk(0x0)
  62         , fv0fmdntrk(0x0)
  63         , fv0fmdspdntrk(0x0)
  64         , fv0fmdspdtpcntrk(0x0)
  65
  66         , fhStatsFlow(0x0)
  67 {
  68         //
  69         // standard constructor
  70         //
  71         // slot in TaskSE must start from 1
  72         DefineOutput(1, TList::Class());
  73
  74         // initialize gap information
  75         for (Int_t iGap = 0; iGap < kMax; ++iGap) {
  76                 fGapInformation[iGap] = 0;
  77         }
  78 }
  79
  80
  81 //------------------------------------------------------------------------------
  82 AliAnalysisTaskCDex::~AliAnalysisTaskCDex()
  83 {
  84         //
  85         // destructor
  86         //
  87
  88         if (fHist) {
  89                 fHist->Clear();
  90                 delete fHist;
  91                 fHist = 0x0;
  92         }
  93
  94         if (fTracks) {
  95                 delete fTracks;
  96                 fTracks = 0x0;
  97         }
  98 }
  99
 100
 101 //------------------------------------------------------------------------------
 102 void AliAnalysisTaskCDex::UserCreateOutputObjects()
 103 {
 104         //
 105         // createOutputObjects
 106         //
 107
 108         //= TList for Histograms =====================================================
 109         fHist = new TList;
 110         fHist->SetOwner(); // ensures that the histograms are all deleted on exit!
 111
 112         //= MULTIPLICITY PER GAP CONDITION =
 113         fv0ntrk = new TH2D("b00_v0ntrk", ";number of tracks;gap condition",
 114                            80, 0., 80., 4, 1., 5.);
 115         //x: ntrk; y: V0
 116         fHist->Add(fv0ntrk);
 117
 118         fv0fmdntrk = new TH2D("b01_v0fmdntrk", ";number of tracks;gap condition",
 119                                       80, 0., 80., 4, 1., 5.);
 120         //x: ntrk; y: V0FMD
 121         fHist->Add(fv0fmdntrk);
 122
 123         fv0fmdspdntrk =
 124                 new TH2D("b02_v0fmdspdntrk", ";number of tracks;gap condition",
 125                          80, 0., 80., 4, 1., 5.);
 126         //x: ntrk; y: V0FMDSPD
 127         fHist->Add(fv0fmdspdntrk);
 128
 129         fv0fmdspdtpcntrk =
 130                 new TH2D("b03_v0fmdspdtpcntrk", ";number of tracks;gap condition",
 131                          80, 0., 80., 4, 1., 5.);
 132         //x: ntrk; y: V0FMDSPDTPC
 133         fHist->Add(fv0fmdspdtpcntrk);
 134
 135
 136         //= STATISTICS FLOW =
 137         fhStatsFlow = AliCDMesonBaseStripped::GetHistStatsFlow();
 138         fHist->Add(fhStatsFlow);
 139 }
 140
 141
 142 //------------------------------------------------------------------------------
 143 void AliAnalysisTaskCDex::UserExec(Option_t *)
 144 {
 145         //
 146         // executed for every event which passed the physics selection
 147         //
 148         // in order to select only correct minimum bias events,
 149         // SetCollisionCandidates(AliVEvent::kMB) should be used
 150
 151         fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinTotalInput); // stats flow
 152
 153         //= INPUT DATA SANITY TESTS ==================================================
 154         if(!CheckInput()) {
 155                 PostOutputs();
 156                 return;
 157         }
 158
 159         fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinGoodInput); // stats flow
 160
 161         //= EVENT SELECTION ==========================================================
 162         Bool_t eventIsValid = (fDoAOD) ?
 163                 AliCDMesonUtilsStripped::CutEvent(fAODEvent) :
 164                 AliCDMesonUtilsStripped::CutEvent(fESDEvent);
 165         if (!eventIsValid) {
 166                 PostOutputs();
 167                 return;
 168         }
 169
 170         fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinEventsAfterCuts); // stats flow
 171
 172         //= PILE UP ==================================================================
 173         const Bool_t isPileup = (fDoAOD) ?
 174                 fAODEvent->IsPileupFromSPD(2, 0.8, 3., 2., 5.) :
 175                 fESDEvent->IsPileupFromSPD(2, 0.8, 3., 2., 5.);
 176         // using only 2 instead of three contributors
 177
 178         if (isPileup) {
 179                 PostOutputs();
 180                 return;
 181         }
 182
 183         fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinEventsWithOutPileUp);
 184
 185         //= GAP ======================================================================
 186         // determine the complete gap configuration (for all gap-tagging detectors)
 187         if (!DetermineGap()) {
 188                 PostOutputs();
 189                 return;
 190         }
 191
 192         //= VERTEX COINCIDENCE AND POSITION ==========================================
 193         AnalyzeVtx();
 194         if (!(abs(fVtxZ) < 4.)) { // vertex from tracks within +/-4cm
 195                 //PostOutputs();
 196                 //return;
 197         }
 198
 199         //= TRACK CUTS ===============================================================
 200         fTracks->ProcessEvent(fAODEvent, fESDEvent, kTRUE);
 201         // apply cuts (including soft)
 202         DoMultiplicityStudy(); // fill corresponding histograms
 203
 204         // is multiplicity within the desired range of  2 to 3?
 205         Int_t nch = fTracks->GetTracks();
 206         Int_t ncombined = fTracks->GetCombinedTracks();
 207
 208         //============================================================================
 209         //=== USER ANALYSIS CODE =====================================================
 210         //============================================================================
 211         for (Int_t iTrack = 0; iTrack < ncombined; ++iTrack) { // including soft
 212                 AliVTrack* trk = fTracks->GetTrack(iTrack);
 213                 trk->GetID(); // prevent warning...
 214         }
 215         for (Int_t iTrack = 0; iTrack < nch; ++iTrack) { // excluding soft tracks
 216                 AliVTrack* trk = fTracks->GetTrack(iTrack);
 217                 trk->GetID(); // prevent warning...
 218         }
 219
 220         if (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBaseStripped::kBinDG) {
 221                 // event is a full double gap event
 222         }
 223
 224         //============================================================================
 225         PostOutputs();
 226 }
 227
 228
 229 //------------------------------------------------------------------------------
 230 void AliAnalysisTaskCDex::PostOutputs()
 231 {
 232         //
 233         // PostData
 234         //
 235         // this function is main of use with multiple output containers
 236
 237         PostData(1, fHist);
 238 }
 239
 240
 241 //------------------------------------------------------------------------------
 242 Bool_t AliAnalysisTaskCDex::CheckInput()
 243 {
 244         //
 245         // general protection of the task against malicious input data
 246         //
 247         if (const AliESDInputHandler *esdH =
 248             dynamic_cast<AliESDInputHandler*>(fInputHandler)){
 249                 fESDEvent = esdH->GetEvent();
 250         }
 251         else if (const AliAODInputHandler *aodH =
 252                  dynamic_cast<AliAODInputHandler*>(fInputHandler)){
 253                 fAODEvent = aodH->GetEvent();
 254                 fDoAOD = kTRUE;
 255         }
 256         fPIDResponse = (AliPIDResponse*)fInputHandler->GetPIDResponse();
 257
 258         if(!fESDEvent && !fAODEvent){
 259                 printf("AliAnalysisTaskex - No valid event\n");
 260                 return kFALSE;
 261         }
 262
 263         if(!fPIDResponse){
 264                 printf("AliAnalysisTaskex -  No PIDd\n");
 265                 // PID is fixed to unknown
 266                 //return kFALSE;
 267         }
 268
 269         if(fDoAOD && fAODEvent && fabs(fAODEvent->GetMagneticField())<1){
 270                 printf("AliAnalysisTaskex - strange Bfield! %f\n",
 271                        fAODEvent->GetMagneticField());
 272                 return kFALSE;
 273         }
 274         else if((!fDoAOD) && fESDEvent && fabs(fESDEvent->GetMagneticField())<1){
 275                 printf("AliAnalysisTaskex - strange Bfield! %f\n",
 276                        fESDEvent->GetMagneticField());
 277                 return kFALSE;
 278         }
 279
 280         Int_t tmprun = 0;
 281         if (fDoAOD && fAODEvent) {
 282                 tmprun = fAODEvent->GetRunNumber();
 283         }
 284         else if (fESDEvent) {
 285                 tmprun = fESDEvent->GetRunNumber();
 286         }
 287
 288         if(fRun!=tmprun){
 289                 fRun = tmprun;
 290                 AliCDMesonUtilsStripped::SPDLoadGeom(fRun);
 291         }
 292
 293         return kTRUE;
 294 }
 295
 296
 297 //------------------------------------------------------------------------------
 298 Bool_t AliAnalysisTaskCDex::DetermineGap()
 299 {
 300         // determines the gap configuration for all gap tagging detectors based on the
 301         // data set which is available
 302         //
 303
 304         if (fDoAOD) {
 305                 AliAODHandler* aodHandler =
 306                         (AliAODHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 307                 TTree *aodTree = aodHandler->GetTree();
 308                 aodTree->SetBranchAddress("gapCondition", &fCurrentGapCondition);
 309                 aodTree->GetEvent(Entry()); // seems to be needed! (loads current event)
 310                 if (!fCurrentGapCondition) {
 311                         fCurrentGapCondition = 0xfffe;
 312                         puts("AliAnalysisTaskCDex - ");
 313                         puts("error while gap condition determination using AODs\n");
 314                         return kFALSE;
 315                 }
 316         }
 317         else {
 318                 // gap determination from detector information
 319                 fCurrentGapCondition = AliCDMesonUtilsStripped::GetGapConfig(fESDEvent);
 320
 321                 // gap determination from preprocessed detector information
 322                 /*
 323                 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 324                 TTree* esdTree = am->GetInputEventHandler()->GetTree(); // get ESD tree
 325                 if (esdTree) {
 326                         esdTree->SetBranchAddress("gapCondition", &fCurrentGapCondition);
 327                         esdTree->GetEvent(Entry()); // seems to be needed! (loads current event)
 328                 }
 329                 */
 330
 331                 if (!fCurrentGapCondition) {
 332                         fCurrentGapCondition = 0xfffe;
 333                         puts("AliAnalysisTaskCDex - ");
 334                         puts("error while gap condition determination using ESDs\n");
 335                         return kFALSE;
 336                 }
 337         }
 338
 339         // disentagle the contributions to the gap conditions of different "tightness"
 340         fGapInformation[kV0] =
 341                 AliCDMesonBaseStripped::GetGapBin("V0", fCurrentGapCondition);
 342         fGapInformation[kV0FMD] =
 343                 AliCDMesonBaseStripped::GetGapBin("V0FMD", fCurrentGapCondition);
 344         fGapInformation[kV0FMDSPD] =
 345                 AliCDMesonBaseStripped::GetGapBin("V0FMDSPD", fCurrentGapCondition);
 346         fGapInformation[kV0FMDSPDTPC] =
 347                 AliCDMesonBaseStripped::GetGapBin("V0FMDSPDTPC", fCurrentGapCondition);
 348         fGapInformation[kFMD] =
 349                 AliCDMesonBaseStripped::GetGapBin("FMD",fCurrentGapCondition);
 350         fGapInformation[kSPD] =
 351                 AliCDMesonBaseStripped::GetGapBin("SPD",fCurrentGapCondition);
 352         fGapInformation[kTPC] =
 353                 AliCDMesonBaseStripped::GetGapBin("TPC",fCurrentGapCondition);
 354
 355         return kTRUE;
 356 }
 357
 358
 359 //------------------------------------------------------------------------------
 360 void AliAnalysisTaskCDex::DoMultiplicityStudy()
 361 {
 362         // stores the multiplicity distributions for different gap conditions and
 363         // adds some information to the statsFlow histogram
 364         //
 365
 366         // retrieve values from the track object
 367         Int_t ntrk0 = fTracks->GetTracksBeforeCuts(); // number of tracks before cuts
 368         //Int_t nch = fTracks->GetTracks(); // number of good ITS-TPC primaries
 369         Int_t ncombined = fTracks->GetCombinedTracks(); // number ITSTPC and ITS only
 370         Int_t nITSpureSA = fTracks->GetITSpureSACount(); // number ITS standalone
 371
 372         // determine the residual tracks / tracklets
 373         fResidualTracks = ntrk0 - ncombined - nITSpureSA;
 374         fResidualTracklets = fTracks->GetRemainingTrackletsCentralBarrel();
 375
 376         // multiplicity distributions for different gaps
 377         fv0ntrk->Fill(ncombined, fGapInformation[kV0]);
 378         fv0fmdntrk->Fill(ncombined, fGapInformation[kV0FMD]);
 379         fv0fmdspdntrk->Fill(ncombined, fGapInformation[kV0FMDSPD]);
 380         fv0fmdspdtpcntrk->Fill(ncombined, fGapInformation[kV0FMDSPDTPC]);
 381
 382         if (fGapInformation[kV0] == AliCDMesonBaseStripped::kBinDG) {
 383                 fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinv0Gap);
 384         }
 385         if (fGapInformation[kV0FMD] == AliCDMesonBaseStripped::kBinDG) {
 386                 fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinv0fmdGap);
 387         }
 388         if (fGapInformation[kV0FMDSPD] == AliCDMesonBaseStripped::kBinDG) {
 389                 fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinv0fmdspdGap);
 390         }
 391         if (fGapInformation[kV0FMDSPDTPC] == AliCDMesonBaseStripped::kBinDG) {
 392                 fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinv0fmdspdtpcGap);
 393         }
 394
 395         // event cleanliness
 396         if (fResidualTracks == 0) {
 397                 fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinResidualTracks);
 398         }
 399         if (fResidualTracklets == 0) {
 400                 fhStatsFlow->Fill(AliCDMesonBaseStripped::kBinResidualTracklets);
 401         }
 402 }
 403
 404
 405 //--------------------------------------------------------------------------
 406 void AliAnalysisTaskCDex::AnalyzeVtx()
 407 {
 408         // calculates the distance between the vertex obtain from tracks and the
 409         // vertex obtain from spd tracklets
 410         // stores the z position of the primary vertex from tracks
 411
 412         fVtxDst = 0.; // reset distance
 413
 414         // retrieve the pointers of the current primary vertices
 415         AliVVertex* trackVtx = (fDoAOD) ?
 416                 (AliVVertex*)fAODEvent->GetPrimaryVertex() :
 417                 (AliVVertex*)fESDEvent->GetPrimaryVertexTracks();
 418         AliVVertex* spdVtx = (fDoAOD) ?
 419                 (AliVVertex*)fAODEvent->GetPrimaryVertexSPD() :
 420                 (AliVVertex*)fESDEvent->GetPrimaryVertexSPD();
 421
 422         fVtxZ = trackVtx->GetZ(); // store the vertex z position
 423
 424         if (fDoAOD && (trackVtx == spdVtx)) { // no primary track vertex in the AOD
 425                 fVtxDst = -5; // set arbitrary distance (counted in underflow bin!)
 426         }
 427         else { // do proper calculation of the geometrical distance
 428                 fVtxDst += (trackVtx->GetX()-spdVtx->GetX())
 429                         * (trackVtx->GetX()-spdVtx->GetX());
 430                 fVtxDst += (trackVtx->GetY()-spdVtx->GetY())
 431                         * (trackVtx->GetY()-spdVtx->GetY());
 432                 fVtxDst += (trackVtx->GetZ()-spdVtx->GetZ())
 433                         * (trackVtx->GetZ()-spdVtx->GetZ());
 434                 fVtxDst = TMath::Sqrt(fVtxDst);
 435         }
 436 }