T0/AliT0MultiplicityTask.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TList.h"
   6 #include "TClonesArray.h"
   7 #include "TCanvas.h"
   8 #include "TStyle.h"
   9 #include "TLine.h"
  10 #include "TText.h"
  11 #include "TFile.h"
  12 #include "TBenchmark.h"
  13
  14 #include "AliAnalysisTask.h"
  15 #include "AliAnalysisManager.h"
  16
  17 #include "AliLog.h"
  18 #include "AliVEvent.h"
  19 #include "AliESDEvent.h"
  20 #include "AliMCEvent.h"
  21 #include "AliESDtrack.h"
  22 #include "AliStack.h"
  23 #include "AliTrackReference.h"
  24
  25 #include "AliT0MultiplicityTask.h"
  26
  27 ClassImp(AliT0MultiplicityTask)
  28
  29
  30 AliT0MultiplicityTask::AliT0MultiplicityTask()
  31   : AliAnalysisTaskSE("AliT0MultiplicityTask"),
  32     fListOfHistos(0),
  33     fOrA(0),
  34     fOrC(0),
  35     fMean(0),
  36     fVertex(0),
  37     fTime(0),
  38     fAmp(0),
  39     fTotalMult(0),
  40     fMultRecTotal(0),
  41     fMultRecRealA(0),
  42     fMultRecRealC(0),
  43     fPrim(0),
  44     fRef(0),
  45     fRec(0)
  46 {
  47   // Default constructor
  48   AliInfo("Default constructor AliT0MultiplicityTask");
  49   // Define input and output slots here
  50   // Input slot #0 works with a TChain
  51   DefineInput(0, TChain::Class());
  52   // Output slot #1 TList
  53   DefineOutput(1, TList::Class());
  54 }
  55
  56
  57 AliT0MultiplicityTask::AliT0MultiplicityTask(const char* name)
  58   : AliAnalysisTaskSE(name),
  59     fListOfHistos(0),
  60     fOrA(0),
  61     fOrC(0),
  62     fMean(0),
  63     fVertex(0),
  64     fTime(0),
  65     fAmp(0),
  66     fTotalMult(0),
  67     fMultRecTotal(0),
  68     fMultRecRealA(0),
  69     fMultRecRealC(0),
  70     fPrim(0),
  71     fRef(0),
  72     fRec(0)
  73 {
  74   // Constructor
  75   // Define input and output slots here
  76   // Input slot #0 works with a TChain
  77   DefineInput(0, TChain::Class());
  78   // Output slot #1 TList
  79   DefineOutput(1, TList::Class());
  80 }
  81
  82
  83 void AliT0MultiplicityTask::UserCreateOutputObjects()
  84 {
  85   // Create histograms
  86   // Called once
  87   // Create output container
  88   fListOfHistos = new TList();
  89
  90
  91   fOrA = new TH1F("hOrA"," T0 A", 100,2450,2700);
  92   fOrC = new TH1F("hOrC"," T0 C", 100,2450,2700);
  93   fMean= new TH1F("hMean"," T0 ",100, 12000, 13000);
  94   fVertex = new TH1F("hVertex","Z position of vertex",   100,-30, 30);
  95
  96   fTime = new TH1F("Time", " Amp vs Time",100, 12000, 13000);
  97   fAmp = new TH1F("fAmp", " Amplitude", 100, 0, 200);
  98
  99   //  fHighMult  = new TH1F("fHighMult", " events with high amp ",100, 0, 10);
 100
 101   fTotalMult = new TH1F("fTotalMult","total multiplicity",500,0,5000);
 102
 103   fMultRecRealA = new TH2F("fMultRecRealA"," ",100,0.,200,100,0,200);
 104   fMultRecRealC = new TH2F("fMultRecRealC"," ",100,0,200,100,0,200);
 105   fMultRecTotal = new TH2F("fMultRecTotal"," ",100,0,200,100,0,200);
 106
 107   fPrim = new TH1F("fPrim", " primary",100, 0, 200);
 108   fRef  = new TH1F("fRef", " from TR ",100, 0, 200);
 109   fRec  = new TH1F("fRec", " in ESD ",100, 0, 200);
 110
 111   fListOfHistos->Add(fOrA);
 112   fListOfHistos->Add(fOrC);
 113   fListOfHistos->Add(fMean);
 114   fListOfHistos->Add(fVertex);
 115   fListOfHistos->Add(fAmp);
 116   fListOfHistos->Add(fTime);
 117   fListOfHistos->Add(fTotalMult);
 118   fListOfHistos->Add(fMultRecRealA);
 119   fListOfHistos->Add(fMultRecRealC);
 120   fListOfHistos->Add(fMultRecTotal);
 121   fListOfHistos->Add(fPrim);
 122   fListOfHistos->Add(fRef);
 123   fListOfHistos->Add(fRec);
 124 }
 125
 126
 127 void AliT0MultiplicityTask::UserExec(Option_t *)
 128 {
 129   // Main loop
 130   // Called for each event
 131
 132
 133   // MC information
 134   AliMCEvent* mcEvent = MCEvent();
 135   if (!mcEvent) {
 136     Printf("ERROR: Could not retrieve MC event");
 137     return;
 138   }
 139
 140
 141   Int_t primaryA=0;
 142   Int_t primaryC=0;
 143   Int_t numPrim = 0;
 144   Int_t refT0=0;
 145
 146   AliStack* stack = mcEvent->Stack();
 147   // printf("AliT0MultiplicityAnalysis: Number of tracks on stack %5d\n", stack->GetNtrack());
 148   Int_t nTracks  = mcEvent->GetNumberOfTracks();
 149   for (Int_t ipart=0; ipart<nTracks; ipart++)
 150     {
 151       //    TParticle* particle = stack->Particle(ipart);
 152       AliMCParticle* track = mcEvent->GetTrack(ipart);
 153       if (!track) {
 154         Printf("ERROR: Could not receive track %d (mc loop)", ipart);
 155         continue;
 156       }
 157       Int_t label = track->GetLabel();
 158       if(stack->IsPhysicalPrimary(label) == kFALSE)
 159         continue;
 160
 161       numPrim++;
 162       Double_t eta=track->Eta();
 163
 164       if(eta<-2.97 && eta>-3.28)  primaryC++;
 165       if (eta >4.61 && eta<4.92)  primaryA++;
 166
 167        // Loop over Track References
 168        AliTrackReference* trackRef = 0;
 169
 170       for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
 171         trackRef = track->GetTrackReference(iTrackRef);
 172         if(trackRef) {
 173           Int_t detectorId = trackRef->DetectorId();
 174           if (detectorId == 7)  refT0++;
 175         }
 176       }
 177     }
 178   fTotalMult->Fill(numPrim);
 179   fPrim->Fill(primaryC+primaryA);
 180   fRef->Fill(refT0);
 181   //     printf (" tracks %i primaries %d \n",nTracks, numPrim);
 182
 183
 184   // ESD information
 185   AliVEvent* event = InputEvent();
 186   if (!event) {
 187     Printf("ERROR: Could not retrieve event");
 188     return;
 189   }
 190
 191   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 192
 193   Bool_t eventTriggered = false;
 194   ULong64_t triggerMask = esd->GetTriggerMask();
 195   // definitions from p-p.cfg
 196   ULong64_t spdFO = (1 << 14);
 197   ULong64_t v0left = (1 << 11);
 198   ULong64_t v0right = (1 << 12);
 199
 200   if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
 201     eventTriggered == true;
 202
 203     //= AliPWG0Helper::IsEventTriggered(esd, AliPWG0Helper::kMB1);
 204   //  printf("!!!!! eventTriggered %i",eventTriggered);
 205   // if(!eventTriggered) return;
 206    Double_t besttimeC=99999.;
 207   Double_t besttimeA=99999.;
 208
 209 //  Float_t coefA = 0.891466;
 210 //  Float_t coefC = 0.922772;
 211   Float_t coefA = 1;
 212   Float_t coefC = 1;
 213   Float_t sumampA=0;
 214   Float_t sumampC=0;
 215   Int_t highA=0, highC=0;
 216   const Double_t *amp = esd->GetT0amplitude();
 217   const Double_t *time = esd->GetT0time();
 218   Float_t vertex = esd->GetT0zVertex();
 219   if(vertex<999) fVertex->Fill(vertex);
 220   Float_t   mean = esd->GetT0();
 221   fMean->Fill(mean);
 222
 223   for (Int_t i=0; i<12; i++) {
 224     sumampC += amp[i];
 225     fTime->Fill(amp[i],time[i]);
 226     fAmp->Fill(amp[i]);
 227     if(time[i]<besttimeC && time[i]>0) besttimeC=time[i]; //timeC
 228   }
 229   fOrC->Fill(besttimeC);
 230
 231   for (Int_t i=12; i<24; i++){
 232     sumampA += amp[i];
 233     fAmp->Fill(amp[i]);
 234     fTime->Fill(time[i]);
 235     if(time[i]<besttimeA && time[i]>0) besttimeA=time[i]; //timeC
 236   }
 237
 238   fOrA->Fill(besttimeA);
 239
 240   fMultRecRealA ->Fill(primaryA, sumampA*coefA);
 241   fMultRecRealC ->Fill(primaryC, sumampC*coefC);
 242   fMultRecTotal->Fill(numPrim, sumampA*coefA +  sumampC*coefC);
 243   fRec->Fill(Int_t(sumampA+sumampC));
 244
 245   // Post output data.
 246   PostData(1, fListOfHistos);
 247 }
 248
 249
 250 void AliT0MultiplicityTask::Terminate(Option_t *)
 251 {
 252   // Draw result to the screen
 253   // Called once at the end of the query
 254   //  printf(" AliT0MultiplicityTask::Terminate ");
 255
 256   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
 257   if (!fListOfHistos) {
 258     Printf("ERROR: fListOfHistos not available");
 259     return;
 260   }
 261   //  printf(" before cast ");
 262   fOrA = dynamic_cast<TH1F*>(fListOfHistos->At(0));
 263   fOrC = dynamic_cast<TH1F*>(fListOfHistos->At(1));
 264   fMean  = dynamic_cast<TH1F*>(fListOfHistos->At(2));
 265   fVertex  = dynamic_cast<TH1F*>(fListOfHistos->At(3));
 266   fAmp = dynamic_cast<TH1F*>(fListOfHistos->At(4));
 267   fTime = dynamic_cast<TH1F*>(fListOfHistos->At(5));
 268   fTotalMult = dynamic_cast<TH1F*>(fListOfHistos->At(6));
 269   fMultRecTotal = dynamic_cast<TH2F*>(fListOfHistos->At(7));
 270   fMultRecRealA = dynamic_cast<TH2F*>(fListOfHistos->At(8));
 271   fMultRecRealC= dynamic_cast<TH2F*>(fListOfHistos->At(9));
 272   fPrim= dynamic_cast<TH1F*>(fListOfHistos->At(10));
 273   fRef= dynamic_cast<TH1F*>(fListOfHistos->At(11));
 274   fRec= dynamic_cast<TH1F*>(fListOfHistos->At(12));
 275
 276   TFile fc("MultHist.root\n", "RECREATE");
 277   //  printf(" File MultHist.root recreated\n");
 278   fOrC->Write();
 279   //  printf("1\n");
 280   fOrA->Write();
 281   //  printf("2\n");
 282   fMean->Write();
 283   //  printf("3\n");
 284   fVertex->Write();
 285   //  printf("4\n");
 286   fAmp->Write();
 287   //  printf("5\n");
 288   fTime->Write();
 289  fTotalMult->Write();
 290  //  printf("6\n");
 291  fMultRecTotal->Write();
 292  //  printf("7\n");
 293  fMultRecRealA->Write();
 294  // printf("8\n");
 295   fMultRecRealC->Write();
 296   // printf("9\n");
 297   fPrim->Write();
 298   fRef->Write();
 299   fRec->Write();
 300
 301   fc.Close();
 302  printf(" fc.Close()\n");
 303
 304 }