CORRFW/test/AliCFRsnTask.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //-----------------------------------------------------------------------
  18 // Author : R. Vernet, Consorzio Cometa - Catania (it)
  19 //-----------------------------------------------------------------------
  20
  21
  22 #include "AliCFRsnTask.h"
  23 #include "TCanvas.h"
  24 #include "AliStack.h"
  25 #include "TParticle.h"
  26 #include "TH1I.h"
  27 #include "TChain.h"
  28 #include "AliMCEvent.h"
  29 #include "AliESDEvent.h"
  30 #include "AliCFManager.h"
  31 #include "AliCFCutBase.h"
  32 #include "AliCFContainer.h"
  33 #include "AliESDtrack.h"
  34 #include "AliLog.h"
  35 #include "AliRsnDaughter.h"
  36 #include "AliCFPair.h"
  37 #include "AliRsnMCInfo.h"
  38 #include "AliRsnPairParticle.h"
  39 #include "AliAODEvent.h"
  40 #include "AliRsnReader.h"
  41
  42 //__________________________________________________________________________
  43 AliCFRsnTask::AliCFRsnTask() :
  44   AliAnalysisTaskSE(),
  45   fRsnPDG(0),
  46   fCFManager(0x0),
  47   fHistEventsProcessed(0x0)
  48 {
  49   //
  50   //Default ctor
  51   //
  52 }
  53 //___________________________________________________________________________
  54 AliCFRsnTask::AliCFRsnTask(const Char_t* name) :
  55   AliAnalysisTaskSE(name),
  56   fRsnPDG(0),
  57   fCFManager(0x0),
  58   fHistEventsProcessed(0x0)
  59 {
  60   //
  61   // Constructor. Initialization of Inputs and Outputs
  62   //
  63   Info("AliCFRsnTask","Calling Constructor");
  64   /*
  65     DefineInput(0) and DefineOutput(0)
  66     are taken care of by AliAnalysisTaskSE constructor
  67   */
  68   DefineOutput(1,TH1I::Class());
  69   DefineOutput(2,AliCFContainer::Class());
  70   //   DefineOutput(3,TList::Class());
  71 }
  72
  73 //___________________________________________________________________________
  74 AliCFRsnTask& AliCFRsnTask::operator=(const AliCFRsnTask& c)
  75 {
  76   //
  77   // Assignment operator
  78   //
  79   if (this!=&c) {
  80     AliAnalysisTaskSE::operator=(c) ;
  81     fRsnPDG     = c.fRsnPDG;
  82     fCFManager  = c.fCFManager;
  83     fHistEventsProcessed = c.fHistEventsProcessed;
  84   }
  85   return *this;
  86 }
  87
  88 //___________________________________________________________________________
  89 AliCFRsnTask::AliCFRsnTask(const AliCFRsnTask& c) :
  90   AliAnalysisTaskSE(c),
  91   fRsnPDG(c.fRsnPDG),
  92   fCFManager(c.fCFManager),
  93   fHistEventsProcessed(c.fHistEventsProcessed)
  94 {
  95   //
  96   // Copy Constructor
  97   //
  98 }
  99
 100 //___________________________________________________________________________
 101 AliCFRsnTask::~AliCFRsnTask() {
 102   //
 103   //destructor
 104   //
 105   Info("~AliCFRsnTask","Calling Destructor");
 106   if (fCFManager)           delete fCFManager ;
 107   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 108 }
 109
 110 //_________________________________________________
 111 void AliCFRsnTask::UserExec(Option_t *)
 112 {
 113   //
 114   // Main loop function
 115   //
 116   Info("UserExec","") ;
 117
 118   if (!fInputEvent) {
 119     Error("UserExec","NO EVENT FOUND!");
 120     return;
 121   }
 122
 123   if (!fMCEvent) Error("UserExec","NO MC INFO FOUND!");
 124   fCFManager->SetEventInfo(fMCEvent);
 125
 126   AliStack*   stack   = fMCEvent->Stack();
 127
 128   Bool_t isESDEvent = strcmp(fInputEvent->ClassName(),"AliESDEvent") == 0 ? kTRUE : kFALSE ;
 129   Bool_t isAODEvent = strcmp(fInputEvent->ClassName(),"AliAODEvent") == 0 ? kTRUE : kFALSE ;
 130
 131   // MC-event selection
 132   Double_t containerInput[2] ;
 133
 134   //loop on the MC event
 135   Info("UserExec","Looping on MC event");
 136   for (Int_t ipart=0; ipart<stack->GetNprimary(); ipart++) {
 137     AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);
 138
 139     //check the MC-level cuts
 140     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 141     containerInput[0] = mcPart->Pt();
 142     containerInput[1] = mcPart->Y() ;
 143     //fill the container for Gen-level selection
 144     fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
 145
 146     //check the Acceptance-level cuts
 147     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
 148     //fill the container for Acceptance-level selection
 149     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
 150   }
 151
 152
 153   //Now go to rec level
 154   Info("UserExec","Looping on %s",fInputEvent->ClassName());
 155
 156   // Loop on negative tracks
 157   for (Int_t iTrack1 = 0; iTrack1<fInputEvent->GetNumberOfTracks(); iTrack1++) {
 158     AliVParticle* track1 = fInputEvent->GetTrack(iTrack1);
 159     //track1 is negative
 160     if (track1->Charge()>=0) continue;
 161     Int_t label1 = track1->GetLabel();
 162     if (label1<0) continue;
 163
 164     //Loop on positive tracks
 165     for (Int_t iTrack2 = 0; iTrack2<fInputEvent->GetNumberOfTracks(); iTrack2++) {
 166       AliVParticle* track2 = fInputEvent->GetTrack(iTrack2);
 167       //track2 is positive
 168       if (track2->Charge()<=0) continue;
 169       Int_t label2 = track2->GetLabel();
 170       if (label2<0) continue;
 171
 172       AliRsnDaughter daughter1;
 173       AliRsnDaughter daughter2;
 174
 175       AliRsnReader reader ; //needed for conversion ESD/AOD->AliRsnDaughter
 176
 177       if (isESDEvent) {
 178         if (!reader.ConvertTrack(&daughter1,(AliESDtrack*)track1)) Error("UserExec","Could not convert ESD track") ;
 179         if (!reader.ConvertTrack(&daughter2,(AliESDtrack*)track2)) Error("UserExec","Could not convert ESD track") ;
 180       }
 181       else if (isAODEvent) {
 182         if (!reader.ConvertTrack(&daughter1,(AliAODTrack*)track1)) Error("UserExec","Could not convert AOD track") ;
 183         if (!reader.ConvertTrack(&daughter2,(AliAODTrack*)track2)) Error("UserExec","Could not convert AOD track") ;
 184       }
 185       else {
 186         Error("UserExec","Error: input data file is not an ESD nor an AOD");
 187         return;
 188       }
 189
 190       AliCFPair pair(track1,track2); // This object is used for cuts
 191                                      // (to be replaced when AliRsnPairParticle
 192                                      // inherits from AliVParticle)
 193
 194       //Set MC PDG information to resonance daughters
 195       TParticle *part1 = stack->Particle(label1);
 196       daughter1.InitMCInfo(part1);
 197       daughter1.GetMCInfo()->SetPDG(part1->GetPdgCode());
 198       daughter1.SetM(part1->GetCalcMass()); // assign true mass
 199
 200       Int_t mother1 = part1->GetFirstMother();
 201       daughter1.GetMCInfo()->SetMother(mother1);
 202       if (mother1 >= 0) {
 203         TParticle *mum = stack->Particle(mother1);
 204         daughter1.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
 205       }
 206
 207       TParticle *part2 = stack->Particle(label2);
 208       daughter2.InitMCInfo(part2);
 209       daughter2.GetMCInfo()->SetPDG(part2->GetPdgCode());
 210       daughter2.SetM(part2->GetCalcMass()); // assign true mass
 211
 212       Int_t mother2 = part2->GetFirstMother();
 213       daughter2.GetMCInfo()->SetMother(mother2);
 214       if (mother2 >= 0) {
 215         TParticle *mum = stack->Particle(mother2);
 216         daughter2.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
 217       }
 218
 219       //make a mother resonance from the 2 candidate daughters
 220       AliRsnPairParticle rsn;
 221       rsn.SetPair(&daughter1,&daughter2);
 222
 223       //check if true resonance
 224       if (!rsn.IsTruePair(fRsnPDG)) continue;
 225       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,&pair)) continue;
 226
 227       //check if associated MC resonance passes the cuts
 228       Int_t motherLabel = rsn.GetDaughter(0)->GetMCInfo()->Mother() ;
 229       if (motherLabel<0) continue ;
 230
 231       AliMCParticle* mcRsn = fMCEvent->GetTrack(motherLabel);
 232       if (!mcRsn) continue;
 233       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcRsn)) continue;
 234
 235       // fill the container
 236       Double_t rsnEnergy = rsn.GetDaughter(0)->E()  + rsn.GetDaughter(1)->E()  ;
 237       Double_t rsnPz     = rsn.GetDaughter(0)->Pz() + rsn.GetDaughter(1)->Pz() ;
 238
 239       containerInput[0] = rsn.GetPt() ;
 240       containerInput[1] = GetRapidity(rsnEnergy,rsnPz);
 241       fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
 242
 243       if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,&pair)) continue ;
 244       fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
 245     }
 246   }
 247
 248   fHistEventsProcessed->Fill(0);
 249   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 250   PostData(1,fHistEventsProcessed) ;
 251   PostData(2,fCFManager->GetParticleContainer()) ;
 252
 253   //   TList * list = new TList();
 254   //   fCFManager->AddQAHistosToList(list);
 255   //   PostData(2,list) ;
 256 }
 257
 258
 259 //___________________________________________________________________________
 260 void AliCFRsnTask::Terminate(Option_t*)
 261 {
 262   // The Terminate() function is the last function to be called during
 263   // a query. It always runs on the client, it can be used to present
 264   // the results graphically or save the results to file.
 265
 266   Info("Terminate","");
 267   AliAnalysisTaskSE::Terminate();
 268
 269
 270   //draw some example plots....
 271
 272   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 273
 274   TH1D* h00 =   cont->ShowProjection(0,0) ;
 275   TH1D* h01 =   cont->ShowProjection(0,1) ;
 276   TH1D* h02 =   cont->ShowProjection(0,2) ;
 277   TH1D* h03 =   cont->ShowProjection(0,3) ;
 278
 279   TH1D* h10 =   cont->ShowProjection(1,0) ;
 280   TH1D* h11 =   cont->ShowProjection(1,1) ;
 281   TH1D* h12 =   cont->ShowProjection(1,2) ;
 282   TH1D* h13 =   cont->ShowProjection(1,3) ;
 283
 284   Double_t max1 = h00->GetMaximum();
 285   Double_t max2 = h10->GetMaximum();
 286
 287   h00->GetYaxis()->SetRangeUser(0,max1*1.2);
 288   h01->GetYaxis()->SetRangeUser(0,max1*1.2);
 289   h02->GetYaxis()->SetRangeUser(0,max1*1.2);
 290   h03->GetYaxis()->SetRangeUser(0,max1*1.2);
 291
 292   h10->GetYaxis()->SetRangeUser(0,max2*1.2);
 293   h11->GetYaxis()->SetRangeUser(0,max2*1.2);
 294   h12->GetYaxis()->SetRangeUser(0,max2*1.2);
 295   h13->GetYaxis()->SetRangeUser(0,max2*1.2);
 296
 297   h00->SetMarkerStyle(23) ;
 298   h01->SetMarkerStyle(24) ;
 299   h02->SetMarkerStyle(25) ;
 300   h03->SetMarkerStyle(26) ;
 301
 302   h10->SetMarkerStyle(23) ;
 303   h11->SetMarkerStyle(24) ;
 304   h12->SetMarkerStyle(25) ;
 305   h13->SetMarkerStyle(26) ;
 306
 307   TCanvas * c =new TCanvas("c","",1400,800);
 308   c->Divide(4,2);
 309
 310   c->cd(1);
 311   h00->Draw("p");
 312   c->cd(2);
 313   h01->Draw("p");
 314   c->cd(3);
 315   h02->Draw("p");
 316   c->cd(4);
 317   h03->Draw("p");
 318   c->cd(5);
 319   h10->Draw("p");
 320   c->cd(6);
 321   h11->Draw("p");
 322   c->cd(7);
 323   h12->Draw("p");
 324   c->cd(8);
 325   h13->Draw("p");
 326
 327   c->SaveAs("plots.eps");
 328 }
 329
 330 //___________________________________________________________________________
 331 void AliCFRsnTask::UserCreateOutputObjects() {
 332   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 333   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 334   //
 335   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 336
 337   //slot #1
 338   OpenFile(1);
 339   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
 340 }
 341
 342 //___________________________________________________________________________
 343 Double_t AliCFRsnTask::GetRapidity(Double_t energy, Double_t pz) {
 344   if (energy == pz || energy == -pz) {
 345     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 346     return 999;
 347   }
 348   if (energy < pz) {
 349     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
 350     return 999;
 351   }
 352   Double_t y = 0.5 * log((energy + pz) / (energy - pz));
 353   return y;
 354 }
 355
 356