CORRFW/test/AliCFRsnTask.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 // Author : R. Vernet, Consorzio Cometa - Catania (it)
  18 //-----------------------------------------------------------------------
  19
  20
  21 #include <TROOT.h>
  22 #include <TInterpreter.h>
  23
  24 #include "AliCFRsnTask.h"
  25 #include "TCanvas.h"
  26 #include "AliStack.h"
  27 #include "TParticle.h"
  28 #include "TH1I.h"
  29 #include "TChain.h"
  30 #include "AliMCEventHandler.h"
  31 #include "AliMCEvent.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliESDEvent.h"
  34 #include "AliCFManager.h"
  35 #include "AliCFCutBase.h"
  36 #include "AliCFContainer.h"
  37 #include "TChain.h"
  38 #include "AliESDtrack.h"
  39 #include "AliLog.h"
  40 #include "AliRsnDaughter.h"
  41 #include "AliCFPair.h"
  42 #include <memory>
  43
  44 //__________________________________________________________________________
  45 AliCFRsnTask::AliCFRsnTask() :
  46   fRsnPDG(313),
  47   fChain(0x0),
  48   fESD(0x0),
  49   fCFManager(0x0),
  50   fHistEventsProcessed(0x0)
  51 {
  52   //Defual ctor
  53 }
  54 //___________________________________________________________________________
  55 AliCFRsnTask::AliCFRsnTask(const Char_t* name) :
  56   AliAnalysisTask(name,"AliCFRsnTask"),
  57   fRsnPDG(313),
  58   fChain(0x0),
  59   fESD(0x0),
  60   fCFManager(0x0),
  61   fHistEventsProcessed(0x0)
  62 {
  63   //
  64   // Constructor. Initialization of Inputs and Outputs
  65   //
  66   Info("AliCFRsnTask","Calling Constructor");
  67   DefineInput (0,TChain::Class());
  68   DefineOutput(0,TH1I::Class());
  69   DefineOutput(1,AliCFContainer::Class());
  70   //   DefineOutput(2,TList::Class());
  71 }
  72
  73 //___________________________________________________________________________
  74 AliCFRsnTask& AliCFRsnTask::operator=(const AliCFRsnTask& c)
  75 {
  76   //
  77   // Assignment operator
  78   //
  79   if (this!=&c) {
  80     AliAnalysisTask::operator=(c) ;
  81     fRsnPDG     = c.fRsnPDG;
  82     fChain      = c.fChain;
  83     fESD        = c.fESD;
  84     fCFManager  = c.fCFManager;
  85     fHistEventsProcessed = c.fHistEventsProcessed;
  86   }
  87   return *this;
  88 }
  89
  90 //___________________________________________________________________________
  91 AliCFRsnTask::AliCFRsnTask(const AliCFRsnTask& c) :
  92   AliAnalysisTask(c),
  93   fRsnPDG(c.fRsnPDG),
  94   fChain(c.fChain),
  95   fESD(c.fESD),
  96   fCFManager(c.fCFManager),
  97   fHistEventsProcessed(c.fHistEventsProcessed)
  98 {
  99   //
 100   // Copy Constructor
 101   //
 102 }
 103
 104 //___________________________________________________________________________
 105 AliCFRsnTask::~AliCFRsnTask() {
 106   //
 107   //destructor
 108   //
 109   Info("~AliCFRsnTask","Calling Destructor");
 110   if (fChain)               delete fChain ;
 111   if (fESD)                 delete fESD ;
 112   if (fCFManager)           delete fCFManager ;
 113   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 114 }
 115
 116 //___________________________________________________________________________
 117
 118 void AliCFRsnTask::Init()
 119 {
 120
 121 }
 122 //_________________________________________________
 123 void AliCFRsnTask::Exec(Option_t *)
 124 {
 125   //
 126   // Main loop function
 127   //
 128   Info("Exec","") ;
 129   // Get the mc truth
 130   AliMCEventHandler* mcTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 131
 132   if (!mcTruth) Error("Exec","NO MC INFO FOUND... EXITING");
 133
 134   // transform possible old AliESD into AliESDEvent
 135
 136   if (fESD->GetAliESDOld()) fESD->CopyFromOldESD(); //transition to new ESD format
 137
 138   //pass the MC evt handler to the cuts that need it
 139
 140   fCFManager->SetEventInfo(mcTruth);
 141
 142   // Get the MC event
 143   AliMCEvent* mcEvent = mcTruth->MCEvent();
 144   AliStack*   stack   = mcEvent->Stack();
 145
 146   // MC-event selection
 147   Double_t containerInput[2] ;
 148
 149   //loop on the MC event
 150   Info("Exec","Looping on MC event");
 151   for (Int_t ipart=0; ipart<stack->GetNprimary(); ipart++) {
 152     AliMCParticle *mcPart  = mcEvent->GetTrack(ipart);
 153
 154     //check the MC-level cuts
 155     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 156     containerInput[0] = mcPart->Pt();
 157     containerInput[1] = mcPart->Eta() ;
 158     //fill the container for Gen-level selection
 159     fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
 160
 161     //check the Acceptance-level cuts
 162     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
 163     //fill the container for Acceptance-level selection
 164     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
 165   }
 166
 167
 168   //Now go to rec level
 169   Info("Exec","Looping on ESD event");
 170
 171   //SET THE ESD AS EVENT INFO IN RECONSTRUCTION CUTS
 172   TObjArray* fCutsReco = fCFManager->GetParticleCutsList(AliCFManager::kPartRecCuts);
 173   TObjArray* fCutsPID = fCFManager->GetParticleCutsList(AliCFManager::kPartSelCuts);
 174   TObjArrayIter iter1(fCutsReco);
 175   TObjArrayIter iter2(fCutsPID);
 176   AliCFCutBase *cut = 0;
 177   while ( (cut = (AliCFCutBase*)iter1.Next()) ) {
 178     cut->SetEvtInfo(fESD);
 179   }
 180   while ( (cut = (AliCFCutBase*)iter2.Next()) ) {
 181     cut->SetEvtInfo(fESD);
 182   }
 183
 184   for (Int_t iTrack1 = 0; iTrack1<fESD->GetNumberOfTracks(); iTrack1++) {
 185     AliESDtrack* esdTrack1 = fESD->GetTrack(iTrack1);
 186     //track1 is negative
 187     if (esdTrack1->Charge()>=0) continue;
 188     Int_t esdLabel1 = esdTrack1->GetLabel();
 189     if (esdLabel1<0) continue;
 190
 191     for (Int_t iTrack2 = 0; iTrack2<fESD->GetNumberOfTracks(); iTrack2++) {
 192       AliESDtrack* esdTrack2 = fESD->GetTrack(iTrack2);
 193       //track2 is positive
 194       if (esdTrack2->Charge()<=0) continue;
 195       Int_t esdLabel2 = esdTrack2->GetLabel();
 196       if (esdLabel2<0) continue;
 197
 198       // copy ESDtrack data into RsnDaughter (and make Bayesian PID)
 199       //if problem with copy constructor, use the following special command to destroy the pointer when exiting scope
 200       //std::auto_ptr<AliRsnDaughter> track1(AliRsnDaughter::Adopt(esdTrack1,iTrack1));
 201
 202       AliRsnDaughter* tmp1=AliRsnDaughter::Adopt(esdTrack1,iTrack1);
 203       AliRsnDaughter* tmp2=AliRsnDaughter::Adopt(esdTrack2,iTrack2);
 204       AliRsnDaughter track1(*tmp1);
 205       AliRsnDaughter track2(*tmp2);
 206       delete tmp1;
 207       delete tmp2;
 208
 209       TParticle *part1 = stack->Particle(esdLabel1);
 210       track1.SetTruePDG(part1->GetPdgCode());
 211       Int_t mother1 = part1->GetFirstMother();
 212       track1.SetMother(mother1);
 213       if (mother1 >= 0) {
 214         TParticle *mum = stack->Particle(mother1);
 215         track1.SetMotherPDG(mum->GetPdgCode());
 216       }
 217
 218       TParticle *part2 = stack->Particle(esdLabel2);
 219       track2.SetTruePDG(part2->GetPdgCode());
 220       Int_t mother2 = part2->GetFirstMother();
 221       track2.SetMother(mother2);
 222       if (mother2 >= 0) {
 223         TParticle *mum = stack->Particle(mother2);
 224         track2.SetMotherPDG(mum->GetPdgCode());
 225       }
 226
 227       //make a mother resonance from the 2 candidate daughters
 228       AliRsnDaughter rsn = AliRsnDaughter::Sum(track1,track2) ;
 229       AliCFPair pair(esdTrack1,esdTrack2);
 230
 231       //check if true resonance
 232       if (rsn.GetMotherPDG() != fRsnPDG) continue;
 233       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,&pair)) continue;
 234
 235       //check if associated MC resonance passes the cuts
 236       Int_t motherLabel=rsn.GetLabel();
 237       if (motherLabel<0) continue ;
 238       AliMCParticle* mcRsn = mcEvent->GetTrack(motherLabel);
 239       if (!mcRsn) continue;
 240       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcRsn)) continue;
 241
 242       //fill the container
 243       containerInput[0] = rsn.GetPt() ;
 244       containerInput[1] = GetRapidity(rsn.GetEnergy(),rsn.GetPz());
 245       fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
 246
 247       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,&pair)) continue ;
 248       fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
 249     }
 250   }
 251
 252   fHistEventsProcessed->Fill(0);
 253   PostData(0,fHistEventsProcessed) ;
 254   PostData(1,fCFManager->GetParticleContainer()) ;
 255
 256   //   TList * list = new TList();
 257   //   fCFManager->AddQAHistosToList(list);
 258   //   PostData(2,list) ;
 259 }
 260
 261
 262 //___________________________________________________________________________
 263 void AliCFRsnTask::Terminate(Option_t*)
 264 {
 265   // The Terminate() function is the last function to be called during
 266   // a query. It always runs on the client, it can be used to present
 267   // the results graphically or save the results to file.
 268
 269   Info("Terminate","");
 270   AliAnalysisTask::Terminate();
 271
 272
 273   Double_t max1 = fCFManager->GetParticleContainer()->ShowProjection(0,0)->GetMaximum();
 274   Double_t max2 = fCFManager->GetParticleContainer()->ShowProjection(1,0)->GetMaximum();
 275
 276   fCFManager->GetParticleContainer()->ShowProjection(0,0)->GetYaxis()->SetRangeUser(0,max1*1.2);
 277   fCFManager->GetParticleContainer()->ShowProjection(0,1)->GetYaxis()->SetRangeUser(0,max1*1.2);
 278   fCFManager->GetParticleContainer()->ShowProjection(0,2)->GetYaxis()->SetRangeUser(0,max1*1.2);
 279   fCFManager->GetParticleContainer()->ShowProjection(0,3)->GetYaxis()->SetRangeUser(0,max1*1.2);
 280
 281   fCFManager->GetParticleContainer()->ShowProjection(1,0)->GetYaxis()->SetRangeUser(0,max2*1.2);
 282   fCFManager->GetParticleContainer()->ShowProjection(1,1)->GetYaxis()->SetRangeUser(0,max2*1.2);
 283   fCFManager->GetParticleContainer()->ShowProjection(1,2)->GetYaxis()->SetRangeUser(0,max2*1.2);
 284   fCFManager->GetParticleContainer()->ShowProjection(1,3)->GetYaxis()->SetRangeUser(0,max2*1.2);
 285
 286   fCFManager->GetParticleContainer()->ShowProjection(0,0)->SetMarkerStyle(23) ;
 287   fCFManager->GetParticleContainer()->ShowProjection(0,1)->SetMarkerStyle(24) ;
 288   fCFManager->GetParticleContainer()->ShowProjection(0,2)->SetMarkerStyle(25) ;
 289   fCFManager->GetParticleContainer()->ShowProjection(0,3)->SetMarkerStyle(26) ;
 290
 291   fCFManager->GetParticleContainer()->ShowProjection(1,0)->SetMarkerStyle(23) ;
 292   fCFManager->GetParticleContainer()->ShowProjection(1,1)->SetMarkerStyle(24) ;
 293   fCFManager->GetParticleContainer()->ShowProjection(1,2)->SetMarkerStyle(25) ;
 294   fCFManager->GetParticleContainer()->ShowProjection(1,3)->SetMarkerStyle(26) ;
 295
 296   TCanvas * c =new TCanvas("c","",1400,800);
 297   c->Divide(4,2);
 298
 299   c->cd(1);
 300   fCFManager->GetParticleContainer()->ShowProjection(0,0)->Draw("p");
 301   c->cd(2);
 302   fCFManager->GetParticleContainer()->ShowProjection(0,1)->Draw("p");
 303   c->cd(3);
 304   fCFManager->GetParticleContainer()->ShowProjection(0,2)->Draw("p");
 305   c->cd(4);
 306   fCFManager->GetParticleContainer()->ShowProjection(0,3)->Draw("p");
 307   c->cd(5);
 308   fCFManager->GetParticleContainer()->ShowProjection(1,0)->Draw("p");
 309   c->cd(6);
 310   fCFManager->GetParticleContainer()->ShowProjection(1,1)->Draw("p");
 311   c->cd(7);
 312   fCFManager->GetParticleContainer()->ShowProjection(1,2)->Draw("p");
 313   c->cd(8);
 314   fCFManager->GetParticleContainer()->ShowProjection(1,3)->Draw("p");
 315
 316   c->SaveAs("plots.eps");
 317
 318   delete fHistEventsProcessed ;
 319 }
 320
 321 //___________________________________________________________________________
 322 void AliCFRsnTask::ConnectInputData(Option_t *) {
 323   //
 324   // Initialize branches.
 325   //
 326   Info("ConnectInputData","ConnectInputData of task %s\n",GetName());
 327
 328   fChain = (TChain*)GetInputData(0);
 329   fChain->SetBranchStatus("*FMD*",0);
 330   fChain->SetBranchStatus("*CaloClusters*",0);
 331   fESD = new AliESDEvent();
 332   fESD->ReadFromTree(fChain);
 333 }
 334
 335 //___________________________________________________________________________
 336 void AliCFRsnTask::CreateOutputObjects() {
 337   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 338   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 339   //
 340   Info("CreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 341
 342   //slot #0
 343   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
 344 }
 345
 346 //___________________________________________________________________________
 347 Double_t AliCFRsnTask::GetRapidity(Double_t energy, Double_t pz) {
 348   if (energy == pz || energy == -pz) {
 349     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 350     return 999;
 351   }
 352   if (energy < pz) {
 353     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 354     return 999;
 355   }
 356   Double_t y = 0.5 * log((energy + pz) / (energy - pz));
 357   return y;
 358 }
 359
 360