CORRFW/test/AliCFV0Task.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 on AliEn (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 ALICFV0TASK_CXX
  28 #define ALICFV0TASK_CXX
  29
  30 #include "AliCFV0Task.h"
  31 #include "TCanvas.h"
  32 #include "AliStack.h"
  33 #include "TParticle.h"
  34 #include "TH1I.h"
  35 #include "AliMCEvent.h"
  36 #include "AliCFManager.h"
  37 #include "AliCFContainer.h"
  38 #include "AliESDtrack.h"
  39 #include "AliLog.h"
  40 #include "AliESDv0.h"
  41 #include "AliV0vertexer.h"
  42 #include "AliCFPair.h"
  43 #include "AliESDEvent.h"
  44 #include "TChain.h"
  45 #include "AliCFParticleGenCuts.h"
  46
  47 //__________________________________________________________________________
  48 AliCFV0Task::AliCFV0Task() :
  49   AliAnalysisTaskSE(),
  50   fRebuildV0s(0),
  51   fV0PDG(310),
  52   fCFManager(0x0),
  53   fHistEventsProcessed(0x0)
  54 {
  55   //
  56   //Default ctor
  57   //
  58 }
  59 //___________________________________________________________________________
  60 AliCFV0Task::AliCFV0Task(const Char_t* name) :
  61   AliAnalysisTaskSE(name),
  62   fRebuildV0s(0),
  63   fV0PDG(310),
  64   fCFManager(0x0),
  65   fHistEventsProcessed(0x0)
  66 {
  67   //
  68   // Constructor. Initialization of Inputs and Outputs
  69   //
  70   Info("AliCFV0Task","Calling Constructor");
  71   /*
  72     DefineInput(0) and DefineOutput(0)
  73     are taken care of by AliAnalysisTaskSE constructor
  74   */
  75   DefineOutput(1,TH1I::Class());
  76   DefineOutput(2,AliCFContainer::Class());
  77 //   DefineOutput(3,TList::Class());
  78 }
  79
  80 //___________________________________________________________________________
  81 AliCFV0Task& AliCFV0Task::operator=(const AliCFV0Task& c)
  82 {
  83   //
  84   // Assignment operator
  85   //
  86   if (this!=&c) {
  87     AliAnalysisTaskSE::operator=(c) ;
  88     fRebuildV0s = c.fRebuildV0s;
  89     fV0PDG      = c.fV0PDG;
  90     fCFManager  = c.fCFManager;
  91     fHistEventsProcessed = c.fHistEventsProcessed;
  92   }
  93   return *this;
  94 }
  95
  96 //___________________________________________________________________________
  97 AliCFV0Task::AliCFV0Task(const AliCFV0Task& c) :
  98   AliAnalysisTaskSE(c),
  99   fRebuildV0s(c.fRebuildV0s),
 100   fV0PDG(c.fV0PDG),
 101   fCFManager(c.fCFManager),
 102   fHistEventsProcessed(c.fHistEventsProcessed)
 103 {
 104   //
 105   // Copy Constructor
 106   //
 107 }
 108
 109 //___________________________________________________________________________
 110 AliCFV0Task::~AliCFV0Task() {
 111   //
 112   //destructor
 113   //
 114   Info("~AliCFV0Task","Calling Destructor");
 115   if (fCFManager)           delete fCFManager ;
 116   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 117 }
 118
 119 //_________________________________________________
 120 void AliCFV0Task::UserExec(Option_t *)
 121 {
 122   //
 123   // Main loop function
 124   //
 125   Info("UserExec","") ;
 126
 127   AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 128   if (!fESD) {
 129     Error("UserExec","NO ESD FOUND!");
 130     return;
 131   }
 132
 133   if (!fMCEvent) Error("UserExec","NO MC INFO FOUND!");
 134   fCFManager->SetEventInfo(fMCEvent);
 135
 136   // MC-event selection
 137   Double_t containerInput[2] ;
 138
 139   //loop on the MC event
 140   Info("UserExec","Looping on MC event");
 141   for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 142     AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);
 143
 144     //check the MC-level cuts
 145     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 146     containerInput[0] = mcPart->Pt();
 147     containerInput[1] = mcPart->Eta() ;
 148     //fill the container for Gen-level selection
 149     fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
 150
 151     //check the Acceptance-level cuts
 152     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
 153     //fill the container for Acceptance-level selection
 154     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
 155   }
 156
 157
 158   //Now go to rec level
 159   Info("UserExec","Looping on ESD event");
 160
 161   if (fRebuildV0s) RebuildV0s(fESD) ;
 162   Info("UserExec","There are %d V0s in event",fESD->GetNumberOfV0s());
 163   for (Int_t iV0 = 0; iV0<fESD->GetNumberOfV0s(); iV0++) {
 164
 165     AliESDv0* esdV0 = fESD->GetV0(iV0);
 166
 167     //check if mother reconstructed V0 can be associated to a MC V0
 168     Int_t labMCV0 = IsMcV0(esdV0) ;
 169     if (labMCV0 <0) continue;
 170
 171     esdV0->ChangeMassHypothesis(fV0PDG); //important to do that before entering the cut check
 172
 173     AliCFPair pair(esdV0,fESD);
 174     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,&pair)) continue;
 175
 176     //check if associated MC v0 passes the cuts
 177     AliMCParticle* mcV0 = fMCEvent->GetTrack(labMCV0);
 178     if (!mcV0) continue;
 179
 180     Info("UserExec","is v0 primary : %d",AliCFParticleGenCuts::IsPrimary(mcV0,fMCEvent->Stack()));
 181
 182     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcV0)) continue;
 183
 184     //fill the container
 185     Double_t mom[3];
 186     esdV0->GetPxPyPz(mom[0],mom[1],mom[2]);
 187     Double_t pt2 = mom[0]*mom[0]+mom[1]*mom[1] ;
 188     Double_t pt  = TMath::Sqrt(pt2);
 189     Double_t energy  = TMath::Sqrt(pt2 + mom[2]*mom[2] + TMath::Power(TDatabasePDG::Instance()->GetParticle(fV0PDG)->Mass(),2));
 190
 191     containerInput[0] = pt ;
 192     containerInput[1] = GetRapidity(energy,mom[2]);
 193     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
 194     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,&pair)) continue ;
 195     fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
 196   }
 197
 198   fHistEventsProcessed->Fill(0);
 199   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 200   PostData(1,fHistEventsProcessed) ;
 201   PostData(2,fCFManager->GetParticleContainer()) ;
 202 }
 203
 204
 205 //___________________________________________________________________________
 206 void AliCFV0Task::Terminate(Option_t*)
 207 {
 208   // The Terminate() function is the last function to be called during
 209   // a query. It always runs on the client, it can be used to present
 210   // the results graphically or save the results to file.
 211
 212   Info("Terminate","");
 213   AliAnalysisTaskSE::Terminate();
 214
 215
 216   //draw some example plots....
 217
 218   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 219
 220   TH1D* h00 =   cont->ShowProjection(0,0) ;
 221   TH1D* h01 =   cont->ShowProjection(0,1) ;
 222   TH1D* h02 =   cont->ShowProjection(0,2) ;
 223   TH1D* h03 =   cont->ShowProjection(0,3) ;
 224
 225   TH1D* h10 =   cont->ShowProjection(1,0) ;
 226   TH1D* h11 =   cont->ShowProjection(1,1) ;
 227   TH1D* h12 =   cont->ShowProjection(1,2) ;
 228   TH1D* h13 =   cont->ShowProjection(1,3) ;
 229
 230   Double_t max1 = h00->GetMaximum();
 231   Double_t max2 = h10->GetMaximum();
 232
 233   h00->GetYaxis()->SetRangeUser(0,max1*1.2);
 234   h01->GetYaxis()->SetRangeUser(0,max1*1.2);
 235   h02->GetYaxis()->SetRangeUser(0,max1*1.2);
 236   h03->GetYaxis()->SetRangeUser(0,max1*1.2);
 237
 238   h10->GetYaxis()->SetRangeUser(0,max2*1.2);
 239   h11->GetYaxis()->SetRangeUser(0,max2*1.2);
 240   h12->GetYaxis()->SetRangeUser(0,max2*1.2);
 241   h13->GetYaxis()->SetRangeUser(0,max2*1.2);
 242
 243   h00->SetMarkerStyle(23) ;
 244   h01->SetMarkerStyle(24) ;
 245   h02->SetMarkerStyle(25) ;
 246   h03->SetMarkerStyle(26) ;
 247
 248   h10->SetMarkerStyle(23) ;
 249   h11->SetMarkerStyle(24) ;
 250   h12->SetMarkerStyle(25) ;
 251   h13->SetMarkerStyle(26) ;
 252
 253   TCanvas * c =new TCanvas("c","",1400,800);
 254   c->Divide(4,2);
 255
 256   c->cd(1);
 257   h00->Draw("p");
 258   c->cd(2);
 259   h01->Draw("p");
 260   c->cd(3);
 261   h02->Draw("p");
 262   c->cd(4);
 263   h03->Draw("p");
 264   c->cd(5);
 265   h10->Draw("p");
 266   c->cd(6);
 267   h11->Draw("p");
 268   c->cd(7);
 269   h12->Draw("p");
 270   c->cd(8);
 271   h13->Draw("p");
 272
 273   c->SaveAs("plots.eps");
 274 }
 275
 276 //___________________________________________________________________________
 277 void AliCFV0Task::UserCreateOutputObjects() {
 278   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 279   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 280   //
 281   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 282
 283   //slot #1
 284   OpenFile(1);
 285   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
 286 }
 287
 288 //___________________________________________________________________________
 289 Int_t AliCFV0Task::GetV0Label(UInt_t lab1, UInt_t lab2) const {
 290   //
 291   // returns the label of the V0, given the labels of the 2 daughter tracks
 292   // returns -1 if the V0 is fake
 293   //
 294
 295   AliStack* stack = fMCEvent->Stack();
 296   TParticle* part1 = stack->Particle(lab1) ;
 297   TParticle* part2 = stack->Particle(lab2) ;
 298
 299   Int_t part1MotherLab=part1->GetFirstMother();
 300   Int_t part2MotherLab=part2->GetFirstMother();
 301
 302   if (part1MotherLab==-1 || part2MotherLab==-1) return -1 ;
 303   if (part1MotherLab != part2MotherLab )        return -1 ;
 304
 305   if (stack->Particle(part1MotherLab)->GetPdgCode() != fV0PDG ) return -1 ;
 306
 307   switch (fV0PDG) {
 308   case kK0Short :
 309     if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()==  211)  ||
 310          (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
 311     break ;
 312   case kLambda0 :
 313     if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()== 2212)  ||
 314          (part1->GetPdgCode()==  2212 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
 315     break ;
 316   case kLambda0Bar :
 317     if ( (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -2212) ||
 318          (part1->GetPdgCode()== -2212 && part2->GetPdgCode()==   211)  )   return part1MotherLab ;
 319     break ;
 320   default :
 321     return -1;
 322     break ;
 323   }
 324
 325   return -1 ;
 326 }
 327
 328 //___________________________________________________________________________
 329 Double_t AliCFV0Task::GetRapidity(Double_t energy, Double_t pz) {
 330   if (energy == pz || energy == -pz) {
 331     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 332     return 999;
 333   }
 334   if (energy < pz) {
 335     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 336     return 999;
 337   }
 338   Double_t y = 0.5 * log((energy + pz) / (energy - pz));
 339   return y;
 340 }
 341
 342 //___________________________________________________________________________
 343 Int_t AliCFV0Task::IsMcV0(AliESDv0* v0) const {
 344
 345   AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 346
 347   Int_t nindex = v0->GetNindex();
 348   Int_t pindex = v0->GetPindex();
 349   AliESDtrack *nTrack = fESD->GetTrack(nindex) ;
 350   AliESDtrack *pTrack = fESD->GetTrack(pindex) ;
 351
 352   if (!nTrack || !pTrack) return -1 ;
 353
 354   Int_t nlab  = nTrack->GetLabel() ;
 355   Int_t plab  = pTrack->GetLabel() ;
 356
 357   if (nlab <0 || plab <0) return -1 ;
 358
 359   return GetV0Label((UInt_t)nlab,(UInt_t)plab) ;
 360 }
 361
 362
 363 //___________________________________________________________________________
 364 void AliCFV0Task::RebuildV0s(AliESDEvent* fESD) {
 365
 366   fESD->ResetV0s();
 367
 368   //These are pp cuts : to change if Pb+Pb !
 369   Double_t cuts[]={33,  // max. allowed chi2
 370                    0.01,// min. allowed negative daughter's impact parameter
 371                    0.01,// min. allowed positive daughter's impact parameter
 372                    0.5,// max. allowed DCA between the daughter tracks
 373                    0.98,// max. allowed cosine of V0's pointing angle
 374                    0.2,  // min. radius of the fiducial volume
 375                    100.   // max. radius of the fiducial volume
 376   };
 377   //   AliV0vertexer* v0Vertexer = new AliV0vertexer(cuts);
 378   //   v0Vertexer->SetVertex(primVertexPos);
 379   AliV0vertexer* v0Vertexer = new AliV0vertexer();
 380   v0Vertexer->SetCuts(cuts) ;
 381   v0Vertexer->Tracks2V0vertices(fESD);
 382   delete v0Vertexer ;
 383 }
 384
 385 #endif