PWG/muon/AliCFMuonSingleTask1.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 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------
  19 // Example of task (running locally, on AliEn and CAF),
  20 // which provides standard way of calculating acceptance and efficiency
  21 // between different steps of the procedure.
  22 // The ouptut of the task is a AliCFContainer from which the efficiencies
  23 // can be calculated
  24 //-----------------------------------------------------------------------
  25 // Author : R. Vernet, Consorzio Cometa - Catania (it)
  26 //-----------------------------------------------------------------------
  27 // Modification done by X. Lopez - LPC Clermont (fr)
  28 //-----------------------------------------------------------------------
  29
  30 #ifndef ALICFMUONSINGLETASK1_CXX
  31 #define ALICFMUONSINGLETASK1_CXX
  32
  33 #include "AliCFMuonSingleTask1.h"
  34 #include "AliHeader.h"
  35 #include "AliESDHeader.h"
  36 #include "AliStack.h"
  37 #include "TParticle.h"
  38 #include "TLorentzVector.h"
  39 #include "TH1I.h"
  40 #include "AliMCEvent.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDEvent.h"
  43 #include "AliAODEvent.h"
  44 #include "AliCFManager.h"
  45 #include "AliCFCutBase.h"
  46 #include "AliCFContainer.h"
  47 #include "TChain.h"
  48 #include "AliESDtrack.h"
  49 #include "AliLog.h"
  50 #include "AliESDMuonTrack.h"
  51 #include "AliESDtrack.h"
  52 #include "AliESDInputHandler.h"
  53 #include "TCanvas.h"
  54
  55 ClassImp(AliCFMuonSingleTask1)
  56
  57 //__________________________________________________________________________
  58 AliCFMuonSingleTask1::AliCFMuonSingleTask1() :
  59   fReadAODData(0),
  60   fCFManager(0x0),
  61   fQAHistList(0x0),
  62   fHistEventsProcessed(0x0),
  63   fNevt(0),
  64   fIsMC(kFALSE)
  65 {
  66   //
  67   //Default ctor
  68   //
  69 }
  70 //___________________________________________________________________________
  71 AliCFMuonSingleTask1::AliCFMuonSingleTask1(const Char_t* name) :
  72   AliAnalysisTaskSE(name),
  73   fReadAODData(0),
  74   fCFManager(0x0),
  75   fQAHistList(0x0),
  76   fHistEventsProcessed(0x0),
  77   fNevt(0),
  78   fIsMC(kFALSE)
  79 {
  80   //
  81   // Constructor. Initialization of Inputs and Outputs
  82   //
  83   Info("AliCFMuonSingleTask1","Calling Constructor");
  84
  85   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
  86
  87   DefineOutput(1,TH1I::Class());
  88   DefineOutput(2,AliCFContainer::Class());
  89
  90 }
  91
  92 //___________________________________________________________________________
  93 AliCFMuonSingleTask1& AliCFMuonSingleTask1::operator=(const AliCFMuonSingleTask1& c)
  94 {
  95   //
  96   // Assignment operator
  97   //
  98   if (this!=&c) {
  99     AliAnalysisTaskSE::operator=(c) ;
 100     fReadAODData = c.fReadAODData ;
 101     fCFManager  = c.fCFManager;
 102     fQAHistList = c.fQAHistList ;
 103     fHistEventsProcessed = c.fHistEventsProcessed;
 104     fNevt = c.fNevt ;
 105   }
 106   return *this;
 107 }
 108
 109 //___________________________________________________________________________
 110 AliCFMuonSingleTask1::AliCFMuonSingleTask1(const AliCFMuonSingleTask1& c) :
 111   AliAnalysisTaskSE(c),
 112   fReadAODData(c.fReadAODData),
 113   fCFManager(c.fCFManager),
 114   fQAHistList(c.fQAHistList),
 115   fHistEventsProcessed(c.fHistEventsProcessed),
 116   fNevt(c.fNevt),
 117   fIsMC(kFALSE)
 118 {
 119   //
 120   // Copy Constructor
 121   //
 122 }
 123
 124 //___________________________________________________________________________
 125 AliCFMuonSingleTask1::~AliCFMuonSingleTask1() {
 126   //
 127   //destructor
 128   //
 129   Info("~AliCFMuonSingleTask1","Calling Destructor");
 130   if (fCFManager)           delete fCFManager ;
 131   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 132   if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
 133 }
 134
 135 //_________________________________________________
 136 void AliCFMuonSingleTask1::UserExec(Option_t *)
 137 {
 138   //
 139   // Main loop function
 140   //
 141
 142   Info("UserExec"," ") ;
 143   if (!fMCEvent && fIsMC) {
 144     Error("UserExec","NO MC EVENT FOUND!");
 145     return;
 146   }
 147
 148   fNevt++;
 149   fCFManager->SetEventInfo(fMCEvent);
 150   Double_t containerInput[18] ;
 151
 152 ////////
 153 //// MC
 154 ////////
 155
 156 // loop on the MC part
 157
 158   if(fIsMC){
 159
 160       for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 161
 162           AliMCParticle *mcPart  = (AliMCParticle*) fMCEvent->GetTrack(ipart);
 163           TParticle *part = mcPart->Particle();
 164
 165           // rapidity and Pt cuts (default -4<y<-2.5 et 0<pt<20)
 166           if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
 167
 168           // Selection of mu
 169           Int_t pdg=part->GetPdgCode();
 170           if(TMath::Abs(pdg)==13){
 171
 172               Float_t emc = part->Energy();
 173               Float_t pzmc = part->Pz();
 174               Float_t rapmc = Rap(emc,pzmc);
 175               Float_t phimc = part->Phi();
 176               phimc = Phideg(phimc);
 177               Float_t ptmc = part->Pt();
 178               Float_t pmc = part->P();
 179               Float_t zmc = part->Vz();
 180               Float_t etamc = part->Eta();
 181               Float_t chargemc=0;
 182               if(pdg==13) chargemc=-1;
 183               if(pdg==-13) chargemc=1;
 184               Float_t thetamc;
 185               thetamc = TMath::Abs(-TMath::Pi()+TMath::ASin(ptmc/pmc)); // convention for EVGEN
 186               Float_t ymc = part->Vy();
 187               Float_t xmc = part->Vx();
 188               Float_t rmc = TMath::Sqrt(xmc*xmc+ymc*ymc);
 189
 190               containerInput[0] = etamc ;
 191               containerInput[1] = rapmc ;
 192               containerInput[2] = phimc ;
 193               containerInput[3] = ptmc ;
 194               containerInput[4] = pmc ;
 195               containerInput[14] = zmc ;
 196               containerInput[16] = chargemc ;
 197               containerInput[17] = thetamc ;
 198
 199               // 10 var calculated only for ESD .i.e set at 1 in MC step
 200               for(Int_t i=5; i<=13; i++){
 201                   containerInput[i] = 1;
 202               }
 203
 204               containerInput[7] = 3;
 205               containerInput[10] = rmc ; // radius at production
 206               containerInput[15] = 175. ;  // rabsmc
 207
 208               // fill the container at the first step
 209               fCFManager->GetParticleContainer()->Fill(containerInput,0);
 210           }
 211       }
 212   }
 213
 214
 215 ////////
 216 //// ESD
 217 ////////
 218
 219   Int_t trig1MuL = 1 ,trig1MuH = 2, trigUSL = 4, trigUSH = 8, trigLSL = 16, trigLSH = 32;
 220   Int_t trig = 0;
 221
 222   AliESDEvent *fESD;
 223   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
 224       (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 225   if ( ! esdH) {
 226     AliError("Cannot get input event handler");
 227     return;
 228   }
 229   fESD = esdH->GetEvent();
 230
 231 // trigger
 232
 233   if (fESD->GetTriggerMask() & trig1MuL) trig = 1;
 234   if (fESD->GetTriggerMask() & trig1MuH) trig = 2;
 235   if (fESD->GetTriggerMask() & trigUSL)  trig = 3;
 236   if (fESD->GetTriggerMask() & trigUSH)  trig = 4;
 237   if (fESD->GetTriggerMask() & trigLSL)  trig = 5;
 238   if (fESD->GetTriggerMask() & trigLSH)  trig = 6;
 239
 240 // vertex
 241
 242   Float_t vx = -200 , vy = -200 , vz = -200, vt=-200;
 243
 244   Double_t Ncont = 0.;
 245   AliESDVertex* vertex = (AliESDVertex*) fESD->GetVertex();
 246   if (vertex->GetNContributors()) {
 247       vz = vertex->GetZ();
 248       vy = vertex->GetY();
 249       vx = vertex->GetX();
 250       vt=TMath::Sqrt(vx*vx+vy*vy);
 251       Ncont=vertex->GetNContributors();
 252   }
 253
 254 // tracks
 255
 256   Int_t mult1 = fESD->GetNumberOfMuonTracks() ;
 257
 258     for (Int_t j = 0; j<mult1; j++) {
 259         AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
 260
 261         if (!mu1->ContainTrackerData()) { mu1=0; continue; }
 262
 263         Float_t charger = mu1->Charge();
 264         Float_t er = mu1->E();
 265         Float_t pzr = mu1->Pz();
 266         Float_t rapr=Rap(er,pzr);
 267         Float_t phir = mu1->Phi();
 268         phir = Phideg(phir);
 269         Float_t ptr = mu1->Pt();
 270         Float_t pr = mu1->P();
 271         Float_t hit = mu1->GetNHit();
 272         Float_t chi2 = mu1->GetChi2() / (2.0 * hit - 5);
 273         Int_t matchtrig = mu1->GetMatchTrigger();
 274         Float_t chi2match = mu1->GetChi2MatchTrigger();
 275         Float_t dcar = mu1->GetDCA();
 276         Float_t zr = mu1->GetZ();
 277         Float_t etar = mu1->Eta();
 278         Float_t rabs = mu1->GetRAtAbsorberEnd();
 279         Float_t thetar ;
 280         thetar = TMath::Abs(-TMath::Pi()+TMath::ASin(ptr/pr));
 281
 282 // rapidity and Pt cuts (default -4<y<-2.5 et 0<pt<20)
 283         if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu1)) continue;
 284
 285         Float_t thetaAbs = 180.*(1.-TMath::ATan(rabs/505.)/TMath::Pi());
 286         if (thetaAbs<171 || thetaAbs>178) continue;
 287
 288         containerInput[0] = etar ;
 289         containerInput[1] = rapr ;
 290         containerInput[2] = phir ;
 291         containerInput[3] = ptr ;
 292         containerInput[4] = pr ;
 293         containerInput[5] = hit ;
 294         containerInput[6] = chi2 ;
 295         containerInput[7] = matchtrig ;
 296         containerInput[8] = chi2match ;
 297         containerInput[9] = Ncont ;
 298         containerInput[10] = vt ;
 299         containerInput[11] = vz ;
 300         containerInput[12] = trig ;
 301         containerInput[13] = dcar ;
 302         containerInput[14] = zr ;
 303         containerInput[15] = thetaAbs;
 304         containerInput[16] = charger ;
 305         containerInput[17] = thetar ;
 306
 307 // fill the container at the second step (for simu, first for data)
 308         if(fIsMC){
 309             fCFManager->GetParticleContainer()->Fill(containerInput,1);
 310         }
 311         else{
 312             fCFManager->GetParticleContainer()->Fill(containerInput,0);
 313         }
 314
 315
 316     }      // mu Loop
 317
 318 //  ----------
 319   fHistEventsProcessed->Fill(0);
 320   PostData(1,fHistEventsProcessed) ;
 321   PostData(2,fCFManager->GetParticleContainer()) ;
 322
 323 }
 324 //________________________________________________________________________
 325 Float_t AliCFMuonSingleTask1::Rap(Float_t e, Float_t pz)
 326 {
 327 // calculate rapidity
 328     Float_t rap;
 329     if(e!=pz){
 330         rap = 0.5*TMath::Log((e+pz)/(e-pz));
 331         return rap;
 332     }
 333     else{
 334         rap = -200;
 335         return rap;
 336     }
 337 }
 338 //________________________________________________________________________
 339 Float_t AliCFMuonSingleTask1::Phideg(Float_t phi)
 340 {
 341 // calculate Phi in range [-180,180]
 342     Float_t phideg;
 343
 344         phideg = phi-TMath::Pi();
 345         phideg = phideg*57.296;
 346         return phideg;
 347 }
 348 //________________________________________________________________________
 349 void AliCFMuonSingleTask1::Terminate(Option_t *)
 350 {
 351   // project pt (var[3]) from the two steps MC(0) and ESD(1)
 352
 353     AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
 354     if ( ! cont ) {
 355       AliError("Cannot find container in file");
 356       return;
 357     }
 358
 359     TH1D *kpt = cont->ShowProjection(16,0);
 360     TH1D *rpt = cont->ShowProjection(16,1);
 361
 362     TCanvas *c1 = new TCanvas("AliCFMuonSingleTask1"," MC & ESD",10,10,510,510);
 363     c1->Divide(1,2);
 364     c1->cd(1);
 365     kpt->Draw("HIST");
 366     c1->cd(2);
 367     rpt->Draw("HIST");
 368
 369 }
 370 //________________________________________________________________________
 371
 372 #endif