PWGLF/FORWARD/photons/AliAnalysisTaskPMD.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 "AliESDVZERO.h"
  13 #include "AliAnalysisTaskPMD.h"
  14
  15 // AnalysisTask For PMD
  16 // Authors: Sudipan De, Subhash Singha
  17
  18 ClassImp(AliAnalysisTaskPMD)
  19
  20 //________________________________________________________________________
  21 AliAnalysisTaskPMD::AliAnalysisTaskPMD(const char *name)
  22 : AliAnalysisTaskSE(name),
  23   fESD(0),
  24   fOutputList(0),
  25   fHistTotEvent(0),
  26   fHistTotEventAfterPhySel(0),
  27   fHistTotEventAfterVtx(0),
  28   fHistVtxZ(0),
  29   fHistXYPre(0),
  30   fHistEta(0),
  31   fHistEta1(0),
  32   fHistMultMeas(0),
  33   fHistMultMeas1(0)
  34 {
  35   for(Int_t i=0; i<10; i++){
  36     fHistMultMeasEtaBinA[i] = 0;
  37     fHistMultMeasEtaBinA1[i] = 0;
  38   }
  39   // Constructor
  40
  41   // Define input and output slots here
  42   // Input slot #0 works with a TChain
  43   DefineInput(0, TChain::Class());
  44   // Output slot #0 id reserved by the base class for AOD
  45   // Output slot #1 writes into a TH1 container
  46   DefineOutput(1, TList::Class());
  47 }
  48
  49 //________________________________________________________________________
  50 void AliAnalysisTaskPMD::UserCreateOutputObjects()
  51 {
  52   // Create histograms
  53   // Called once
  54
  55   fOutputList = new TList();
  56   Int_t kNbinsMultA = 50;  Float_t XminMultA = -0.5; Float_t XmaxMultA = 49.5;
  57   Int_t kNbinsXY = 200; Float_t XminXY = -100.0; Float_t XmaxXY  = 100.0;
  58   Int_t kNBinsEvent = 10; Float_t XminEvent = 0; Float_t XmaxEvent = 10;
  59   Int_t kNBinsEta = 10; Float_t XminEta = 2.1; Float_t XmaxEta = 4.1;
  60
  61   fHistTotEvent = new TH1F("TotEvent","TotEvent",
  62                            kNBinsEvent,XminEvent,XmaxEvent);
  63   fOutputList->Add(fHistTotEvent);
  64   fHistTotEventAfterPhySel = new TH1F("TotEventAfterPhySel","TotEventAfterPhySel",
  65                                       kNBinsEvent,XminEvent,XmaxEvent);
  66   fOutputList->Add(fHistTotEventAfterPhySel);
  67   fHistTotEventAfterVtx = new TH1F("TotEventAfterVtx","TotEventAfterVtx",
  68                                    kNBinsEvent,XminEvent,XmaxEvent);
  69   fOutputList->Add(fHistTotEventAfterVtx);
  70   fHistVtxZ = new TH1F("VtxZ","VtxZ",100,-50,50);
  71   fOutputList->Add(fHistVtxZ);
  72   fHistXYPre = new TH2F("XYPre","XYPre",kNbinsXY,XminXY,XmaxXY,
  73                         kNbinsXY,XminXY,XmaxXY);
  74   fOutputList->Add(fHistXYPre);
  75   fHistEta = new TH1F ("Eta","Eta",kNBinsEta,XminEta,XmaxEta);
  76   fOutputList->Add(fHistEta);
  77   fHistEta1 = new TH1F ("Eta1","Eta1",kNBinsEta,XminEta,XmaxEta);
  78   fOutputList->Add(fHistEta1);
  79   fHistMultMeas = new TH1F("MultM","MultM",100,-0.5,99.5);
  80   fOutputList->Add(fHistMultMeas);
  81   fHistMultMeas1 = new TH1F("MultM1","MultM1",100,-0.5,99.5);
  82   fOutputList->Add(fHistMultMeas1);
  83
  84   Char_t nameM[256], nameM1[256];
  85   for(Int_t i=0; i<10; i++){
  86     sprintf(nameM,"MultM_EtaBin%d",i+1);
  87     fHistMultMeasEtaBinA[i] =
  88       new TH1F(nameM,nameM,kNbinsMultA,
  89                XminMultA,XmaxMultA);
  90     fOutputList->Add(fHistMultMeasEtaBinA[i]);
  91     sprintf(nameM1,"MultM1_EtaBin%d",i+1);
  92     fHistMultMeasEtaBinA1[i] =
  93       new TH1F(nameM1,nameM1,kNbinsMultA,
  94                XminMultA,XmaxMultA);
  95     fOutputList->Add(fHistMultMeasEtaBinA1[i]);
  96   }//i loop
  97
  98 }
  99
 100 //________________________________________________________________________
 101 void AliAnalysisTaskPMD::UserExec(Option_t *)
 102 {
 103   // Main loop
 104   // Called for each event
 105   Float_t MipCut1 = 432;//MPV=72=> 6*MPV=432
 106   Float_t MipCut2 = 648;//MPV=72=> 9*MPV=648
 107   Float_t etacls=0., theta=0., rdist=0.;
 108   Int_t PhotonCls = 0;//# of photon measured within 2.3 to 3.9 #eta
 109   Int_t PhotonCls1 = 0;//# of photon measured within 2.3 to 3.9 #eta
 110   Int_t PhotonClsAEtaBin[10] = {0};//# of photon measured with diff Eta bin
 111   Int_t PhotonClsAEtaBin1[10] = {0};//# of photon measured with diff Eta bin
 112   // Post output data.
 113   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 114   if (!fESD) {
 115     printf("ERROR: fESD not available\n");
 116     return;
 117   }
 118   fHistTotEvent->Fill(5);
 119   //event selection
 120   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 121   if (! isSelected) return;
 122   fHistTotEventAfterPhySel->Fill(5);
 123   //printf("There are %d tracks in this event\n", fESD->GetNumberOfPmdTracks());
 124   //vertex selection
 125   const AliESDVertex *vertex = fESD->GetPrimaryVertex();
 126   Float_t Vz = vertex->GetZ();
 127   Bool_t zVerStatus = vertex->GetStatus();
 128   if(zVerStatus){
 129     fHistVtxZ->Fill(Vz);
 130     if(TMath::Abs(Vz)<10){
 131       fHistTotEventAfterVtx->Fill(5);
 132       Int_t ptracks = fESD->GetNumberOfPmdTracks();
 133       //PMDtrack Loop
 134       for(Int_t kk=0;kk<ptracks;kk++){
 135         AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(kk);
 136         Int_t   det   = pmdtr->GetDetector();
 137         Float_t clsX  = pmdtr->GetClusterX();
 138         Float_t clsY  = pmdtr->GetClusterY();
 139         Float_t clsZ  = pmdtr->GetClusterZ();
 140         clsZ -= Vz;
 141         Float_t ncell = pmdtr->GetClusterCells();
 142         Float_t adc   = pmdtr->GetClusterADC();
 143         //Int_t smn = pmdtr->GetSmn();
 144         //Calculation of #eta
 145         rdist = TMath::Sqrt(clsX*clsX + clsY*clsY);
 146         if(clsZ!=0) theta = TMath::ATan2(rdist,clsZ);
 147         etacls  = -TMath::Log(TMath::Tan(0.5*theta));
 148
 149         if(det==0 && adc>72)fHistXYPre->Fill(clsX,clsY);
 150         if(det==0 && adc>MipCut1 && ncell>2){
 151           if(etacls > 2.1 && etacls<= 2.3)PhotonClsAEtaBin[0]++;
 152           if(etacls > 2.3 && etacls<= 2.5)PhotonClsAEtaBin[1]++;
 153           if(etacls > 2.5 && etacls<= 2.7)PhotonClsAEtaBin[2]++;
 154           if(etacls > 2.7 && etacls<= 2.9)PhotonClsAEtaBin[3]++;
 155           if(etacls > 2.9 && etacls<= 3.1)PhotonClsAEtaBin[4]++;
 156           if(etacls > 3.1 && etacls<= 3.3)PhotonClsAEtaBin[5]++;
 157           if(etacls > 3.3 && etacls<= 3.5)PhotonClsAEtaBin[6]++;
 158           if(etacls > 3.5 && etacls<= 3.7)PhotonClsAEtaBin[7]++;
 159           if(etacls > 3.7 && etacls<= 3.9)PhotonClsAEtaBin[8]++;
 160           if(etacls > 3.9 && etacls<= 4.1)PhotonClsAEtaBin[9]++;
 161           if(etacls > 2.3 && etacls<= 3.9)PhotonCls++;
 162           fHistEta->Fill(etacls);
 163         }//if MipCut1
 164         if(det==0 && adc>MipCut2 && ncell>2){
 165           if(etacls > 2.1 && etacls<= 2.3)PhotonClsAEtaBin1[0]++;
 166           if(etacls > 2.3 && etacls<= 2.5)PhotonClsAEtaBin1[1]++;
 167           if(etacls > 2.5 && etacls<= 2.7)PhotonClsAEtaBin1[2]++;
 168           if(etacls > 2.7 && etacls<= 2.9)PhotonClsAEtaBin1[3]++;
 169           if(etacls > 2.9 && etacls<= 3.1)PhotonClsAEtaBin1[4]++;
 170           if(etacls > 3.1 && etacls<= 3.3)PhotonClsAEtaBin1[5]++;
 171           if(etacls > 3.3 && etacls<= 3.5)PhotonClsAEtaBin1[6]++;
 172           if(etacls > 3.5 && etacls<= 3.7)PhotonClsAEtaBin1[7]++;
 173           if(etacls > 3.7 && etacls<= 3.9)PhotonClsAEtaBin1[8]++;
 174           if(etacls > 3.9 && etacls<= 4.1)PhotonClsAEtaBin1[9]++;
 175           if(etacls > 2.3 && etacls<= 3.9)PhotonCls1++;
 176           fHistEta1->Fill(etacls);
 177         }//if Mipcut2
 178       } //PMDtrack loop
 179       for(Int_t i=0; i<10; i++){
 180         fHistMultMeasEtaBinA[i]->Fill(PhotonClsAEtaBin[i]);
 181         fHistMultMeasEtaBinA1[i]->Fill(PhotonClsAEtaBin1[i]);
 182       }//i loop
 183       fHistMultMeas->Fill(PhotonCls);
 184       fHistMultMeas1->Fill(PhotonCls1);
 185     }// Vz loop
 186   }// Vzcut loop
 187   PostData(1, fOutputList);
 188 }
 189
 190 //________________________________________________________________________
 191 void AliAnalysisTaskPMD::Terminate(Option_t *)
 192 {
 193   // Draw result to the screen
 194   // Called once at the end of the query
 195
 196   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 197   if (!fOutputList) {
 198     printf("ERROR: Output list not available\n");
 199     return;
 200   }
 201 }