PWGLF/SPECTRA/XtAnalysis/AliXtAnalysis.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2014, 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  *     AliXtAnalysis:
  18  * This class constructs inclusive and isolated (based on charged tracks)
  19  * invariant spectra. Isolated case uses a different efficiency where the
  20  * contamination from cases where non-isolated particle appears as an
  21  * isolated one due to finite detector efficiency is taken into account.
  22  *
  23  * contact: Sami Räsänen
  24  *          University of Jyväskylä, Finland
  25  *          sami.s.rasanen@jyu.fi
  26 **************************************************************************/
  27
  28 // general + root classes
  29 #include <TList.h>
  30 #include <TChain.h>
  31 #include <TObjArray.h>
  32
  33 // AliRoot classes
  34 #include <AliAnalysisManager.h>
  35 #include <AliInputEventHandler.h>
  36 #include <AliAnalysisUtils.h>
  37 #include <AliAODTrack.h>
  38 #include <AliAODEvent.h>
  39
  40 // Jyväskylä classes
  41 #include <AliJCard.h>  // this class at PWGCF / JCORRAN
  42 #include "AliJXtHistos.h"
  43 //#include "AliIsolatedEfficiency.h"
  44
  45 // This analysis
  46 #include "AliXtAnalysis.h"
  47
  48 ClassImp(AliXtAnalysis)
  49
  50 //________________________________________________________________________
  51 AliXtAnalysis::AliXtAnalysis()
  52       : AliAnalysisTaskSE(),
  53         fOutputList(0x0),
  54         fAnaUtils(0x0),
  55         fCard(0x0),
  56         fHistDir(0x0),
  57         fHistos(0x0),
  58         fhEvents(0x0),
  59         //fEfficiency(0x0),
  60         fChargedList(0x0),
  61         fIsolatedChargedList(0x0),
  62         fCentBin(-1),
  63         fZBin(-1),
  64         fZVert(999999.),
  65         fevt(0),
  66         fDebugMode(0)
  67 {
  68 // Constructor
  69 }
  70 //________________________________________________________________________
  71 AliXtAnalysis::AliXtAnalysis(const char *name, const char * cardname)
  72       : AliAnalysisTaskSE(name),
  73         fOutputList(0x0),
  74         fAnaUtils(0x0),
  75         fCard(0x0),
  76         fHistDir(0x0),
  77         fHistos(0x0),
  78         fhEvents(0x0),
  79         //fEfficiency(0x0),
  80         fChargedList(0x0),
  81         fIsolatedChargedList(0x0),
  82         fCentBin(-1),
  83         fZBin(-1),
  84         fZVert(999999.),
  85         fevt(0),
  86         fDebugMode(0)
  87 {
  88         // All parameters of the analysis
  89         fCard = new AliJCard(cardname);
  90
  91         cout << "debug: card created to address " << fCard << endl;
  92
  93         fCard->PrintOut();
  94
  95         fAnaUtils = new AliAnalysisUtils();
  96         fAnaUtils->SetUseOutOfBunchPileUp( kTRUE );
  97         fAnaUtils->SetUseSPDCutInMultBins( kTRUE);
  98
  99         // Define input and output slots here
 100         // Input slot #0 works with a TChain
 101         DefineInput(0, TChain::Class());
 102         // Output slot #0 writes into a TH1 container
 103         DefineOutput(1, TList::Class());
 104         // JHistos into TDirectory
 105         DefineOutput(2, TDirectory::Class());
 106 }
 107 //________________________________________________________________________
 108 AliXtAnalysis::AliXtAnalysis(const AliXtAnalysis& a)
 109       : AliAnalysisTaskSE(a.GetName()),
 110         fOutputList(a.fOutputList),
 111         fAnaUtils(a.fAnaUtils),
 112         fCard(a.fCard),
 113         fHistDir(a.fHistDir),
 114         fHistos(a.fHistos),
 115         fhEvents(a.fhEvents),
 116         //fEfficiency(a.fEfficiency),
 117         fChargedList(a.fChargedList),
 118         fIsolatedChargedList(a.fIsolatedChargedList),
 119         fCentBin(-1),
 120         fZBin(-1),
 121         fZVert(999999.),
 122         fevt(0),
 123         fDebugMode(0)
 124 {
 125     //copy constructor
 126     fhEvents = (TH1D*) a.fhEvents->Clone(a.fhEvents->GetName());
 127 }
 128 //________________________________________________________________________
 129 AliXtAnalysis& AliXtAnalysis::operator = (const AliXtAnalysis& ap){
 130         // assignment operator
 131         this->~AliXtAnalysis();
 132         new(this) AliXtAnalysis(ap);
 133         return *this;
 134 }
 135 //________________________________________________________________________
 136 AliXtAnalysis::~AliXtAnalysis() {
 137     // destructor
 138     delete fOutputList;
 139     delete [] fAnaUtils;
 140     delete [] fhEvents;
 141     delete fHistos;
 142     //delete fEfficiency;
 143     delete fCard;
 144     delete fChargedList;
 145     delete fIsolatedChargedList;
 146 }
 147 //________________________________________________________________________
 148 void AliXtAnalysis::UserCreateOutputObjects(){
 149     // Create histograms
 150     // Called once
 151
 152     cout<<"\n=============== CARD =============="<<endl;
 153     fCard->PrintOut();
 154     cout<<"===================================\n"<<endl;
 155
 156     fhEvents = new TH1D("hEvents","events passing cuts", 11, -0.5, 10.5 );
 157     fhEvents->SetXTitle( "0 - all, 1 - pileup, 2 - kMB selected, 3 - has vertex, 4 - good vertex, 5 - MB + good vertex, 6 - MB + good vertex + centrality" );
 158     fhEvents->Sumw2();
 159
 160     fOutputList = new TList();
 161     fOutputList->SetOwner(kTRUE);
 162     //for any MC track filled with MC pt
 163     fOutputList->Add(fhEvents);
 164
 165     fChargedList = new TObjArray;
 166     fChargedList->SetOwner(kTRUE);
 167     fIsolatedChargedList = new TObjArray;
 168
 169     PostData(1, fOutputList);
 170
 171     bool orignalTH1AdddirectoryStatus=TH1::AddDirectoryStatus();
 172     TH1::AddDirectory(kTRUE);
 173     if( !orignalTH1AdddirectoryStatus ) cout<<"DEBUG : TH1::AddDirectory is turned on"<<endl;
 174     //TFile * file2 = OpenFile(2);
 175
 176     fHistDir=gDirectory;
 177     fHistos = new AliJXtHistos(fCard);
 178     fHistos->CreateXtHistos();
 179
 180     PostData(2, fHistDir);
 181     TH1::AddDirectory( orignalTH1AdddirectoryStatus );
 182     cout<<"DEBUG : TH1::AddDirectory get orignal Value = "<<( orignalTH1AdddirectoryStatus?"True":"False" )<<endl;
 183
 184     //fEfficiency = new AliJEfficiency();
 185
 186     fevt = 0;
 187 }
 188 //________________________________________________________________________
 189 void AliXtAnalysis::UserExec(Option_t *) {
 190
 191     // Main loop - called for each event
 192
 193     fevt++;
 194     if(fevt % 10000 == 0) cout << "Number of event scanned = "<< fevt << endl;
 195
 196     Int_t filterBit = AliAODTrack::kTrkGlobal; //fCard->Get("TrackFilterBit");
 197
 198     AliVEvent *event = InputEvent();
 199     if(!event) return;
 200
 201     // check if the event was triggered or not and vertex
 202     if( !IsGoodEvent( event ) ) return; // fZBin is set there
 203     fhEvents->Fill( 5 );
 204     if(fDebugMode > 0) cout << "zvtx = " << fZVert << endl;
 205
 206     AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(event);
 207     if(!aodEvent) return;
 208
 209     // centrality - in case of pPb and PbPb, dictated by DataType in the card
 210     double fcent;
 211     if( fCard->Get("DataType") > 0.5 ){
 212         AliCentrality *cent = event->GetCentrality();
 213         if( ! cent ) return;
 214         fcent = cent->GetCentralityPercentile("V0M");
 215         fCentBin = fCard->GetBin(kCentrType, fcent);;
 216         //cout <<"Centrality = "<< fcent <<"\t"<< fCentBin << endl;
 217     }else{
 218         fcent = 0.0;
 219         fCentBin = 0;
 220     }
 221
 222     if(fCentBin<0) return;
 223     fhEvents->Fill( 6 );
 224     Int_t nt = aodEvent->GetNumberOfTracks();
 225
 226     fHistos->FillCentralityHistos(fcent, fCentBin);
 227
 228     // clear them up for every event
 229     fChargedList->Clear();
 230     fIsolatedChargedList->Clear();
 231
 232     // Is there a better way than hard code this into card???
 233     double sqrts = fCard->Get("Sqrts");
 234
 235     // Isolation method: 0 = relative, 1 = absolute
 236     double isolMethod = fCard->Get("IsolationMethod");
 237
 238     // Minimum isolation pT that is considered
 239     double minIsolationPt = fCard->Get("MinimumIsolatedPt");
 240
 241     // Dummy numbers; either of these is updated based on selected method (other is void)
 242     double isolationFraction  = 99999.;
 243     double isolationThreshold = 99999.;
 244
 245     if( isolMethod > 0.5 ){
 246         isolationThreshold = fCard->Get("IsolationThreshold");  // this is fixed threshold
 247     }else{
 248         isolationFraction = fCard->Get("IsolationFraction");  // fraction in pT, threshold updated later in a loop
 249     }
 250
 251     // Isolation radius
 252     double isolR = fCard->Get("IsolationRadius");
 253
 254     // Get good tracks to the list
 255     for(Int_t it = 0; it < nt; it++) {
 256         AliAODTrack *track = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(it));
 257         if( !track ) continue;
 258         if( !track->TestFilterBit(filterBit) ) continue;
 259         double eta = track->Eta();
 260         if( !fCard->IsInEtaRange( eta ) ) continue; // Do not use fiducial cut here
 261         AliAODTrack *acceptedTrack = new AliAODTrack( *track );
 262         double pT = acceptedTrack->Pt();
 263         double xT = 2.*pT/sqrts;
 264         double phi = acceptedTrack->Phi();
 265
 266         // TODO:
 267         double effCorr = 1.0;
 268
 269         // Fill histos
 270         fHistos->FillInclusiveHistograms(pT, xT, eta, phi, effCorr, fCentBin);
 271
 272         // Add an accepted track into list of inclusive charged particles
 273         fChargedList->Add( acceptedTrack );
 274     }
 275
 276     // Check isolation of the found good tracks
 277     Int_t nAcc = fChargedList->GetEntriesFast();
 278
 279     // Get eta range for fiducial cut
 280     double etaRange = fCard->Get("EtaRange");
 281
 282     for( Int_t it = 0; it < nAcc; it++ ){
 283         AliAODTrack *track = (AliAODTrack*)fChargedList->At(it);
 284         // Fiducial eta cut for isolated particles
 285         double eta = track->Eta();
 286         if( abs(eta) > etaRange - isolR ) continue;
 287
 288         // To check the isolation with cone, use TLorentzVector -routines
 289         TLorentzVector lvTrack = TLorentzVector( track->Px(), track->Py(), track->Pz(), track->P() );
 290
 291         // Set isolation threshold, if relative isolation
 292         double pT = lvTrack.Pt();
 293         if( pT < minIsolationPt ) continue;
 294         if( isolMethod < 0.5 ) isolationThreshold = pT * isolationFraction; // here: relative isolation
 295
 296         // Isolation check
 297         double sum = 0.0;
 298         for( Int_t ia = 0; ia < nAcc; ia++ ){
 299             if( ia == it ) continue; // Do not count the particle itself
 300             AliAODTrack *assocTrack = (AliAODTrack*)fChargedList->At(ia);
 301             TLorentzVector lvAssoc = TLorentzVector( assocTrack->Px(), assocTrack->Py(), assocTrack->Pz(), assocTrack->P() );
 302             if( lvTrack.DeltaR( lvAssoc ) < isolR ) sum += lvAssoc.Pt();
 303         }
 304
 305         // Cone activity
 306         fHistos->FillInclusiveConeActivities(pT,sum);  // efficiency correction?
 307
 308         // If the particle was isolated, fill histograms and the list
 309         if( sum < isolationThreshold ){
 310
 311             double xT = 2.*pT/sqrts;
 312
 313             if( fDebugMode > 0 ){
 314                 cout << "DEBUG: isolated particle found " << endl;
 315                 cout << "fCentBin = " << fCentBin << ", pT = " << pT << ", xT = " << xT << " and eta = " << eta << endl;
 316             }
 317
 318             fHistos->FillIsolatedConeActivities(pT, sum);  // efficiency correction?
 319
 320             // TODO:
 321             double effCorr = 1.0;
 322
 323             // Fill histos
 324             fHistos->FillIsolatedHistograms(pT, xT, eta, track->Phi(), effCorr, fCentBin);
 325
 326             // Add tracks to the isolated list
 327             AliAODTrack *acceptedIsolatedTrack = new AliAODTrack( *track );
 328             fIsolatedChargedList->Add( acceptedIsolatedTrack );
 329         }
 330     }
 331
 332     // For further analysis, there is now lists of all charged and isolated charged available, will be usefull.
 333
 334     return;  // END
 335 }
 336 //________________________________________________________________________
 337 void AliXtAnalysis::Terminate(Option_t *)
 338 {
 339     cout<<"Successfully finished"<<endl;
 340 }
 341 //________________________________________________________________________
 342 bool AliXtAnalysis::IsGoodEvent(AliVEvent *event) {
 343
 344     // This function checks that the event has the requested trigger
 345     // and that the z-vertex is in given range
 346
 347
 348     fhEvents->Fill( 0 );
 349
 350     if(fAnaUtils->IsPileUpEvent(event)) {
 351         return kFALSE;
 352     } else {
 353         Bool_t triggeredEventMB = kFALSE; //init
 354
 355         fhEvents->Fill( 1 );
 356
 357         Bool_t triggerkMB = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & ( AliVEvent::kMB );
 358
 359         if( triggerkMB ){
 360             triggeredEventMB = kTRUE;  //event triggered as minimum bias
 361             fhEvents->Fill( 2 );
 362         }
 363         //--------------------------------------------------------------
 364         // check reconstructed vertex
 365         int ncontributors = 0;
 366         Bool_t goodRecVertex = kFALSE;
 367         const AliVVertex *vtx = event->GetPrimaryVertex();
 368         if(vtx){
 369             fhEvents->Fill( 3 );
 370             fZVert = vtx->GetZ();
 371             fHistos->FillRawVertexHisto(fZVert);
 372             ncontributors = vtx->GetNContributors();
 373             if(ncontributors > 0){
 374                 if(fCard->VertInZRange(fZVert)) {
 375                     goodRecVertex = kTRUE;
 376                     fhEvents->Fill( 4 );
 377                     fHistos->FillAcceptedVertexHisto(fZVert, fCentBin);
 378                     fZBin  = fCard->GetBin(kZVertType, fZVert);
 379                 }
 380             }
 381         }
 382         return triggerkMB && goodRecVertex;
 383     }
 384     //---------------------------------
 385 }