FMD/analysis/AliFMDAnalysisTaskSharing.cxx

   1
   2 #include <TROOT.h>
   3 #include <TSystem.h>
   4 #include <TInterpreter.h>
   5 #include <TChain.h>
   6 #include <TFile.h>
   7 #include <TList.h>
   8 #include <iostream>
   9 #include <TMath.h>
  10 #include "AliFMDDebug.h"
  11 #include "AliFMDAnalysisTaskSharing.h"
  12 #include "AliAnalysisManager.h"
  13 #include "AliESDFMD.h"
  14 #include "AliFMDGeometry.h"
  15 #include "AliMCEventHandler.h"
  16 #include "AliStack.h"
  17 #include "AliESDVertex.h"
  18 #include "AliMultiplicity.h"
  19 #include "AliFMDAnaParameters.h"
  20 #include "AliFMDParameters.h"
  21
  22 ClassImp(AliFMDAnalysisTaskSharing)
  23
  24 //_____________________________________________________________________
  25 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing()
  26 : fDebug(0),
  27   fESD(0x0),
  28   foutputESDFMD(),
  29   fEnergy(0),
  30   fNstrips(0),
  31   fSharedThis(kFALSE),
  32   fSharedPrev(kFALSE),
  33   fDiagList(),
  34   fStandalone(kTRUE),
  35   fEsdVertex(0),
  36   fStatus(kTRUE)
  37 {
  38   // Default constructor
  39   DefineInput (0, AliESDEvent::Class());
  40   DefineOutput(0, AliESDFMD::Class());
  41   DefineOutput(1, AliESDVertex::Class());
  42   DefineOutput(2, AliESDEvent::Class());
  43   DefineOutput(3, TList::Class());
  44 }
  45 //_____________________________________________________________________
  46 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing(const char* name, Bool_t SE):
  47     AliAnalysisTask(name, "AnalysisTaskFMD"),
  48     fDebug(0),
  49     fESD(0x0),
  50     foutputESDFMD(),
  51     fEnergy(0),
  52     fNstrips(0),
  53     fSharedThis(kFALSE),
  54     fSharedPrev(kFALSE),
  55     fDiagList(),
  56     fStandalone(kTRUE),
  57     fEsdVertex(0),
  58     fStatus(kTRUE)
  59 {
  60   fStandalone = SE;
  61   if(fStandalone) {
  62     DefineInput (0, AliESDEvent::Class());
  63     DefineOutput(0, AliESDFMD::Class());
  64     DefineOutput(1, AliESDVertex::Class());
  65     DefineOutput(2, AliESDEvent::Class());
  66     DefineOutput(3, TList::Class());
  67   }
  68 }
  69 //_____________________________________________________________________
  70 void AliFMDAnalysisTaskSharing::CreateOutputObjects()
  71 {
  72   if(!foutputESDFMD)
  73     foutputESDFMD = new AliESDFMD();
  74
  75   if(!fEsdVertex)
  76     fEsdVertex    = new AliESDVertex();
  77   //Diagnostics
  78   fDiagList.SetName("Sharing diagnostics");
  79   for(Int_t det = 1; det<=3; det++) {
  80     Int_t nRings = (det==1 ? 1 : 2);
  81
  82     for(Int_t iring = 0;iring<nRings; iring++) {
  83       Char_t ringChar = (iring == 0 ? 'I' : 'O');
  84       TH1F* hEdist        = new TH1F(Form("Edist_before_sharing_FMD%d%c", det, ringChar),
  85                                      Form("Edist_before_sharing_FMD%d%c", det, ringChar),
  86                                      1000,0,25);
  87       TH1F* hEdist_after  = new TH1F(Form("Edist_after_sharing_FMD%d%c", det, ringChar),
  88                                      Form("Edist_after_sharing_FMD%d%c", det, ringChar),
  89                                      1000,0,25);
  90
  91
  92       TH1F* hNstripsHit    = new TH1F(Form("N_strips_hit_FMD%d%c",det,ringChar),
  93                                      Form("N_strips_hit_FMD%d%c",det,ringChar),
  94                                      25,0,25);
  95       fDiagList.Add(hEdist);
  96       fDiagList.Add(hEdist_after);
  97       fDiagList.Add(hNstripsHit);
  98
  99     }
 100   }
 101 }
 102 //_____________________________________________________________________
 103 void AliFMDAnalysisTaskSharing::ConnectInputData(Option_t */*option*/)
 104 {
 105   if(fStandalone)
 106     fESD = (AliESDEvent*)GetInputData(0);
 107 }
 108 //_____________________________________________________________________
 109 void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
 110 {
 111
 112   AliESD* old = fESD->GetAliESDOld();
 113   if (old) {
 114     fESD->CopyFromOldESD();
 115   }
 116
 117   foutputESDFMD->Clear();
 118
 119   Double_t vertex[3];
 120   GetVertex(vertex);
 121   fEsdVertex->SetXYZ(vertex);
 122
 123   if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0) {
 124
 125     fStatus = kFALSE;
 126     return;
 127   }
 128   else
 129     fStatus = kTRUE;
 130   const AliMultiplicity* testmult = fESD->GetMultiplicity();
 131
 132   Int_t nTrackLets = testmult->GetNumberOfTracklets();
 133
 134   if(nTrackLets < 1000) foutputESDFMD->SetUniqueID(kTRUE);
 135   else foutputESDFMD->SetUniqueID(kFALSE);
 136
 137   AliESDFMD* fmd = fESD->GetFMDData();
 138
 139   if (!fmd) return;
 140   Int_t nHits = 0;
 141   for(UShort_t det=1;det<=3;det++) {
 142     Int_t nRings = (det==1 ? 1 : 2);
 143     for (UShort_t ir = 0; ir < nRings; ir++) {
 144       Char_t   ring = (ir == 0 ? 'I' : 'O');
 145       UShort_t nsec = (ir == 0 ? 20  : 40);
 146       UShort_t nstr = (ir == 0 ? 512 : 256);
 147
 148       TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
 149
 150       for(UShort_t sec =0; sec < nsec;  sec++) {
 151         fSharedThis      = kFALSE;
 152         fSharedPrev      = kFALSE;
 153         fEnergy = 0;
 154         fNstrips = 0;
 155         for(UShort_t strip = 0; strip < nstr; strip++) {
 156           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
 157           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
 158
 159           if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
 160
 161           Double_t eta  = fmd->Eta(det,ring,sec,strip);//EtaFromStrip(det,ring,sec,strip,vertex[2]);
 162           //std::cout<<EtaFromStrip(det,ring,sec,strip,vertex[2]) <<"    "<<fmd->Eta(det,ring,sec,strip)<<std::endl;
 163
 164           hEdist->Fill(mult);
 165           if(fmd->IsAngleCorrected())
 166             mult = mult/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip)));
 167           Float_t Eprev = 0;
 168           Float_t Enext = 0;
 169           if(strip != 0)
 170             if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
 171               Eprev = fmd->Multiplicity(det,ring,sec,strip-1);
 172               if(fmd->IsAngleCorrected())
 173                 Eprev = Eprev/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
 174             }
 175           if(strip != nstr - 1)
 176             if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
 177               Enext = fmd->Multiplicity(det,ring,sec,strip+1);
 178               if(fmd->IsAngleCorrected())
 179                 Enext = Enext/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
 180             }
 181
 182           Float_t merged_energy = GetMultiplicityOfStrip(mult,eta,Eprev,Enext,det,ring,sec,strip);
 183
 184           if(merged_energy > 0 )
 185             nHits++;
 186           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,merged_energy);
 187           foutputESDFMD->SetEta(det,ring,sec,strip,eta);
 188
 189         }
 190       }
 191     }
 192   }
 193
 194   if(fStandalone) {
 195     PostData(0, foutputESDFMD);
 196     PostData(1, fEsdVertex);
 197     PostData(2, fESD);
 198     PostData(3, &fDiagList);
 199   }
 200 }
 201 //_____________________________________________________________________
 202 Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
 203                                                           Float_t eta,
 204                                                           Float_t Eprev,
 205                                                           Float_t Enext,
 206                                                           UShort_t   det,
 207                                                           Char_t  ring,
 208                                                           UShort_t /*sec*/,
 209                                                           UShort_t /*strip*/) {
 210   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 211
 212   Float_t merged_energy = 0;
 213   Float_t nParticles = 0;
 214   Float_t cutLow  = 0.2;
 215   Float_t cutHigh = pars->GetMPV(det,ring,eta) -1*pars->GetSigma(det,ring,eta);
 216   // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
 217   Float_t Etotal  = mult;
 218
 219   //if(mult > 5)
 220   //  std::cout<<mult<<"    "<<det<<"    "<<ring<<"   "<<sec<<"    "<<strip<<std::endl;
 221
 222   if(foutputESDFMD->GetUniqueID() == kTRUE) {
 223
 224     if(mult > cutLow ) {
 225       fEnergy = fEnergy + mult;
 226       fNstrips++;
 227     }
 228   if((Enext <0.01 && fEnergy >0) || fNstrips >2 ) {
 229
 230
 231     //if((fEnergy*TMath::Cos(Eta2Theta(eta))) > cutPart || fNstrips > 1) {
 232       nParticles = 1;
 233       merged_energy = fEnergy*TMath::Cos(Eta2Theta(eta));
 234       TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
 235       hEdist->Fill(fEnergy);
 236       TH1F* hNstrips = (TH1F*)fDiagList.FindObject(Form("N_strips_hit_FMD%d%c",det,ring));
 237       hNstrips->Fill(fNstrips);
 238       //  std::cout<<Form("Merged signals %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d",Eprev, mult, Enext, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,strip,sec,ring,det )<<std::endl;
 239
 240       // }
 241     // else
 242     //std::cout<<Form("NO HIT  for  %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",Eprev, mult, Enext, fEnergy/TMath::Cos(Eta2Theta(eta)),fEnergy,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
 243
 244     fEnergy  = 0;
 245     fNstrips = 0;
 246     return merged_energy;
 247   }
 248
 249   return 0;
 250
 251   }
 252   else {
 253
 254   if(fSharedThis) {
 255     fSharedThis      = kFALSE;
 256     fSharedPrev      = kTRUE;
 257     return 0.;
 258   }
 259
 260   if(mult < cutLow) {
 261     fSharedThis      = kFALSE;
 262     fSharedPrev      = kFALSE;
 263     return 0;
 264   }
 265
 266   if(Eprev > cutLow && Eprev < cutHigh && !fSharedPrev ) {
 267     Etotal += Eprev;
 268   }
 269
 270   if(Enext > cutLow && Enext < cutHigh ) {
 271     Etotal += Enext;
 272     fSharedThis      = kTRUE;
 273   }
 274   TH1F* hEdist = (TH1F*)fDiagList.FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
 275   hEdist->Fill(Etotal);
 276
 277   Etotal = Etotal*TMath::Cos(Eta2Theta(eta));
 278   if(Etotal > 0) {
 279     merged_energy = Etotal;
 280     fSharedPrev      = kTRUE;
 281     //if(det == 3 && ring =='I')
 282     //  std::cout<<Form("Merged signals %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det )<<std::endl;
 283   }
 284     else{// if(Etotal > 0) {
 285       //if(det == 3 && ring =='I')
 286       //        std::cout<<Form("NO HIT  for  %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
 287     fSharedThis      = kFALSE;
 288     fSharedPrev      = kFALSE;
 289   }
 290   // merged_energy = mult;
 291
 292   return merged_energy;
 293   }
 294 }
 295 //_____________________________________________________________________
 296 void AliFMDAnalysisTaskSharing::GetVertex(Double_t* vertexXYZ)
 297 {
 298   const AliESDVertex* vertex = 0;
 299   vertex = fESD->GetPrimaryVertex();
 300   if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
 301     vertex = fESD->GetPrimaryVertexSPD();
 302   if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
 303     vertex = fESD->GetPrimaryVertexTPC();
 304   if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
 305     vertex = fESD->GetVertex();
 306   if (vertex && (vertexXYZ[0] != 0 || vertexXYZ[1] != 0 || vertexXYZ[2] != 0)) {
 307     vertex->GetXYZ(vertexXYZ);
 308     //std::cout<<vertex->GetName()<<"   "<< vertex->GetTitle() <<"   "<< vertex->GetZv()<<std::endl;
 309     return;
 310   }
 311   else if (fESD->GetESDTZERO()) {
 312     vertexXYZ[0] = 0;
 313     vertexXYZ[1] = 0;
 314     vertexXYZ[2] = fESD->GetT0zVertex();
 315
 316     return;
 317   }
 318
 319   return;
 320
 321 }
 322 //_____________________________________________________________________
 323 Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) {
 324
 325   Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
 326
 327   if(eta < 0)
 328     theta = theta-TMath::Pi();
 329
 330   std::cout<<"From eta2Theta: "<<theta<<"   "<<eta<<std::endl;
 331   return theta;
 332
 333
 334
 335 }
 336 //_____________________________________________________________________
 337 Double_t AliFMDAnalysisTaskSharing::EtaFromStrip(UShort_t det,
 338                                                 Char_t ring,
 339                                                 UShort_t sector,
 340                                                 UShort_t strip,
 341                                                 Double_t zvtx)
 342 {
 343
 344   AliFMDGeometry* geo = AliFMDGeometry::Instance();
 345
 346   Double_t x,y,z;
 347   geo->Detector2XYZ(det,ring,sector,strip,x,y,z);
 348
 349   Double_t r = TMath::Sqrt(x*x+y*y);
 350
 351   Double_t z_real      = z-zvtx;
 352   Double_t theta       = TMath::ATan2(r,z_real);
 353   // std::cout<<"From EtaFromStrip "<<theta<<std::endl;
 354   Double_t eta         =  -1*TMath::Log(TMath::Tan(0.5*theta));
 355
 356   // std::cout<<det<<"   "<<ring<<"   "<<sector<<"   "<<strip<<"   "<<r<<"    "<<z_real<<"   "<<theta<<"    "<<eta<<std::endl;
 357
 358   return eta;
 359 }
 360 //_____________________________________________________________________
 361 //
 362 // EOF
 363 //