PWGLF/FORWARD/photons/AliAnalysisTaskPMDSim.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TCanvas.h"
   6 #include "AliAnalysisTask.h"
   7 #include "AliAnalysisManager.h"
   8 #include "AliESDEvent.h"
   9 #include "AliESDInputHandler.h"
  10 #include "AliESDPmdTrack.h"
  11 #include "AliESDVertex.h"
  12 #include "AliAnalysisTaskPMDSim.h"
  13 #include <TParticle.h>
  14 #include <AliLog.h>
  15 #include <AliStack.h>
  16 #include <AliHeader.h>
  17 #include <AliGenEventHeader.h>
  18 #include <AliHeader.h>
  19 #include <AliGenPythiaEventHeader.h>
  20 #include <AliMCEventHandler.h>
  21 #include <AliMCEvent.h>
  22 #include "AliTriggerAnalysis.h"
  23 #include "AliPhysicsSelection.h"
  24
  25 // AnalysisTask For PMD
  26 // Authors: Sudipan De, Subhash Singha
  27
  28 ClassImp(AliAnalysisTaskPMDSim)
  29
  30 //________________________________________________________________________
  31 AliAnalysisTaskPMDSim::AliAnalysisTaskPMDSim(const char *name)
  32   : AliAnalysisTaskSE(name),
  33   fESD(0),
  34   fOutputList(0),
  35   fHistTotEvent(0),
  36   fHistTotEventAfterPhySel(0),
  37   fHistTotEventAfterVtx(0),
  38   fVtxZ(0),
  39   fHistXYPre(0),
  40   fHistEtaPhM(0),
  41   fHistEtaPhM1(0),
  42   fHistEtaT(0),
  43   fMultMeasured(0),
  44   fMultMeasured1(0),
  45   fMultTrue(0),
  46   fMultCorr(0),
  47   fMultCorr1(0)
  48 {
  49   for(Int_t i=0; i<10; i++){
  50     fHistMultMeasEtaBinA[i] = 0;
  51     fHistMultMeasEtaBinA1[i] = 0;
  52     fHistMultTrueEtaBinA[i] = 0;
  53     fHistMultCorrEtaBinA[i] = 0;
  54     fHistMultCorrEtaBinA1[i] = 0;
  55   }
  56
  57   // Constructor
  58
  59   // Define input and output slots here
  60   // Input slot #0 works with a TChain
  61   DefineInput(0, TChain::Class());
  62   // Output slot #0 id reserved by the base class for AOD
  63   // Output slot #1 writes into a TH1 container
  64   DefineOutput(1, TList::Class());
  65 }
  66
  67 //________________________________________________________________________
  68 void AliAnalysisTaskPMDSim::UserCreateOutputObjects()
  69 {
  70   // Create histograms
  71   // Called once
  72
  73   fOutputList = new TList();
  74
  75   Int_t kNbinsMult = 100;  Float_t XminMult = -0.5;
  76   Float_t XmaxMult = 99.5;//mult
  77   Int_t kNbinsMultA = 50;  Float_t XminMultA = -0.5;
  78   Float_t XmaxMultA = 49.5;//mult
  79   Int_t kNbinEtaPMDCov = 10; Float_t XminEtaPMDCov = 2.1; Float_t XmaxEtaPMDCov = 4.1;//etaPMD
  80   Int_t kNBinsEvent = 10; Float_t XminEvent = 0; Float_t XmaxEvent = 10;
  81   Int_t kNbinsXY = 200; Float_t XminXY = -100.0; Float_t XmaxXY  = 100.0;
  82   fHistTotEvent = new TH1F("TotEvent","TotEvent",
  83                            kNBinsEvent,XminEvent,XmaxEvent);
  84   fOutputList->Add(fHistTotEvent);
  85   fHistTotEventAfterPhySel = new TH1F("TotEventAfterPhySel","TotEventAfterPhySel",
  86                                       kNBinsEvent,XminEvent,XmaxEvent);
  87   fOutputList->Add(fHistTotEventAfterPhySel);
  88   fHistTotEventAfterVtx = new TH1F("TotEventAfterVtx","TotEventAfterVtx",
  89                                   kNBinsEvent,XminEvent,XmaxEvent);
  90   fOutputList->Add(fHistTotEventAfterVtx);
  91   fVtxZ     = new TH1F("VtxZ","VtxZ",200,-50.0,50.0);
  92   fOutputList->Add(fVtxZ);
  93   fHistXYPre = new TH2F("XYPre","XYPre",kNbinsXY,XminXY,XmaxXY,
  94                         kNbinsXY,XminXY,XmaxXY);
  95   fOutputList->Add(fHistXYPre);
  96   fHistEtaPhM = new TH1F("fHistEtaPhM","fHistEtaPhM",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
  97   fOutputList->Add(fHistEtaPhM);
  98   fHistEtaPhM1 = new TH1F("fHistEtaPhM1","fHistEtaPhM1",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
  99   fOutputList->Add(fHistEtaPhM1);
 100   fHistEtaT = new TH1F("fHistEtaT","fHistEtaT",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
 101   fOutputList->Add(fHistEtaT);
 102   fMultMeasured = new TH1F("MultM","MultM",kNbinsMult,XminMult,XmaxMult);
 103   fOutputList->Add(fMultMeasured);
 104   fMultMeasured1 = new TH1F("MultM1","MultM1",kNbinsMult,XminMult,XmaxMult);
 105   fOutputList->Add(fMultMeasured1);
 106   fMultTrue = new TH1F("MultT","MultT",kNbinsMult,XminMult,XmaxMult);
 107   fOutputList->Add(fMultTrue);
 108   fMultCorr = new TH2F("MultCorr","MultCorr",kNbinsMult,XminMult,XmaxMult,
 109                        kNbinsMult,XminMult,XmaxMult);
 110   fOutputList->Add(fMultCorr);
 111   fMultCorr1 = new TH2F("MultCorr1","MultCorr1",kNbinsMult,XminMult,XmaxMult,
 112                        kNbinsMult,XminMult,XmaxMult);
 113   fOutputList->Add(fMultCorr1);
 114
 115
 116   Char_t nameT[256], nameM[256], nameCorr[256],nameM1[256], nameCorr1[256];
 117   for(Int_t i=0; i<10; i++){
 118     sprintf(nameM,"MultM_EtaBin%d",i+1);
 119     sprintf(nameM1,"MultM1_EtaBin%d",i+1);
 120     sprintf(nameT,"MultT_EtaBin%d",i+1);
 121     sprintf(nameCorr,"MultCorr_EtaBin%d",i+1);
 122     sprintf(nameCorr1,"MultCorr1_EtaBin%d",i+1);
 123     fHistMultMeasEtaBinA[i] =
 124       new TH1F(nameM,nameM,kNbinsMultA,XminMultA,XmaxMultA);
 125     fHistMultMeasEtaBinA1[i] =
 126       new TH1F(nameM1,nameM1,kNbinsMultA,XminMultA,XmaxMultA);
 127     fHistMultTrueEtaBinA[i] =
 128       new TH1F(nameT,nameT,kNbinsMultA,XminMultA,XmaxMultA);
 129     fHistMultCorrEtaBinA[i] =
 130       new TH2F(nameCorr,nameCorr,kNbinsMultA,XminMultA,XmaxMultA,
 131                kNbinsMultA,XminMultA,XmaxMultA);
 132     fHistMultCorrEtaBinA1[i] =
 133       new TH2F(nameCorr1,nameCorr1,kNbinsMultA,XminMultA,XmaxMultA,
 134                kNbinsMultA,XminMultA,XmaxMultA);
 135     fOutputList->Add(fHistMultMeasEtaBinA[i]);
 136     fOutputList->Add(fHistMultMeasEtaBinA1[i]);
 137     fOutputList->Add(fHistMultTrueEtaBinA[i]);
 138     fOutputList->Add(fHistMultCorrEtaBinA[i]);
 139     fOutputList->Add(fHistMultCorrEtaBinA1[i]);
 140   }//i loop
 141   //fOutputList->Add(fHistPt);
 142 }
 143
 144 //________________________________________________________________________
 145 void AliAnalysisTaskPMDSim::UserExec(Option_t *)
 146 {
 147   // Main loop
 148   // Called for each event
 149   Float_t MipCut1 = 432;//MPV=432=> 6*MPV=432
 150   Float_t MipCut2 = 648;//MPV=432=> 6*MPV=648
 151   Float_t etacls=0., theta=0., rdist=0.;
 152   Int_t PhotonCls = 0;
 153   Int_t PhotonCls1 = 0;
 154   Int_t PhotonClsAEtaBin[10] = {0};//#of photon measured with diff Eta bin
 155   Int_t PhotonClsAEtaBin1[10] = {0};//#of photon measured with diff Eta bin
 156   Int_t PhotonTrueAEtaBin[10] = {0};//#of photon Incident with diff Eta bin
 157   // Post output data.
 158   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 159   if (!fESD) {
 160     printf("ERROR: fESD not available\n");
 161     return;
 162   }
 163   fHistTotEvent->Fill(5);
 164   //event selection
 165   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 166   if (! isSelected) return;
 167   fHistTotEventAfterPhySel->Fill(5);
 168   //Vertex selection
 169   const AliESDVertex *vertex = fESD->GetPrimaryVertex();
 170   Float_t Vz = vertex->GetZ();
 171   // Float_t Vx = vertex->GetX();
 172   //Float_t Vy = vertex->GetY();
 173   Bool_t zVerStatus = vertex->GetStatus();
 174   if(zVerStatus){
 175     fVtxZ->Fill(Vz);
 176     if(TMath::Abs(Vz)<10){
 177       fHistTotEventAfterVtx->Fill(5);
 178       Int_t ptracks = fESD->GetNumberOfPmdTracks();
 179       for(Int_t kk=0;kk<ptracks;kk++){
 180         AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(kk);
 181         Int_t   det   = pmdtr->GetDetector();
 182         Float_t clsX  = pmdtr->GetClusterX();
 183         Float_t clsY  = pmdtr->GetClusterY();
 184         Float_t clsZ  = pmdtr->GetClusterZ();
 185         clsZ -= Vz;
 186         Float_t ncell = pmdtr->GetClusterCells();
 187         Float_t adc   = pmdtr->GetClusterADC();
 188         //Float_t pid   = pmdtr->GetClusterPID();
 189         //Float_t isotag = pmdtr->GetClusterSigmaX();
 190         //Int_t trno = pmdtr->GetClusterTrackNo();
 191         //calculation of #eta
 192         rdist = TMath::Sqrt(clsX*clsX + clsY*clsY);
 193         if(clsZ!=0) theta = TMath::ATan2(rdist,clsZ);
 194         etacls   = -TMath::Log(TMath::Tan(0.5*theta));
 195
 196         if(det==0 && adc>0)fHistXYPre->Fill(clsX,clsY);
 197         if(det==0 && adc>MipCut1 && ncell>2){
 198           fHistEtaPhM->Fill(etacls);
 199           if(etacls > 2.3 && etacls <= 3.9)  PhotonCls++;
 200           if(etacls > 2.1 && etacls <= 2.3)  PhotonClsAEtaBin[0]++;
 201           if(etacls > 2.3 && etacls <= 2.5)  PhotonClsAEtaBin[1]++;
 202           if(etacls > 2.5 && etacls <= 2.7)  PhotonClsAEtaBin[2]++;
 203           if(etacls > 2.7 && etacls <= 2.9)  PhotonClsAEtaBin[3]++;
 204           if(etacls > 2.9 && etacls <= 3.1)  PhotonClsAEtaBin[4]++;
 205           if(etacls > 3.1 && etacls <= 3.3)  PhotonClsAEtaBin[5]++;
 206           if(etacls > 3.3 && etacls <= 3.5)  PhotonClsAEtaBin[6]++;
 207           if(etacls > 3.5 && etacls <= 3.7)  PhotonClsAEtaBin[7]++;
 208           if(etacls > 3.7 && etacls <= 3.9)  PhotonClsAEtaBin[8]++;
 209           if(etacls > 3.9 && etacls <= 4.1)  PhotonClsAEtaBin[9]++;
 210         }//if MipCut1
 211         if( det==0 && adc>MipCut2 && ncell>2){
 212           fHistEtaPhM1->Fill(etacls);
 213           if(etacls > 2.3 && etacls <= 3.9)  PhotonCls1++;
 214           if(etacls > 2.1 && etacls <= 2.3)  PhotonClsAEtaBin1[0]++;
 215           if(etacls > 2.3 && etacls <= 2.5)  PhotonClsAEtaBin1[1]++;
 216           if(etacls > 2.5 && etacls <= 2.7)  PhotonClsAEtaBin1[2]++;
 217           if(etacls > 2.7 && etacls <= 2.9)  PhotonClsAEtaBin1[3]++;
 218           if(etacls > 2.9 && etacls <= 3.1)  PhotonClsAEtaBin1[4]++;
 219           if(etacls > 3.1 && etacls <= 3.3)  PhotonClsAEtaBin1[5]++;
 220           if(etacls > 3.3 && etacls <= 3.5)  PhotonClsAEtaBin1[6]++;
 221           if(etacls > 3.5 && etacls <= 3.7)  PhotonClsAEtaBin1[7]++;
 222           if(etacls > 3.7 && etacls <= 3.9)  PhotonClsAEtaBin1[8]++;
 223           if(etacls > 3.9 && etacls <= 4.1)  PhotonClsAEtaBin1[9]++;
 224         }//if MipCut2
 225       } //track loop
 226       //reading MC info.
 227       AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*>
 228         (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 229       if (!eventHandler) {
 230         Printf("ERROR: Could not retrieve MC event handler");
 231         return;
 232       }
 233       AliMCEvent* mcEvent = eventHandler->MCEvent();
 234       if (!mcEvent) {
 235         Printf("ERROR: Could not retrieve MC event");
 236         return;
 237       }
 238       AliStack* stack = mcEvent->Stack();
 239       if (!stack)
 240         {
 241           AliDebug(AliLog::kError, "Stack not available");
 242           return;
 243         }
 244       if(stack){
 245         Int_t nPrim  = stack->GetNprimary();
 246         Int_t PhotonTrue = 0;
 247         for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 248           {
 249             TParticle *MPart = stack->Particle(iMc);
 250             Int_t mpart  = MPart->GetPdgCode();
 251             Float_t eta   = MPart->Eta();
 252             if(mpart == 22){
 253               fHistEtaT->Fill(eta);
 254               if(eta > 2.3 && eta <= 3.9)PhotonTrue++;
 255               if(eta > 2.1 && eta <= 2.3)PhotonTrueAEtaBin[0]++;
 256               if(eta > 2.3 && eta <= 2.5)PhotonTrueAEtaBin[1]++;
 257               if(eta > 2.5 && eta <= 2.7)PhotonTrueAEtaBin[2]++;
 258               if(eta > 2.7 && eta <= 2.9)PhotonTrueAEtaBin[3]++;
 259               if(eta > 2.9 && eta <= 3.1)PhotonTrueAEtaBin[4]++;
 260               if(eta > 3.1 && eta <= 3.3)PhotonTrueAEtaBin[5]++;
 261               if(eta > 3.3 && eta <= 3.5)PhotonTrueAEtaBin[6]++;
 262               if(eta > 3.5 && eta <= 3.7)PhotonTrueAEtaBin[7]++;
 263               if(eta > 3.7 && eta <= 3.9)PhotonTrueAEtaBin[8]++;
 264               if(eta > 3.9 && eta <= 4.1)PhotonTrueAEtaBin[9]++;
 265             }//mpart
 266           }//track loop
 267         for(Int_t i=0; i<10; i++){
 268           fHistMultMeasEtaBinA[i]->Fill(PhotonClsAEtaBin[i]);
 269           fHistMultMeasEtaBinA1[i]->Fill(PhotonClsAEtaBin1[i]);
 270           fHistMultTrueEtaBinA[i]->Fill(PhotonTrueAEtaBin[i]);
 271           fHistMultCorrEtaBinA[i]->Fill(PhotonTrueAEtaBin[i],PhotonClsAEtaBin[i]);
 272           fHistMultCorrEtaBinA1[i]->Fill(PhotonTrueAEtaBin[i],PhotonClsAEtaBin1[i]);
 273         }//i loop
 274         fMultMeasured->Fill(PhotonCls);
 275         fMultMeasured1->Fill(PhotonCls1);
 276         fMultTrue->Fill(PhotonTrue);
 277         fMultCorr->Fill(PhotonTrue,PhotonCls);
 278         fMultCorr1->Fill(PhotonTrue,PhotonCls1);
 279       }//if stack
 280     }//if Vz!= 0
 281   }//Vz loop
 282   PostData(1, fOutputList);
 283 }//UserExec()
 284
 285 //_______________________________________________________________________
 286 void AliAnalysisTaskPMDSim::Terminate(Option_t *)
 287 {
 288   // Draw result to the screen
 289   // Called once at the end of the query
 290   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 291   if (!fOutputList) {
 292     printf("ERROR: Output list not available\n");
 293     return;
 294   }
 295 }//Terminate