PWG2/AliProtonCorrectionTask.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 //-----------------------------------------------------------------------
  17 // Example of task (running locally, on AliEn and CAF),
  18 // which provides standard way of calculating acceptance and efficiency
  19 // between different steps of the procedure.
  20 // The ouptut of the task is a AliCFContainer from which the efficiencies
  21 // can be calculated
  22 //-----------------------------------------------------------------------
  23 // Author : R. Vernet, Consorzio Cometa - Catania (it)
  24 //-----------------------------------------------------------------------
  25
  26
  27 #ifndef ALIPROTONCORRECTIONTASK_CXX
  28 #define ALIPROTONCORRECTIONTASK_CXX
  29
  30 #include "AliProtonCorrectionTask.h"
  31 #include "TCanvas.h"
  32 #include "AliStack.h"
  33 #include "TParticle.h"
  34 #include "TH1I.h"
  35 #include "AliMCEvent.h"
  36 #include "AliAnalysisManager.h"
  37 #include "AliESDEvent.h"
  38 #include "AliAODEvent.h"
  39 #include "AliCFManager.h"
  40 #include "AliCFCutBase.h"
  41 #include "AliCFContainer.h"
  42 #include "TChain.h"
  43 #include "AliESDtrack.h"
  44 #include "AliLog.h"
  45 //#include "AliTPCtrack.h"
  46
  47 ClassImp(AliProtonCorrectionTask)
  48
  49 //__________________________________________________________________________
  50 AliProtonCorrectionTask::AliProtonCorrectionTask() :
  51   fReadTPCTracks(0),
  52   fReadAODData(0),
  53   fCFManagerProtons(0x0),
  54   fCFManagerAntiProtons(0x0),
  55   fQAHistList(0x0),
  56   fHistEventsProcessed(0x0)
  57 {
  58   //
  59   //Default ctor
  60   //
  61 }
  62 //___________________________________________________________________________
  63 AliProtonCorrectionTask::AliProtonCorrectionTask(const Char_t* name) :
  64   AliAnalysisTaskSE(name),
  65   fReadTPCTracks(0),
  66   fReadAODData(0),
  67   fCFManagerProtons(0x0),
  68   fCFManagerAntiProtons(0x0),
  69   fQAHistList(0x0),
  70   fHistEventsProcessed(0x0) {
  71   //
  72   // Constructor. Initialization of Inputs and Outputs
  73   //
  74   Info("AliProtonCorrectionTask","Calling Constructor");
  75
  76   /*
  77     DefineInput(0) and DefineOutput(0)
  78     are taken care of by AliAnalysisTaskSE constructor
  79   */
  80   DefineOutput(1,TH1I::Class());
  81   DefineOutput(2,AliCFContainer::Class());
  82   DefineOutput(3,AliCFContainer::Class());
  83   DefineOutput(4,TList::Class());
  84 }
  85
  86 //___________________________________________________________________________
  87 AliProtonCorrectionTask& AliProtonCorrectionTask::operator=(const AliProtonCorrectionTask& c) {
  88   //
  89   // Assignment operator
  90   //
  91   if (this!=&c) {
  92     AliAnalysisTaskSE::operator=(c) ;
  93     fReadTPCTracks = c.fReadTPCTracks ;
  94     fReadAODData = c.fReadAODData ;
  95     fCFManagerProtons  = c.fCFManagerProtons;
  96     fCFManagerAntiProtons  = c.fCFManagerAntiProtons;
  97     fQAHistList = c.fQAHistList ;
  98     fHistEventsProcessed = c.fHistEventsProcessed;
  99   }
 100   return *this;
 101 }
 102
 103 //___________________________________________________________________________
 104 AliProtonCorrectionTask::AliProtonCorrectionTask(const AliProtonCorrectionTask& c) :
 105   AliAnalysisTaskSE(c),
 106   fReadTPCTracks(c.fReadTPCTracks),
 107   fReadAODData(c.fReadAODData),
 108   fCFManagerProtons(c.fCFManagerProtons),
 109   fCFManagerAntiProtons(c.fCFManagerAntiProtons),
 110   fQAHistList(c.fQAHistList),
 111   fHistEventsProcessed(c.fHistEventsProcessed) {
 112   //
 113   // Copy Constructor
 114   //
 115 }
 116
 117 //___________________________________________________________________________
 118 AliProtonCorrectionTask::~AliProtonCorrectionTask() {
 119   //
 120   //destructor
 121   //
 122   Info("~AliProtonCorrectionTask","Calling Destructor");
 123   if (fCFManagerProtons) delete fCFManagerProtons ;
 124   if (fCFManagerAntiProtons) delete fCFManagerAntiProtons ;
 125   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 126   if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
 127 }
 128
 129 //_________________________________________________
 130 void AliProtonCorrectionTask::UserExec(Option_t *) {
 131   //
 132   // Main loop function
 133   //
 134   Info("UserExec","") ;
 135
 136   AliVEvent*    fEvent = fInputEvent ;
 137   AliVParticle* track ;
 138
 139   if (!fEvent) {
 140     Error("UserExec","NO EVENT FOUND!");
 141     return;
 142   }
 143
 144   if (!fMCEvent) Error("UserExec","NO MC INFO FOUND");
 145
 146   //pass the MC evt handler to the cuts that need it
 147   fCFManagerProtons->SetEventInfo(fMCEvent);
 148   fCFManagerAntiProtons->SetEventInfo(fMCEvent);
 149
 150   // MC-event selection
 151   Double_t containerInput[2] ;
 152
 153   //__________________________________________________________//
 154   //loop on the MC event
 155   for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 156     AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);
 157
 158     //check the MC-level cuts
 159     if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 160
 161     containerInput[0] = Rapidity(mcPart->Px(),
 162                                  mcPart->Py(),
 163                                  mcPart->Pz());
 164     containerInput[1] = (Float_t)mcPart->Pt();
 165     //fill the container for Gen-level selection
 166     fCFManagerProtons->GetParticleContainer()->Fill(containerInput,kStepGenerated);
 167
 168     //check the Acceptance-level cuts
 169     if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
 170     //fill the container for Acceptance-level selection
 171     fCFManagerProtons->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
 172   }//loop over MC particles (protons)
 173   //__________________________________________________________//
 174   //loop on the MC event
 175   for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 176     AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);
 177
 178     //check the MC-level cuts
 179     if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 180
 181     containerInput[0] = Rapidity(mcPart->Px(),
 182                                  mcPart->Py(),
 183                                  mcPart->Pz());
 184     containerInput[1] = (Float_t)mcPart->Pt();
 185     //fill the container for Gen-level selection
 186     fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepGenerated);
 187
 188     //check the Acceptance-level cuts
 189     if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
 190     //fill the container for Acceptance-level selection
 191     fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
 192   }//loop over MC particles (antiprotons)
 193
 194   //__________________________________________________________//
 195   //ESD track loop
 196   for (Int_t iTrack = 0; iTrack<fEvent->GetNumberOfTracks(); iTrack++) {
 197     track = fEvent->GetTrack(iTrack);
 198
 199     if (fReadTPCTracks) {
 200       if (fReadAODData) {
 201         Error("UserExec","TPC-only tracks are not supported with AOD");
 202         return ;
 203       }
 204       AliESDtrack* esdTrack    = (AliESDtrack*) track;
 205       AliESDtrack* esdTrackTPC = new AliESDtrack();
 206       if (!esdTrack->FillTPCOnlyTrack(*esdTrackTPC)) {
 207         Error("UserExec","Could not retrieve TPC info");
 208         continue;
 209       }
 210       track = esdTrackTPC ;
 211     }
 212
 213     if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartRecCuts,track)) continue;
 214
 215     // is track associated to particle ?
 216     Int_t label = track->GetLabel();
 217     if (label<0) continue;
 218     AliMCParticle *mcPart  = fMCEvent->GetTrack(label);
 219
 220     // check if this track was part of the signal
 221     if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 222
 223     //fill the container
 224     containerInput[0] = Rapidity(track->Px(),
 225                                  track->Py(),
 226                                  track->Pz());
 227     containerInput[1] = track->Pt();
 228     fCFManagerProtons->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
 229
 230     if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartSelCuts,track)) continue ;
 231     fCFManagerProtons->GetParticleContainer()->Fill(containerInput,kStepSelected);
 232
 233     if (fReadTPCTracks) delete track;
 234   }
 235   //__________________________________________________________//
 236   //ESD track loop
 237   for (Int_t iTrack = 0; iTrack<fEvent->GetNumberOfTracks(); iTrack++) {
 238     track = fEvent->GetTrack(iTrack);
 239
 240     if (fReadTPCTracks) {
 241       if (fReadAODData) {
 242         Error("UserExec","TPC-only tracks are not supported with AOD");
 243         return ;
 244       }
 245       AliESDtrack* esdTrack    = (AliESDtrack*) track;
 246       AliESDtrack* esdTrackTPC = new AliESDtrack();
 247       if (!esdTrack->FillTPCOnlyTrack(*esdTrackTPC)) {
 248         Error("UserExec","Could not retrieve TPC info");
 249         continue;
 250       }
 251       track = esdTrackTPC ;
 252     }
 253
 254     if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartRecCuts,track)) continue;
 255
 256     // is track associated to particle ?
 257     Int_t label = track->GetLabel();
 258     if (label<0) continue;
 259     AliMCParticle *mcPart  = fMCEvent->GetTrack(label);
 260
 261     // check if this track was part of the signal
 262     if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 263
 264     //fill the container
 265     containerInput[0] = Rapidity(track->Px(),
 266                                  track->Py(),
 267                                  track->Pz());
 268     containerInput[1] = track->Pt();
 269     fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
 270
 271     if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartSelCuts,track)) continue ;
 272     fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepSelected);
 273
 274     if (fReadTPCTracks) delete track;
 275   }
 276
 277   fHistEventsProcessed->Fill(0);
 278
 279   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 280   PostData(1,fHistEventsProcessed) ;
 281   PostData(2,fCFManagerProtons->GetParticleContainer()) ;
 282   PostData(3,fCFManagerAntiProtons->GetParticleContainer()) ;
 283   PostData(4,fQAHistList) ;
 284 }
 285
 286
 287 //___________________________________________________________________________
 288 void AliProtonCorrectionTask::Terminate(Option_t*) {
 289   // The Terminate() function is the last function to be called during
 290   // a query. It always runs on the client, it can be used to present
 291   // the results graphically or save the results to file.
 292
 293   Info("Terminate","");
 294   AliAnalysisTaskSE::Terminate();
 295
 296   //draw some example plots....
 297
 298   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 299
 300   TH1D* h00 =   cont->ShowProjection(0,0) ;
 301   TH1D* h01 =   cont->ShowProjection(0,1) ;
 302   TH1D* h02 =   cont->ShowProjection(0,2) ;
 303   TH1D* h03 =   cont->ShowProjection(0,3) ;
 304
 305   TH1D* h10 =   cont->ShowProjection(1,0) ;
 306   TH1D* h11 =   cont->ShowProjection(1,1) ;
 307   TH1D* h12 =   cont->ShowProjection(1,2) ;
 308   TH1D* h13 =   cont->ShowProjection(1,3) ;
 309
 310   TCanvas * c =new TCanvas("c","",1400,800);
 311   c->Divide(4,2);
 312
 313   c->cd(1); h00->Draw("p");
 314   c->cd(2); h01->Draw("p");
 315   c->cd(3); h02->Draw("p");
 316   c->cd(4); h03->Draw("p");
 317   c->cd(5); h10->Draw("p");
 318   c->cd(6); h11->Draw("p");
 319   c->cd(7); h12->Draw("p");
 320   c->cd(8); h13->Draw("p");
 321
 322   c->SaveAs("plots.eps");
 323 }
 324
 325
 326 //___________________________________________________________________________
 327 void AliProtonCorrectionTask::UserCreateOutputObjects() {
 328   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 329   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 330   //
 331   Info("CreateOutputObjects","CreateOutputObjects of task %s", GetName());
 332
 333   //slot #1
 334   OpenFile(1);
 335   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
 336
 337 //   OpenFile(2);
 338 //   OpenFile(3);
 339 }
 340
 341 //___________________________________________________________________________
 342 Double_t AliProtonCorrectionTask::Rapidity(Double_t Px, Double_t Py, Double_t Pz) {
 343   //returns the rapidity of the (anti)proton
 344   Double_t fMass = 9.38270000000000048e-01;
 345
 346   Double_t P = TMath::Sqrt(TMath::Power(Px,2) +
 347                            TMath::Power(Py,2) +
 348                            TMath::Power(Pz,2));
 349   Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
 350   Double_t y = -999;
 351   if(energy != Pz)
 352     y = 0.5*TMath::Log((energy + Pz)/(energy - Pz));
 353
 354   return y;
 355 }
 356
 357 #endif