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