MUON/AliMUONClusterInput.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 #include <TClonesArray.h>
  19 #include <TMinuit.h>
  20
  21 #include "AliRun.h"
  22 #include "AliMUON.h"
  23 #include "AliMUONChamber.h"
  24 #include "AliMUONConstants.h"
  25 #include "AliMUONClusterInput.h"
  26 #include "AliMUONMathieson.h"
  27 #include "AliMUONRawCluster.h"
  28 #include "AliMUONDigit.h"
  29 #include "AliLog.h"
  30
  31 ClassImp(AliMUONClusterInput)
  32
  33 AliMUONClusterInput* AliMUONClusterInput::fgClusterInput = 0;
  34 TMinuit* AliMUONClusterInput::fgMinuit = 0;
  35 AliMUONMathieson* AliMUONClusterInput::fgMathieson = 0;
  36
  37 AliMUONClusterInput::AliMUONClusterInput()
  38   : TObject(),
  39     fCluster(0),
  40     fChargeCorrel(1.) // in case not defined
  41
  42 {
  43   fDigits[0]=0;
  44   fDigits[1]=0;
  45   fSegmentation[0]=0;
  46   fSegmentation[1]=0;
  47 }
  48
  49 AliMUONClusterInput* AliMUONClusterInput::Instance()
  50 {
  51 // return pointer to the singleton instance
  52     if (fgClusterInput == 0) {
  53         fgClusterInput = new AliMUONClusterInput();
  54         fgMinuit = new TMinuit(8);
  55     }
  56
  57     return fgClusterInput;
  58 }
  59
  60 AliMUONClusterInput::~AliMUONClusterInput()
  61 {
  62 // Destructor
  63     delete fgMinuit;
  64     delete fgMathieson;
  65 }
  66
  67 AliMUONClusterInput::AliMUONClusterInput(const AliMUONClusterInput& clusterInput):TObject(clusterInput)
  68 {
  69 // Protected copy constructor
  70
  71   AliFatal("Not implemented.");
  72 }
  73
  74 void AliMUONClusterInput::SetDigits(Int_t chamber, TClonesArray* dig1, TClonesArray* dig2)
  75 {
  76 // Set pointer to digits with corresponding segmentations and responses (two cathode planes)
  77     fChamber=chamber;
  78     fDigits[0]=dig1;
  79     fDigits[1]=dig2;
  80     fNDigits[0]=dig1->GetEntriesFast();
  81     fNDigits[1]=dig2->GetEntriesFast();
  82
  83     AliMUON *pMUON;
  84     AliMUONChamber* iChamber;
  85
  86     pMUON = (AliMUON*) gAlice->GetModule("MUON");
  87     if ((fSegmentationType = pMUON->WhichSegmentation()) != 1)
  88       AliFatal("Wrong segmentation type");
  89
  90     iChamber =  &(pMUON->Chamber(chamber));
  91     fgMathieson = new AliMUONMathieson();
  92
  93     fSegmentation[0]=iChamber->SegmentationModel(1);
  94     fSegmentation[1]=iChamber->SegmentationModel(2);
  95
  96     fNseg = 2;
  97     if (chamber < AliMUONConstants::NTrackingCh()) {
  98       if (chamber > 1 ) {
  99         fgMathieson->SetPitch(AliMUONConstants::Pitch());
 100         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3());
 101         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3());
 102         fChargeCorrel = AliMUONConstants::ChargeCorrel();
 103       } else {
 104         fgMathieson->SetPitch(AliMUONConstants::PitchSt1());
 105         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3St1());
 106         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3St1());
 107         fChargeCorrel = AliMUONConstants::ChargeCorrelSt1();
 108       }
 109     }
 110 }
 111
 112 void AliMUONClusterInput::SetDigits(Int_t chamber, TClonesArray* dig)
 113 {
 114 // Set pointer to digits with corresponding segmentations and responses (one cathode plane)
 115     fDigits[0]=dig;
 116     AliMUON *pMUON;
 117     AliMUONChamber* iChamber;
 118
 119     pMUON = (AliMUON*) gAlice->GetModule("MUON");
 120     if ((fSegmentationType = pMUON->WhichSegmentation()) != 1)
 121       AliFatal("Wrong segmentation type");
 122
 123     iChamber =  &(pMUON->Chamber(chamber));
 124     fSegmentation[0]=iChamber->SegmentationModel(1);
 125
 126     fNseg=1;
 127 }
 128 void AliMUONClusterInput::SetDigits(Int_t chamber, Int_t idDE, TClonesArray* dig1, TClonesArray* dig2)
 129 {
 130   // Set pointer to digits with corresponding segmentations and responses (two cathode planes)
 131     fChamber = chamber;
 132     fDetElemId = idDE;
 133     fDigits[0]  = dig1;
 134     fDigits[1]  = dig2;
 135     fNDigits[0] = dig1->GetEntriesFast();
 136     fNDigits[1] = dig2->GetEntriesFast();
 137
 138     AliMUON *pMUON;
 139     AliMUONChamber* iChamber;
 140
 141     pMUON = (AliMUON*) gAlice->GetModule("MUON");
 142     iChamber =  &(pMUON->Chamber(chamber));
 143
 144     fgMathieson = new AliMUONMathieson();
 145     if ((fSegmentationType = pMUON->WhichSegmentation()) != 2)
 146       AliFatal("Wrong segmentation type");
 147
 148     fSegmentation2[0]=iChamber->SegmentationModel2(1);
 149     fSegmentation2[1]=iChamber->SegmentationModel2(2);
 150
 151     fNseg = 2;
 152     if (chamber < AliMUONConstants::NTrackingCh()) {
 153       if (chamber > 1 ) {
 154         fgMathieson->SetPitch(AliMUONConstants::Pitch());
 155         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3());
 156         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3());
 157         fChargeCorrel = AliMUONConstants::ChargeCorrel();
 158       } else {
 159         fgMathieson->SetPitch(AliMUONConstants::PitchSt1());
 160         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3St1());
 161         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3St1());
 162         fChargeCorrel = AliMUONConstants::ChargeCorrelSt1();
 163       }
 164     }
 165 }
 166
 167 void AliMUONClusterInput::SetDigits(Int_t chamber, Int_t idDE, TClonesArray* dig)
 168 {
 169 // Set pointer to digits with corresponding segmentations and responses (one cathode plane)
 170
 171     fChamber = chamber;
 172     fDetElemId = idDE;
 173     fDigits[0] = dig;
 174
 175     AliMUON *pMUON;
 176     AliMUONChamber* iChamber;
 177
 178     pMUON = (AliMUON*) gAlice->GetModule("MUON");
 179     iChamber =  &(pMUON->Chamber(chamber));
 180     if ((fSegmentationType = pMUON->WhichSegmentation()) != 2)
 181       AliFatal("Wrong segmentation type");
 182
 183     fSegmentation2[0]=iChamber->SegmentationModel2(1);
 184
 185     fNseg=1;
 186 }
 187
 188 void  AliMUONClusterInput::SetCluster(AliMUONRawCluster* cluster)
 189 {
 190 // Set the current cluster
 191     //PH printf("\n %p \n", cluster);
 192     fCluster=cluster;
 193     Float_t qtot;
 194     Int_t   i, cath, ix, iy;
 195     AliMUONDigit* digit;
 196     fNmul[0]=cluster->GetMultiplicity(0);
 197     fNmul[1]=cluster->GetMultiplicity(1);
 198     //PH printf("\n %p %p ", fDigits[0], fDigits[1]);
 199
 200     for (cath=0; cath<2; cath++) {
 201         qtot=0;
 202         for (i=0; i<fNmul[cath]; i++) {
 203             // pointer to digit
 204             digit =(AliMUONDigit*)
 205                 (fDigits[cath]->UncheckedAt(cluster->GetIndex(i,cath)));
 206             // pad coordinates
 207             ix = digit->PadX();
 208             iy = digit->PadY();
 209             // pad charge
 210             fCharge[i][cath] = digit->Signal();
 211             // pad centre coordinates
 212 //          fSegmentation[cath]->GetPadCxy(ix, iy, x, y);
 213             // globals kUsed in fitting functions
 214             fix[i][cath]=ix;
 215             fiy[i][cath]=iy;
 216             // total charge per cluster
 217             qtot+=fCharge[i][cath];
 218             // Current z
 219             Float_t xc, yc;
 220             if (fSegmentationType == 1)
 221               fSegmentation[cath]->GetPadC(ix,iy,xc,yc,fZ);
 222             else
 223               fSegmentation2[cath]->GetPadC(fDetElemId,ix,iy,xc,yc,fZ);
 224
 225         } // loop over cluster digits
 226         fQtot[cath]=qtot;
 227         fChargeTot[cath]=Int_t(qtot);
 228     }  // loop over cathodes
 229 }
 230
 231
 232
 233 Float_t AliMUONClusterInput::DiscrChargeS1(Int_t i,Double_t *par)
 234 {
 235 // Compute the charge on first cathod only.
 236 return DiscrChargeCombiS1(i,par,0);
 237 }
 238
 239 Float_t AliMUONClusterInput::DiscrChargeCombiS1(Int_t i,Double_t *par, Int_t cath)
 240 {
 241 // par[0]    x-position of cluster
 242 // par[1]    y-position of cluster
 243
 244     Float_t q1;
 245     if (fSegmentationType == 1) {
 246
 247       fSegmentation[cath]->SetPad(fix[i][cath], fiy[i][cath]);
 248       //  First Cluster
 249       fSegmentation[cath]->SetHit(par[0],par[1],fZ);
 250       q1 = fgMathieson->IntXY(fSegmentation[cath]);
 251
 252     } else {
 253
 254       fSegmentation2[cath]->SetPad(fDetElemId, fix[i][cath], fiy[i][cath]);
 255       //  First Cluster
 256       fSegmentation2[cath]->SetHit(fDetElemId, par[0],par[1],fZ);
 257       q1 = fgMathieson->IntXY(fDetElemId, fSegmentation2[cath]);
 258    }
 259
 260    Float_t value = fQtot[cath]*q1;
 261    return value;
 262 }
 263
 264
 265 Float_t AliMUONClusterInput::DiscrChargeS2(Int_t i,Double_t *par)
 266 {
 267 // par[0]    x-position of first  cluster
 268 // par[1]    y-position of first  cluster
 269 // par[2]    x-position of second cluster
 270 // par[3]    y-position of second cluster
 271 // par[4]    charge fraction of first  cluster
 272 // 1-par[4]  charge fraction of second cluster
 273
 274   Float_t q1, q2;
 275
 276   if (fSegmentationType == 1) {
 277
 278     fSegmentation[0]->SetPad(fix[i][0], fiy[i][0]);
 279     //  First Cluster
 280     fSegmentation[0]->SetHit(par[0],par[1],fZ);
 281     q1 = fgMathieson->IntXY(fSegmentation[0]);
 282
 283     //  Second Cluster
 284     fSegmentation[0]->SetHit(par[2],par[3],fZ);
 285     q2 = fgMathieson->IntXY(fSegmentation[0]);
 286
 287   } else {
 288
 289     fSegmentation2[0]->SetPad(fDetElemId, fix[i][0], fiy[i][0]);
 290     //  First Cluster
 291     fSegmentation2[0]->SetHit(fDetElemId, par[0],par[1],fZ);
 292     q1 = fgMathieson->IntXY(fSegmentation[0]);
 293
 294     //  Second Cluster
 295     fSegmentation2[0]->SetHit(fDetElemId,par[2],par[3],fZ);
 296     q2 = fgMathieson->IntXY(fDetElemId, fSegmentation2[0]);
 297   }
 298
 299   Float_t value = fQtot[0]*(par[4]*q1+(1.-par[4])*q2);
 300   return value;
 301 }
 302
 303 Float_t AliMUONClusterInput::DiscrChargeCombiS2(Int_t i,Double_t *par, Int_t cath)
 304 {
 305 // par[0]    x-position of first  cluster
 306 // par[1]    y-position of first  cluster
 307 // par[2]    x-position of second cluster
 308 // par[3]    y-position of second cluster
 309 // par[4]    charge fraction of first  cluster - first cathode
 310 // 1-par[4]  charge fraction of second cluster
 311 // par[5]    charge fraction of first  cluster - second cathode
 312
 313   Float_t q1, q2;
 314
 315   if (fSegmentationType == 1) {
 316
 317     fSegmentation[cath]->SetPad(fix[i][cath], fiy[i][cath]);
 318     //  First Cluster
 319     fSegmentation[cath]->SetHit(par[0],par[1],fZ);
 320     q1 = fgMathieson->IntXY(fSegmentation[cath]);
 321
 322     //  Second Cluster
 323     fSegmentation[cath]->SetHit(par[2],par[3],fZ);
 324     q2 = fgMathieson->IntXY(fSegmentation[cath]);
 325
 326   } else {
 327     fSegmentation2[cath]->SetPad(fDetElemId,fix[i][cath], fiy[i][cath]);
 328     //  First Cluster
 329     fSegmentation2[cath]->SetHit(fDetElemId,par[0],par[1],fZ);
 330     q1 = fgMathieson->IntXY(fDetElemId, fSegmentation2[cath]);
 331
 332     //  Second Cluster
 333     fSegmentation2[cath]->SetHit(fDetElemId,par[2],par[3],fZ);
 334     q2 = fgMathieson->IntXY(fDetElemId, fSegmentation2[cath]);
 335   }
 336   Float_t value;
 337   if (cath==0) {
 338     value = fQtot[0]*(par[4]*q1+(1.-par[4])*q2);
 339   } else {
 340     value = fQtot[1]*(par[5]*q1+(1.-par[5])*q2);
 341   }
 342   return value;
 343 }
 344
 345 AliMUONClusterInput& AliMUONClusterInput
 346 ::operator = (const AliMUONClusterInput& rhs)
 347 {
 348 // Protected assignement operator
 349
 350   if (this == &rhs) return *this;
 351
 352   AliFatal("Not implemented.");
 353
 354   return *this;
 355 }