PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateCorrection.cxx

   1 #include "AliFMDAnalysisTaskGenerateCorrection.h"
   2 #include "AliESDEvent.h"
   3 #include "iostream"
   4 #include "AliESDFMD.h"
   5 #include "TH2F.h"
   6 #include "AliTrackReference.h"
   7 #include "AliStack.h"
   8 #include "AliFMDAnaParameters.h"
   9 #include "AliFMDStripIndex.h"
  10 #include "AliStack.h"
  11 #include "AliMCParticle.h"
  12 #include "AliMCEvent.h"
  13 //#include "AliFMDGeometry.h"
  14 #include "TArray.h"
  15 #include "AliGenEventHeader.h"
  16 #include "AliHeader.h"
  17 #include "AliFMDAnaCalibBackgroundCorrection.h"
  18 #include "AliFMDAnaCalibEventSelectionEfficiency.h"
  19 //#include "AliCDBManager.h"
  20 //#include "AliCDBId.h"
  21 //#include "AliCDBMetaData.h"
  22 #include "TSystem.h"
  23 #include "TROOT.h"
  24 #include "TAxis.h"
  25 ClassImp(AliFMDAnalysisTaskGenerateCorrection)
  26
  27 //_____________________________________________________________________
  28 AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection():
  29 AliAnalysisTaskSE(),
  30   fListOfHits(),
  31   fListOfPrimaries(),
  32   fListOfCorrection(),
  33   fVertexBins(),
  34   fLastTrackByStrip(0),
  35   fHitsByStrip(0),
  36   fZvtxCut(10),
  37   fNvtxBins(10),
  38   fNbinsEta(200),
  39   fBackground(0),
  40   fEventSelectionEff(0)
  41 {
  42   // Default constructor
  43 }
  44 //_____________________________________________________________________
  45 AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection(const char* name):
  46   AliAnalysisTaskSE(name),
  47   fListOfHits(),
  48   fListOfPrimaries(),
  49   fListOfCorrection(),
  50   fVertexBins(),
  51   fLastTrackByStrip(0),
  52   fHitsByStrip(0),
  53   fZvtxCut(10),
  54   fNvtxBins(10),
  55   fNbinsEta(200),
  56   fBackground(0),
  57   fEventSelectionEff(0)
  58 {
  59
  60   DefineOutput(1, TList::Class());
  61   DefineOutput(2, TList::Class());
  62   DefineOutput(3, TH1F::Class());
  63   DefineOutput(4, TList::Class());
  64 }
  65 //_____________________________________________________________________
  66 void AliFMDAnalysisTaskGenerateCorrection::UserCreateOutputObjects()
  67 {
  68 // Create the output containers
  69 //
  70
  71   std::cout<<"Creating output objects"<<std::endl;
  72   for(Int_t iring = 0; iring<2;iring++) {
  73     Char_t ringChar = (iring == 0 ? 'I' : 'O');
  74     Int_t nSec = (iring == 1 ? 40 : 20);
  75     for(Int_t v=0; v<fNvtxBins;v++) {
  76
  77       TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
  78                                       Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
  79                                       fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
  80       hPrimary->Sumw2();
  81       fListOfPrimaries.Add(hPrimary);
  82     }
  83   }
  84
  85
  86
  87   for(Int_t det =1; det<=3;det++) {
  88     Int_t nRings = (det==1 ? 1 : 2);
  89     for(Int_t ring = 0;ring<nRings;ring++) {
  90       Int_t nSec = (ring == 1 ? 40 : 20);
  91       Char_t ringChar = (ring == 0 ? 'I' : 'O');
  92       TH1F* doubleHits = new TH1F(Form("DoubleHits_FMD%d%c",det,ringChar),
  93                                   Form("DoubleHits_FMD%d%c",det,ringChar),fNbinsEta, -6,6);
  94       TH1F* allHits = new TH1F(Form("allHits_FMD%d%c",det,ringChar),
  95                                Form("allHits_FMD%d%c",det,ringChar), fNbinsEta, -6,6);
  96
  97       doubleHits->Sumw2();
  98       allHits->Sumw2();
  99       fListOfHits.Add(allHits);
 100       fListOfHits.Add(doubleHits);
 101
 102       for(Int_t v=0; v<fNvtxBins;v++) {
 103         TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 104                                Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 105                                fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
 106         hHits->Sumw2();
 107         fListOfHits.Add(hHits);
 108
 109       }
 110     }
 111   }
 112
 113   TH1F* hEventsSelected  = new TH1F("EventsSelected","EventsSelected",fNvtxBins,0,fNvtxBins);
 114   TH1F* hEventsAll    = new TH1F("EventsAll","EventsAll",fNvtxBins,0,fNvtxBins);
 115   TH1F* hEventsSelectedVtx  = new TH1F("EventsSelectedVtx","EventsSelectedVtx",fNvtxBins,0,fNvtxBins);
 116   TH1F* hEventsSelectedTrigger  = new TH1F("EventsSelectedTrigger","EventsSelectedTrigger",fNvtxBins,0,fNvtxBins);
 117
 118
 119
 120   //  TH1F* hTriggered    = new TH1F("Triggered","Triggered",fNvtxBins,0,fNvtxBins);
 121   // TH1F* hTriggeredAll = new TH1F("TriggeredAll","TriggeredAll",fNvtxBins,0,fNvtxBins);
 122   hEventsSelected->Sumw2();
 123   hEventsAll->Sumw2();
 124   fListOfHits.Add(hEventsSelected);
 125   fListOfHits.Add(hEventsSelectedVtx);
 126   fListOfHits.Add(hEventsSelectedTrigger);
 127   fListOfPrimaries.Add(hEventsAll);
 128
 129   // fListOfHits.Add(hTriggered);
 130   //  fListOfPrimaries.Add(hTriggeredAll);
 131
 132   fVertexBins.SetName("VertexBins");
 133   fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 134
 135 }
 136 //_____________________________________________________________________
 137 void AliFMDAnalysisTaskGenerateCorrection::Init()
 138 {
 139   fLastTrackByStrip.Reset(-1);
 140
 141
 142 }
 143 //_____________________________________________________________________
 144 void AliFMDAnalysisTaskGenerateCorrection::UserExec(Option_t */*option*/)
 145 {
 146
 147   fLastTrackByStrip.Reset(-1);
 148   fHitsByStrip.Reset(0);
 149   AliMCEvent* mcevent = MCEvent();
 150
 151   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 152
 153
 154   AliESDEvent* esdevent = (AliESDEvent*)InputEvent();
 155   Double_t esdvertex[3];
 156   Bool_t vtxStatus =  pars->GetVertex(esdevent,esdvertex);
 157
 158   AliMCParticle* particle = 0;
 159   AliStack* stack = mcevent->Stack();
 160
 161   UShort_t det,sec,strip;
 162   Char_t   ring;
 163
 164   Int_t nTracks                = mcevent->GetNumberOfTracks();
 165   AliHeader* header            = mcevent->Header();
 166   AliGenEventHeader* genHeader = header->GenEventHeader();
 167
 168   TArrayF vertex;
 169   genHeader->PrimaryVertex(vertex);
 170
 171   if(TMath::Abs(vertex.At(2)) > fZvtxCut)
 172     return;
 173
 174   Double_t delta           = 2*fZvtxCut/fNvtxBins;
 175   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 176   Int_t    vertexBin       = (Int_t)vertexBinDouble;
 177
 178   // Vertex determination correction
 179   TH1F* hEventsSelected           = (TH1F*)fListOfHits.FindObject("EventsSelected");
 180   TH1F* hEventsSelectedVtx        = (TH1F*)fListOfHits.FindObject("EventsSelectedVtx");
 181   TH1F* hEventsSelectedTrigger    = (TH1F*)fListOfHits.FindObject("EventsSelectedTrigger");
 182   TH1F* hEventsAll                = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
 183
 184   // TH1F* hTriggered      = (TH1F*)fListOfHits.FindObject("Triggered");
 185   //  TH1F* hTriggeredAll   = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
 186
 187   Bool_t vtxFound = kTRUE;
 188   if(!vtxStatus)
 189     vtxFound = kFALSE;
 190
 191   Bool_t isTriggered = pars->IsEventTriggered(esdevent);
 192
 193   if(vtxFound && isTriggered) hEventsSelected->Fill(vertexBin);
 194
 195   if(vtxFound) hEventsSelectedVtx->Fill(vertexBin);
 196   if(isTriggered) hEventsSelectedTrigger->Fill(vertexBin);
 197
 198   hEventsAll->Fill(vertexBin);
 199
 200   //  if(!vtxFound || !isTriggered) return;
 201
 202   for(Int_t i = 0 ;i<nTracks;i++) {
 203     particle = (AliMCParticle*) mcevent->GetTrack(i);
 204
 205     if(!particle)
 206       continue;
 207
 208     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
 209
 210
 211       TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
 212       TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin));
 213       hPrimaryInner->Fill(particle->Eta(),particle->Phi());
 214       hPrimaryOuter->Fill(particle->Eta(),particle->Phi());
 215     }
 216
 217     for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
 218
 219       AliTrackReference* ref = particle->GetTrackReference(j);
 220
 221       if(ref->DetectorId() != AliTrackReference::kFMD)
 222         continue;
 223       AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
 224       Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
 225       if(particle->Charge() != 0 && i != thisStripTrack ) {
 226
 227         Float_t phi = pars->GetPhiFromSector(det,ring,sec);
 228         Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));
 229
 230         TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
 231         hHits->Fill(eta,phi);
 232         Float_t nstrips = (ring =='O' ? 256 : 512);
 233         fHitsByStrip(det,ring,sec,strip) +=1;
 234         TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
 235         TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
 236
 237         if(fHitsByStrip(det,ring,sec,strip) == 1)
 238           allHits->Fill(eta);
 239
 240         doubleHits->Fill(eta);
 241
 242         fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
 243         if(strip >0)
 244           fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
 245         if(strip < (nstrips - 1))
 246           fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
 247       }
 248     }
 249
 250   }
 251
 252
 253   PostData(1, &fListOfHits);
 254   PostData(2, &fListOfPrimaries);
 255   PostData(3, &fVertexBins);
 256 }
 257 //_____________________________________________________________________
 258 void AliFMDAnalysisTaskGenerateCorrection::Terminate(Option_t */*option*/)
 259 {
 260   /*  TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
 261   TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
 262
 263   doubleHits->Divide(allHits);
 264   GenerateCorrection();
 265   PostData(1, &fListOfHits);
 266   PostData(4, &fListOfCorrection);*/
 267
 268 }
 269 //_____________________________________________________________________
 270 void AliFMDAnalysisTaskGenerateCorrection::GenerateCorrection() {
 271
 272   fBackground         = new AliFMDAnaCalibBackgroundCorrection();
 273   fEventSelectionEff  = new AliFMDAnaCalibEventSelectionEfficiency();
 274
 275   //TH1F* hTriggered      = (TH1F*)fListOfHits.FindObject("Triggered");
 276   //TH1F* hTriggeredAll   = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
 277
 278   TH1F* hEventsSelected           = (TH1F*)fListOfHits.FindObject("EventsSelected");
 279   TH1F* hEventsSelectedVtx        = (TH1F*)fListOfHits.FindObject("EventsSelectedVtx");
 280   TH1F* hEventsSelectedTrigger    = (TH1F*)fListOfHits.FindObject("EventsSelectedTrigger");
 281   TH1F* hEventsAll                = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
 282
 283   //  hEventsAll->Divide(hEventsAll,hEventsSelected,1,1,"B");
 284   hEventsSelectedVtx->Divide(hEventsAll);
 285   hEventsSelectedTrigger->Divide(hEventsAll);
 286
 287   for(Int_t i = 1; i<=hEventsSelected->GetNbinsX(); i++) {
 288     if(hEventsSelected->GetBinContent(i) == 0 )
 289       continue;
 290     Float_t b    = hEventsSelected->GetBinContent(i);
 291     Float_t db   = hEventsSelected->GetBinError(i);
 292     Float_t sum  = hEventsAll->GetBinContent(i);
 293     Float_t dsum = hEventsAll->GetBinError(i);
 294     Float_t a    = sum-b;
 295     Float_t da   = TMath::Sqrt(TMath::Power(db,2) + TMath::Power(dsum,2));
 296
 297     Float_t cor  = sum / b;
 298     Float_t ecor = TMath::Sqrt(TMath::Power(da,2) + TMath::Power(a/(b*db),2)) / b;
 299
 300     hEventsAll->SetBinContent(i,cor);
 301     hEventsAll->SetBinError(i,ecor);
 302
 303   }
 304
 305   fEventSelectionEff->SetCorrection(hEventsAll);
 306
 307   for(Int_t det= 1; det <=3; det++) {
 308     Int_t nRings = (det==1 ? 1 : 2);
 309
 310     for(Int_t iring = 0; iring<nRings; iring++) {
 311       Char_t ring = (iring == 0 ? 'I' : 'O');
 312       TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
 313       TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
 314       allHits->Divide(doubleHits);
 315
 316       fBackground->SetDoubleHitCorrection(det,ring,allHits);
 317
 318       for(Int_t vertexBin=0;vertexBin<fNvtxBins  ;vertexBin++) {
 319         TH2F* hHits          = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
 320         TH2F* hPrimary  = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",ring,vertexBin));
 321         TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
 322         hCorrection->Divide(hPrimary);
 323         /*for(Int_t i = 1; i<=hCorrection->GetNbinsX();i++)  {
 324           for(Int_t j=1; j<=hCorrection->GetNbinsY();j++) {
 325             if(hCorrection()->GetBinContent(i,j) == 0)
 326               continue;
 327             Float_t a = h
 328
 329
 330           }
 331         }
 332         */
 333
 334         hCorrection->SetTitle(hCorrection->GetName());
 335         fListOfCorrection.Add(hCorrection);
 336         fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
 337
 338
 339       }
 340
 341     }
 342   }
 343   TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 344   fBackground->SetRefAxis(&refAxis);
 345
 346 }
 347 //_____________________________________________________________________
 348 void AliFMDAnalysisTaskGenerateCorrection::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /*runNo*/) {
 349
 350   TFile infile(filename);
 351   TH1F* hVertex = (TH1F*)infile.Get("VertexBins");
 352   fZvtxCut = hVertex->GetXaxis()->GetXmax();
 353   fNvtxBins = hVertex->GetXaxis()->GetNbins();
 354   fVertexBins.SetName("VertexBins");
 355   fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
 356
 357   TList* listOfHits = (TList*)infile.Get("Hits");
 358   TList* listOfPrim = (TList*)infile.Get("Primaries");
 359
 360   TH1F* hEventsSelected           = (TH1F*)listOfHits->FindObject("EventsSelected");
 361   TH1F* hEventsSelectedVtx        = (TH1F*)listOfHits->FindObject("EventsSelectedVtx");
 362   TH1F* hEventsSelectedTrigger    = (TH1F*)listOfHits->FindObject("EventsSelectedTrigger");
 363   TH1F* hEventsAll                = (TH1F*)listOfPrim->FindObject("EventsAll");
 364
 365   fListOfHits.Add(hEventsSelected);
 366   fListOfHits.Add(hEventsSelectedVtx);
 367   fListOfHits.Add(hEventsSelectedTrigger);
 368   fListOfPrimaries.Add(hEventsAll);
 369
 370   for(Int_t det =1; det<=3;det++)
 371       {
 372         Int_t nRings = (det==1 ? 1 : 2);
 373         for(Int_t ring = 0;ring<nRings;ring++)
 374           {
 375             Char_t ringChar = (ring == 0 ? 'I' : 'O');
 376             TH1F* allHits = (TH1F*)listOfHits->FindObject(Form("allHits_FMD%d%c",det,ringChar));
 377             TH1F* doubleHits = (TH1F*)listOfHits->FindObject(Form("DoubleHits_FMD%d%c",det,ringChar));
 378             fListOfHits.Add(allHits);
 379             fListOfHits.Add(doubleHits);
 380             for(Int_t v=0; v<fNvtxBins;v++)
 381               {
 382
 383                 TH2F* hHits          = (TH2F*)listOfHits->FindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v));
 384                 fListOfHits.Add(hHits);
 385               }
 386           }
 387       }
 388   for(Int_t iring = 0; iring<2;iring++) {
 389     Char_t ringChar = (iring == 0 ? 'I' : 'O');
 390     for(Int_t v=0; v<fNvtxBins;v++) {
 391
 392       TH2F* hPrimary       = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
 393       fListOfPrimaries.Add(hPrimary);
 394
 395     }
 396   }
 397   GenerateCorrection();
 398
 399   TFile fout("backgroundFromFile.root","recreate");
 400   fListOfHits.Write();
 401   fListOfPrimaries.Write();
 402   fListOfCorrection.Write();
 403   fVertexBins.Write();
 404
 405   fout.Close();
 406   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 407   if(storeInOCDB) {
 408     TFile fbg(pars->GetPath(pars->GetBackgroundID()),"RECREATE");
 409     fBackground->Write(AliFMDAnaParameters::GetBackgroundID());
 410     fbg.Close();
 411     TFile feselect(pars->GetPath(pars->GetEventSelectionEffID()),"RECREATE");
 412     fEventSelectionEff->Write(AliFMDAnaParameters::GetEventSelectionEffID());
 413     feselect.Close();
 414
 415   }
 416
 417
 418
 419
 420 }
 421 //_____________________________________________________________________
 422 //
 423 // EOF
 424 //